Generate a yearly sun graph for a location. See: http://adrn.github.io/posts/2016/07/yearly-sun-graph/

sun_plot.py

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