Solar position calculation

[dopplershift]

Including some kind of solar position calculator would be useful to a variety of calculations, including wet bulb globe temperate (WBGT #1374 ) and plotting theoretical solar radiation curves for meteograms (#208 ).

I'm dumping a bunch of open tabs I have from a brief investigation, where of course there are many ways to go about this:

• Rough calculation from ESRL
• Algorithm documentation from NREL
• ESRL spreadsheet calculators -- may be useful for validation
• BAMS article on the Sun4Cast system
• Paper on solar position algorithm in the journal on Solar Energy

[wholmgren]

pvlib python's solarposition.py (wrappers plus low precision algorithms) and spa.py (the NREL algorithm implemented with numba support) modules are very well tested. Here are the API docs. BSD-3, so copy away.

[dopplershift]

@wholmgren Awesome! It never occurred to me to search for this already existing, but I should have checked out pvlib. Thanks for pointing us that way.

[danieldjewell]

I'm assuming you're looking for things like the azimuth/elevation of the sun from a given Lat/Lon/elevation?

If so, I'd highly recommend Skyfield /// add'l homepage

It uses SPICE kernels from the JPL for ephem data - and it's super simple:

import skyfield.api as sapi

#Read about BSP SPICE Kernels
#Automatically downloads from the JPL/NASA
#see: https://naif.jpl.nasa.gov/naif/tutorials.html
ephem = sapi.load('de421.bsp')
timescale = sapi.load.timescale()

tsnow = timescale.now()

earth = ephem['earth']
sun = ephem['sun']

waikiki=sapi.Topos(latitude_degrees=21.281194, longitude_degrees=-157.837911, elevation_m=1)

me_at_waikiki = earth + waikiki
#(yes, making a location on earth is really as easy as adding the earth to the Topos() object!)

sunpos = me_at_waikiki.at(tsnow).observe(sun)

sunpos.apparent().altaz()
#or any number of things...

Check out the documentation -- it really explains it well! (There's also astropy and poliastro but they are probably a bit overboard for your application 😁)

[Schamschula]

I stumbled upon this issue researching a matplotlib deprecation issue in the Meteogram example https://unidata.github.io/MetPy/latest/examples/meteogram_metpy.html

The weewx weather station software uses the ephem Python package. Works great.

[sgdecker]

I've used Astropy for this purpose in the past.

[dopplershift]

astropy would be a "nuclear-powered sledgehammer" for what MetPy would use it for. 😉 PyEphem's docs say to use skyfield, so skyfield is what I figure we would depend on. It's pure Python, already on conda-forge, and I have good faith in the project's lead dev (Brandon Rhodes). So I'd lean that direction or borrowing pvlib's code.

[wholmgren]

FWIW we created a pvlib asv benchmarker in the time since this issue first came up. The solar position functions are relatively well covered, but pvlib does not wrap skyfield so we can't say anything about that comparison.

One thing to keep in mind is that pvlib's functions, like many other libraries, are tailored to time series at a specific location. Occasionally people from the weather/climate community will come to pvlib and want to run simulations over grids. See this SO post for example. That means they have to hack at the lower level functions or they need to loop over each grid point.

I think metpy could do a great service by reimplementing a good enough algorithm in a vectorized way that supported both spatial and temporal dimensions. A working definition of good enough:

pvlib analytical, NAM/GFS algorithm, WRF ~3.5 algorithm < good enough <? NREL algorithm, WRF ~3.6+ < astropy, pyephem, skyfield

(~approximate for WRF versions because I don't remember exactly when the more accurate solar position algorithm was copied from the WRF Solar fork to the main WRF code base)

pvlib.solarposition.ephemeris might be the sweet spot, and it can probably work with array latitude/longitude with minor changes, but unfortunately it has the least-solid reference. It compares favorably with the "better" algorithms for all locations and times that I've used, but I'd recommend documenting the differences at more pole-ward locations before relying on it.

[deeplycloudy]

+1 for paying attention to vectorizing as suggested by @wholmgren.

The (related) example below, adapted from the skyfield docs, calculates time since solar noon. A user might need many of these time offests for arbitrarily located weather stations.

import datetime as dt
from skyfield import almanac
import skyfield.api as skyapi

def local_solar_noon_offset(utc, longitude, latitude=0.0, apparent=True):
    """ Calculate the time offset from local solar noon for a position on earth.
    
    Returns utc - local_noon_utc.
    
    utc: Python datetime, in UTC
    longitude: east longitude in degrees.
    latitude (optional, default 0.0): north latitude in degrees.
    
    Note that UTC is not exactly the same as GMT=UT1; they are different by about 1 second.
    """
    utc_tz = utc.replace(tzinfo=skyapi.utc)
    ts = skyapi.load.timescale()
    midnight = utc_tz.replace(hour=0, minute=0, second=0, microsecond=0)
    next_midnight = midnight + dt.timedelta(days=1)
    t0 = ts.from_datetime(midnight)
    t1 = ts.from_datetime(next_midnight)
    
    location = skyapi.wgs84.latlon(latitude * skyapi.N, longitude * skyapi.E)

    eph = skyapi.load('de421.bsp')
    f = almanac.meridian_transits(eph, eph['Sun'], location)
    times, events = almanac.find_discrete(t0, t1, f)
    # Select transits instead of antitransits.
    times = times[events == 1]

    if len(times) > 0:
        t_noon = times[0].utc_datetime()    
        return utc_tz - t_noon
    else:
        return None

local_noon_offset = local_solar_noon_offset(dt.datetime(2021, 11, 17, 19, 0, 0), -101.5)
print(local_noon_offset)

which gives 0:28:57.123186.

Naively calculating the offset for 360 longitudes takes 12 seconds:

%%timeit 
lons = np.arange(-180, 180, 1)
offsets = [local_solar_noon_offset(dt.datetime(2021, 11, 17, 19, 0, 0), lon)
               for lon in lons]