pvlib
Function to output exact time, in seconds or microseconds or nanoseconds if possible, of sun reaching specific elevation.
Pvlib team, thank you for your amazing work.
First of all, I'm not a developer, I code as a hobby.
Regarding my request, I wanted to implement a function to calculate the time when the sun reaches a certain elevation, I used iterations to do so, but it takes a lot of time to output the time, and even then, the time it outputs is not accurate, it is close to the exact time, but still.
It would be amazing if there were a function that could do this directly without taking up so much time.
My code:
import datetime
import pvlib
import pytz
import requests
response = requests.get('https://ipinfo.io')
data = response.json()
lat, long = map(float, data['loc'].split(','))
timezone_str = data['timezone']
timezone_object = pytz.timezone(timezone_str)
current_time = datetime.datetime.now(timezone_object)
current_time = current_time.replace(hour=0, minute=0, second=0, microsecond=0)
angle = int(input('angle: '))
for h in range(24):
for m in range(60):
current_time = current_time.replace(hour=h, minute=m)
solpos = pvlib.solarposition.get_solarposition(current_time, lat, long)
elev = solpos.elevation.values[0]
if abs(angle-elev) < 0.1:
print(elev, '\t', current_time)
A solar position algorithm does this: solar azimuth, solar elevation = function (datetime given (longitude, latitude))
You want the inverse: datetime = function (elevation given (longitude, latitude))
This thread concludes that the solar azimuth would need to be specified somehow. If it is assumed that the ordinal date is known, then it seems possible to arrive at an inverse that would give both times in that date, since the ordinal date corresponds to a specific sun path from the earth's viewpoint.
I think you will get significantly faster results using SciPy Optimize instead of brute force. For example, minimize_scalar. I think you have the following challenges to consider:
- You need to convert from datetimes to something continuously variable like seconds since the common epoch (aka UNIX time).
- Everyday will have two different times when the elevation is reached, once in the morning and again in the evening.
- There are days when the desired elevation may never be reached.
Here's some code that covers the entire year in about 1 minute:
from calendar import timegm
from datetime import datetime, timedelta
from time import gmtime
import pandas as pd
import pytz
from pvlib.solarposition import get_solarposition
import numpy as np
from scipy.optimize import minimize_scalar, root_scalar
import logging
logging.basicConfig()
LOGGER = logging.getLogger()
LOGGER.setLevel(logging.DEBUG)
def sp_from_epoch(epochtime, lat, lon, tz):
time = pd.Timestamp(*gmtime(epochtime)[:6], tz=tz)
times = pd.DatetimeIndex([time]) # pvlib expects an array
return get_solarposition(times, lat, lon)
def objective(epochtime, lat, lon, tz, elev):
df = sp_from_epoch(epochtime, lat, lon, tz)
#LOGGER.debug(df)
residual = df.apparent_elevation.values[0] - elev
return residual*residual
# CONSTANTS
SP_NAT = {
'apparent_zenith': np.nan,
'zenith': np.nan,
'apparent_elevation': np.nan,
'elevation': np.nan,
'azimuth': np.nan,
'equation_of_time': np.nan}
SP_NAT = pd.DataFrame(SP_NAT, index=pd.DatetimeIndex(['NaT']))
NDAYS = 365
TOL = 1.0 # return NaN if resnorm greater than TOL
def get_elevation_time(day_of_year, latitude, longitude, timezone, elevation, year):
"""
Get time when solar positon is at specified elevation.
Parameters
----------
day_of_year : int
day of year starting at zero for JAN1
latitude : float
latitude in degrees, measured from the equator, positive is north
longitude : float
longitude in degrees, measured from prime meridian, positive is east
timezone : pytz.timezone
Hint: Use 'Etc/GMT+X' format where -X is hours from GMT
elevation : float
specify desired elevation
year : int
IDK what will happen during leap year
Returns
-------
Dictionary of dataframes of solar positions corresponding to desired apparent
elevation for morning and evening.
"""
args = (latitude, longitude, timezone, elevation)
delta = datetime(year, 1, 1) + timedelta(days=day_of_year)
# morning
a = timegm(delta.timetuple())
b = timegm((delta + timedelta(hours=12)).timetuple())
#LOGGER.debug(a)
#LOGGER.debug(b)
res = minimize_scalar(objective, args=args, bounds=(a, b))
if res.fun < TOL:
morning = sp_from_epoch(res.x, *args[:3])
else:
morning = SP_NAT
# evening
c = timegm((delta + timedelta(hours=24)).timetuple())
res = minimize_scalar(objective, args=args, bounds=(b, c))
if res.fun < TOL:
evening = sp_from_epoch(res.x, *args[:3])
else:
evening = SP_NAT
return {'morning': morning, 'evening': evening}
def get_elevation_times(latitude, longitude, timezone, elevation, year):
"""
Get times for a full year when solar positon is at specified elevation.
Parameters
----------
latitude : float
latitude in degrees, measured from the equator, positive is north
longitude : float
longitude in degrees, measured from prime meridian, positive is east
timezone : pytz.timezone
Hint: Use 'Etc/GMT+X' format where -X is hours from GMT
elevation : float
specify desired elevation
year : int
IDK what will happen during leap year
Returns
-------
Dictionary of dataframes of solar positions corresponding to desired apparent
elevation for mornings and evenings.
"""
mornings = []
evenings = []
args = (latitude, longitude, timezone, elevation, year)
for d in range(NDAYS):
res = get_elevation_time(d, *args)
mornings.append(res['morning'])
evenings.append(res['evening'])
mornings = pd.concat(mornings)
evenings = pd.concat(evenings)
return {'mornings': mornings, 'evenings': evenings}
Here's a test example:
# test drive
LAT = 38.5
LON = -122.0
TZ = pytz.timezone('Etc/GMT+8')
ELEV = 15 # arcdegrees
YEAR = 2023
res = get_elevation_times(LAT, LON, TZ, ELEV, YEAR)
res['mornings']
res['evenings']
The test output for mornings:
and for evenings:
I've attached the jupyter notebook and html version in a zip file pvlib_gh1981.zip using the following packages:
- pvlib: 0.9.2
- pandas: 2.0.3
- numpy: 1.24.4
- scipy: 1.11.1
- python: 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]
I think you can increase the speed even more if you vectorize, perhaps using root_scalar and the Newton method, but IDK.
Note that pvlib has a function that uses scipy.optimize.brentq with pyephem: pvlib.solarposition.calc_time.
