Open dopplershift opened 3 years ago
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.
@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.
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 :grin:)
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.
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.
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.
+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]
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: