acgc.solar

Module to calculate solar zenith angle, solar declination, and equation of time Results should be accurate to < 0.1 degree, but other modules should be used for high-precision calculations.

C.D. Holmes 9 Nov 2018

View Source

  1#!/usr/local/bin/env python3
  2''' Module to calculate solar zenith angle, solar declination, and equation of time 
  3Results should be accurate to < 0.1 degree, but other modules should be used for 
  4high-precision calculations.
  5
  6C.D. Holmes 9 Nov 2018
  7'''
  8
  9import warnings
 10import numpy as np
 11import pandas as pd
 12
 13pi180 = np.pi / 180
 14
 15def solar_azimuth_angle( lat, lon, datetimeUTC ):
 16    '''Solar azimuth angle (degrees) for a latitude, longitude, date and time
 17    
 18    SAA is degrees clockwise from north.
 19    
 20    Parameters
 21    ----------
 22    lat : float or ndarray
 23        latitude in degrees
 24    lon : float or ndarray
 25        longitudes in degrees
 26    datetimeUTC : pandas.Timestamp, datetime, or str
 27        date and time in UTC
 28
 29    Returns
 30    -------
 31    saa : float or ndarray
 32        solar azimuth angle in degrees (clockwise from north)
 33    '''
 34    # Convert to pandas Timestamp, if needed
 35    datetimeUTC = _to_timestamp(datetimeUTC)
 36
 37    # Solar declination, degrees
 38    dec = solar_declination( datetimeUTC )
 39
 40    # Hour angle, degrees
 41    Ha = solar_hour_angle( lon, datetimeUTC )
 42
 43    # Solar zenith angle, degrees
 44    # Use true sza, without refraction
 45    zen = sza( lat, lon, datetimeUTC, refraction=False )
 46
 47    # Solar azimuth angle, degrees
 48    saa = np.arcsin( -np.sin( Ha*pi180 ) * np.cos( dec*pi180 ) /
 49            np.sin( zen*pi180 ) ) / pi180
 50
 51    # Change range [-180,180] to [0,360]
 52    return np.mod( saa+360, 360 )
 53
 54def solar_elevation_angle( lat, lon, alt, datetimeUTC,
 55                       refraction=False, temperature=10., pressure=101325. ):
 56    '''Solar elevation angle (degrees) above the horizon
 57
 58    The altitude parameter should be the vertical distance 
 59    above the surrounding terrain that defines the horizon,
 60    not necessarily the altitude above sea level or the altitude above ground level.
 61    For example, on a mountain peak that is 4000 m above sea level and 
 62    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
 63    For an observer on the plateau, the relevant altitude is 0 m.
 64
 65    See documentation for `solar_zenith_angle` and `horizon_zenith_angle`.
 66
 67    Parameters
 68    ----------
 69    lat : float or ndarray
 70        latitude in degrees
 71    lon : float or ndarray
 72        longitudes in degrees
 73    alt : float or ndarray
 74        altitude above surrounding terrain that defines the horizon, meters
 75    datetimeUTC : pandas.Timestamp, datetime, or str
 76        date and time in UTC
 77    refraction : bool, optional (default=False)
 78        specifies whether to account for atmospheric refraction
 79    temperature : float or ndarray, optional (default=10)
 80        surface atmospheric temperature (Celsius), only used for refraction calculation
 81    pressure : float or ndarray, optional (default=101325)
 82        surface atmospheric pressure (Pa), only used for refraction calculation
 83    
 84    Returns
 85    -------
 86    sea : float or ndarray
 87        solar elevation angle in degrees at the designated locations and times
 88        If refraction=False, this is the true solar elevation angle
 89        If refraction=True, this is the apparent solar elevation angle
 90    
 91    '''
 92
 93    if refraction and np.any(alt):
 94        warnings.warn( 'Atmospheric refraction is calculated for surface conditions, '
 95                    + 'but an altitude above the surface was specified',
 96                     category=UserWarning,
 97                     stacklevel=2 )
 98
 99    sea = horizon_zenith_angle( alt ) \
100         - solar_zenith_angle( lat, lon, datetimeUTC, refraction, temperature, pressure )
101
102    return sea
103
104def solar_zenith_angle( lat, lon, datetimeUTC, 
105                       refraction=False, temperature=10., pressure=101325. ):
106    '''Solar zenith angle (degrees) for a given latitude, longitude, date and time.
107    
108    Accounts for equation of time and (optionally) for atmospheric refraction.
109    Altitude of the observer is not accounted for, which can be important when the sun 
110    is near the horizon. 
111    
112    Results are accurate to tenths of a degree, except where altitude is important
113    (< 20 degrees solar elevation)
114
115    Parameters
116    ----------
117    lat : float or ndarray
118        latitude in degrees
119    lon : float or ndarray
120        longitudes in degrees
121    datetimeUTC : pandas.Timestamp, datetime, or str
122        date and time in UTC
123    refraction : bool, optional (default=False)
124        specifies whether to account for atmospheric refraction
125    temperature : float or ndarray, optional (default=10)
126        surface atmospheric temperature (Celsius), only used for refraction calculation
127    pressure : float or ndarray, optional (default=101325)
128        surface atmospheric pressure (Pa), only used for refraction calculation
129    
130    Returns
131    -------
132    sza : float or ndarray
133        solar zenith angle in degrees at the designated locations and times
134        If refraction=False, this is the true solar zenith angle
135        If refraction=True, this is the apparent solar zenith angle
136    '''
137    # Convert to pandas Timestamp, if needed
138    datetimeUTC = _to_timestamp(datetimeUTC)
139
140    # Solar declination, degrees
141    dec = solar_declination( datetimeUTC )
142
143    # Hour angle, degrees
144    Ha = solar_hour_angle( lon, datetimeUTC )
145
146    # True solar zenith angle, radians
147    sza = np.arccos( np.sin(lat*pi180) * np.sin(dec*pi180) + \
148          np.cos(lat*pi180) * np.cos(dec*pi180) * np.cos(Ha*pi180) )
149
150    # Convert radians -> degrees
151    sza /= pi180
152
153    if refraction:
154        # Subtract refraction angle (degrees) from zenith angle.
155        # SZA is always smaller due to refraction.
156        sza -= refraction_angle( 90-sza, pressure, temperature )
157
158    return sza
159
160def equation_of_time( date ):
161    '''Equation of time (degrees) for specified date
162    
163    Implements the "alternative equation" from Wikipedia, derived from
164    https://web.archive.org/web/20120323231813/http://www.green-life-innovators.org/tiki-index.php?page=The%2BLatitude%2Band%2BLongitude%2Bof%2Bthe%2BSun%2Bby%2BDavid%2BWilliams
165    Results checked against NOAA solar calculator and agree within 10 seconds.
166    
167    Note: Leap years are not accounted for.
168
169    Argument
170    --------
171    date : pandas.Timestamp, date, or datetime
172        date for calculation
173
174    Returns
175    -------
176    eot : float
177        equation of time in degrees on the specified date
178    '''
179    # Convert to pandas Timestamp, if needed
180    date = _to_timestamp(date)
181
182    # Equation of time, accounts for the solar day differing slightly from 24 hr
183    try:
184        doy = date.dt.dayofyear
185    except AttributeError:
186        doy = date.dayofyear
187    W = 360 / 365.24
188    A = W * (doy+10)
189    B = A + 1.914 * np.sin( W * (doy-2) * pi180 )
190    C = ( A - np.arctan2( np.tan(B*pi180), np.cos(23.44*pi180) ) / pi180 ) / 180
191
192    # Equation of time in minutes of an hour
193    eotmin = 720 * ( C - np.round(C) )
194
195    # Equation of time, minutes -> degrees
196    eot = eotmin / 60 * 360 / 24
197
198    return eot
199
200def horizon_zenith_angle( alt ):
201    '''Angle from the zenith to the horizon
202    
203    The horizon is the locii of points where a line from the 
204    observation location to the ellipsoid is tangent to the ellipsoid surface.
205    
206    The altitude parameter should be the vertical distance 
207    above the surrounding terrain that defines the horizon,
208    not necessarily the altitude above sea level or the altitude above ground level.
209    For example, on a mountain peak that is 4000 m above sea level and 
210    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
211    For an observer on the plateau, the relevant altitude is 0 m.
212
213    The implementation below assumes a spherical Earth.
214    Results using the WGS84 ellipsoid (see commented code below)
215    differ from the spherical case by << 1°. Terrain,
216    which is neglected here, has a larger effect on the horizon
217    location, so the simpler spherical calculation is appropriate. 
218
219    Parameters
220    ----------
221    lat : float or ndarray
222        latitude in degrees
223    alt : float or ndarray
224        altitude above surrounding terrain that defines the horizon, meters
225        
226    Returns
227    -------
228    hza : float or ndarray
229        horizon zenith angle in degrees
230    '''
231
232    # WGS84 ellipsoid parameters
233    # semi-major radius, m
234    r_earth = 6378137.0
235    # ellipsoidal flattening, unitless
236    f = 1/298.257223563
237
238    # Horizon zenith angle, degrees (spherical earth)
239    hza = 180 - np.arcsin( r_earth / ( r_earth + alt ) ) / pi180
240
241    ## Ellipsoidal Earth
242    # # Eccentricity of ellipsoid
243    # ecc = f * (2-f)
244    # # Local (i.e. prime vertical) radius of curvature at latitude
245    # N = r_earth / np.sqrt( 1 - ecc**2 * np.sin(lat*pi180)**2 )
246    # # Horizon zenith angle, degrees
247    # hza = 180 - np.arcsin( N / (N+alt) ) / pi180
248
249    return hza
250
251def solar_declination( date ):
252    '''Calculate solar declination (degrees) for specified date
253    
254    Implements Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling
255    
256    Argument
257    --------
258    date : pandas.Timestamp, date, datetime, or str
259        date for calculation
260
261    Returns
262    -------
263    dec : float
264        solar declination in degrees at the specified date
265    '''
266
267    # Convert to pandas Timestamp, if needed
268    date = _to_timestamp(date)
269
270     # Number of days since beginning of 2000
271    NJD = ( date - np.datetime64('2000-01-01') )
272    try:
273        NJD = NJD.dt.days
274    except AttributeError:
275        NJD = NJD.days
276
277    # Obliquity, degrees
278    ob = 23.439 - 4e-7 * NJD
279
280    # Parameters for ecliptic, degrees
281    gm = 357.528 + 0.9856003 * NJD
282    lm = 280.460 + 0.9856474 * NJD
283
284    # Ecliptic longitude of sun, degrees
285    ec = lm + 1.915 * np.sin( gm * pi180 ) + 0.020 * np.sin( 2 * gm * pi180 )
286
287    #Solar declination, degrees
288    dec = np.arcsin( np.sin( ob * pi180 ) * np.sin( ec * pi180 ) ) / pi180
289
290    return dec
291
292def solar_hour_angle( lon, datetimeUTC ):
293    '''Solar hour angle (degrees) for specified longitude, date and time
294
295    Hour angle is the angular displacement of the sun from the local meridian.
296    It is zero at local noon, negative in the morning, and positive is afternoon.
297    
298    Parameters
299    ----------
300    lon : float
301        longitude in degrees east
302    datetimeUTC : pandas.Timestamp or datetime
303        date and time for calculation, must be UTC
304    
305    Returns
306    -------
307    ha : float
308        hour angle in degrees at the specified location and time
309    '''
310
311    # Convert to pandas Timestamp, if needed
312    datetimeUTC = _to_timestamp(datetimeUTC)
313
314    # Hour angle for mean solar time.
315    # Actual solar position has a small offset given by the equation of time (below)
316    try: 
317        # Treat as xarray.DataArray or pandas.Series
318        Ha = lon + 15 * ( datetimeUTC.dt.hour + 
319                          datetimeUTC.dt.minute / 60 + 
320                          datetimeUTC.dt.second / 3600 - 12 )
321    except AttributeError:
322        Ha = lon + 15 * ( datetimeUTC.hour + 
323                          datetimeUTC.minute / 60 + 
324                          datetimeUTC.second / 3600 - 12 )
325
326    # Add equation of time to the hour angle, degrees
327    Ha += equation_of_time( datetimeUTC )
328
329    return Ha
330
331def refraction_angle( true_elevation_angle, pressure=101325., temperature_celsius=10. ):
332    '''Atmospheric refraction angle for light passing through Earth's atmosphere
333
334    The apparent locations in the sky of objects outsides Earth's atmosphere 
335    differs from their true locations due to atmospheric refraction. 
336    (e.g. sun and moon when they rise and set)
337    The apparent elevation of an object above the horizon is
338    apparent elevation angle = (true elevation angle) + (refraction angle)
339    
340    The equations here are from Saemundsson/Bennett, whose calculations use
341    a typical vertical profile of atmospheric density (i.e. temperature and pressure).
342    The profiles can be rescaled to a particular surface temperature and pressure
343    to approximately account for varying atmospheric conditions.
344    Accurate refraction calculations should use fully specified vertical profile
345    of temperature and pressure, which cannot be done here.
346
347    Parameters
348    ----------
349    true_elevation_angle : float
350        degrees above horizon of sun or other object
351    pressure : float (default=101325)
352        surface atmospheric pressure (Pa)
353    temperature_celsius : float (default=10)
354        surface atmospheric temperature (C)
355
356    Returns
357    -------
358    angle : float
359        refraction angle in degrees. Value is zero when apparent elevation is below horizon
360    '''
361    # Refraction angle, arcminutes
362    R = 1.02 / np.tan( ( true_elevation_angle + 10.3 / (true_elevation_angle + 5.11) ) * pi180 )
363    # Account for temperature and pressure, arcminutes
364    R = R * pressure / 101325 * 283 / ( 273 + temperature_celsius )
365    # Convert arcminutes -> degrees
366    R /= 60
367
368    # Result must be positive
369    R = np.maximum(R,0)
370
371    # Refraction defined only when the apparent elevation angle is positive
372    # Set refraction to zero when the apparent elevation is below horizon
373    refraction_angle = np.where( true_elevation_angle + R <= 0, 0, R)
374
375    return refraction_angle
376
377def _to_timestamp(time_in):
378    '''Convert input to Pandas Timestamp or Series of datetime64 
379    
380    Arguments
381    ---------
382    time_in : datetime-like or array
383        time to be converted
384
385    Returns
386    -------
387    time_out : pandas.Timestamp or pandas.Series of datetime64
388    '''
389    if hasattr(time_in,'dt'):
390        time_out = time_in
391    elif isinstance(time_in, pd.DatetimeIndex ):
392        time_out = pd.Series(time_in)
393    else:
394        try:
395            # Convert list of times
396            time_out = pd.Series(pd.DatetimeIndex(time_in))
397        except TypeError:
398            time_out = pd.Timestamp(time_in)
399
400    return time_out
401
402# Aliases for functions
403sza = solar_zenith_angle
404saa = solar_azimuth_angle
405sea = solar_elevation_angle
406# Additional aliases for backwards compatibility
407equationOfTime = equation_of_time
408solarDeclination = solar_declination

pi180 = 0.017453292519943295

def solar_azimuth_angle(lat, lon, datetimeUTC): View Source

16def solar_azimuth_angle( lat, lon, datetimeUTC ):
17    '''Solar azimuth angle (degrees) for a latitude, longitude, date and time
18    
19    SAA is degrees clockwise from north.
20    
21    Parameters
22    ----------
23    lat : float or ndarray
24        latitude in degrees
25    lon : float or ndarray
26        longitudes in degrees
27    datetimeUTC : pandas.Timestamp, datetime, or str
28        date and time in UTC
29
30    Returns
31    -------
32    saa : float or ndarray
33        solar azimuth angle in degrees (clockwise from north)
34    '''
35    # Convert to pandas Timestamp, if needed
36    datetimeUTC = _to_timestamp(datetimeUTC)
37
38    # Solar declination, degrees
39    dec = solar_declination( datetimeUTC )
40
41    # Hour angle, degrees
42    Ha = solar_hour_angle( lon, datetimeUTC )
43
44    # Solar zenith angle, degrees
45    # Use true sza, without refraction
46    zen = sza( lat, lon, datetimeUTC, refraction=False )
47
48    # Solar azimuth angle, degrees
49    saa = np.arcsin( -np.sin( Ha*pi180 ) * np.cos( dec*pi180 ) /
50            np.sin( zen*pi180 ) ) / pi180
51
52    # Change range [-180,180] to [0,360]
53    return np.mod( saa+360, 360 )

Solar azimuth angle (degrees) for a latitude, longitude, date and time

SAA is degrees clockwise from north.

Parameters

lat (float or ndarray): latitude in degrees
lon (float or ndarray): longitudes in degrees
datetimeUTC (pandas.Timestamp, datetime, or str): date and time in UTC

Returns

saa (float or ndarray): solar azimuth angle in degrees (clockwise from north)

