acgc.solar

Module for calculating solar position and TOA insolation

The functions here are are vectorized and generally broadcast over xarray dimensions, making this program faster than PySolar and pvlib. Calculations here use orbital parameters from the NOAA Solar Position Calculator, following Jean Meeus Astronomical Algorithms, unless the "fast" keyword is used. These accurate calculations are suitable for years -2000 to +3000. The "fast" calculations have lower accuracy orbital parameters and coarser approximation for the equation of time. All calculations here use geocentric solar position, neglecting the parallax effect of viewing the sun from different points on Earth (i.e. topocentric vs. geocentric in NREL SPA algorithm).

Accuracy: The NREL SPA algorithm in pvlib is used as an accurate reference. The maximum error in solar zenith angle computed here is 0.02° over 1900-2100. The maximum error in overall solar angular position is 0.022°. Large apparent differences in azimuth alone can occur when the sun is near zenith or nadir, where a small angular displacement results in a large azimuthal change.

The "fast" calculations have typical errors of ~0.2°.

NOAA Solar Calculator https://gml.noaa.gov/grad/solcalc/calcdetails.html

  1#!/usr/local/bin/env python3
  2'''Module for calculating solar position and TOA insolation
  3
  4The functions here are are vectorized and generally broadcast over xarray dimensions,
  5making this program faster than PySolar and pvlib. Calculations here use orbital parameters
  6from the NOAA Solar Position Calculator, following Jean Meeus Astronomical Algorithms, unless
  7the "fast" keyword is used. These accurate calculations are suitable for years -2000 to +3000.
  8The "fast" calculations have lower accuracy orbital parameters and coarser approximation 
  9for the equation of time. All calculations here use geocentric solar position, neglecting the 
 10parallax effect of viewing the sun from different points on Earth (i.e. topocentric vs. 
 11geocentric in NREL SPA algorithm).
 12
 13Accuracy:
 14The NREL SPA algorithm in pvlib is used as an accurate reference.
 15The maximum error in solar zenith angle computed here is 0.02° over 1900-2100.
 16The maximum error in overall solar angular position is 0.022°.
 17Large apparent differences in azimuth alone can occur when the sun is near zenith or nadir,
 18where a small angular displacement results in a large azimuthal change.
 19
 20The "fast" calculations have typical errors of ~0.2°.
 21
 22NOAA Solar Calculator
 23https://gml.noaa.gov/grad/solcalc/calcdetails.html
 24'''
 25
 26from collections import namedtuple
 27import warnings
 28import numpy as np
 29import pandas as pd
 30import xarray as xr
 31
 32def insolation_toa( lat, lon, datetime, solar_pos=None, **kwargs ):
 33    '''Insolation at top of the atmosphere, accounting for solar zenith angle
 34    
 35    Parameters
 36    ----------
 37    lat : float or ndarray
 38        latitude in degrees
 39    lon : float or ndarray
 40        longitudes in degrees
 41    datetime : datetime-like or str
 42        date and time. Include time zone or UTC will be assumed
 43    solar_pos : `SolarPositionResults`, optional
 44        solar position parameters from a prior call to `solar_position`
 45    **kwargs passed to `solar_zenith_angle`
 46    
 47    Returns
 48    -------
 49    Insolation : float
 50        radiation flux density accounting for solar zenith angle, W/m2
 51    '''
 52
 53    if solar_pos is None:
 54        solar_pos = solar_position( datetime )
 55
 56    S = solar_constant(datetime, solar_pos=solar_pos)
 57    sza = solar_zenith_angle(lat, lon, datetime, **kwargs, solar_pos=solar_pos )
 58
 59    return S * np.cos(sza)
 60
 61def solar_azimuth_angle( lat, lon, datetime, solar_pos=None ):
 62    '''Solar azimuth angle (degrees) for a latitude, longitude, date and time
 63    
 64    SAA is degrees clockwise from north.
 65    
 66    Parameters
 67    ----------
 68    lat : float or ndarray
 69        latitude in degrees
 70    lon : float or ndarray
 71        longitudes in degrees
 72    datetime : datetime-like or str
 73        date and time. Include time zone or UTC will be assumed
 74    solar_pos : `SolarPositionResults`, optional
 75        solar position parameters from a prior call to `solar_position`
 76
 77    Returns
 78    -------
 79    saa : float or ndarray
 80        solar azimuth angle in degrees (clockwise from north)
 81    '''
 82    # Convert to pandas Timestamp, if needed
 83    datetime = _to_datetime(datetime)
 84
 85    # Subsolar point, latitude longitude, degrees
 86    solar_lat = solar_latitude( datetime, solar_pos=solar_pos )
 87    solar_lon = solar_longitude( datetime, solar_pos=solar_pos )
 88
 89    # Vector pointing toward sun
 90    x = np.cos( solar_lat * pi180 ) * np.sin( (solar_lon - lon) * pi180 )
 91    y = np.cos( lat*pi180 ) * np.sin( solar_lat*pi180 ) \
 92        - np.sin( lat*pi180 ) * np.cos( solar_lat*pi180 ) \
 93            * np.cos( (solar_lon - lon) * pi180 )
 94
 95    # Azimuth angle from north, degrees
 96    saa = np.arctan2( x, y ) / pi180
 97
 98    # Change range [-180,180] to [0,360]
 99    return np.mod( saa+360, 360 )
100
101def solar_elevation_angle( lat, lon, datetime, alt=0,
102                        refraction=False, temperature=10., pressure=101325.,
103                        solar_pos=None ):
104    '''Solar elevation angle (degrees) above the horizon
105
106    The altitude parameter should be the vertical distance 
107    above the surrounding terrain that defines the horizon,
108    not necessarily the altitude above sea level or the altitude above ground level.
109    For example, on a mountain peak that is 4000 m above sea level and 
110    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
111    For an observer on the plateau, the relevant altitude is 0 m.
112
113    See documentation for `solar_zenith_angle` and `horizon_zenith_angle`.
114
115    Parameters
116    ----------
117    lat : float or ndarray
118        latitude in degrees
119    lon : float or ndarray
120        longitudes in degrees
121    datetime : datetime-like or str
122        date and time. Include time zone or UTC will be assumed
123    alt : float or ndarray (default=0)
124        altitude above surrounding terrain that defines the horizon, meters
125    refraction : bool (default=False)
126        specifies whether to account for atmospheric refraction
127    temperature : float or ndarray (default=10)
128        surface atmospheric temperature (Celsius), only used for refraction calculation
129    pressure : float or ndarray (default=101325)
130        surface atmospheric pressure (Pa), only used for refraction calculation
131    solar_pos : `SolarPositionResults`, optional
132        solar position parameters from a prior call to `solar_position`
133    
134    Returns
135    -------
136    sea : float or ndarray
137        solar elevation angle in degrees at the designated locations and times
138        If refraction=False, this is the true solar elevation angle
139        If refraction=True, this is the apparent solar elevation angle
140    
141    '''
142
143    if refraction and np.any(alt):
144        warnings.warn( 'Atmospheric refraction is calculated for surface conditions, '
145                    + 'but an altitude above the surface was specified',
146                     category=UserWarning,
147                     stacklevel=2 )
148
149    sea = horizon_zenith_angle( alt ) \
150         - solar_zenith_angle( lat, lon, datetime, refraction, temperature, pressure, 
151                               solar_pos=solar_pos )
152
153    return sea
154
155def solar_zenith_angle( lat, lon, datetime,
156                        refraction=False, temperature=10., pressure=101325.,
157                        solar_pos=None ):
158    '''Solar zenith angle (degrees) for a given latitude, longitude, date and time.
159    
160    Accounts for equation of time and (optionally) for atmospheric refraction.
161    Altitude of the observer is not accounted for, which can be important when the sun 
162    is near the horizon. 
163    
164    Results are accurate to tenths of a degree, except where altitude is important
165    (< 20 degrees solar elevation)
166
167    Parameters
168    ----------
169    lat : float or ndarray
170        latitude in degrees
171    lon : float or ndarray
172        longitudes in degrees
173    datetime : datetime-like or str
174        date and time. Include time zone or UTC will be assumed
175    refraction : bool (default=False)
176        specifies whether to account for atmospheric refraction
177    temperature : float or ndarray (default=10)
178        surface atmospheric temperature (Celsius), only used for refraction calculation
179    pressure : float or ndarray (default=101325)
180        surface atmospheric pressure (Pa), only used for refraction calculation
181    solar_pos : `SolarPositionResults`, optional
182        solar position parameters from a prior call to `solar_position`
183    
184    Returns
185    -------
186    sza : float or ndarray
187        solar zenith angle in degrees at the designated locations and times
188        If refraction=False, this is the true solar zenith angle
189        If refraction=True, this is the apparent solar zenith angle
190    '''
191    # Convert to pandas Timestamp, if needed
192    datetime = _to_datetime(datetime)
193
194    # Solar declination, degrees
195    if solar_pos is None:
196        dec = solar_declination( datetime )
197    else:
198        dec = solar_pos.declination
199
200    # Hour angle, degrees
201    Ha = solar_hour_angle( lon, datetime, solar_pos )
202
203    # True solar zenith angle, radians
204    sza = np.arccos( np.sin(lat*pi180) * np.sin(dec*pi180) + \
205          np.cos(lat*pi180) * np.cos(dec*pi180) * np.cos(Ha*pi180) )
206
207    # Convert radians -> degrees
208    sza /= pi180
209
210    if refraction:
211        # Subtract refraction angle (degrees) from zenith angle.
212        # SZA is always smaller due to refraction.
213        sza -= refraction_angle( 90-sza, pressure, temperature )
214
215    return sza
216
217def sunrise_time( *args, **kwargs ):
218    '''Compute sunrise time
219    
220    See `sun_times` for Parameters.'''
221    result = sun_times( *args, **kwargs )
222    return result[0]
223
224def sunset_time( *args, **kwargs ):
225    '''Compute sunset time
226    
227    See `sun_times` for Parameters.'''
228    result = sun_times( *args, **kwargs )
229    return result[1]
230
231def day_length( *args, **kwargs ):
232    '''Compute length of daylight
233    
234    See `sun_times` for Parameters.'''
235    result = sun_times( *args, **kwargs )
236    return result[2]
237
238def solar_noon( *args, **kwargs ):
239    '''Compute time of solar noon (meridian transit)
240    
241    See `sun_times` for Parameters.'''
242    result = sun_times( *args, **kwargs )
243    return result[3]
244
245def solar_constant( datetime, solar_pos=None ):
246    '''Compute solar constant for specific date or dates
247    
248    Parameters
249    ----------
250    datetime : datetime-like
251    solar_pos : `SolarPositionResults`, optional
252        solar position parameters from a prior call to `solar_position`
253
254    Returns
255    -------
256    S : float
257        Solar direct beam radiation flux density, W/m2
258    '''
259    if solar_pos is None:
260        solar_pos = solar_position( datetime )
261    S = 1361/solar_pos.distance**2
262
263    return S
264
265def solar_declination( datetime, fast=False, solar_pos=None ):
266    '''Calculate solar declination (degrees) for specified date
267        
268    Parameters
269    ----------
270    datetime : datetime-like or str
271        date and time. Include time zone or UTC will be assumed
272    fast : bool (default=False)
273        Specifies using a faster but less accurate calculation
274    solar_pos : `SolarPositionResults`, optional
275        solar position parameters from a prior call to `solar_position`
276
277    Returns
278    -------
279    dec : float
280        solar declination in degrees at the specified date
281    '''
282    # Convert to pandas Timestamp, if needed
283    datetime = _to_datetime(datetime)
284
285    # Select the accurate or fast calculation
286    accurate = not fast
287
288    if solar_pos is not None:
289        dec = solar_pos.declination
290
291    else:
292        if accurate:
293
294            # Solar declination, degrees
295            dec, junk, junk, junk, junk = solar_position( datetime )
296
297        else:
298            # The fast method implements
299            # Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling
300
301            # Number of days since beginning of 2000
302            NJD = datetime - np.datetime64('2000-01-01')
303            try:
304                NJD = NJD.dt.days
305            except AttributeError:
306                NJD = NJD.days
307
308            # Obliquity, degrees
309            ob = 23.439 - 4e-7 * NJD
310
311            # Parameters for ecliptic, degrees
312            gm = 357.528 + 0.9856003 * NJD
313            lm = 280.460 + 0.9856474 * NJD
314
315            # Ecliptic longitude of sun, degrees
316            ec = lm + 1.915 * np.sin( gm * pi180 ) + 0.020 * np.sin( 2 * gm * pi180 )
317
318            #Solar declination, degrees
319            dec = np.arcsin( np.sin( ob * pi180 ) * np.sin( ec * pi180 ) ) / pi180
320
321    return dec
322
323def solar_latitude( *args, **kwargs ):
324    '''Latitude of the subsolar point
325    
326    Parameters
327    ----------
328    datetime : datetime-like or str
329        date and time. Include time zone or UTC will be assumed
330    fast : bool
331        see `solar_declination`
332    solar_pos : `SolarPositionResults`, optional
333        solar position parameters from a prior call to `solar_position`
334        
335    Returns
336    -------
337    latitude : float
338        degrees of latitude
339    '''
340    return solar_declination( *args, **kwargs )
341
342def solar_longitude( datetime, solar_pos=None ):
343    '''Longitude of the subsolar point, degrees
344    
345    Parameters
346    ----------
347    datetimeUTC : datetime-like or str
348        date and time. Include time zone or UTC will be assumed
349    solar_pos : `SolarPositionResults`, optional
350        solar position parameters from a prior call to `solar_position`
351    
352    Returns
353    -------
354    longitude : float
355        degrees of longitude
356    '''
357    # Convert to pandas Timestamp, if needed
358    datetimeUTC, tz_in = _to_datetime_utc(datetime)
359
360    # Longitude of subsolar point, degrees
361    # Equation of time will be added below
362    try:
363        # Treat as xarray.DataArray or pandas.Series
364        solar_lon = - 15 * ( datetimeUTC.dt.hour +
365                          datetimeUTC.dt.minute / 60 +
366                          datetimeUTC.dt.second / 3600 - 12 )
367    except AttributeError:
368        solar_lon = - 15 * ( datetimeUTC.hour +
369                          datetimeUTC.minute / 60 +
370                          datetimeUTC.second / 3600 - 12 )
371
372    # Add equation of time to the solar longitude, degrees
373    if solar_pos is None:
374        eot = equation_of_time( datetimeUTC, degrees=True )
375    else:
376        eot = solar_pos.equation_of_time
377    solar_lon -= eot
378
379    return solar_lon
380
381def solar_hour_angle( lon, datetime, solar_pos=None ):
382    '''Solar hour angle (degrees) for specified longitude, date and time
383
384    Hour angle is the angular displacement of the sun from the local meridian.
385    It is zero at local noon, negative in the morning, and positive is afternoon.
386    
387    Parameters
388    ----------
389    lon : float
390        longitude in degrees east
391    datetimeUTC : datetime-like or str
392        date and time. Include time zone or UTC will be assumed
393    solar_pos : `SolarPositionResults`, optional
394        solar position parameters from a prior call to `solar_position`
395    
396    Returns
397    -------
398    ha : float
399        hour angle in degrees at the specified location and time
400    '''
401
402    # Subsolar longitude, degrees
403    solar_lon = solar_longitude(datetime, solar_pos )
404
405    # Hour angle, degrees
406    Ha = lon - solar_lon
407
408    return Ha
409
410def equation_of_time( datetime, degrees=False, fast=False ):
411    '''Equation of time for specified date
412    
413    Accounts for the solar day being slightly different from 24 hours
414
415    Parameters
416    ----------
417    datetime : datetime-like or str
418        date and time. Include time zone or UTC will be assumed
419    degrees : bool (default=False)
420        If True, then return value in compass degrees
421        If False, then return value in minutes of an hour
422    fast : bool (default=False)
423        specifies whether to use a faster, but less accurate calculation
424        
425    Returns
426    -------
427    eot : float
428        equation of time on the specified date, degrees or minutes
429    '''
430    # Convert to pandas Timestamp, if needed
431    datetime = _to_datetime(datetime)
432
433    # Determine whether to use the fast or accurate calculation
434    accurate = not fast
435
436    if accurate:
437
438        # Equation of time, minutes
439        junk, junk, eot, junk, junk = solar_position( datetime )
440
441    else:
442        # Implements the "alternative equation" from Wikipedia, derived from
443        # 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
444        # When compared to the NREL SPA algorithm, differences reach are up to about 0.5 minute.
445        # Note: Leap years are not accounted for.
446
447        # Equation of time, accounts for the solar day differing slightly from 24 hr
448        try:
449            doy = datetime.dt.dayofyear
450        except AttributeError:
451            doy = datetime.dayofyear
452        W = 360 / 365.24
453        A = W * (doy+10)
454        B = A + 1.914 * np.sin( W * (doy-2) * pi180 )
455        C = ( A - np.arctan2( np.tan(B*pi180), np.cos(23.44*pi180) ) / pi180 ) / 180
456
457        # Equation of time in minutes of an hour (1440 minutes per day)
458        eot = 720 * ( C - np.round(C) )
459
460    # Equation of time, minutes -> degrees (360 degrees per day)
461    if degrees:
462        eot = eot / 60 * 360 / 24
463
464    return eot
465
466def solar_position( datetime ):
467    '''Compute position of sun (declination, right ascension, equation of time, distance)) on specified date
468    
469    Calculations follow the NOAA solar calculator spreadsheet
470    Applicable to years 1900-2100.
471
472    Parameters
473    ----------
474    date : datetime-like or str
475        date and time. Include time zone or UTC will be assumed
476    
477    Returns
478    -------
479    result : `SolarPositionResults`
480        declination, right ascension, equation of time, Earth-sun distance, datetime
481    '''
482    # Ensure time is Timestamp in UTC
483    datetimeUTC, tz_in = _to_datetime_utc(datetime)
484
485    # Raise warning if any dates are outside date range
486    # recommended for orbital parameters used here
487    if np.logical_or( np.any( datetimeUTC < np.datetime64('1900-01-01') ),
488                      np.any( datetimeUTC > np.datetime64('2100-01-01') ) ):
489        warnings.warn('Solar position accuracy declines for dates outside 1900-2100', \
490                      RuntimeWarning )
491
492    # Number of days since 1 Jan 2000
493    NJD = datetimeUTC - np.datetime64('2000-01-01')
494    try:
495        NJD = NJD.dt.days \
496            + NJD.dt.seconds / 86400.
497    except AttributeError:
498        NJD = NJD.days \
499            + NJD.seconds / 86400.
500
501    # Julian day (since 12:00UTC 1 Jan 4713 BCE)
502    NJD += 2451544.50
503
504    # Julian century
505    JC = (NJD-2451545)/36525
506
507    # Earth orbital eccentricity, unitless
508    ec = 0.016708634 - JC*( 0.000042037 + 0.0000001267*JC )
509
510    # Earth mean orbital obliquity, degree
511    mean_ob = 23 + ( 26 + ( (21.448
512                             - JC * (46.815
513                                    + JC * (0.00059 - JC * 0.001813) ) ) )/60 )/60
514
515    # Earth true orbital obliquity, corrected for nutation, degree 
516    ob = mean_ob + 0.00256 * np.cos( (125.04 - 1934.136*JC ) * pi180 )
517
518    # Sun Mean ecliptic longitude, degree
519    mean_ec_lon = np.mod( 280.46646 + JC*( 36000.76983 + JC*0.0003032 ), 360 )
520
521    # Sun Mean anomaly, degree
522    mean_anom = 357.52911 + JC*( 35999.05029 - 0.0001537*JC )
523
524    # Sun Equation of center, degree
525    eq_center = np.sin(mean_anom*pi180) * (1.914602 - JC*( 0.004817 + 0.000014*JC )) \
526                    + np.sin(2*mean_anom*pi180) * (0.019993 - 0.000101*JC) \
527                    + np.sin(3*mean_anom*pi180) * 0.000289
528
529    # Sun True ecliptic longitude, degrees
530    true_ec_lon = mean_ec_lon + eq_center
531
532    # Sun True anomaly, degree
533    true_anom = mean_anom + eq_center
534
535    # Earth-Sun distance, AU
536    distance = (1.000001018 * (1-ec**2) ) / (1 + ec * np.cos( true_anom * pi180 ))
537
538    # Sun Apparent ecliptic longitude, corrected for nutation, degrees
539    ec_lon = true_ec_lon - 0.00569 - 0.00478 * np.sin( (125.04 - 1934.136*JC ) * pi180)
540
541    # Sun Right ascension, deg
542    right_ascension = np.arctan2( np.cos(ob*pi180) * np.sin(ec_lon*pi180),
543                                  np.cos(ec_lon*pi180) ) / pi180
544
545    # Sun Declination, deg
546    declination = np.arcsin( np.sin(ob*pi180) * np.sin(ec_lon*pi180) ) / pi180
547
548    # var y
549    vary = np.tan( ob/2 * pi180 )**2
550
551    # Equation of time, minutes
552    eot = vary * np.sin( 2 * mean_ec_lon * pi180) \
553        - 2 * ec * np.sin( mean_anom * pi180 ) \
554        + 4 * ec * vary * np.sin( mean_anom * pi180 ) * np.cos( 2 * mean_ec_lon * pi180) \
555        - 0.5 * vary**2 * np.sin( 4 * mean_ec_lon * pi180) \
556        - 1.25 * ec**2 * np.sin( 2 * mean_anom * pi180)
557    eot = eot * 4 / pi180
558
559    # Bundle results
560    result = SolarPositionResults(declination,
561                                  right_ascension,
562                                  eot,
563                                  distance,
564                                  datetimeUTC)
565
566    return result
567
568def sun_times( lat, lon, datetime, tz_out=None, sza_sunrise=90.833,
569              fast=False, solar_pos=None ):
570    '''Compute times of sunrise, sunset, solar noon, and day length
571    
572    Common options for solar zenith angle at sunrise
573    1. 90.833 for first edge of sun rising, typical (0.567°) refraction (default)
574    2. 90.267 for first edge of sun rising, no refraction
575    3. 90 degrees for center of sun rising, no refraction
576
577    Note: horizon_zenith_angle can be used to compute a more accurate horizon location
578    for sites at elevation.
579    
580    Parameters
581    ----------
582    lat : float or ndarray
583        latitude in degrees
584    lon : float or ndarray
585        longitudes in degrees
586    datetime : datetime-like or str
587        datetime, provide a time zone or UTC will be assumed 
588    tz_out : str, pytz.timezone, datetime.tzinfo, optional
589        timezone to be used for output times. 
590        If None is provided, then result will be in same time zone as input or UTC
591    sza_sunrise : float (default=90.833)
592        Solar zenith angle at which sunrise and sunset are calculated, degrees
593    fast : bool (default=False)
594        Select a faster but less accurate calculation
595    solar_pos : `SolarPositionResults`, optional
596        solar position parameters from a prior call to `solar_position`
597
598    Returns
599    -------
600    result : `SunTimesResults`
601        times of sunrise, sunset, solar noon, and day length
602    '''
603    # Convert to pandas Timestamp in UTC, if needed
604    datetimeUTC, tz_in = _to_datetime_utc(datetime)
605
606    # If no output timezone is specified, use the input time zone
607    if tz_out is None:
608        tz_out = tz_in
609
610    # Select fast or accurate calculation
611    accurate = not fast
612
613    # Solar declination (degrees) and equation of time (minutes)
614    if solar_pos is not None:
615        dec = solar_pos.declination
616        eot = solar_pos.eot
617    else:
618        if accurate:
619            dec, junk, eot, junk, junk = solar_position( datetimeUTC )
620        else:
621            dec = solar_declination( datetimeUTC )
622            eot = equation_of_time( datetimeUTC )
623
624    # Sunrise hour angle, degree
625    # Degrees east of the local meridian where sun rises
626    ha_sunrise = np.arccos( np.cos(sza_sunrise*pi180) /
627                           (np.cos(lat*pi180) * np.cos(dec*pi180))
628                           - np.tan(lat*pi180)*np.tan(dec*pi180) ) / pi180
629
630    # Solar noon, local standard time, day fraction
631    solar_noon = (720 - 4*lon - eot ) / 1440
632
633    # Sunrise and sunset, local standard time, day fraction
634    t_sunrise = solar_noon - 4 * ha_sunrise / 1440
635    t_sunset  = solar_noon + 4 * ha_sunrise / 1440
636
637    # Midnight UTC
638    # datetimeUTC is in UTC but time-zone-naive
639    try:
640        # Series time objects
641        dateUTC = datetimeUTC.dt.floor('D')
642    except AttributeError:
643        # Scalar time objects
644        dateUTC = datetimeUTC.floor('D')
645
646    # Convert day fraction -> date time
647    solar_noon = dateUTC + solar_noon * pd.Timedelta( 1, 'day' )
648    t_sunrise  = dateUTC + t_sunrise  * pd.Timedelta( 1, 'day' )
649    t_sunset   = dateUTC + t_sunset   * pd.Timedelta( 1, 'day' )
650
651    # Convert to output timezone, if any is provided
652    if tz_out is not None:
653        if isinstance(solar_noon,(xr.DataArray,np.ndarray)) or \
654            isinstance(t_sunrise,(xr.DataArray,np.ndarray)):
655            # These types don't localize tz, but we can add offset to the tz-naive time
656            if hasattr(datetimeUTC,'tz_localize'):
657                # For scalar datetime, there is a single time offset, which we can add
658                utcoffset = np.timedelta64( datetimeUTC.tz_localize(tz_out).utcoffset() )
659                solar_noon += utcoffset
660                t_sunrise  += utcoffset
661                t_sunset   += utcoffset
662            else:
663                # For Series datetime, there are potentially multiple offsets. We can only add one
664                unique_datetimeUTC = pd.DatetimeIndex(np.unique(datetimeUTC))
665                unique_utcoffsets = np.unique( unique_datetimeUTC.tz_localize('UTC') \
666                                        - unique_datetimeUTC.tz_localize(tz_out) )
667                if len(unique_utcoffsets)==1:
668                    utcoffset = unique_utcoffsets[0]
669                    solar_noon += utcoffset
670                    t_sunrise  += utcoffset
671                    t_sunset   += utcoffset
672                else:
673                    # We might be able to handle multiple offsets if we can
674                    # determine which dimension of `solar_noon` and `t_sunrise`
675                    # is the time dimension.
676                    raise ValueError('Multiple timezone offsets not supported. '
677                                     +'Request output in UTC or reduce number of input times.')
678        else:
679            try:
680                # Series time objects
681                solar_noon = solar_noon.dt.tz_localize('UTC')\
682                                    .dt.tz_convert(tz_out)
683                t_sunrise  = t_sunrise.dt.tz_localize('UTC')\
684                                    .dt.tz_convert(tz_out)
685                t_sunset   = t_sunset.dt.tz_localize('UTC')\
686                                    .dt.tz_convert(tz_out)
687            except AttributeError:
688                # Scale time objects
689                solar_noon = solar_noon.tz_localize('UTC')\
690                                    .tz_convert(tz_out)
691                t_sunrise  = t_sunrise.tz_localize('UTC')\
692                                    .tz_convert(tz_out)
693                t_sunset   = t_sunset.tz_localize('UTC')\
694                                    .tz_convert(tz_out)
695
696    # Sunlight duration, minutes
697    day_length = 8 * ha_sunrise * pd.Timedelta(1, 'minute')
698
699    result = SunTimesResults(t_sunrise, t_sunset, day_length, solar_noon, datetimeUTC)
700
701    return result
702
703def horizon_zenith_angle( alt ):
704    '''Angle from the zenith to the horizon
705    
706    The horizon is the locii of points where a line from the 
707    observation location to the ellipsoid is tangent to the ellipsoid surface.
708    
709    The altitude parameter should be the vertical distance 
710    above the surrounding terrain that defines the horizon,
711    not necessarily the altitude above sea level or the altitude above ground level.
712    For example, on a mountain peak that is 4000 m above sea level and 
713    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
714    For an observer on the plateau, the relevant altitude is 0 m.
715
716    The implementation below assumes a spherical Earth.
717    Results using the WGS84 ellipsoid (see commented code below)
718    differ from the spherical case by << 1°. Terrain,
719    which is neglected here, has a larger effect on the horizon
720    location, so the simpler spherical calculation is appropriate. 
721
722    Parameters
723    ----------
724    lat : float or ndarray
725        latitude in degrees
726    alt : float or ndarray
727        altitude above surrounding terrain that defines the horizon, meters
728        
729    Returns
730    -------
731    hza : float or ndarray
732        horizon zenith angle in degrees
733    '''
734
735    # WGS84 ellipsoid parameters
736    # semi-major radius, m
737    r_earth = 6378137.0
738    # ellipsoidal flattening, unitless
739    f = 1/298.257223563
740
741    # Horizon zenith angle, degrees (spherical earth)
742    hza = 180 - np.arcsin( r_earth / ( r_earth + alt ) ) / pi180
743
744    ## Ellipsoidal Earth
745    # # Eccentricity of ellipsoid
746    # ecc = f * (2-f)
747    # # Local (i.e. prime vertical) radius of curvature at latitude
748    # N = r_earth / np.sqrt( 1 - ecc**2 * np.sin(lat*pi180)**2 )
749    # # Horizon zenith angle, degrees
750    # hza = 180 - np.arcsin( N / (N+alt) ) / pi180
751
752    return hza
753
754def refraction_angle( true_elevation_angle, pressure=101325., temperature_celsius=10. ):
755    '''Atmospheric refraction angle for light passing through Earth's atmosphere
756
757    The apparent locations in the sky of objects outsides Earth's atmosphere 
758    differs from their true locations due to atmospheric refraction. 
759    (e.g. sun and moon when they rise and set)
760    The apparent elevation of an object above the horizon is
761    apparent elevation angle = (true elevation angle) + (refraction angle)
762    
763    The equations here are from Saemundsson/Bennett, whose calculations use
764    a typical vertical profile of atmospheric density (i.e. temperature and pressure).
765    The profiles can be rescaled to a particular surface temperature and pressure
766    to approximately account for varying atmospheric conditions.
767    Accurate refraction calculations should use fully specified vertical profile
768    of temperature and pressure, which cannot be done here.
769
770    Parameters
771    ----------
772    true_elevation_angle : float
773        degrees above horizon of sun or other object
774    pressure : float (default=101325)
775        surface atmospheric pressure (Pa)
776    temperature_celsius : float (default=10)
777        surface atmospheric temperature (C)
778
779    Returns
780    -------
781    angle : float
782        refraction angle in degrees. Value is zero when apparent elevation is below horizon
783    '''
784    # Refraction angle, arcminutes
785    R = 1.02 / np.tan( ( true_elevation_angle + 10.3 / (true_elevation_angle + 5.11) ) * pi180 )
786    # Account for temperature and pressure, arcminutes
787    R = R * pressure / 101325 * 283 / ( 273 + temperature_celsius )
788    # Convert arcminutes -> degrees
789    R /= 60
790
791    # Result must be positive
792    R = np.maximum(R,0)
793
794    # Refraction defined only when the apparent elevation angle is positive
795    # Set refraction to zero when the apparent elevation is below horizon
796    refraction_angle = np.where( true_elevation_angle + R <= 0, 0, R)
797
798    return refraction_angle
799
800def _to_datetime(time_in):
801    '''Convert to Timestamp, Series of datetime64, or DataArray of datetime64 
802    
803    Parameters
804    ----------
805    time_in : datetime-like or array
806        time to be converted
807
808    Returns
809    -------
810    time_out : pandas.Timestamp, pandas.Series of datetime64, DataArray of datetime64
811    '''
812    if hasattr(time_in,'dt'):
813        time_out = time_in
814    # elif isinstance(time_in, pd.DatetimeIndex ):
815    #     # Unnecessary; DatetimeIndex will work fine in the else cases
816    #     time_out = pd.Series(time_in)
817    #     tz = time_out.dt.tz
818    else:
819        try:
820            # Convert list of times
821            time_out = pd.Series(pd.DatetimeIndex(time_in))
822        except TypeError:
823            # Single datetime or str
824            time_out = pd.Timestamp(time_in)
825
826    return time_out
827
828def _to_datetime_utc( datetime_in ):
829    '''Convert to Timestamp, Series of datetime64, or DataArray of datetime64 and in UTC
830    
831    Parameters
832    ----------
833    datetime_in : datetime-like
834        date and time to be converted
835        
836    Returns
837    -------
838    datetimeUTC : Timestamp, Series, or DataArray
839        date and time in UTC but tz-naive
840    tz_in : datetime.timezone
841        timezone of datetime_in
842    '''
843
844    # Ensure input is a timestamp
845    datetime_in = _to_datetime( datetime_in )
846
847    # Convert to UTC, then strip timezone
848    try:
849        if hasattr(datetime_in,'dt'):
850            # Pandas Series objects
851            tz_in = datetime_in.dt.tz
852            datetimeUTC = datetime_in.dt.tz_convert('UTC').dt.tz_localize(None)
853        else:
854            # Scalar time objects
855            tz_in = datetime_in.tzinfo
856            datetimeUTC = datetime_in.tz_convert('UTC').tz_localize(None)
857    except (TypeError, AttributeError):
858        # tz-naive time objects: Timestamp, Series, (all) DataArrays
859        # No timezone info, so assume it is already UTC
860        warnings.warn('Time does not contain timezone. UTC will be assumed',RuntimeWarning)
861        datetimeUTC = datetime_in
862        tz_in = None
863
864    return datetimeUTC, tz_in
865
866SolarPositionResults = namedtuple('SolarPositionResults',
867                    'declination right_ascension equation_of_time distance datetimeUTC')
868'''Namedtuple containing results of `solar_position`
869
870All values are geocentric, strictly accurate for the center of the Earth, not a point on 
871Earth's surface. The parallax angle from Earth's center to surface is 4e-5 degrees.
872
873Attributes
874----------
875declination : float
876    position of the sun relative to Earth's equatorial plane, degrees
877right_ascension : float
878    position of the sun along Earth's equatorial plane relative to the vernal equinox, degrees
879equation_of_time : float
880    equation of time (minutes) between mean solar time and true solar time.
881    Divide by 4 minutes per degree to obtain equation of time in degrees.
882distance : float
883    Earth-sun distance in AU (1 AU = 1.495978707e11 m)
884datetimeUTC : Timestamp, Series, or DataArray of numpy.datetime64
885    date and time input for calculations, UTC
886'''
887
888SunTimesResults = namedtuple('SunTimesResults',
889                    'sunrise sunset day_length solar_noon datetimeUTC')
890'''Namedtuple containing results of `sun_times`
891
892Attributes
893----------
894sunrise : pandas.DatetimeIndex
895    sunrise time, UTC if not specified otherwise
896sunset : pandas.DatetimeIndex
897    sunset time, UTC if not specified otherwise
898day_length : pandas.Timedelta
899    duration of daylight
900solar_noon : pandas.DatetimeIndex
901    time of meridian transit, UTC if not specified otherwise
902datetimeUTC : Timestamp, Series, or DataArray of numpy.datetime64
903    date and time input for calculations, UTC
904'''
905
906pi180 = np.pi / 180
907'''Constant $\pi/180$'''
908
909# Aliases for functions
910sza = solar_zenith_angle
911'''Alias for `solar_zenith_angle`'''
912saa = solar_azimuth_angle
913'''Alias for `solar_azimuth_angle`'''
914sea = solar_elevation_angle
915'''Alias for `solar_elevation_angle`'''
916
917# Deprecated aliases
918equationOfTime = equation_of_time
919'''Alias for `equation_of_time`'''
920solarDeclination = solar_declination
921'''Alias for `solar_declination`'''
def insolation_toa(lat, lon, datetime, solar_pos=None, **kwargs): 
33def insolation_toa( lat, lon, datetime, solar_pos=None, **kwargs ):
34    '''Insolation at top of the atmosphere, accounting for solar zenith angle
35    
36    Parameters
37    ----------
38    lat : float or ndarray
39        latitude in degrees
40    lon : float or ndarray
41        longitudes in degrees
42    datetime : datetime-like or str
43        date and time. Include time zone or UTC will be assumed
44    solar_pos : `SolarPositionResults`, optional
45        solar position parameters from a prior call to `solar_position`
46    **kwargs passed to `solar_zenith_angle`
47    
48    Returns
49    -------
50    Insolation : float
51        radiation flux density accounting for solar zenith angle, W/m2
52    '''
53
54    if solar_pos is None:
55        solar_pos = solar_position( datetime )
56
57    S = solar_constant(datetime, solar_pos=solar_pos)
58    sza = solar_zenith_angle(lat, lon, datetime, **kwargs, solar_pos=solar_pos )
59
60    return S * np.cos(sza)

