Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Generate a yearly sun graph for a location. See: http://adrn.github.io/posts/2016/07/yearly-sun-graph/
from datetime import datetime
# Third-party
import astropy.coordinates as coord
import astropy.units as u
import astropy.time as t
import matplotlib as mpl
import matplotlib.pyplot as pl
import matplotlib.dates as mdates
import numpy as np
import pytz
def sun_info_plot(earth_location, tz, start=None, day_grid=None,
n_hour_grid=512, fig=None):
"""
Visualize how solar noon, solar midnight, and twilight (civil,
nautical, and astronomical) vary over a range of dates at a given
location.
Parameters
----------
earth_location : `~astropy.coordinates.EarthLocation`
The location to produce the info for (e.g., latitude, longitude)
as an Astropy `~astropy.coordinates.EarthLocation` object.
tz : `~pytz.tzinfo.BaseTzInfo`
The timezone of the location.
start : `~astropy.time.Time`, optional
day_grid : `~astropy.units.Quantity`, optional
n_hour_grid : int, optional
Returns
-------
fig : `matplotlib.figure.Figure`
"""
if start is None and day_grid is None: # default
year = datetime.now().year
start = t.Time('{}-01-01'.format(year), format='iso', scale='utc')
day_grid = start + np.arange(0,365+1).astype(int) * u.day
elif start is not None and day_grid is not None:
day_grid = start + day_grid
elif start is not None and day_grid is None:
day_grid = start + np.arange(0,365+1).astype(int) * u.day
elif start is None and day_grid is not None:
year = datetime.now().year
start = t.Time('{}-01-01'.format(year), format='iso', scale='utc')
day_grid = start + day_grid
else: # should never reach here
raise ValueError("How did I get here?")
twil_grid = np.zeros((len(day_grid), n_hour_grid), dtype=int)
_solar_noon = []
_solar_midnight = []
for i,day in enumerate(day_grid):
utc_offset = tz.utcoffset(day.datetime).total_seconds()*u.second
hour_grid_loc = day + np.linspace(0, 24-1E-7, n_hour_grid)*u.hour - utc_offset
loc_hr = [d.hour + d.minute/60. + d.second/3600.
for d in (hour_grid_loc + utc_offset).datetime]*u.hour
# get position of the Sun at all times during this day
sun = coord.get_sun(hour_grid_loc)
sun_altaz = sun.transform_to(coord.AltAz(location=loc))
# solar noon and midnight
_solar_noon.append(loc_hr[sun_altaz.alt.argmax()]) # Sun at max altitude
_solar_midnight.append(loc_hr[sun_altaz.alt.argmin()]) # Sun at min altitude
# civil, nautxical, astronomical twilights
for key in _twil.keys():
idx = (sun_altaz.alt >= _twil[key][0]) & (sun_altaz.alt < _twil[key][1])
twil_grid[i,idx] = _twil_name_map[key]
idx = (sun_altaz.alt < _twil['astronomical'][0])
twil_grid[i,idx] = _twil_name_map['night']
# convert lists of Quantity objects to Quantity arrays
solar_noon = u.Quantity(_solar_noon)
solar_midnight = u.Quantity(_solar_midnight)
# -------------------------------------------------------------------
# Plotting
#
if fig is None:
fig,ax = pl.subplots(1,1,figsize=(15,8))
else:
ax = fig.axes[0]
# matplotlib date trickery: see http://matplotlib.org/api/dates_api.html
xlim = mdates.date2num([day_grid.datetime.min(), day_grid.datetime.max()]).tolist()
ylim = [loc_hr.value.min(), loc_hr.value.max()]
# use imshow to visualize the stages of daylight
ax.imshow(twil_grid.T, origin='bottom', aspect='auto', interpolation='nearest',
cmap='bone_r', extent=xlim+ylim)
# don't connect discontinuities with a line, which happen when solar midnight shifts
# from just before 00:00 to just after.
idx = np.where(np.abs(np.diff(solar_midnight.to(u.hour).value)) >= 5)[0]+1
solar_midnight = np.insert(solar_midnight, idx, np.nan)
x_midnight = np.insert(mdates.date2num(day_grid.datetime), idx, np.nan)
ax.plot(day_grid.datetime, solar_noon, color='#fecc5c', marker=None, linewidth=2)
ax.plot(x_midnight, solar_midnight, color='#2b83ba', marker=None, linewidth=2)
# assign date locator / formatter to the x-axis to get proper labels
months = mdates.MonthLocator()
if day_grid[0].datetime.year == day_grid[-1].datetime.year:
date_fmt = mdates.DateFormatter('%b')
xlabel = "{:d}".format(day_grid[0].datetime.year)
else:
date_fmt = mdates.DateFormatter('%b-%Y')
xlabel = ""
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(date_fmt)
ax.yaxis.set_ticks(np.arange(0,24+2,2))
ax.tick_params(axis='both', colors='#cccccc')
[lbl.set_color("k") for lbl in ax.get_xticklabels() + ax.get_yticklabels()];
ax.set_ylim(-0.1,24.1)
ax.set_xlabel(xlabel)
fig.tight_layout()
return fig
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.