add sunrise / sunset to Location

[dacoex]

In reference to changes in #93:

add the following functions to :py:class:Location``:

get_sunrise and get_sunset and deferred from these get_dayduration.

This helps filtering values, e.g. set night values to np.nan.

[wholmgren]

Probably a good idea, though it's not likely to happen as part of that PR. Are you aware of the get_sun_rise_set_transit function? It returns both sunrise and sunset. I'm less interested in get_dayduration since I think it's less likely to be used and it can easily be calculated if needed (though not totally opposed to it).

All that being said, I've found that filtering using zenith or elevation angle is usually a better approach.

[cwhanse]

I'm in need of this functionality, and am looking for opinions about implementation. I'd prefer methods attached to a PVsystem or Location object that mimic pyephem, e.g., Location.get_sunrise(date), Location.get_next_sunrise(datetime).

Currently get_sun_rise_set_transit runs solarposition.spa and returns sunrise, sunset and transit. #316 describes some odd behavior of get_sun_rise_set_transit and ideas about fixing it. Alternatively, we could request sunrise etc. from solarposition.spa (which is what get_sun_rise_set_transit is doing) or use solarposition.pyphem by adding a call to the appropriate method of the Sun object.

I'd like to avoid a workflow that calculates ephemeris twice (once to get solar position, and again to get sunrise and sunset), which is what happens with the current get_sun_rise_set_transit.

Thoughts?

[mikofski]

couple of ideas

1. when I fixed the timezone bug in spa, I changed the SPA function code to SPA_ALL which calculates sunrise, sun-transit, and sunset times, so if it's not implemented in the numpy version of spa, you could use the NREL C version

2. I think sunrise/set is a perfect use for the analytical zenith function, just find where it equals zero, either using a solver, or by solving it analytically:

   no.pi/2 = np.arccos(
       np.cos(declination) * np.cos(latitude) * np.cos(hour_angle) +
       np.sin(declination) * np.sin(latitude)
   )
   0 = (np.cos(declination) * np.cos(latitude) * np.cos(hour_angle) +
        np.sin(declination) * np.sin(latitude))
   

   so

   np.cos(hour_angle) = ((0 - np.sin(declination)*np.sin(latitude)) /
                         np.cos(declination)*np.cos(latitude))
   

or

np.cos(hour_angle) = -np.tan(declination)*np.tan(latitude)

Is that right?

[wholmgren]

Location still seems like an appropriate place for this functionality.

Sorry if I'm picking at a typo, but so we're all on the same page...

• there is no solarposition.spa.
• pvlib.solarposition.spa_python() and pvlib.solarposition.get_sun_rise_set_transit() each import pvlib.spa by calling pvlib.solarposition._spa_python_import(). This allows for optional numba support that would not be possible using the standard imports at the top of pvlib.solarposition.
• pvlib.solarposition.get_sun_rise_set_transit is a wrapper around pvlib.spa.transit_sunrise_sunset.

I'd like to avoid a workflow that calculates ephemeris twice (once to get solar position, and again to get sunrise and sunset), which is what happens with the current get_sun_rise_set_transit.

I'm not sure how avoid two different functions because it seems that they do two different things. A function to calculate sunrise/sunset should accept a date and return a collection of datetimes. A function to calculate solar position should accept a datetime and return a collection of angles. I'm having a hard time thinking outside of this box.

In any case, the get_sun_rise_set_transit function could certainly be modified to use either pyephem or pvlib.spa or an analytical solution (controllable by a kwarg). We'd want to make the API consistent for all of the methods.

[mikofski]

Silly me, sunset/rise is at zenith of pi/2 not zero, comments above corrected

[markcampanelli]

Assuming it holds water, I'm curious how much more performant the above analytical solution might be. Perhaps it could be added as one of the how options to get_sun_rise_set_transit() as part of the proposed addition the to the Location API?

[mikofski]

i guess it seems to work, if I compare it to the online spa calculator I get nearly the correct sunrise/sunset times:

from datetime import datetime, timedelta
import numpy as np
import pytz
import pandas as pd
from pvlib.solarposition import (
    equation_of_time_pvcdrom, declination_spencer71)

dt = pd.DatetimeIndex(start='2018-01-01 0:00:00',
                      end='2018-12-31 23:59:59',
                      freq='H').tz_localize('Etc/GMT+7')