def solar_elevation_angle( lat, lon, alt, datetimeUTC, refraction=False, temperature=10.0, pressure=101325.0): View Source

 55def solar_elevation_angle( lat, lon, alt, datetimeUTC,
 56                       refraction=False, temperature=10., pressure=101325. ):
 57    '''Solar elevation angle (degrees) above the horizon
 58
 59    The altitude parameter should be the vertical distance 
 60    above the surrounding terrain that defines the horizon,
 61    not necessarily the altitude above sea level or the altitude above ground level.
 62    For example, on a mountain peak that is 4000 m above sea level and 
 63    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
 64    For an observer on the plateau, the relevant altitude is 0 m.
 65
 66    See documentation for `solar_zenith_angle` and `horizon_zenith_angle`.
 67
 68    Parameters
 69    ----------
 70    lat : float or ndarray
 71        latitude in degrees
 72    lon : float or ndarray
 73        longitudes in degrees
 74    alt : float or ndarray
 75        altitude above surrounding terrain that defines the horizon, meters
 76    datetimeUTC : pandas.Timestamp, datetime, or str
 77        date and time in UTC
 78    refraction : bool, optional (default=False)
 79        specifies whether to account for atmospheric refraction
 80    temperature : float or ndarray, optional (default=10)
 81        surface atmospheric temperature (Celsius), only used for refraction calculation
 82    pressure : float or ndarray, optional (default=101325)
 83        surface atmospheric pressure (Pa), only used for refraction calculation
 84    
 85    Returns
 86    -------
 87    sea : float or ndarray
 88        solar elevation angle in degrees at the designated locations and times
 89        If refraction=False, this is the true solar elevation angle
 90        If refraction=True, this is the apparent solar elevation angle
 91    
 92    '''
 93
 94    if refraction and np.any(alt):
 95        warnings.warn( 'Atmospheric refraction is calculated for surface conditions, '
 96                    + 'but an altitude above the surface was specified',
 97                     category=UserWarning,
 98                     stacklevel=2 )
 99
100    sea = horizon_zenith_angle( alt ) \
101         - solar_zenith_angle( lat, lon, datetimeUTC, refraction, temperature, pressure )
102
103    return sea

Solar elevation angle (degrees) above the horizon

The altitude parameter should be the vertical distance above the surrounding terrain that defines the horizon, not necessarily the altitude above sea level or the altitude above ground level. For example, on a mountain peak that is 4000 m above sea level and 1500 m above the surrounding plateau, the relevant altitude is 1500 m. For an observer on the plateau, the relevant altitude is 0 m.

See documentation for solar_zenith_angle and horizon_zenith_angle.

Parameters

lat (float or ndarray): latitude in degrees
lon (float or ndarray): longitudes in degrees
alt (float or ndarray): altitude above surrounding terrain that defines the horizon, meters
datetimeUTC (pandas.Timestamp, datetime, or str): date and time in UTC
refraction (bool, optional (default=False)): specifies whether to account for atmospheric refraction
temperature (float or ndarray, optional (default=10)): surface atmospheric temperature (Celsius), only used for refraction calculation
pressure (float or ndarray, optional (default=101325)): surface atmospheric pressure (Pa), only used for refraction calculation

Returns

sea (float or ndarray): solar elevation angle in degrees at the designated locations and times If refraction=False, this is the true solar elevation angle If refraction=True, this is the apparent solar elevation angle

def solar_zenith_angle( lat, lon, datetimeUTC, refraction=False, temperature=10.0, pressure=101325.0): View Source

105def solar_zenith_angle( lat, lon, datetimeUTC, 
106                       refraction=False, temperature=10., pressure=101325. ):
107    '''Solar zenith angle (degrees) for a given latitude, longitude, date and time.
108    
109    Accounts for equation of time and (optionally) for atmospheric refraction.
110    Altitude of the observer is not accounted for, which can be important when the sun 
111    is near the horizon. 
112    
113    Results are accurate to tenths of a degree, except where altitude is important
114    (< 20 degrees solar elevation)
115
116    Parameters
117    ----------
118    lat : float or ndarray
119        latitude in degrees
120    lon : float or ndarray
121        longitudes in degrees
122    datetimeUTC : pandas.Timestamp, datetime, or str
123        date and time in UTC
124    refraction : bool, optional (default=False)
125        specifies whether to account for atmospheric refraction
126    temperature : float or ndarray, optional (default=10)
127        surface atmospheric temperature (Celsius), only used for refraction calculation
128    pressure : float or ndarray, optional (default=101325)
129        surface atmospheric pressure (Pa), only used for refraction calculation
130    
131    Returns
132    -------
133    sza : float or ndarray
134        solar zenith angle in degrees at the designated locations and times
135        If refraction=False, this is the true solar zenith angle
136        If refraction=True, this is the apparent solar zenith angle
137    '''
138    # Convert to pandas Timestamp, if needed
139    datetimeUTC = _to_timestamp(datetimeUTC)
140
141    # Solar declination, degrees
142    dec = solar_declination( datetimeUTC )
143
144    # Hour angle, degrees
145    Ha = solar_hour_angle( lon, datetimeUTC )
146
147    # True solar zenith angle, radians
148    sza = np.arccos( np.sin(lat*pi180) * np.sin(dec*pi180) + \
149          np.cos(lat*pi180) * np.cos(dec*pi180) * np.cos(Ha*pi180) )
150
151    # Convert radians -> degrees
152    sza /= pi180
153
154    if refraction:
155        # Subtract refraction angle (degrees) from zenith angle.
156        # SZA is always smaller due to refraction.
157        sza -= refraction_angle( 90-sza, pressure, temperature )
158
159    return sza

Solar zenith angle (degrees) for a given latitude, longitude, date and time.

Accounts for equation of time and (optionally) for atmospheric refraction. Altitude of the observer is not accounted for, which can be important when the sun is near the horizon.

