Molly Peeples (@astronomolly) recently tweeted some really cool plots from timeanddate.com that visualize local solar properties over a year from a location (e.g., twilights, solar noon, etc.), e.g.:

My first thought was: can we make those plots with Astropy? This post shows one way to make figures like this in Python using Astropy and Matplotlib.

First, some imports:

In [1]:
```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
pl.style.use('apw-notebook')
%matplotlib inline
import pytz
```

This code relies on a brand new feature of the development version of Astropy that uses the Google maps API to allow you to get an EarthLocation object given any arbitrary address. For example, a fully-specified street address:

In [2]:
```loc = coord.EarthLocation.of_address('1002 Holy Grail Court, St. Louis, MO')
loc.longitude, loc.latitude
```
Out[2]:
`(<Longitude -90.3064519 deg>, <Latitude 38.514261999999995 deg>)`

A city name:

In [3]:
```loc = coord.EarthLocation.of_address('Parrot, ME')
loc.longitude, loc.latitude
```
Out[3]:
`(<Longitude -69.5539404 deg>, <Latitude 45.2550505 deg>)`

Or even a point of interest:

In [4]:
```loc = coord.EarthLocation.of_address('Camelot Castle Hotel, Tintagel, UK')
# (it's only a model)
loc.longitude, loc.latitude
```
Out[4]:
`(<Longitude -4.7554539 deg>, <Latitude 50.668217 deg>)`

We'll use this feature to get an `EarthLocation` (e.g., latitude and longitude) for any city we'd like to make a sun graph for.

Next we define some more objects to define the various types of twilight. The first dictionary contains the sun altitude specifications for civil, nautical, and astronomical twilight (the range or minimum to maximum sun altitude for each). The next dictionary maps the names of the various stages of sunlight to integers that are ranked with the brightness. This will be used in plotting later on.

In [5]:
```_twil = {
'civil': [-6., 0]*u.degree,
'nautical': [-12.,-6]*u.degree,
'astronomical': [-18., -12]*u.degree
}

_twil_name_map = {
'day': 0,
'civil': 1,
'nautical': 2,
'astronomical': 3,
'night': 4
}
```

Now we can write a function to make a sun plot for a given location! At minimum, the function will need to accept an `EarthLocation`, which contains the latitude/longitude information, and a `pytz` timezone information object. We'll write the function so that, by default, it produces a plot for the entire current year by computing the solar properties for every day. But, we should let the user have some control over this -- for example, if we want lower resolution over only three months of the year starting from May 15, it should be possible to specify a custom grid of days and a start date. The full function is somewhat long, so I will break it up into chunks and explain through each chunk in succession. Skip down to the full function cell to see the whole thing.

The first few lines of the plot function look at the inputs that specify the date grid over which we will compute the solar properties. The two relevant arguments are `start` and `day_grid`, and the possible combinations and specifications are:

• `start` is not specified and `day_grid` is not specified: this is the default and will produce a plot with solar properties for every day in a full year starting on January 1 of the current year.
• `start` is specified and `day_grid` is specified: if both are specified, we simply need to add the grid of days to the input start `astropy.time.Time` object.
• `start` is specified and `day_grid` is not specified: if the start datetime is specified but the grid of days isn't, it uses a grid of 1 year with a spacing of 1 day.
• `start` is not specified but `day_grid` is specified: assumes the start is January 1 of the current year.
In [6]:
```def sun_info_plot1(earth_location, tz, start=None, day_grid=None,
n_hour_grid=512, fig=None):

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?")

# ...
```

Now we will loop through the days for which we will compute the solar properties, convert the local times to UTC, get the sun's position throughout the day and use our definitions of twilight to figure out at what times the conditions are met. We start by defining a 2D array: the 0th axis represents the grid of days throughout the year, and the 1st axis represents a grid of hours throughout each day when we will compute the solar properties.

We now loop over each day in the year. For each day, we generate a grid of hours and convert these local times to UTC times. We then use `astropy.coordinates.get_sun` to get the position of the sun at each of these times, and transform these positions to the Altitude-Azimuth system using the input location on the Earth. Solar noon and midnight are computed as extrema of the altitude. For each type of twilight, we record the times at which the twilight constraints (defined above) are met by storing an integer in the 2D array `twil_grid` (0 corresponding to daytime, 1 for civil twilight, 2 for nautical twilight, etc. as defined above). Finally, when the sun is below -18º, we set it to night time.

In [7]:
```def sun_info_plot2(earth_location, tz, start=None, day_grid=None,
n_hour_grid=512, fig=None):

# ... previous code ...

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, nautical, 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)
```

The final block of code actually generates the plot using `matplotlib`. This part is actually pretty easy because of how we set up the 2D array `twil_grid`. In anticipation of making the plot using `matplotlib.pyplot.imshow`, we filled `twil_grid` so that the values at each position in the array represent the brightness using our map from twilight name to an integer value. `imshow` then uses a colormap to turn these different numbers into different colors. The rest of the code is just to fancify various aspects of the plot.

In [8]:
```def sun_info_plot3(earth_location, tz, start=None, day_grid=None,
n_hour_grid=512, fig=None):

# ... previous code ...

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()
```

Finally, the function in its entirety:

In [10]:
```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
```

Now let's try this out on a few locations! (By the way: it'd be awesome if there were a way in Python to go from an Earth location to a timezone directly, but as far as I can tell there isn't! Instead, we have to specify the timezone name explicitly...kind of annoying).

In [10]:
```address = 'Baltimore, MD'
tz = pytz.timezone('US/Eastern')
fig = sun_info_plot(loc, tz)
```
In [ ]:
```address = 'Trondheim, Norway'
tz = pytz.timezone('Europe/Oslo')
fig = sun_info_plot(loc, tz)
```address = 'McMurdo, Antarctica'