In this example, we’ll show how to generate a sun path diagram for any given day for a preselected location. By using the skyfield library, we can obtain extremely accurate predictions of the sun’s position in the sky, and in contrast to many other libraries, we can use the same method to predict the position of other celestial bodies.
#!/usr/bin/env python3 from skyfield import api from skyfield import almanac from datetime import datetime from datetime import timedelta from matplotlib import pyplot as plt import pytz import dateutil.parser from collections import namedtuple from UliEngineering.Utils.Date import * import matplotlib.ticker as mtick # Matplotlib formatter def format_degrees(value, pos=None): return f'{value:.0f} °' ts = api.load.timescale() ephem = api.load_file('de421.bsp') sun = ephem["Sun"] earth = ephem["Earth"] # Data types AltAz = namedtuple("AltAz", ["altitude", "azimuth", "distance"]) TimeAltAz = namedtuple("TimeAltAz", ["time", "altitude", "azimuth", "distance"]) def sun_altaz(planet_at_location, ts): # Compute the sun position as seen from the observer at <location> sun_pos = planet_at_location.at(ts).observe(sun).apparent() # Compute apparent altitude & azimuth for the sun's position altitude, azimuth, distance = sun_pos.altaz() return AltAz(altitude, azimuth, distance) def sun_altaz_for_day(location, year, month, day, tz): earth_at_location = (earth + location) minutes = list(yield_minutes_on_day(year=year, month=month, day=day, tz=tz)) skyfield_minutes = ts.from_datetimes(minutes) minutely_altaz = [sun_altaz(earth_at_location, ts) for ts in skyfield_minutes] # Extract components for plotting minutely_alt = [altaz.altitude.degrees for altaz in minutely_altaz] minutely_az = [altaz.azimuth.degrees for altaz in minutely_altaz] minutely_distance = [altaz.distance for altaz in minutely_altaz] return TimeAltAz(minutes, minutely_alt, minutely_az, minutely_distance) # Random location near Munich location = api.Topos('48.324777 N', '11.405610 E', elevation_m=519) # Compute Alt/Az for two different days jun = sun_altaz_for_day(location, 2022, 6, 15, pytz.timezone("Europe/Berlin")) dec = sun_altaz_for_day(location, 2022, 12, 15, pytz.timezone("Europe/Berlin")) # Plot! plt.gca().yaxis.set_major_formatter(mtick.FuncFormatter(format_degrees)) plt.gca().xaxis.set_major_formatter(mtick.FuncFormatter(format_degrees)) plt.plot(jun.azimuth, jun.altitude, label="15th of June") plt.plot(dec.azimuth, dec.altitude, label="15th of December") plt.ylim([0, None]) # Do not show sun angles below the horizon plt.ylabel("Sun altitude") plt.xlabel("Sun azimuth") plt.title("Apparent sun position at our office") plt.grid() plt.gca().legend() plt.gcf().set_size_inches(10,5) plt.savefig("Sun path.svg")