Results are accurate to tenths of a degree, except where altitude is important (< 20 degrees solar elevation)

Parameters

lat (float or ndarray): latitude in degrees
lon (float or ndarray): longitudes in degrees
datetimeUTC (pandas.Timestamp, datetime, or str): date and time in UTC
refraction (bool, optional (default=False)): specifies whether to account for atmospheric refraction
temperature (float or ndarray, optional (default=10)): surface atmospheric temperature (Celsius), only used for refraction calculation
pressure (float or ndarray, optional (default=101325)): surface atmospheric pressure (Pa), only used for refraction calculation

Returns

sza (float or ndarray): solar zenith angle in degrees at the designated locations and times If refraction=False, this is the true solar zenith angle If refraction=True, this is the apparent solar zenith angle

def equation_of_time(date): View Source

161def equation_of_time( date ):
162    '''Equation of time (degrees) for specified date
163    
164    Implements the "alternative equation" from Wikipedia, derived from
165    https://web.archive.org/web/20120323231813/http://www.green-life-innovators.org/tiki-index.php?page=The%2BLatitude%2Band%2BLongitude%2Bof%2Bthe%2BSun%2Bby%2BDavid%2BWilliams
166    Results checked against NOAA solar calculator and agree within 10 seconds.
167    
168    Note: Leap years are not accounted for.
169
170    Argument
171    --------
172    date : pandas.Timestamp, date, or datetime
173        date for calculation
174
175    Returns
176    -------
177    eot : float
178        equation of time in degrees on the specified date
179    '''
180    # Convert to pandas Timestamp, if needed
181    date = _to_timestamp(date)
182
183    # Equation of time, accounts for the solar day differing slightly from 24 hr
184    try:
185        doy = date.dt.dayofyear
186    except AttributeError:
187        doy = date.dayofyear
188    W = 360 / 365.24
189    A = W * (doy+10)
190    B = A + 1.914 * np.sin( W * (doy-2) * pi180 )
191    C = ( A - np.arctan2( np.tan(B*pi180), np.cos(23.44*pi180) ) / pi180 ) / 180
192
193    # Equation of time in minutes of an hour
194    eotmin = 720 * ( C - np.round(C) )
195
196    # Equation of time, minutes -> degrees
197    eot = eotmin / 60 * 360 / 24
198
199    return eot

Equation of time (degrees) for specified date

Implements the "alternative equation" from Wikipedia, derived from https://web.archive.org/web/20120323231813/http://www.green-life-innovators.org/tiki-index.php?page=The%2BLatitude%2Band%2BLongitude%2Bof%2Bthe%2BSun%2Bby%2BDavid%2BWilliams Results checked against NOAA solar calculator and agree within 10 seconds.

Note: Leap years are not accounted for.

Argument

date : pandas.Timestamp, date, or datetime date for calculation

Returns

eot (float): equation of time in degrees on the specified date

def horizon_zenith_angle(alt): View Source

201def horizon_zenith_angle( alt ):
202    '''Angle from the zenith to the horizon
203    
204    The horizon is the locii of points where a line from the 
205    observation location to the ellipsoid is tangent to the ellipsoid surface.
206    
207    The altitude parameter should be the vertical distance 
208    above the surrounding terrain that defines the horizon,
209    not necessarily the altitude above sea level or the altitude above ground level.
210    For example, on a mountain peak that is 4000 m above sea level and 
211    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
212    For an observer on the plateau, the relevant altitude is 0 m.
213
214    The implementation below assumes a spherical Earth.
215    Results using the WGS84 ellipsoid (see commented code below)
216    differ from the spherical case by << 1°. Terrain,
217    which is neglected here, has a larger effect on the horizon
218    location, so the simpler spherical calculation is appropriate. 
219
220    Parameters
221    ----------
222    lat : float or ndarray
223        latitude in degrees
224    alt : float or ndarray
225        altitude above surrounding terrain that defines the horizon, meters
226        
227    Returns
228    -------
229    hza : float or ndarray
230        horizon zenith angle in degrees
231    '''
232
233    # WGS84 ellipsoid parameters
234    # semi-major radius, m
235    r_earth = 6378137.0
236    # ellipsoidal flattening, unitless
237    f = 1/298.257223563
238
239    # Horizon zenith angle, degrees (spherical earth)
240    hza = 180 - np.arcsin( r_earth / ( r_earth + alt ) ) / pi180
241
242    ## Ellipsoidal Earth
243    # # Eccentricity of ellipsoid
244    # ecc = f * (2-f)
245    # # Local (i.e. prime vertical) radius of curvature at latitude
246    # N = r_earth / np.sqrt( 1 - ecc**2 * np.sin(lat*pi180)**2 )
247    # # Horizon zenith angle, degrees
248    # hza = 180 - np.arcsin( N / (N+alt) ) / pi180
249
250    return hza

Angle from the zenith to the horizon

The horizon is the locii of points where a line from the observation location to the ellipsoid is tangent to the ellipsoid surface.

The implementation below assumes a spherical Earth. Results using the WGS84 ellipsoid (see commented code below) differ from the spherical case by << 1°. Terrain, which is neglected here, has a larger effect on the horizon location, so the simpler spherical calculation is appropriate.

Parameters

lat (float or ndarray): latitude in degrees
alt (float or ndarray): altitude above surrounding terrain that defines the horizon, meters

Returns

hza (float or ndarray): horizon zenith angle in degrees

def solar_declination(date): View Source

