diff --git a/ABV.py b/ABV.py index d41a7fc..1413e9c 100755 --- a/ABV.py +++ b/ABV.py @@ -12,195 +12,157 @@ import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import CubicSpline +from matplotlib import rcParams # constants -KP_THRESHOLD = 5 # threshold for kp index to highlight -CACHE_DIR = './cache' # directory to store cached data -SUNSPOT_FILE = f'{CACHE_DIR}/SN_y_tot_V2.0.csv' # path to the sunspot csv file +KP_THRESHOLD = 5 +CACHE_DIR = './cache' +SUNSPOT_FILE = f'{CACHE_DIR}/SN_y_tot_V2.0.csv' -# ensure the cache directory exists -if not os.path.exists(CACHE_DIR): - os.makedirs(CACHE_DIR) +# ensure cache directory exists +os.makedirs(CACHE_DIR, exist_ok=True) -# function to handle logging, respects --quiet flag def log(message, quiet=False): if not quiet: print(message) -# load kp index data from a local json file or fetch from api if not cached def load_or_fetch_kp_index(year, quiet=False): cache_file = f'{CACHE_DIR}/kp_index_{year}.json' - # check if the year is in the future current_year = datetime.now().year if year > current_year: log(f"Year {year} not supported by data or not current", quiet) sys.exit(1) - # check if cache file exists if os.path.exists(cache_file): log(f"Loading Kp index data for {year} from cache", quiet) with open(cache_file, 'r') as f: data = json.load(f) times = [datetime.strptime(t, '%Y-%m-%dT%H:%M:%SZ') for t in data['times']] - kp_index = data['kp_index'] - return times, kp_index + kp_values = data['kp_index'] + return times, kp_values else: try: - # fetch kp index data from the api log(f"Fetching Kp index data for {year} from the API", quiet) - start_date = f'{year}-01-01' - end_date = f'{year}-12-31' - - # get kp index data using the external script - times, kp_index, status = getKpindex(start_date, end_date, 'Kp', 'def') - - # convert time strings to datetime objects + start_date, end_date = f'{year}-01-01', f'{year}-12-31' + times, kp_values, status = getKpindex(start_date, end_date, 'Kp', 'def') times = [datetime.strptime(t, '%Y-%m-%dT%H:%M:%SZ') for t in times] - # save the data to the cache file with open(cache_file, 'w') as f: - json.dump({'times': [t.strftime('%Y-%m-%dT%H:%M:%SZ') for t in times], 'kp_index': kp_index}, f) + json.dump({'times': [t.strftime('%Y-%m-%dT%H:%M:%SZ') for t in times], 'kp_index': kp_values}, f) - return times, kp_index - except NameError: - log(f"Year {year} not supported by data or not current", quiet) + return times, kp_values + except Exception: + log(f"Error fetching Kp index data for {year}.", quiet) sys.exit(1) -# load sunspot data for a range of years from the csv file def load_sunspot_data(start_year, end_year, quiet=False): if not os.path.exists(SUNSPOT_FILE): log(f"Sunspot data file not found: {SUNSPOT_FILE}", quiet) return None sunspot_data = {} - # read the csv file without skipping the header (since it's crude and headerless) with open(SUNSPOT_FILE, 'r') as csvfile: reader = csv.reader(csvfile, delimiter=';') - for row in reader: - row_year = float(row[0]) # use float to capture the `.5` - int_year = int(row_year) # convert to integer to ignore the `.5` - if start_year <= int_year <= end_year: - sunspot_value = float(row[1]) # second column is the sunspot value - sunspot_data[int_year] = sunspot_value + try: + row_year = int(float(row[0])) # conversions to ignore the .5 float in the data + if start_year <= row_year <= end_year: + sunspot_data[row_year] = float(row[1]) # second column is the sunspot value + except (ValueError, IndexError): + continue return sunspot_data -# plot kp index for a specific year with local cache support and overlay sunspot data def plot_kp_index_for_year(year, quiet=False, clean=False): - # load or fetch kp index data - times, kp_index = load_or_fetch_kp_index(year, quiet) + times, kp_values = load_or_fetch_kp_index(year, quiet) - # debugging: print the retrieved data length log(f"Retrieved {len(times)} entries for year {year}", quiet) - # load sunspot data for the current year, the previous 5 years, and the next year start_year = year - 5 - prev_year = year - 1 # fetch the previous year's value - next_year = year + 1 # fetch the next year's value - sunspot_data = load_sunspot_data(start_year, next_year, quiet) + sunspot_data = load_sunspot_data(start_year, year + 1, quiet) - # handle missing sunspot data if not sunspot_data: log("Sunspot data not available, skipping sunspot plot.", quiet) - sunspot_data = {} # set to an empty dictionary to avoid future errors + sunspot_data = {} # set to an empty dictionary to avoid future errors # find the last available year for sunspot data last_available_year = max(sunspot_data.keys(), default=None) + sunspot_values = [sunspot_data.get(yr) for yr in range(start_year, last_available_year + 1)] if last_available_year else None + if last_available_year and last_available_year < year: log(f"Sunspot data only available up to {last_available_year}", quiet) - else: - last_available_year = year # use the current year if data exists for it - # get the sunspot values in the range [start_year, last_available_year] - sunspot_years = [yr for yr in range(start_year, last_available_year + 1) if yr in sunspot_data] # only use years with data - sunspot_values = [sunspot_data.get(yr) for yr in sunspot_years] # skip years without sunspot data + plt.figure(figsize=(20, 10)) # log sunspot data for the previous year (if available) - if prev_year in sunspot_data: - log(f"Previous year sunspot: {sunspot_data[prev_year]}", quiet) + if year - 1 in sunspot_data: + log(f"Sunspots before {year}: {sunspot_data[year - 1]}", quiet) + + plt.plot(times, kp_values, color='#dddddd', linewidth=0.5, zorder=1) # output sunspot value for logging (only if available for the current year) if year in sunspot_data: - log(f"Sunspot value for {year}: {sunspot_data[year]}", quiet) - - # handle next year's sunspot value with fallback to the current year if missing - if last_available_year == year: - next_year_sunspot = sunspot_data.get(next_year, sunspot_data.get(year, None)) - if next_year_sunspot: - log(f"Next year sunspot: {next_year_sunspot}", quiet) - - # use cubic spline interpolation for a smooth transition between sunspot values (if data exists) - if sunspot_years and sunspot_values: - cubic_spline = CubicSpline(np.linspace(0, len(times) - 1, len(sunspot_values)), sunspot_values) - sunspot_curve_y = cubic_spline(np.linspace(0, len(times) - 1, len(times))) - else: - sunspot_curve_y = None # no sunspot data to plot - - # prepare a figure with increased size for better visibility - plt.figure(figsize=(20, 10)) # larger figure size + log(f"Sunspots for {year}: {sunspot_data[year]}", quiet) - # plot all kp index values (line) in light gray - plt.plot(times, kp_index, color='#dddddd', linewidth=0.5) # base line in light gray + for i, value in enumerate(kp_values): + if value > KP_THRESHOLD: + plt.axvline(x=times[i], color='#1F51FF', linewidth=1, zorder=3) + plt.scatter([times[i]], [value], color='r', marker='o', linewidth=2, zorder=5) - # highlight values above threshold in blue and mark them with red dashes - for i in range(len(kp_index)): - if kp_index[i] > KP_THRESHOLD: - plt.axvline(x=times[i], color='b', linewidth=3, zorder=3) # draw vertical line for values > threshold - # add horizontal markers for the kp index level - plt.scatter([times[i]], [kp_index[i]], color='r', marker='o', zorder=5) # red markers + if year + 1 in sunspot_data or year in sunspot_data: + next_year_sunspot = sunspot_data.get(year + 1, sunspot_data.get(year, None)) + if next_year_sunspot: + log(f"Sunspots after {year}: {next_year_sunspot}", quiet) - # keep the y-axis fixed at 0-9 as requested plt.ylim(0, 9) - # overlay the green line smoothly transitioning from 5 years ago to the input year (if sunspot data exists) - if sunspot_curve_y is not None: + # draw the green band smoothly transitioning from 5 years ago to the input year (if sunspot data exists) + if sunspot_values: + cubic_spline = CubicSpline(np.linspace(0, len(times) - 1, len(sunspot_values)), sunspot_values) + sunspot_curve_y = cubic_spline(np.linspace(0, len(times) - 1, len(times))) scaled_curve_y = 9 * (sunspot_curve_y - min(sunspot_curve_y)) / (max(sunspot_curve_y) - min(sunspot_curve_y)) - plt.plot(times, scaled_curve_y, color='g', linewidth=2, linestyle='-', zorder=2, label=f"Sunspot activity trend between {start_year} and {last_available_year}") + plt.plot(times, scaled_curve_y, color='#AFE1AF', linewidth=15, label=f"Sunspots trend", zorder=1) - # create custom legend handles for kp index, sunspot, and events - base_line_handle = plt.Line2D([], [], color='#dddddd', linewidth=2, label="Magnetosphere disturbance") - line_handle = plt.Line2D([0], [0], color='b', linewidth=2, label=f'Good visibility days in {year}') - sunspot_handle = plt.Line2D([0], [0], color='g', linewidth=2, label=f'Sunspot activity trend between {start_year} and {last_available_year}' if sunspot_curve_y is not None else '') - event_handle = plt.Line2D([0], [0], color='r', marker='o', linestyle='None', markersize=5, label='Aurora events') + # set sketch parameters only for the yellow zig-zag equinox lines + rcParams['path.sketch'] = (50, 250, 1) # (scale, length of wiggle, randomness) - more amplitude, smoother - plt.title(f'Aurora borealis visibility throughout {year}') # dynamic year in title - plt.ylabel(r'$\it{K}_p$ index') # k italic, p subscript, index lowercase - plt.grid(axis='y') # enable only horizontal grid lines + # plot zig-zag yellow lines for march 21st and september 21st + highlight_yellow = '#FFFAA0' + for date_str in [f"{year}-03-21", f"{year}-09-21"]: + equinox_date = datetime.strptime(date_str, '%Y-%m-%d') + if equinox_date in times: + plt.axvline(x=equinox_date, color=highlight_yellow, linewidth=15, label='Equinox season start', zorder=1) - # format the x-axis to show only the month and day, hide the year - date_format = mdates.DateFormatter('%b %d') # e.g., 'jan 01' - plt.gca().xaxis.set_major_formatter(date_format) + # restore normal line settings for everything else + rcParams['path.sketch'] = None - # set x-ticks for every day in the year with padding - plt.xticks(times, rotation=90, fontsize=4) # keep all dates but rotate for visibility - plt.gca().xaxis.set_major_locator(plt.MaxNLocator(nbins=365, integer=True)) # allow for all ticks with spacing + plt.title(f'Aurora borealis visibility throughout {year}') + plt.ylabel(r'$\it{K}_p$ index') # k italic, p subscript, index lowercase + plt.grid(axis='y') - # adjust x-axis limits to add space around the plot - plt.xlim([times[0] - timedelta(days=1), times[-1] + timedelta(days=1)]) # adding padding on the x-axis + plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %d')) # format the x-axis to show only the month and day, hide the year + plt.xticks(times, rotation=90, fontsize=4) # keep all dates but rotate for visibility + plt.gca().xaxis.set_major_locator(plt.MaxNLocator(nbins=365, integer=True)) # allow for all ticks with spacing + plt.xlim([times[0] - timedelta(days=1), times[-1] + timedelta(days=1)]) # adding padding on the x-axis - # hide the x-axis label - plt.xlabel('') - - # if --clean option is set, skip adding the legend if not clean: - # add a custom legend with a white background, including the green line for sunspot activity - legend_handles = [sunspot_handle, base_line_handle, line_handle, event_handle] # adjust order of handles - legend = plt.legend(handles=legend_handles, loc='upper left', facecolor='#ffffff', framealpha=0.75) - legend.get_frame().set_edgecolor('#000000') # add a border to the legend box + legend_handles = [ + plt.Line2D([0], [0], color='#dddddd', linewidth=5, label="Magnetosphere disturbance"), + plt.Line2D([0], [0], color='#AFE1AF', linewidth=5, label=f'Sunspots trend between {start_year}-today'), + plt.Line2D([0], [0], color=highlight_yellow, linewidth=5, label='Equinox season start'), + plt.Line2D([0], [0], color='#1F51FF', linewidth=5, label=f'Good visibility days this year'), + plt.Line2D([0], [0], color='r', marker='o', linestyle='None', markersize=8, label='Aurora events') + ] + plt.legend(handles=legend_handles, loc='upper left', facecolor='#ffffff', framealpha=0.75) - # save the plot as a png file without showing it file_name = f'kp_index_{year}.png' plt.savefig(file_name, format='png', dpi=300, bbox_inches='tight') - - # close the plot to avoid display plt.close() log(f"Plot saved as {file_name}", quiet) -# print the help message for script usage def print_help(): help_message = """ Usage: python script.py [--quiet] [--clean] [--help] @@ -214,20 +176,14 @@ def print_help(): print(help_message) if __name__ == '__main__': - # check for the --quiet, --clean, and --help flags + if '--help' in sys.argv or len(sys.argv) == 1: + print_help() + sys.exit(0) + quiet = '--quiet' in sys.argv clean = '--clean' in sys.argv - show_help = '--help' in sys.argv - - # remove --quiet, --clean, and --help from arguments args = [arg for arg in sys.argv if arg not in ['--quiet', '--clean', '--help']] - # show help message if --help is provided or no arguments are given - if show_help or len(args) == 1: - print_help() - sys.exit(0) - - # check if the user provided a year argument if len(args) != 2: print("Usage: python script.py [--quiet] [--clean] [--help]") sys.exit(1) diff --git a/README.markdown b/README.markdown index 074b78d..4cdcf59 100644 --- a/README.markdown +++ b/README.markdown @@ -5,7 +5,7 @@ - [Running it](#running-it) * [Setup](#setup) * [Plot the aurora events, this is it!](#plot-the-aurora-events-this-is-it) - + [Caveats](#caveats) + + [Notes and caveats](#notes-and-caveats) * [Tips](#tips) - [Know more](#know-more) - [Cached data](#cached-data) @@ -13,7 +13,7 @@ -ABV will graph aurora borealis visibility throughout years, with the planetary [_K_p-index](https://en.wikipedia.org/wiki/K-index) plotted and with decent thresholds. It should also show the trend line of sunspots visible, which correlates with a better chance of seeing auroras in the upcoming year (or not). Just you wait, the current 3 years period is fantastic! +ABV will graph aurora borealis visibility throughout years, with the planetary [_K_p-index](https://en.wikipedia.org/wiki/K-index) plotted and with decent thresholds. It should also show the trend line of sunspots visible, which correlates with a better chance of seeing auroras in the upcoming year (or not). Just you wait, the current 3 years period is looking fantastic... This was created to show to my wife the best time of the year and annual average occurrences of auroras in the northern hemisphere, because she didn't believe me when I told her about the _K_p-index and multi year cycles of correlated sunspots etc. She wanted to see an aurora but nobody could tell her the best time nor place to go with minimal guarantee to see one. @@ -21,7 +21,7 @@ This was created to show to my wife the best time of the year and annual average ![Graph example for 2024](example.png "Graph example for 2024") -I'm a complete ignorant and armchair nerd about this, so take it all with a sol-sized grain of salt. +I'm a complete ignorant and armchair nerd about this, so take it all with a sol-sized grain of salt. This is all informative and nice and cool, but that's about it. ## Running it @@ -65,13 +65,14 @@ pip install matplotlib open kp_index_2024.png ``` - -#### Caveats + +#### Notes and caveats - the green trend line of sunspots doesn't match the scale of the rest of the graph, it's just informative - the red aurora events are just the meaningful geomagnetic storms recorded, it doesn't mean you can't see auroras below those dots, you definitely can but they won't be as spectacular -- the solar cycles don't have peaks necessarily, their peak look more like the M letter being that a 3 years cycle can have a "bad" year in between with much less sunspots +- the solar cycles don't have peaks necessarily, their peak looks more like the M letter being that a 3 years cycle can have a "bad" year in between with much less sunspots - once again, there is some correlation between sunspots and geomagnetic storms leading to auroras but such correlation is not strong enough to guarantee sights +- right after the equinox dates in the northern hemisphere is when auroras start to look really beautiful so the period is shown as a reminder ### Tips diff --git a/example.png b/example.png index 23e14cb..950ee09 100644 Binary files a/example.png and b/example.png differ