lat, lon = 39.743, -105.178        
latrad, lonrad = np.radians(lat), np.radians(lon)  # radians

eoq = equation_of_time_pvcdrom(dt.dayofyear)  # minutes
decl = declination_spencer71(dt.dayofyear)  # radians


def get_sunset_angle(decl, latrad):
    return np.arccos(-np.tan(decl) * np.tan(latrad))


# solar noon is at hour angle zero
# so sunrise is just negative of sunset

ss_angrad = get_sunset_angle(decl, latrad) # radians
ss_ang = np.degrees(ss_angrad)

sr_angrad = -ss_angrad # radians
sr_ang = -ss_ang


def local_hours(times, hour_ang, lon, eoq):
    timezone = times.tz.utcoffset(times).total_seconds() / 3600.
    return (hour_ang - lon - eoq/4.)/15. + 12. + timezone


sslh = local_hours(dt, ss_ang, lon, eoq)
srlh = local_hours(dt, sr_ang, lon, eoq)

print(sslh[0], srlh[0])  # analytical hours
#16.69168401501784 7.455555262428694

print(timedelta(hours=sslh[0]) +
    pytz.timezone('Etc/GMT+7').localize(datetime(2018, 1, 1, 0, 0, 0)))
#2018-01-01 16:41:30.062454-07:00

print(timedelta(hours=srlh[0]) +
    pytz.timezone('Etc/GMT+7').localize(datetime(2018, 1, 1, 0, 0, 0)))
#2018-01-01 07:27:19.998945-07:00

#from NREL MIDC SPA at (39.743, -105.178)
#Date,Time,sunrise,transit,sunset
#1/1/2018,0:00:00,7.364246,12.073578,16.785242

from pvlib.spa_c_files.spa_py import spa_calc

result = spa_calc(
    year=2018, month=1, day=1, hour=0, minute=0, second=0,
    time_zone=-7, longitude=-105.178, latitude=39.743,
    elevation=1830.14, pressure=820, temperature=11, delta_t=67)

print(result['sunset'], result['sunrise'])  # SPA hours
#16.785243982117017 7.364247922706669

print(timedelta(hours=16.785242) +
    pytz.timezone('Etc/GMT+7').localize(datetime(2018, 1, 1, 0, 0, 0)))
#2018-01-01 16:47:06.871200-07:00

print(timedelta(hours=7.364246) +
    pytz.timezone('Etc/GMT+7').localize(datetime(2018, 1, 1, 0, 0, 0)))
#2018-01-01 07:21:51.285600-07:00

It's about 6 minute difference 😦

[mikofski]

Maybe I should've tried the other declination or eoq formulas?

[adriesse]

Your formula matches the one in my Vignola book. Maybe the SPA version gives you apparent zenith == 90?

On the rare occasions when I need sunrise/sunset times I usually have a time series of some sort I'm working with, and just interpolate the timestamp for apparent_zenith == 90. But most of the time I'm happy filtering on apparent_zenith > 90.

[cwhanse]

Sorry if I'm picking at a typo, but so we're all on the same page...

You're not, and thanks.

pvlib.solarposition.spa_python() and pvlib.solarposition.get_sun_rise_set_transit() each import pvlib.spa by calling pvlib.solarposition._spa_python_import(). This allows for optional numba support that would not be possible using the standard imports at the top of pvlib.solarposition.

This is important and I had missed this.

pvlib.solarposition.get_sun_rise_set_transit is a wrapper around pvlib.spa.transit_sunrise_sunset.

Right. Now that I look at the spa module it is clear that we won't calculate the ephemeris twice (which I had in mind when I was considering an earlier suggestion that sunrise/sunset are based on zenith, then mistakenly transposed that idea to the spa and pyephem functions).

I'm not sure how avoid two different functions because it seems that they do two different things. A function to calculate sunrise/sunset should accept a date and return a collection of datetimes. A function to calculate solar position should accept a datetime and return a collection of angles. I'm having a hard time thinking outside of this box.

We're in the same box.

[mikofski]

@adriesse I knew the analytical functions don't account for atmospheric refraction, I just didn't realize it would be so large. Visibly the largest difference in zenith to me appears at noon, but it makes sense that refraction would be larger when the angle is greatest, ie sunrise/sunset. Always learning, thanks! @cwhanse glad you can get what you need from spa