252def solar_declination( date ):
253    '''Calculate solar declination (degrees) for specified date
254    
255    Implements Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling
256    
257    Argument
258    --------
259    date : pandas.Timestamp, date, datetime, or str
260        date for calculation
261
262    Returns
263    -------
264    dec : float
265        solar declination in degrees at the specified date
266    '''
267
268    # Convert to pandas Timestamp, if needed
269    date = _to_timestamp(date)
270
271     # Number of days since beginning of 2000
272    NJD = ( date - np.datetime64('2000-01-01') )
273    try:
274        NJD = NJD.dt.days
275    except AttributeError:
276        NJD = NJD.days
277
278    # Obliquity, degrees
279    ob = 23.439 - 4e-7 * NJD
280
281    # Parameters for ecliptic, degrees
282    gm = 357.528 + 0.9856003 * NJD
283    lm = 280.460 + 0.9856474 * NJD
284
285    # Ecliptic longitude of sun, degrees
286    ec = lm + 1.915 * np.sin( gm * pi180 ) + 0.020 * np.sin( 2 * gm * pi180 )
287
288    #Solar declination, degrees
289    dec = np.arcsin( np.sin( ob * pi180 ) * np.sin( ec * pi180 ) ) / pi180
290
291    return dec

Calculate solar declination (degrees) for specified date

Implements Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling

Argument

date : pandas.Timestamp, date, datetime, or str date for calculation

Returns

dec (float): solar declination in degrees at the specified date

def solar_hour_angle(lon, datetimeUTC): View Source

293def solar_hour_angle( lon, datetimeUTC ):
294    '''Solar hour angle (degrees) for specified longitude, date and time
295
296    Hour angle is the angular displacement of the sun from the local meridian.
297    It is zero at local noon, negative in the morning, and positive is afternoon.
298    
299    Parameters
300    ----------
301    lon : float
302        longitude in degrees east
303    datetimeUTC : pandas.Timestamp or datetime
304        date and time for calculation, must be UTC
305    
306    Returns
307    -------
308    ha : float
309        hour angle in degrees at the specified location and time
310    '''
311
312    # Convert to pandas Timestamp, if needed
313    datetimeUTC = _to_timestamp(datetimeUTC)
314
315    # Hour angle for mean solar time.
316    # Actual solar position has a small offset given by the equation of time (below)
317    try: 
318        # Treat as xarray.DataArray or pandas.Series
319        Ha = lon + 15 * ( datetimeUTC.dt.hour + 
320                          datetimeUTC.dt.minute / 60 + 
321                          datetimeUTC.dt.second / 3600 - 12 )
322    except AttributeError:
323        Ha = lon + 15 * ( datetimeUTC.hour + 
324                          datetimeUTC.minute / 60 + 
325                          datetimeUTC.second / 3600 - 12 )
326
327    # Add equation of time to the hour angle, degrees
328    Ha += equation_of_time( datetimeUTC )
329
330    return Ha

Solar hour angle (degrees) for specified longitude, date and time

Hour angle is the angular displacement of the sun from the local meridian. It is zero at local noon, negative in the morning, and positive is afternoon.

Parameters

lon (float): longitude in degrees east
datetimeUTC (pandas.Timestamp or datetime): date and time for calculation, must be UTC

Returns

ha (float): hour angle in degrees at the specified location and time

def refraction_angle(true_elevation_angle, pressure=101325.0, temperature_celsius=10.0): View Source

332def refraction_angle( true_elevation_angle, pressure=101325., temperature_celsius=10. ):
333    '''Atmospheric refraction angle for light passing through Earth's atmosphere
334
335    The apparent locations in the sky of objects outsides Earth's atmosphere 
336    differs from their true locations due to atmospheric refraction. 
337    (e.g. sun and moon when they rise and set)
338    The apparent elevation of an object above the horizon is
339    apparent elevation angle = (true elevation angle) + (refraction angle)
340    
341    The equations here are from Saemundsson/Bennett, whose calculations use
342    a typical vertical profile of atmospheric density (i.e. temperature and pressure).
343    The profiles can be rescaled to a particular surface temperature and pressure
344    to approximately account for varying atmospheric conditions.
345    Accurate refraction calculations should use fully specified vertical profile
346    of temperature and pressure, which cannot be done here.
347
348    Parameters
349    ----------
350    true_elevation_angle : float
351        degrees above horizon of sun or other object
352    pressure : float (default=101325)
353        surface atmospheric pressure (Pa)
354    temperature_celsius : float (default=10)
355        surface atmospheric temperature (C)
356
357    Returns
358    -------
359    angle : float
360        refraction angle in degrees. Value is zero when apparent elevation is below horizon
361    '''
362    # Refraction angle, arcminutes
363    R = 1.02 / np.tan( ( true_elevation_angle + 10.3 / (true_elevation_angle + 5.11) ) * pi180 )
364    # Account for temperature and pressure, arcminutes
365    R = R * pressure / 101325 * 283 / ( 273 + temperature_celsius )
366    # Convert arcminutes -> degrees
367    R /= 60
368
369    # Result must be positive
370    R = np.maximum(R,0)
371
372    # Refraction defined only when the apparent elevation angle is positive
373    # Set refraction to zero when the apparent elevation is below horizon
374    refraction_angle = np.where( true_elevation_angle + R <= 0, 0, R)
375
376    return refraction_angle

Atmospheric refraction angle for light passing through Earth's atmosphere

The apparent locations in the sky of objects outsides Earth's atmosphere differs from their true locations due to atmospheric refraction. (e.g. sun and moon when they rise and set) The apparent elevation of an object above the horizon is apparent elevation angle = (true elevation angle) + (refraction angle)

The equations here are from Saemundsson/Bennett, whose calculations use a typical vertical profile of atmospheric density (i.e. temperature and pressure). The profiles can be rescaled to a particular surface temperature and pressure to approximately account for varying atmospheric conditions. Accurate refraction calculations should use fully specified vertical profile of temperature and pressure, which cannot be done here.

Parameters

