Skip to content

Commit

Permalink
equinox
Browse files Browse the repository at this point in the history
  • Loading branch information
@caiobegotti
caiobegotti committed Oct 18, 2024
1 parent cf1a735 commit 769046c
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 123 deletions.
190 changes: 73 additions & 117 deletions ABV.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <YEAR> [--quiet] [--clean] [--help]
Expand All @@ -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 <YEAR> [--quiet] [--clean] [--help]")
sys.exit(1)
Expand Down
13 changes: 7 additions & 6 deletions README.markdown
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,23 +5,23 @@
- [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)
- [Sources and credits](#sources-and-credits)

<!-- TOC end -->

ABV will graph aurora borealis visibility throughout years, with the planetary [_K_<sub>p</sub>-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_<sub>p</sub>-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_<sub>p</sub>-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.

**There you have it:** pay attention to the roughly 7-years or 11-years long solar cycles suggested by the sunspots trend, plus streak of geomagnetic storms confirmed by _K_<sub>p</sub>-index above certain thresholds.

![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.

<!-- TOC --><a name="running-it"></a>
## Running it
Expand Down Expand Up @@ -65,13 +65,14 @@ pip install matplotlib
open kp_index_2024.png
```

<!-- TOC --><a name="caveats"></a>
#### Caveats
<!-- TOC --><a name="notes-and-caveats"></a>
#### 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

<!-- TOC --><a name="tips"></a>
### Tips
Expand Down
Binary file modified example.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

0 comments on commit 769046c

Please sign in to comment.