Insolation at top of the atmosphere, accounting for solar zenith angle

Parameters
  • lat (float or ndarray): latitude in degrees
  • lon (float or ndarray): longitudes in degrees
  • datetime (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
  • **kwargs passed to solar_zenith_angle
Returns
  • Insolation (float): radiation flux density accounting for solar zenith angle, W/m2
def solar_azimuth_angle(lat, lon, datetime, solar_pos=None): 
 62def solar_azimuth_angle( lat, lon, datetime, solar_pos=None ):
 63    '''Solar azimuth angle (degrees) for a latitude, longitude, date and time
 64    
 65    SAA is degrees clockwise from north.
 66    
 67    Parameters
 68    ----------
 69    lat : float or ndarray
 70        latitude in degrees
 71    lon : float or ndarray
 72        longitudes in degrees
 73    datetime : datetime-like or str
 74        date and time. Include time zone or UTC will be assumed
 75    solar_pos : `SolarPositionResults`, optional
 76        solar position parameters from a prior call to `solar_position`
 77
 78    Returns
 79    -------
 80    saa : float or ndarray
 81        solar azimuth angle in degrees (clockwise from north)
 82    '''
 83    # Convert to pandas Timestamp, if needed
 84    datetime = _to_datetime(datetime)
 85
 86    # Subsolar point, latitude longitude, degrees
 87    solar_lat = solar_latitude( datetime, solar_pos=solar_pos )
 88    solar_lon = solar_longitude( datetime, solar_pos=solar_pos )
 89
 90    # Vector pointing toward sun
 91    x = np.cos( solar_lat * pi180 ) * np.sin( (solar_lon - lon) * pi180 )
 92    y = np.cos( lat*pi180 ) * np.sin( solar_lat*pi180 ) \
 93        - np.sin( lat*pi180 ) * np.cos( solar_lat*pi180 ) \
 94            * np.cos( (solar_lon - lon) * pi180 )
 95
 96    # Azimuth angle from north, degrees
 97    saa = np.arctan2( x, y ) / pi180
 98
 99    # Change range [-180,180] to [0,360]
100    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
  • datetime (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
Returns
  • saa (float or ndarray): solar azimuth angle in degrees (clockwise from north)
def solar_elevation_angle( lat, lon, datetime, alt=0, refraction=False, temperature=10.0, pressure=101325.0, solar_pos=None): 
102def solar_elevation_angle( lat, lon, datetime, alt=0,
103                        refraction=False, temperature=10., pressure=101325.,
104                        solar_pos=None ):
105    '''Solar elevation angle (degrees) above the horizon
106
107    The altitude parameter should be the vertical distance 
108    above the surrounding terrain that defines the horizon,
109    not necessarily the altitude above sea level or the altitude above ground level.
110    For example, on a mountain peak that is 4000 m above sea level and 
111    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
112    For an observer on the plateau, the relevant altitude is 0 m.
113
114    See documentation for `solar_zenith_angle` and `horizon_zenith_angle`.
115
116    Parameters
117    ----------
118    lat : float or ndarray
119        latitude in degrees
120    lon : float or ndarray
121        longitudes in degrees
122    datetime : datetime-like or str
123        date and time. Include time zone or UTC will be assumed
124    alt : float or ndarray (default=0)
125        altitude above surrounding terrain that defines the horizon, meters
126    refraction : bool (default=False)
127        specifies whether to account for atmospheric refraction
128    temperature : float or ndarray (default=10)
129        surface atmospheric temperature (Celsius), only used for refraction calculation
130    pressure : float or ndarray (default=101325)
131        surface atmospheric pressure (Pa), only used for refraction calculation
132    solar_pos : `SolarPositionResults`, optional
133        solar position parameters from a prior call to `solar_position`
134    
135    Returns
136    -------
137    sea : float or ndarray
138        solar elevation angle in degrees at the designated locations and times
139        If refraction=False, this is the true solar elevation angle
140        If refraction=True, this is the apparent solar elevation angle
141    
142    '''
143
144    if refraction and np.any(alt):
145        warnings.warn( 'Atmospheric refraction is calculated for surface conditions, '
146                    + 'but an altitude above the surface was specified',
147                     category=UserWarning,
148                     stacklevel=2 )
149
150    sea = horizon_zenith_angle( alt ) \
151         - solar_zenith_angle( lat, lon, datetime, refraction, temperature, pressure, 
152                               solar_pos=solar_pos )
153
154    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
  • datetime (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • alt (float or ndarray (default=0)): altitude above surrounding terrain that defines the horizon, meters
  • refraction (bool (default=False)): specifies whether to account for atmospheric refraction
  • temperature (float or ndarray (default=10)): surface atmospheric temperature (Celsius), only used for refraction calculation
  • pressure (float or ndarray (default=101325)): surface atmospheric pressure (Pa), only used for refraction calculation
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
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, datetime, refraction=False, temperature=10.0, pressure=101325.0, solar_pos=None): 
156def solar_zenith_angle( lat, lon, datetime,
157                        refraction=False, temperature=10., pressure=101325.,
158                        solar_pos=None ):
159    '''Solar zenith angle (degrees) for a given latitude, longitude, date and time.
160    
161    Accounts for equation of time and (optionally) for atmospheric refraction.
162    Altitude of the observer is not accounted for, which can be important when the sun 
163    is near the horizon. 
164    
165    Results are accurate to tenths of a degree, except where altitude is important
166    (< 20 degrees solar elevation)
167
168    Parameters
169    ----------
170    lat : float or ndarray
171        latitude in degrees
172    lon : float or ndarray
173        longitudes in degrees
174    datetime : datetime-like or str
175        date and time. Include time zone or UTC will be assumed
176    refraction : bool (default=False)
177        specifies whether to account for atmospheric refraction
178    temperature : float or ndarray (default=10)
179        surface atmospheric temperature (Celsius), only used for refraction calculation
180    pressure : float or ndarray (default=101325)
181        surface atmospheric pressure (Pa), only used for refraction calculation
182    solar_pos : `SolarPositionResults`, optional
183        solar position parameters from a prior call to `solar_position`
184    
185    Returns
186    -------
187    sza : float or ndarray
188        solar zenith angle in degrees at the designated locations and times
189        If refraction=False, this is the true solar zenith angle
190        If refraction=True, this is the apparent solar zenith angle
191    '''
192    # Convert to pandas Timestamp, if needed
193    datetime = _to_datetime(datetime)
194
195    # Solar declination, degrees
196    if solar_pos is None:
197        dec = solar_declination( datetime )
198    else:
199        dec = solar_pos.declination
200
201    # Hour angle, degrees
202    Ha = solar_hour_angle( lon, datetime, solar_pos )
203
204    # True solar zenith angle, radians
205    sza = np.arccos( np.sin(lat*pi180) * np.sin(dec*pi180) + \
206          np.cos(lat*pi180) * np.cos(dec*pi180) * np.cos(Ha*pi180) )
207
208    # Convert radians -> degrees
209    sza /= pi180
210
211    if refraction:
212        # Subtract refraction angle (degrees) from zenith angle.
213        # SZA is always smaller due to refraction.
214        sza -= refraction_angle( 90-sza, pressure, temperature )
215
216    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
  • datetime (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • refraction (bool (default=False)): specifies whether to account for atmospheric refraction
  • temperature (float or ndarray (default=10)): surface atmospheric temperature (Celsius), only used for refraction calculation
  • pressure (float or ndarray (default=101325)): surface atmospheric pressure (Pa), only used for refraction calculation
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
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 sunrise_time(*args, **kwargs): 
218def sunrise_time( *args, **kwargs ):
219    '''Compute sunrise time
220    
221    See `sun_times` for Parameters.'''
222    result = sun_times( *args, **kwargs )
223    return result[0]

Compute sunrise time

See sun_times for Parameters.

def sunset_time(*args, **kwargs): 
225def sunset_time( *args, **kwargs ):
226    '''Compute sunset time
227    
228    See `sun_times` for Parameters.'''
229    result = sun_times( *args, **kwargs )
230    return result[1]

Compute sunset time

See sun_times for Parameters.

def day_length(*args, **kwargs): 
232def day_length( *args, **kwargs ):
233    '''Compute length of daylight
234    
235    See `sun_times` for Parameters.'''
236    result = sun_times( *args, **kwargs )
237    return result[2]

Compute length of daylight

See sun_times for Parameters.

def solar_noon(*args, **kwargs): 
239def solar_noon( *args, **kwargs ):
240    '''Compute time of solar noon (meridian transit)
241    
242    See `sun_times` for Parameters.'''
243    result = sun_times( *args, **kwargs )
244    return result[3]

Compute time of solar noon (meridian transit)

See sun_times for Parameters.

def solar_constant(datetime, solar_pos=None): 
246def solar_constant( datetime, solar_pos=None ):
247    '''Compute solar constant for specific date or dates
248    
249    Parameters
250    ----------
251    datetime : datetime-like
252    solar_pos : `SolarPositionResults`, optional
253        solar position parameters from a prior call to `solar_position`
254
255    Returns
256    -------
257    S : float
258        Solar direct beam radiation flux density, W/m2
259    '''
260    if solar_pos is None:
261        solar_pos = solar_position( datetime )
262    S = 1361/solar_pos.distance**2
263
264    return S

Compute solar constant for specific date or dates

Parameters
Returns
  • S (float): Solar direct beam radiation flux density, W/m2
def solar_declination(datetime, fast=False, solar_pos=None): 
266def solar_declination( datetime, fast=False, solar_pos=None ):
267    '''Calculate solar declination (degrees) for specified date
268        
269    Parameters
270    ----------
271    datetime : datetime-like or str
272        date and time. Include time zone or UTC will be assumed
273    fast : bool (default=False)
274        Specifies using a faster but less accurate calculation
275    solar_pos : `SolarPositionResults`, optional
276        solar position parameters from a prior call to `solar_position`
277
278    Returns
279    -------
280    dec : float
281        solar declination in degrees at the specified date
282    '''
283    # Convert to pandas Timestamp, if needed
284    datetime = _to_datetime(datetime)
285
286    # Select the accurate or fast calculation
287    accurate = not fast
288
289    if solar_pos is not None:
290        dec = solar_pos.declination
291
292    else:
293        if accurate:
294
295            # Solar declination, degrees
296            dec, junk, junk, junk, junk = solar_position( datetime )
297
298        else:
299            # The fast method implements
300            # Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling
301
302            # Number of days since beginning of 2000
303            NJD = datetime - np.datetime64('2000-01-01')
304            try:
305                NJD = NJD.dt.days
306            except AttributeError:
307                NJD = NJD.days
308
309            # Obliquity, degrees
310            ob = 23.439 - 4e-7 * NJD
311
312            # Parameters for ecliptic, degrees
313            gm = 357.528 + 0.9856003 * NJD
314            lm = 280.460 + 0.9856474 * NJD
315
316            # Ecliptic longitude of sun, degrees
317            ec = lm + 1.915 * np.sin( gm * pi180 ) + 0.020 * np.sin( 2 * gm * pi180 )
318
319            #Solar declination, degrees
320            dec = np.arcsin( np.sin( ob * pi180 ) * np.sin( ec * pi180 ) ) / pi180
321
322    return dec

Calculate solar declination (degrees) for specified date

Parameters
  • datetime (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • fast (bool (default=False)): Specifies using a faster but less accurate calculation
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
Returns
  • dec (float): solar declination in degrees at the specified date
def solar_latitude(*args, **kwargs): 
324def solar_latitude( *args, **kwargs ):
325    '''Latitude of the subsolar point
326    
327    Parameters
328    ----------
329    datetime : datetime-like or str
330        date and time. Include time zone or UTC will be assumed
331    fast : bool
332        see `solar_declination`
333    solar_pos : `SolarPositionResults`, optional
334        solar position parameters from a prior call to `solar_position`
335        
336    Returns
337    -------
338    latitude : float
339        degrees of latitude
340    '''
341    return solar_declination( *args, **kwargs )

Latitude of the subsolar point

Parameters
Returns
  • latitude (float): degrees of latitude
def solar_longitude(datetime, solar_pos=None): 
343def solar_longitude( datetime, solar_pos=None ):
344    '''Longitude of the subsolar point, degrees
345    
346    Parameters
347    ----------
348    datetimeUTC : datetime-like or str
349        date and time. Include time zone or UTC will be assumed
350    solar_pos : `SolarPositionResults`, optional
351        solar position parameters from a prior call to `solar_position`
352    
353    Returns
354    -------
355    longitude : float
356        degrees of longitude
357    '''
358    # Convert to pandas Timestamp, if needed
359    datetimeUTC, tz_in = _to_datetime_utc(datetime)
360
361    # Longitude of subsolar point, degrees
362    # Equation of time will be added below
363    try:
364        # Treat as xarray.DataArray or pandas.Series
365        solar_lon = - 15 * ( datetimeUTC.dt.hour +
366                          datetimeUTC.dt.minute / 60 +
367                          datetimeUTC.dt.second / 3600 - 12 )
368    except AttributeError:
369        solar_lon = - 15 * ( datetimeUTC.hour +
370                          datetimeUTC.minute / 60 +
371                          datetimeUTC.second / 3600 - 12 )
372
373    # Add equation of time to the solar longitude, degrees
374    if solar_pos is None:
375        eot = equation_of_time( datetimeUTC, degrees=True )
376    else:
377        eot = solar_pos.equation_of_time
378    solar_lon -= eot
379
380    return solar_lon

Longitude of the subsolar point, degrees

Parameters
  • datetimeUTC (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
Returns
  • longitude (float): degrees of longitude
def solar_hour_angle(lon, datetime, solar_pos=None): 
382def solar_hour_angle( lon, datetime, solar_pos=None ):
383    '''Solar hour angle (degrees) for specified longitude, date and time
384
385    Hour angle is the angular displacement of the sun from the local meridian.
386    It is zero at local noon, negative in the morning, and positive is afternoon.
387    
388    Parameters
389    ----------
390    lon : float
391        longitude in degrees east
392    datetimeUTC : datetime-like or str
393        date and time. Include time zone or UTC will be assumed
394    solar_pos : `SolarPositionResults`, optional
395        solar position parameters from a prior call to `solar_position`
396    
397    Returns
398    -------
399    ha : float
400        hour angle in degrees at the specified location and time
401    '''
402
403    # Subsolar longitude, degrees
404    solar_lon = solar_longitude(datetime, solar_pos )
405
406    # Hour angle, degrees
407    Ha = lon - solar_lon
408
409    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 (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
Returns
  • ha (float): hour angle in degrees at the specified location and time
def equation_of_time(datetime, degrees=False, fast=False): 
411def equation_of_time( datetime, degrees=False, fast=False ):
412    '''Equation of time for specified date
413    
414    Accounts for the solar day being slightly different from 24 hours
415
416    Parameters
417    ----------
418    datetime : datetime-like or str
419        date and time. Include time zone or UTC will be assumed
420    degrees : bool (default=False)
421        If True, then return value in compass degrees
422        If False, then return value in minutes of an hour
423    fast : bool (default=False)
424        specifies whether to use a faster, but less accurate calculation
425        
426    Returns
427    -------
428    eot : float
429        equation of time on the specified date, degrees or minutes
430    '''
431    # Convert to pandas Timestamp, if needed
432    datetime = _to_datetime(datetime)
433
434    # Determine whether to use the fast or accurate calculation
435    accurate = not fast
436
437    if accurate:
438
439        # Equation of time, minutes
440        junk, junk, eot, junk, junk = solar_position( datetime )
441
442    else:
443        # Implements the "alternative equation" from Wikipedia, derived from
444        # 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
445        # When compared to the NREL SPA algorithm, differences reach are up to about 0.5 minute.
446        # Note: Leap years are not accounted for.
447
448        # Equation of time, accounts for the solar day differing slightly from 24 hr
449        try:
450            doy = datetime.dt.dayofyear
451        except AttributeError:
452            doy = datetime.dayofyear
453        W = 360 / 365.24
454        A = W * (doy+10)
455        B = A + 1.914 * np.sin( W * (doy-2) * pi180 )
456        C = ( A - np.arctan2( np.tan(B*pi180), np.cos(23.44*pi180) ) / pi180 ) / 180
457
458        # Equation of time in minutes of an hour (1440 minutes per day)
459        eot = 720 * ( C - np.round(C) )
460
461    # Equation of time, minutes -> degrees (360 degrees per day)
462    if degrees:
463        eot = eot / 60 * 360 / 24
464
465    return eot

Equation of time for specified date

Accounts for the solar day being slightly different from 24 hours

Parameters
  • datetime (datetime-like or str): date and time. Include time zone or UTC will be assumed
  • degrees (bool (default=False)): If True, then return value in compass degrees If False, then return value in minutes of an hour
  • fast (bool (default=False)): specifies whether to use a faster, but less accurate calculation
Returns
  • eot (float): equation of time on the specified date, degrees or minutes
def solar_position(datetime): 
467def solar_position( datetime ):
468    '''Compute position of sun (declination, right ascension, equation of time, distance)) on specified date
469    
470    Calculations follow the NOAA solar calculator spreadsheet
471    Applicable to years 1900-2100.
472
473    Parameters
474    ----------
475    date : datetime-like or str
476        date and time. Include time zone or UTC will be assumed
477    
478    Returns
479    -------
480    result : `SolarPositionResults`
481        declination, right ascension, equation of time, Earth-sun distance, datetime
482    '''
483    # Ensure time is Timestamp in UTC
484    datetimeUTC, tz_in = _to_datetime_utc(datetime)
485
486    # Raise warning if any dates are outside date range
487    # recommended for orbital parameters used here
488    if np.logical_or( np.any( datetimeUTC < np.datetime64('1900-01-01') ),
489                      np.any( datetimeUTC > np.datetime64('2100-01-01') ) ):
490        warnings.warn('Solar position accuracy declines for dates outside 1900-2100', \
491                      RuntimeWarning )
492
493    # Number of days since 1 Jan 2000
494    NJD = datetimeUTC - np.datetime64('2000-01-01')
495    try:
496        NJD = NJD.dt.days \
497            + NJD.dt.seconds / 86400.
498    except AttributeError:
499        NJD = NJD.days \
500            + NJD.seconds / 86400.
501
502    # Julian day (since 12:00UTC 1 Jan 4713 BCE)
503    NJD += 2451544.50
504
505    # Julian century
506    JC = (NJD-2451545)/36525
507
508    # Earth orbital eccentricity, unitless
509    ec = 0.016708634 - JC*( 0.000042037 + 0.0000001267*JC )
510
511    # Earth mean orbital obliquity, degree
512    mean_ob = 23 + ( 26 + ( (21.448
513                             - JC * (46.815
514                                    + JC * (0.00059 - JC * 0.001813) ) ) )/60 )/60
515
516    # Earth true orbital obliquity, corrected for nutation, degree 
517    ob = mean_ob + 0.00256 * np.cos( (125.04 - 1934.136*JC ) * pi180 )
518
519    # Sun Mean ecliptic longitude, degree
520    mean_ec_lon = np.mod( 280.46646 + JC*( 36000.76983 + JC*0.0003032 ), 360 )
521
522    # Sun Mean anomaly, degree
523    mean_anom = 357.52911 + JC*( 35999.05029 - 0.0001537*JC )
524
525    # Sun Equation of center, degree
526    eq_center = np.sin(mean_anom*pi180) * (1.914602 - JC*( 0.004817 + 0.000014*JC )) \
527                    + np.sin(2*mean_anom*pi180) * (0.019993 - 0.000101*JC) \
528                    + np.sin(3*mean_anom*pi180) * 0.000289
529
530    # Sun True ecliptic longitude, degrees
531    true_ec_lon = mean_ec_lon + eq_center
532
533    # Sun True anomaly, degree
534    true_anom = mean_anom + eq_center
535
536    # Earth-Sun distance, AU
537    distance = (1.000001018 * (1-ec**2) ) / (1 + ec * np.cos( true_anom * pi180 ))
538
539    # Sun Apparent ecliptic longitude, corrected for nutation, degrees
540    ec_lon = true_ec_lon - 0.00569 - 0.00478 * np.sin( (125.04 - 1934.136*JC ) * pi180)
541
542    # Sun Right ascension, deg
543    right_ascension = np.arctan2( np.cos(ob*pi180) * np.sin(ec_lon*pi180),
544                                  np.cos(ec_lon*pi180) ) / pi180
545
546    # Sun Declination, deg
547    declination = np.arcsin( np.sin(ob*pi180) * np.sin(ec_lon*pi180) ) / pi180
548
549    # var y
550    vary = np.tan( ob/2 * pi180 )**2
551
552    # Equation of time, minutes
553    eot = vary * np.sin( 2 * mean_ec_lon * pi180) \
554        - 2 * ec * np.sin( mean_anom * pi180 ) \
555        + 4 * ec * vary * np.sin( mean_anom * pi180 ) * np.cos( 2 * mean_ec_lon * pi180) \
556        - 0.5 * vary**2 * np.sin( 4 * mean_ec_lon * pi180) \
557        - 1.25 * ec**2 * np.sin( 2 * mean_anom * pi180)
558    eot = eot * 4 / pi180
559
560    # Bundle results
561    result = SolarPositionResults(declination,
562                                  right_ascension,
563                                  eot,
564                                  distance,
565                                  datetimeUTC)
566
567    return result

Compute position of sun (declination, right ascension, equation of time, distance)) on specified date

Calculations follow the NOAA solar calculator spreadsheet Applicable to years 1900-2100.

Parameters
  • date (datetime-like or str): date and time. Include time zone or UTC will be assumed
Returns
  • result (SolarPositionResults): declination, right ascension, equation of time, Earth-sun distance, datetime
def sun_times( lat, lon, datetime, tz_out=None, sza_sunrise=90.833, fast=False, solar_pos=None): 
569def sun_times( lat, lon, datetime, tz_out=None, sza_sunrise=90.833,
570              fast=False, solar_pos=None ):
571    '''Compute times of sunrise, sunset, solar noon, and day length
572    
573    Common options for solar zenith angle at sunrise
574    1. 90.833 for first edge of sun rising, typical (0.567°) refraction (default)
575    2. 90.267 for first edge of sun rising, no refraction
576    3. 90 degrees for center of sun rising, no refraction
577
578    Note: horizon_zenith_angle can be used to compute a more accurate horizon location
579    for sites at elevation.
580    
581    Parameters
582    ----------
583    lat : float or ndarray
584        latitude in degrees
585    lon : float or ndarray
586        longitudes in degrees
587    datetime : datetime-like or str
588        datetime, provide a time zone or UTC will be assumed 
589    tz_out : str, pytz.timezone, datetime.tzinfo, optional
590        timezone to be used for output times. 
591        If None is provided, then result will be in same time zone as input or UTC
592    sza_sunrise : float (default=90.833)
593        Solar zenith angle at which sunrise and sunset are calculated, degrees
594    fast : bool (default=False)
595        Select a faster but less accurate calculation
596    solar_pos : `SolarPositionResults`, optional
597        solar position parameters from a prior call to `solar_position`
598
599    Returns
600    -------
601    result : `SunTimesResults`
602        times of sunrise, sunset, solar noon, and day length
603    '''
604    # Convert to pandas Timestamp in UTC, if needed
605    datetimeUTC, tz_in = _to_datetime_utc(datetime)
606
607    # If no output timezone is specified, use the input time zone
608    if tz_out is None:
609        tz_out = tz_in
610
611    # Select fast or accurate calculation
612    accurate = not fast
613
614    # Solar declination (degrees) and equation of time (minutes)
615    if solar_pos is not None:
616        dec = solar_pos.declination
617        eot = solar_pos.eot
618    else:
619        if accurate:
620            dec, junk, eot, junk, junk = solar_position( datetimeUTC )
621        else:
622            dec = solar_declination( datetimeUTC )
623            eot = equation_of_time( datetimeUTC )
624
625    # Sunrise hour angle, degree
626    # Degrees east of the local meridian where sun rises
627    ha_sunrise = np.arccos( np.cos(sza_sunrise*pi180) /
628                           (np.cos(lat*pi180) * np.cos(dec*pi180))
629                           - np.tan(lat*pi180)*np.tan(dec*pi180) ) / pi180
630
631    # Solar noon, local standard time, day fraction
632    solar_noon = (720 - 4*lon - eot ) / 1440
633
634    # Sunrise and sunset, local standard time, day fraction
635    t_sunrise = solar_noon - 4 * ha_sunrise / 1440
636    t_sunset  = solar_noon + 4 * ha_sunrise / 1440
637
638    # Midnight UTC
639    # datetimeUTC is in UTC but time-zone-naive
640    try:
641        # Series time objects
642        dateUTC = datetimeUTC.dt.floor('D')
643    except AttributeError:
644        # Scalar time objects
645        dateUTC = datetimeUTC.floor('D')
646
647    # Convert day fraction -> date time
648    solar_noon = dateUTC + solar_noon * pd.Timedelta( 1, 'day' )
649    t_sunrise  = dateUTC + t_sunrise  * pd.Timedelta( 1, 'day' )
650    t_sunset   = dateUTC + t_sunset   * pd.Timedelta( 1, 'day' )
651
652    # Convert to output timezone, if any is provided
653    if tz_out is not None:
654        if isinstance(solar_noon,(xr.DataArray,np.ndarray)) or \
655            isinstance(t_sunrise,(xr.DataArray,np.ndarray)):
656            # These types don't localize tz, but we can add offset to the tz-naive time
657            if hasattr(datetimeUTC,'tz_localize'):
658                # For scalar datetime, there is a single time offset, which we can add
659                utcoffset = np.timedelta64( datetimeUTC.tz_localize(tz_out).utcoffset() )
660                solar_noon += utcoffset
661                t_sunrise  += utcoffset
662                t_sunset   += utcoffset
663            else:
664                # For Series datetime, there are potentially multiple offsets. We can only add one
665                unique_datetimeUTC = pd.DatetimeIndex(np.unique(datetimeUTC))
666                unique_utcoffsets = np.unique( unique_datetimeUTC.tz_localize('UTC') \
667                                        - unique_datetimeUTC.tz_localize(tz_out) )
668                if len(unique_utcoffsets)==1:
669                    utcoffset = unique_utcoffsets[0]
670                    solar_noon += utcoffset
671                    t_sunrise  += utcoffset
672                    t_sunset   += utcoffset
673                else:
674                    # We might be able to handle multiple offsets if we can
675                    # determine which dimension of `solar_noon` and `t_sunrise`
676                    # is the time dimension.
677                    raise ValueError('Multiple timezone offsets not supported. '
678                                     +'Request output in UTC or reduce number of input times.')
679        else:
680            try:
681                # Series time objects
682                solar_noon = solar_noon.dt.tz_localize('UTC')\
683                                    .dt.tz_convert(tz_out)
684                t_sunrise  = t_sunrise.dt.tz_localize('UTC')\
685                                    .dt.tz_convert(tz_out)
686                t_sunset   = t_sunset.dt.tz_localize('UTC')\
687                                    .dt.tz_convert(tz_out)
688            except AttributeError:
689                # Scale time objects
690                solar_noon = solar_noon.tz_localize('UTC')\
691                                    .tz_convert(tz_out)
692                t_sunrise  = t_sunrise.tz_localize('UTC')\
693                                    .tz_convert(tz_out)
694                t_sunset   = t_sunset.tz_localize('UTC')\
695                                    .tz_convert(tz_out)
696
697    # Sunlight duration, minutes
698    day_length = 8 * ha_sunrise * pd.Timedelta(1, 'minute')
699
700    result = SunTimesResults(t_sunrise, t_sunset, day_length, solar_noon, datetimeUTC)
701
702    return result

Compute times of sunrise, sunset, solar noon, and day length

Common options for solar zenith angle at sunrise

  1. 90.833 for first edge of sun rising, typical (0.567°) refraction (default)
  2. 90.267 for first edge of sun rising, no refraction
  3. 90 degrees for center of sun rising, no refraction

Note: horizon_zenith_angle can be used to compute a more accurate horizon location for sites at elevation.

Parameters
  • lat (float or ndarray): latitude in degrees
  • lon (float or ndarray): longitudes in degrees
  • datetime (datetime-like or str): datetime, provide a time zone or UTC will be assumed
  • tz_out (str, pytz.timezone, datetime.tzinfo, optional): timezone to be used for output times. If None is provided, then result will be in same time zone as input or UTC
  • sza_sunrise (float (default=90.833)): Solar zenith angle at which sunrise and sunset are calculated, degrees
  • fast (bool (default=False)): Select a faster but less accurate calculation
  • solar_pos (SolarPositionResults, optional): solar position parameters from a prior call to solar_position
Returns
  • result (SunTimesResults): times of sunrise, sunset, solar noon, and day length
def horizon_zenith_angle(alt): 
704def horizon_zenith_angle( alt ):
705    '''Angle from the zenith to the horizon
706    
707    The horizon is the locii of points where a line from the 
708    observation location to the ellipsoid is tangent to the ellipsoid surface.
709    
710    The altitude parameter should be the vertical distance 
711    above the surrounding terrain that defines the horizon,
712    not necessarily the altitude above sea level or the altitude above ground level.
713    For example, on a mountain peak that is 4000 m above sea level and 
714    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
715    For an observer on the plateau, the relevant altitude is 0 m.
716
717    The implementation below assumes a spherical Earth.
718    Results using the WGS84 ellipsoid (see commented code below)
719    differ from the spherical case by << 1°. Terrain,
720    which is neglected here, has a larger effect on the horizon
721    location, so the simpler spherical calculation is appropriate. 
722
723    Parameters
724    ----------
725    lat : float or ndarray
726        latitude in degrees
727    alt : float or ndarray
728        altitude above surrounding terrain that defines the horizon, meters
729        
730    Returns
731    -------
732    hza : float or ndarray
733        horizon zenith angle in degrees
734    '''
735
736    # WGS84 ellipsoid parameters
737    # semi-major radius, m
738    r_earth = 6378137.0
739    # ellipsoidal flattening, unitless
740    f = 1/298.257223563
741
742    # Horizon zenith angle, degrees (spherical earth)
743    hza = 180 - np.arcsin( r_earth / ( r_earth + alt ) ) / pi180
744
745    ## Ellipsoidal Earth
746    # # Eccentricity of ellipsoid
747    # ecc = f * (2-f)
748    # # Local (i.e. prime vertical) radius of curvature at latitude
749    # N = r_earth / np.sqrt( 1 - ecc**2 * np.sin(lat*pi180)**2 )
750    # # Horizon zenith angle, degrees
751    # hza = 180 - np.arcsin( N / (N+alt) ) / pi180
752
753    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 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.

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 refraction_angle(true_elevation_angle, pressure=101325.0, temperature_celsius=10.0): 
755def refraction_angle( true_elevation_angle, pressure=101325., temperature_celsius=10. ):
756    '''Atmospheric refraction angle for light passing through Earth's atmosphere
757
758    The apparent locations in the sky of objects outsides Earth's atmosphere 
759    differs from their true locations due to atmospheric refraction. 
760    (e.g. sun and moon when they rise and set)
761    The apparent elevation of an object above the horizon is
762    apparent elevation angle = (true elevation angle) + (refraction angle)
763    
764    The equations here are from Saemundsson/Bennett, whose calculations use
765    a typical vertical profile of atmospheric density (i.e. temperature and pressure).
766    The profiles can be rescaled to a particular surface temperature and pressure
767    to approximately account for varying atmospheric conditions.
768    Accurate refraction calculations should use fully specified vertical profile
769    of temperature and pressure, which cannot be done here.
770
771    Parameters
772    ----------
773    true_elevation_angle : float
774        degrees above horizon of sun or other object
775    pressure : float (default=101325)
776        surface atmospheric pressure (Pa)
777    temperature_celsius : float (default=10)
778        surface atmospheric temperature (C)
779
780    Returns
781    -------
782    angle : float
783        refraction angle in degrees. Value is zero when apparent elevation is below horizon
784    '''
785    # Refraction angle, arcminutes
786    R = 1.02 / np.tan( ( true_elevation_angle + 10.3 / (true_elevation_angle + 5.11) ) * pi180 )
787    # Account for temperature and pressure, arcminutes
788    R = R * pressure / 101325 * 283 / ( 273 + temperature_celsius )
789    # Convert arcminutes -> degrees
790    R /= 60
791
792    # Result must be positive
793    R = np.maximum(R,0)
794
795    # Refraction defined only when the apparent elevation angle is positive
796    # Set refraction to zero when the apparent elevation is below horizon
797    refraction_angle = np.where( true_elevation_angle + R <= 0, 0, R)
798
799    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
class SolarPositionResults(builtins.tuple):

Namedtuple containing results of solar_position

All values are geocentric, strictly accurate for the center of the Earth, not a point on Earth's surface. The parallax angle from Earth's center to surface is 4e-5 degrees.

Attributes
  • declination (float): position of the sun relative to Earth's equatorial plane, degrees
  • right_ascension (float): position of the sun along Earth's equatorial plane relative to the vernal equinox, degrees
  • equation_of_time (float): equation of time (minutes) between mean solar time and true solar time. Divide by 4 minutes per degree to obtain equation of time in degrees.
  • distance (float): Earth-sun distance in AU (1 AU = 1.495978707e11 m)
  • datetimeUTC (Timestamp, Series, or DataArray of numpy.datetime64): date and time input for calculations, UTC
SolarPositionResults( declination, right_ascension, equation_of_time, distance, datetimeUTC)

Create new instance of SolarPositionResults(declination, right_ascension, equation_of_time, distance, datetimeUTC)

declination

Alias for field number 0

right_ascension

Alias for field number 1

equation_of_time

Alias for field number 2

distance

Alias for field number 3

datetimeUTC

Alias for field number 4

Inherited Members
builtins.tuple
index
count
class SunTimesResults(builtins.tuple):

Namedtuple containing results of sun_times

Attributes
  • sunrise (pandas.DatetimeIndex): sunrise time, UTC if not specified otherwise
  • sunset (pandas.DatetimeIndex): sunset time, UTC if not specified otherwise
  • day_length (pandas.Timedelta): duration of daylight
  • solar_noon (pandas.DatetimeIndex): time of meridian transit, UTC if not specified otherwise
  • datetimeUTC (Timestamp, Series, or DataArray of numpy.datetime64): date and time input for calculations, UTC
SunTimesResults(sunrise, sunset, day_length, solar_noon, datetimeUTC)

Create new instance of SunTimesResults(sunrise, sunset, day_length, solar_noon, datetimeUTC)

sunrise

Alias for field number 0

sunset

Alias for field number 1

day_length

Alias for field number 2

solar_noon

Alias for field number 3

datetimeUTC

Alias for field number 4

Inherited Members
builtins.tuple
index
count
pi180 = 0.017453292519943295

Constant $\pi/180$

def sza( lat, lon, datetime, refraction=False, temperature=10.0, pressure=101325.0, solar_pos=None): 
156def solar_zenith_angle( lat, lon, datetime,
157                        refraction=False, temperature=10., pressure=101325.,
158                        solar_pos=None ):
159    '''Solar zenith angle (degrees) for a given latitude, longitude, date and time.
160    
161    Accounts for equation of time and (optionally) for atmospheric refraction.
162    Altitude of the observer is not accounted for, which can be important when the sun 
163    is near the horizon. 
164    
165    Results are accurate to tenths of a degree, except where altitude is important
166    (< 20 degrees solar elevation)
167
168    Parameters
169    ----------
170    lat : float or ndarray
171        latitude in degrees
172    lon : float or ndarray
173        longitudes in degrees
174    datetime : datetime-like or str
175        date and time. Include time zone or UTC will be assumed
176    refraction : bool (default=False)
177        specifies whether to account for atmospheric refraction
178    temperature : float or ndarray (default=10)
179        surface atmospheric temperature (Celsius), only used for refraction calculation
180    pressure : float or ndarray (default=101325)
181        surface atmospheric pressure (Pa), only used for refraction calculation
182    solar_pos : `SolarPositionResults`, optional
183        solar position parameters from a prior call to `solar_position`
184    
185    Returns
186    -------
187    sza : float or ndarray
188        solar zenith angle in degrees at the designated locations and times
189        If refraction=False, this is the true solar zenith angle
190        If refraction=True, this is the apparent solar zenith angle
191    '''
192    # Convert to pandas Timestamp, if needed
193    datetime = _to_datetime(datetime)
194
195    # Solar declination, degrees
196    if solar_pos is None:
197        dec = solar_declination( datetime )
198    else:
199        dec = solar_pos.declination
200
201    # Hour angle, degrees
202    Ha = solar_hour_angle( lon, datetime, solar_pos )
203
204    # True solar zenith angle, radians
205    sza = np.arccos( np.sin(lat*pi180) * np.sin(dec*pi180) + \
206          np.cos(lat*pi180) * np.cos(dec*pi180) * np.cos(Ha*pi180) )
207
208    # Convert radians -> degrees
209    sza /= pi180
210
211    if refraction:
212        # Subtract refraction angle (degrees) from zenith angle.
213        # SZA is always smaller due to refraction.
214        sza -= refraction_angle( 90-sza, pressure, temperature )
215
216    return sza

Alias for solar_zenith_angle

def saa(lat, lon, datetime, solar_pos=None): 
 62def solar_azimuth_angle( lat, lon, datetime, solar_pos=None ):
 63    '''Solar azimuth angle (degrees) for a latitude, longitude, date and time
 64    
 65    SAA is degrees clockwise from north.
 66    
 67    Parameters
 68    ----------
 69    lat : float or ndarray
 70        latitude in degrees
 71    lon : float or ndarray
 72        longitudes in degrees
 73    datetime : datetime-like or str
 74        date and time. Include time zone or UTC will be assumed
 75    solar_pos : `SolarPositionResults`, optional
 76        solar position parameters from a prior call to `solar_position`
 77
 78    Returns
 79    -------
 80    saa : float or ndarray
 81        solar azimuth angle in degrees (clockwise from north)
 82    '''
 83    # Convert to pandas Timestamp, if needed
 84    datetime = _to_datetime(datetime)
 85
 86    # Subsolar point, latitude longitude, degrees
 87    solar_lat = solar_latitude( datetime, solar_pos=solar_pos )
 88    solar_lon = solar_longitude( datetime, solar_pos=solar_pos )
 89
 90    # Vector pointing toward sun
 91    x = np.cos( solar_lat * pi180 ) * np.sin( (solar_lon - lon) * pi180 )
 92    y = np.cos( lat*pi180 ) * np.sin( solar_lat*pi180 ) \
 93        - np.sin( lat*pi180 ) * np.cos( solar_lat*pi180 ) \
 94            * np.cos( (solar_lon - lon) * pi180 )
 95
 96    # Azimuth angle from north, degrees
 97    saa = np.arctan2( x, y ) / pi180
 98
 99    # Change range [-180,180] to [0,360]
100    return np.mod( saa+360, 360 )

Alias for solar_azimuth_angle

def sea( lat, lon, datetime, alt=0, refraction=False, temperature=10.0, pressure=101325.0, solar_pos=None): 
102def solar_elevation_angle( lat, lon, datetime, alt=0,
103                        refraction=False, temperature=10., pressure=101325.,
104                        solar_pos=None ):
105    '''Solar elevation angle (degrees) above the horizon
106
107    The altitude parameter should be the vertical distance 
108    above the surrounding terrain that defines the horizon,
109    not necessarily the altitude above sea level or the altitude above ground level.
110    For example, on a mountain peak that is 4000 m above sea level and 
111    1500 m above the surrounding plateau, the relevant altitude is 1500 m.
112    For an observer on the plateau, the relevant altitude is 0 m.
113
114    See documentation for `solar_zenith_angle` and `horizon_zenith_angle`.
115
116    Parameters
117    ----------
118    lat : float or ndarray
119        latitude in degrees
120    lon : float or ndarray
121        longitudes in degrees
122    datetime : datetime-like or str
123        date and time. Include time zone or UTC will be assumed
124    alt : float or ndarray (default=0)
125        altitude above surrounding terrain that defines the horizon, meters
126    refraction : bool (default=False)
127        specifies whether to account for atmospheric refraction
128    temperature : float or ndarray (default=10)
129        surface atmospheric temperature (Celsius), only used for refraction calculation
130    pressure : float or ndarray (default=101325)
131        surface atmospheric pressure (Pa), only used for refraction calculation
132    solar_pos : `SolarPositionResults`, optional
133        solar position parameters from a prior call to `solar_position`
134    
135    Returns
136    -------
137    sea : float or ndarray
138        solar elevation angle in degrees at the designated locations and times
139        If refraction=False, this is the true solar elevation angle
140        If refraction=True, this is the apparent solar elevation angle
141    
142    '''
143
144    if refraction and np.any(alt):
145        warnings.warn( 'Atmospheric refraction is calculated for surface conditions, '
146                    + 'but an altitude above the surface was specified',
147                     category=UserWarning,
148                     stacklevel=2 )
149
150    sea = horizon_zenith_angle( alt ) \
151         - solar_zenith_angle( lat, lon, datetime, refraction, temperature, pressure, 
152                               solar_pos=solar_pos )
153
154    return sea

Alias for solar_elevation_angle

def equationOfTime(datetime, degrees=False, fast=False): 
411def equation_of_time( datetime, degrees=False, fast=False ):
412    '''Equation of time for specified date
413    
414    Accounts for the solar day being slightly different from 24 hours
415
416    Parameters
417    ----------
418    datetime : datetime-like or str
419        date and time. Include time zone or UTC will be assumed
420    degrees : bool (default=False)
421        If True, then return value in compass degrees
422        If False, then return value in minutes of an hour
423    fast : bool (default=False)
424        specifies whether to use a faster, but less accurate calculation
425        
426    Returns
427    -------
428    eot : float
429        equation of time on the specified date, degrees or minutes
430    '''
431    # Convert to pandas Timestamp, if needed
432    datetime = _to_datetime(datetime)
433
434    # Determine whether to use the fast or accurate calculation
435    accurate = not fast
436
437    if accurate:
438
439        # Equation of time, minutes
440        junk, junk, eot, junk, junk = solar_position( datetime )
441
442    else:
443        # Implements the "alternative equation" from Wikipedia, derived from
444        # 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
445        # When compared to the NREL SPA algorithm, differences reach are up to about 0.5 minute.
446        # Note: Leap years are not accounted for.
447
448        # Equation of time, accounts for the solar day differing slightly from 24 hr
449        try:
450            doy = datetime.dt.dayofyear
451        except AttributeError:
452            doy = datetime.dayofyear
453        W = 360 / 365.24
454        A = W * (doy+10)
455        B = A + 1.914 * np.sin( W * (doy-2) * pi180 )
456        C = ( A - np.arctan2( np.tan(B*pi180), np.cos(23.44*pi180) ) / pi180 ) / 180
457
458        # Equation of time in minutes of an hour (1440 minutes per day)
459        eot = 720 * ( C - np.round(C) )
460
461    # Equation of time, minutes -> degrees (360 degrees per day)
462    if degrees:
463        eot = eot / 60 * 360 / 24
464
465    return eot

Alias for equation_of_time

def solarDeclination(datetime, fast=False, solar_pos=None): 
266def solar_declination( datetime, fast=False, solar_pos=None ):
267    '''Calculate solar declination (degrees) for specified date
268        
269    Parameters
270    ----------
271    datetime : datetime-like or str
272        date and time. Include time zone or UTC will be assumed
273    fast : bool (default=False)
274        Specifies using a faster but less accurate calculation
275    solar_pos : `SolarPositionResults`, optional
276        solar position parameters from a prior call to `solar_position`
277
278    Returns
279    -------
280    dec : float
281        solar declination in degrees at the specified date
282    '''
283    # Convert to pandas Timestamp, if needed
284    datetime = _to_datetime(datetime)
285
286    # Select the accurate or fast calculation
287    accurate = not fast
288
289    if solar_pos is not None:
290        dec = solar_pos.declination
291
292    else:
293        if accurate:
294
295            # Solar declination, degrees
296            dec, junk, junk, junk, junk = solar_position( datetime )
297
298        else:
299            # The fast method implements
300            # Eq. 9.68-9.72 from M.Z. Jacobson, Fundamentals of Atmospheric Modeling
301
302            # Number of days since beginning of 2000
303            NJD = datetime - np.datetime64('2000-01-01')
304            try:
305                NJD = NJD.dt.days
306            except AttributeError:
307                NJD = NJD.days
308
309            # Obliquity, degrees
310            ob = 23.439 - 4e-7 * NJD
311
312            # Parameters for ecliptic, degrees
313            gm = 357.528 + 0.9856003 * NJD
314            lm = 280.460 + 0.9856474 * NJD
315
316            # Ecliptic longitude of sun, degrees
317            ec = lm + 1.915 * np.sin( gm * pi180 ) + 0.020 * np.sin( 2 * gm * pi180 )
318
319            #Solar declination, degrees
320            dec = np.arcsin( np.sin( ob * pi180 ) * np.sin( ec * pi180 ) ) / pi180
321
322    return dec

Alias for solar_declination