true_elevation_angle (float): degrees above horizon of sun or other object
pressure (float (default=101325)): surface atmospheric pressure (Pa)
temperature_celsius (float (default=10)): surface atmospheric temperature (C)

Returns

angle (float): refraction angle in degrees. Value is zero when apparent elevation is below horizon

def sza( lat, lon, datetimeUTC, refraction=False, temperature=10.0, pressure=101325.0): View Source

105def solar_zenith_angle( lat, lon, datetimeUTC, 
106                       refraction=False, temperature=10., pressure=101325. ):
107    '''Solar zenith angle (degrees) for a given latitude, longitude, date and time.
108    
109    Accounts for equation of time and (optionally) for atmospheric refraction.
110    Altitude of the observer is not accounted for, which can be important when the sun 
111    is near the horizon. 
112    
113    Results are accurate to tenths of a degree, except where altitude is important
114    (< 20 degrees solar elevation)
115
116    Parameters
117    ----------
118    lat : float or ndarray
119        latitude in degrees
120    lon : float or ndarray
121        longitudes in degrees
122    datetimeUTC : pandas.Timestamp, datetime, or str
123        date and time in UTC
124    refraction : bool, optional (default=False)
125        specifies whether to account for atmospheric refraction
126    temperature : float or ndarray, optional (default=10)
127        surface atmospheric temperature (Celsius), only used for refraction calculation
128    pressure : float or ndarray, optional (default=101325)
129        surface atmospheric pressure (Pa), only used for refraction calculation
130    
131    Returns
132    -------
133    sza : float or ndarray
134        solar zenith angle in degrees at the designated locations and times
135        If refraction=False, this is the true solar zenith angle
136        If refraction=True, this is the apparent solar zenith angle
137    '''
138    # Convert to pandas Timestamp, if needed
139    datetimeUTC = _to_timestamp(datetimeUTC)
140
141    # Solar declination, degrees
142    dec = solar_declination( datetimeUTC )
143
144    # Hour angle, degrees
145    Ha = solar_hour_angle( lon, datetimeUTC )
146
147    # True solar zenith angle, radians
148    sza = np.arccos( np.sin(lat*pi180) * np.sin(dec*pi180) + \
149          np.cos(lat*pi180) * np.cos(dec*pi180) * np.cos(Ha*pi180) )
150
151    # Convert radians -> degrees
152    sza /= pi180
153
154    if refraction:
155        # Subtract refraction angle (degrees) from zenith angle.
156        # SZA is always smaller due to refraction.
157        sza -= refraction_angle( 90-sza, pressure, temperature )
158
159    return sza

Solar zenith angle (degrees) for a given latitude, longitude, date and time.

Accounts for equation of time and (optionally) for atmospheric refraction. Altitude of the observer is not accounted for, which can be important when the sun is near the horizon.

Results are accurate to tenths of a degree, except where altitude is important (< 20 degrees solar elevation)

Parameters

lat (float or ndarray): latitude in degrees
lon (float or ndarray): longitudes in degrees
datetimeUTC (pandas.Timestamp, datetime, or str): date and time in UTC
refraction (bool, optional (default=False)): specifies whether to account for atmospheric refraction
temperature (float or ndarray, optional (default=10)): surface atmospheric temperature (Celsius), only used for refraction calculation
pressure (float or ndarray, optional (default=101325)): surface atmospheric pressure (Pa), only used for refraction calculation

Returns

sza (float or ndarray): solar zenith angle in degrees at the designated locations and times If refraction=False, this is the true solar zenith angle If refraction=True, this is the apparent solar zenith angle

def saa(lat, lon, datetimeUTC): View Source

16def solar_azimuth_angle( lat, lon, datetimeUTC ):
17    '''Solar azimuth angle (degrees) for a latitude, longitude, date and time
18    
19    SAA is degrees clockwise from north.
20    
21    Parameters
22    ----------
23    lat : float or ndarray
24        latitude in degrees
25    lon : float or ndarray
26        longitudes in degrees
27    datetimeUTC : pandas.Timestamp, datetime, or str
28        date and time in UTC
29
30    Returns
31    -------
32    saa : float or ndarray
33        solar azimuth angle in degrees (clockwise from north)
34    '''
35    # Convert to pandas Timestamp, if needed
36    datetimeUTC = _to_timestamp(datetimeUTC)
37
38    # Solar declination, degrees
39    dec = solar_declination( datetimeUTC )
40
41    # Hour angle, degrees
42    Ha = solar_hour_angle( lon, datetimeUTC )
43
44    # Solar zenith angle, degrees
45    # Use true sza, without refraction
46    zen = sza( lat, lon, datetimeUTC, refraction=False )
47
48    # Solar azimuth angle, degrees
49    saa = np.arcsin( -np.sin( Ha*pi180 ) * np.cos( dec*pi180 ) /
50            np.sin( zen*pi180 ) ) / pi180
51
52    # Change range [-180,180] to [0,360]
53    return np.mod( saa+360, 360 )

Solar azimuth angle (degrees) for a latitude, longitude, date and time

SAA is degrees clockwise from north.

Parameters

lat (float or ndarray): latitude in degrees
lon (float or ndarray): longitudes in degrees
datetimeUTC (pandas.Timestamp, datetime, or str): date and time in UTC

Returns

saa (float or ndarray): solar azimuth angle in degrees (clockwise from north)

def sea( lat, lon, alt, datetimeUTC, refraction=False, temperature=10.0, pressure=101325.0): View Source

 55def solar_elevation_angle( lat, lon, alt, datetimeUTC,
 56                       refraction=False, temperature=10., pressure=101325. ):
 57    '''Solar elevation angle (degrees) above the horizon
 58
 59    The altitude parameter should be the vertical distance 
 60    above the surrounding terrain that defines the horizon,
 61    not necessarily the altitude above sea level or the altitude above ground level.
 62    For example, on a mountain peak that is 4000 m above sea level and 
 63    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
 64    For an observer on the plateau, the relevant altitude is 0 m.
 65
 66    See documentation for `solar_zenith_angle` and `horizon_zenith_angle`.
 67
 68    Parameters
 69    ----------
 70    lat : float or ndarray
 71        latitude in degrees
 72    lon : float or ndarray
 73        longitudes in degrees
 74    alt : float or ndarray
 75        altitude above surrounding terrain that defines the horizon, meters
 76    datetimeUTC : pandas.Timestamp, datetime, or str
 77        date and time in UTC
 78    refraction : bool, optional (default=False)
 79        specifies whether to account for atmospheric refraction
 80    temperature : float or ndarray, optional (default=10)
 81        surface atmospheric temperature (Celsius), only used for refraction calculation
 82    pressure : float or ndarray, optional (default=101325)
 83        surface atmospheric pressure (Pa), only used for refraction calculation
 84    
 85    Returns
 86    -------
 87    sea : float or ndarray
 88        solar elevation angle in degrees at the designated locations and times
 89        If refraction=False, this is the true solar elevation angle
 90        If refraction=True, this is the apparent solar elevation angle
 91    
 92    '''
 93
 94    if refraction and np.any(alt):
 95        warnings.warn( 'Atmospheric refraction is calculated for surface conditions, '
 96                    + 'but an altitude above the surface was specified',
 97                     category=UserWarning,
 98                     stacklevel=2 )
 99
100    sea = horizon_zenith_angle( alt ) \
101         - solar_zenith_angle( lat, lon, datetimeUTC, refraction, temperature, pressure )
102
103    return sea

Solar elevation angle (degrees) above the horizon

See documentation for solar_zenith_angle and horizon_zenith_angle.

Parameters

lat (float or ndarray): latitude in degrees
lon (float or ndarray): longitudes in degrees
alt (float or ndarray): altitude above surrounding terrain that defines the horizon, meters
datetimeUTC (pandas.Timestamp, datetime, or str): date and time in UTC
refraction (bool, optional (default=False)): specifies whether to account for atmospheric refraction
temperature (float or ndarray, optional (default=10)): surface atmospheric temperature (Celsius), only used for refraction calculation
pressure (float or ndarray, optional (default=101325)): surface atmospheric pressure (Pa), only used for refraction calculation

Returns

sea (float or ndarray): solar elevation angle in degrees at the designated locations and times If refraction=False, this is the true solar elevation angle If refraction=True, this is the apparent solar elevation angle

def equationOfTime(date): View Source

161def equation_of_time( date ):
162    '''Equation of time (degrees) for specified date
163    
164    Implements the "alternative equation" from Wikipedia, derived from
165    https://web.archive.org/web/20120323231813/http://www.green-life-innovators.org/tiki-index.php?page=The%2BLatitude%2Band%2BLongitude%2Bof%2Bthe%2BSun%2Bby%2BDavid%2BWilliams
166    Results checked against NOAA solar calculator and agree within 10 seconds.
167    
168    Note: Leap years are not accounted for.
169
170    Argument
171    --------
172    date : pandas.Timestamp, date, or datetime
173        date for calculation
174
175    Returns
176    -------
177    eot : float
178        equation of time in degrees on the specified date
179    '''
180    # Convert to pandas Timestamp, if needed
181    date = _to_timestamp(date)
182
183    # Equation of time, accounts for the solar day differing slightly from 24 hr
184    try:
185        doy = date.dt.dayofyear
186    except AttributeError:
187        doy = date.dayofyear
188    W = 360 / 365.24
189    A = W * (doy+10)
190    B = A + 1.914 * np.sin( W * (doy-2) * pi180 )
191    C = ( A - np.arctan2( np.tan(B*pi180), np.cos(23.44*pi180) ) / pi180 ) / 180
192
193    # Equation of time in minutes of an hour
194    eotmin = 720 * ( C - np.round(C) )
195
196    # Equation of time, minutes -> degrees
197    eot = eotmin / 60 * 360 / 24
198
199    return eot

Equation of time (degrees) for specified date

Note: Leap years are not accounted for.

Argument

date : pandas.Timestamp, date, or datetime date for calculation

Returns

eot (float): equation of time in degrees on the specified date

def solarDeclination(date): View Source

252def solar_declination( date ):
253    '''Calculate solar declination (degrees) for specified date
254    
255    Implements Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling
256    
257    Argument
258    --------
259    date : pandas.Timestamp, date, datetime, or str
260        date for calculation
261
262    Returns
263    -------
264    dec : float
265        solar declination in degrees at the specified date
266    '''
267
268    # Convert to pandas Timestamp, if needed
269    date = _to_timestamp(date)
270
271     # Number of days since beginning of 2000
272    NJD = ( date - np.datetime64('2000-01-01') )
273    try:
274        NJD = NJD.dt.days
275    except AttributeError:
276        NJD = NJD.days
277
278    # Obliquity, degrees
279    ob = 23.439 - 4e-7 * NJD
280
281    # Parameters for ecliptic, degrees
282    gm = 357.528 + 0.9856003 * NJD
283    lm = 280.460 + 0.9856474 * NJD
284
285    # Ecliptic longitude of sun, degrees
286    ec = lm + 1.915 * np.sin( gm * pi180 ) + 0.020 * np.sin( 2 * gm * pi180 )
287
288    #Solar declination, degrees
289    dec = np.arcsin( np.sin( ob * pi180 ) * np.sin( ec * pi180 ) ) / pi180
290
291    return dec

Calculate solar declination (degrees) for specified date

Implements Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling

Argument

date : pandas.Timestamp, date, datetime, or str date for calculation

Returns

dec (float): solar declination in degrees at the specified date