Add a new and faster way to calculate the solarposition #730

meteoDaniel · 2019-05-21T07:42:08Z

pvlib python pull request guidelines

I am familiar with the contributing guidelines.
Fully tested. Added and/or modified tests to ensure correct behavior for all reasonable inputs. Tests (usually) must pass on the TravisCI and Appveyor testing services.
Updates entries to docs/sphinx/source/api.rst for API changes.
Adds description and name entries in the appropriate docs/sphinx/source/whatsnew file for all changes.
Code quality and style is sufficient. Passes LGTM and SticklerCI checks.
New code is fully documented. Includes sphinx/numpydoc compliant docstrings and comments in the code where necessary.
Pull request is nearly complete and ready for detailed review.

Brief description of the problem and proposed solution (if not already fully described in the issue linked to above):
A new way to calculate solar position. Let us benchmark this against nrel_numba. A first test against the default configuration shows the code runs nearly 100 times faster:
1000 timesteps:
Default: 0:00:00.023037
spencer_mc: 0:00:00.003093
10000 timesteps:
Default: 0:00:00.178758
spencer_mc: 0:00:00.007183
100000 timesteps:
Default: 0:00:01.400652
spencer_mc: 0:00:00.049189

meteoDaniel · 2019-05-21T17:41:03Z

Some issue of stickler-ci are not related to my changes.
Please help me for the other tests.

pvlib/test/test_solarposition.py

cwhanse · 2019-05-21T18:55:18Z

pvlib/solarposition.py

+    return result
+
+
+def datetime2julian(times):


Does this function do the same thing as tools.datetime_to_djd?

No not quite the same but it is similar. Actually datetime_to_djd is referencing to the dublin standard and the function is not able to work with arrays. I will have a look If it is possble to adapt my function that it is working with dublin julian days so that we can use my implementation.

cwhanse

I can't help generally with the stickler complaints. One line error is causing both test errors, I think.

Thanks for the contribution!

meteoDaniel · 2019-05-22T04:32:32Z

@cwhanse I hope i solved all issues for linter and lgtm. Lgtm is still running at the moment.
Can you tell me the difference between apparent and non apparent solar azimuth/zenith/elevation ?
Thanks in advance.

mikofski · 2019-05-22T09:36:12Z

Hi @meteoDaniel , FYI: I pushed a PR with new changes to SolarUtils that allows you to specify either the year or a specific sequence of datetimes, it's very fast

In [1]: from solar_utils import solposAM, get_solposAM
   ...: location = [35.56836, -119.2022, -8.0]
   ...: weather = [1015.62055, 40.0]
   ...: from datetime import datetime, timedelta
   ...: times = [(datetime(2017,1,1,0,0,0)+timedelta(hours=h)).timetuple()[:6] for h in range(1000)]

In [2]: %timeit get_solposAM(location, times, weather)
1.82 ms ± 16.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

On my PC that's about twice as fast as the Spencer method per the 1000-timestep benchmark in this PR. SolarUtils is a very thin Python wrapper around the freely available NREL SOLPOS algorithm source code from Martin Rymes, 1998-2001.

Although SolarUtils runtime scales linearly, while this Spencer method runtime stays fairly low, probably due to NumPy?

method	1000	10,000	100,000
Spencer	3.093	7.183	49.189
SOLPOS	1.82	18.8	210
SPA	23.037	178.758	1400.652

Porting SOLPOS to Python would probably have similar results. The source is freely available, so perhaps if you require even more speed that would be an option.

In [6]: times = [(datetime(2017,1,1,0,0,0)+timedelta(hours=h)).timetuple()[:6] for h in range(10000)]

In [7]: %timeit get_solposAM(location, times, weather)
18.8 ms ± 117 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

One of the stated purpose of pvlib python is to catalog all useful solar algorithms, and IMO this Spencer method is very useful! Thanks so much for contributing this! Let me know if there's anyway I can help.

cwhanse · 2019-05-22T13:52:33Z

Can you tell me the difference between apparent and non apparent solar azimuth/zenith/elevation ?
Thanks in advance.

Generally, the 'apparent' values account for atmospheric refraction, the non-apparent values are geometric.

cwhanse

Test failure was for get_psm3 probably an interruption in server response. I restarted the tests.

wholmgren

Thanks @meteoDaniel.

docs/sphinx/source/whatsnew/v0.7.1.rst

pvlib/solarposition.py

wholmgren · 2019-05-22T14:19:08Z

pvlib/solarposition.py

+def spencer_mc(times, latitude, longitude):
+    """
+    Calculate the solar position using a python implementation of the
+    Spencer (1972) formulation provided by meteocontrol


My preference is a code comment saying who added the code and what their affiliation is. I don't think this belongs in rendered documentation.

Some subtle form of recognition seems fine with me. But if the method is identical to spencer's paper, then the suffix '_mc' might be too much.

pvlib/solarposition.py

wholmgren · 2019-05-22T14:34:21Z

pvlib/tools.py

@@ -425,3 +429,25 @@ def _golden_sect_DataFrame(params, VL, VH, func):
            raise Exception("EXCEPTION:iterations exceeded maximum (50)")

    return func(df, 'V1'), df['V1']
+
+
+def datetime_to_julian(times):


I don't like the duplication of datetime_to_djd. Let's make the old function work with both datetime and DatetimeIndex. I suspect converting a datetime input to a Timestamp will do the job.

pvlib/test/test_tools.py

wholmgren · 2019-05-22T14:37:44Z

pvlib/test/test_tools.py

 import pytest

+import pvlib.tools


tools is already imported on the next line.

Ähh right! ;D

wholmgren · 2019-05-22T14:40:08Z

pvlib/test/test_solarposition.py

@@ -754,8 +754,125 @@ def test_sun_rise_set_transit_geometric(expected_rise_set_spa, golden_mst):
                       atol=np.abs(expected_transit_error).max())


-# put numba tests at end of file to minimize reloading
+def test_spencer_mc():


Please follow this test pattern so that we can be sure the output is correct to within some tolerance:

pvlib-python/pvlib/test/test_solarposition.py

Lines 420 to 455 in d68617c

def test_ephemeris_physical(expected_solpos, golden_mst):

times = pd.date_range(datetime.datetime(2003, 10, 17, 12, 30, 30),

periods=1, freq='D', tz=golden_mst.tz)

ephem_data = solarposition.ephemeris(times, golden_mst.latitude,

golden_mst.longitude,

pressure=82000,

temperature=11)

expected_solpos.index = times

expected_solpos = np.round(expected_solpos, 2)

ephem_data = np.round(ephem_data, 2)

assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])

def test_ephemeris_physical_dst(expected_solpos, golden):

times = pd.date_range(datetime.datetime(2003, 10, 17, 13, 30, 30),

periods=1, freq='D', tz=golden.tz)

ephem_data = solarposition.ephemeris(times, golden.latitude,

golden.longitude, pressure=82000,

temperature=11)

expected_solpos.index = times

expected_solpos = np.round(expected_solpos, 2)

ephem_data = np.round(ephem_data, 2)

assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])

def test_ephemeris_physical_no_tz(expected_solpos, golden_mst):

times = pd.date_range(datetime.datetime(2003, 10, 17, 19, 30, 30),

periods=1, freq='D')

ephem_data = solarposition.ephemeris(times, golden_mst.latitude,

golden_mst.longitude,

pressure=82000,

temperature=11)

expected_solpos.index = times

expected_solpos = np.round(expected_solpos, 2)

ephem_data = np.round(ephem_data, 2)

assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])

Good Idea. I will bring it in the next commit

wholmgren · 2019-05-22T14:41:03Z

pvlib/solarposition.py

+    """
+    Calculate the solar position using a python implementation of the
+    Spencer (1972) formulation provided by meteocontrol
+


I recommend that we add a comment on the relative accuracy and speed of the function. Perhaps it belongs in a Notes section below.

mikofski

a few places where it would be nice to collapse duplicate code with existing functions perhaps with some modifications to make them work better - thanks!

mikofski · 2019-05-22T17:12:59Z

pvlib/solarposition.py

+    eot -= 0.032077 * sin_gamma[0]
+    eot -= 0.014615 * cos_gamma[1]
+    eot -= 0.040849 * sin_gamma[1]
+    eot *= 229.18


I believe this section may be a repeat of equation_of_time_spencer71:

day_angle = _calculate_simple_day_angle(dayofyear) # convert from radians to minutes per day = 24[h/day] * 60[min/h] / 2 / pi eot = (1440.0 / 2 / np.pi) * ( 0.0000075 + 0.001868 * np.cos(day_angle) - 0.032077 * np.sin(day_angle) - 0.014615 * np.cos(2.0 * day_angle) - 0.040849 * np.sin(2.0 * day_angle) ) return eot

I think we should try, if possible, to use DRY vs WET code -- can you pleases consider reusing the existing code functions?

EG: eot = equation_of_time_spencer71(dayofyear)

where I suspect that you can use the Julian dates you've calculated in place of day of year, or adjust them slightly

see _calculate_simple_day_angle

Right. New implementatino uses the equation_of_time_spencer71 function.
Acually I am not really familiar with the Julian Days and Day Of Year usage. I will check it out if there are any ways to use existend code.

mikofski · 2019-05-22T17:16:58Z

pvlib/solarposition.py

+    eot *= 229.18
+
+    # True local time
+    tlt = (day_time + longitude / 15 + eot / 60) % 24 - 12


I suspect this is similar as hour_angle - perhaps there is a way to reuse the existing code, or modify it to suit your needs so we don't have duplicate code?

looking in more detail they are the nearly the same:
tlt = hour_angle / 15.

but IMO this version here, looks a little easier to read, with the use of mod(), so might be nice to replace the existing hour_angle with this, assuming there is better or same performance

I tried hour_angle and it was not working in this code.

But do you agree they are the same equation, that only differ by a factor of 1/15? we have noticed that different formulations have different "offsets" so perhaps adjusting the input may help? Or perhaps leave it, raise an issue to collapse this later, and maybe add a #TODO in the code. If @cwhanse and @wholmgren are okay then so am I

…7.0.0.

mikofski · 2019-05-23T02:38:41Z

@meteoDaniel just add W504,W504 to .stickler.yml. see pep8.org and PyQA. it was a mistake for stickler to raise the warnings

cwhanse · 2019-05-30T14:54:20Z

@meteoDaniel please merge master into this PR, should take care of the test failures for get_psm3

adriesse · 2019-07-10T21:47:37Z

I believe these calculations are identical to the ones in the "analytical" aka "geometric" solar position functions already in pvlib. (#583)

mikofski · 2019-07-10T23:27:37Z

OMG! @adriesse is correct! I guess I didn't review this very well. Might still be worth it, but would be better to use the existing "analytical/geometric" functions or make them better if there's room for improvement.

azimuth

pvlib-python/pvlib/solarposition.py

Lines 1517 to 1524 in 03bb143

    
           # Compute the sun's azimuth angle. 
        
           y = -(np.sin(latitude_radians) * np.sin(elevation) - np.sin(declination)) \ 
        
               / (np.cos(latitude_radians) * np.cos(elevation)) 
        
           azimuth = np.arccos(y) 
        
           # Convert azimuth angle from 0-pi to 0-2pi. 
        
           tlt_filter = 0 <= tlt 
        
           azimuth[tlt_filter] = 2 * np.pi - azimuth[tlt_filter]

is the same as

pvlib-python/pvlib/solarposition.py

Lines 1237 to 1261 in e3635c2

    
           numer = (np.cos(zenith) * np.sin(latitude) - np.sin(declination)) 
        
           denom = (np.sin(zenith) * np.cos(latitude)) 
        
           # cases that would generate new NaN values are safely ignored here 
        
           # since they are dealt with further below 
        
           with np.errstate(invalid='ignore', divide='ignore'): 
        
               cos_azi = numer / denom 
        
           # when zero division occurs, use the limit value of the analytical 
        
           # expression 
        
           cos_azi = \ 
        
               np.where(np.isclose(denom,    0.0, rtol=0.0, atol=1e-8),  1.0, cos_azi) 
        
           # when too many round-ups in floating point math take cos_azi beyond 
        
           # 1.0, use 1.0 
        
           cos_azi = \ 
        
               np.where(np.isclose(cos_azi,  1.0, rtol=0.0, atol=1e-8),  1.0, cos_azi) 
        
           cos_azi = \ 
        
               np.where(np.isclose(cos_azi, -1.0, rtol=0.0, atol=1e-8), -1.0, cos_azi) 
        
           # when NaN values occur in input, ignore and pass to output 
        
           with np.errstate(invalid='ignore'): 
        
               sign_ha = np.sign(hourangle) 
        
           return sign_ha * np.arccos(cos_azi) + np.pi

zenith

pvlib-python/pvlib/solarposition.py

Lines 1490 to 1511 in 03bb143

    
           julians, julian_2000 = datetime_to_julian(times) 
        
           julians_2000 = np.asarray(julians, dtype=np.float) - julian_2000 
        
           latitude_radians = np.radians(latitude) 
        
           day_time = (julians_2000 % 1) * 24 
        
           declination = np.array(declination_spencer71(times.dayofyear)) 
        
           # Equation of time (difference between standard time and solar time). 
        
           eot = np.array(equation_of_time_spencer71(times.dayofyear)) 
        
           # True local time 
        
           tlt = (day_time + longitude / 15 + eot / 60) % 24 - 12 
        
           # Solar hour angle 
        
           ha = np.radians(tlt * 15) 
        
           # Calculate sun elevation. 
        
           sin_sun_elevation = ( 
        
                   np.sin(declination) * np.sin(latitude_radians) + 
        
                   np.cos(declination) * np.cos(latitude_radians) * np.cos(ha) 
        
           )

is the same as

https://github.com/pvlib/pvlib-python/blob/master/pvlib/solarposition.py#L1309-L1312

mikofski

As @adriesse pointed out, you have reproduced the analytical expressions

update your branch from the latest pvlib master
please replace code for zenith calc, julian time, hour angle, declination, etc, with a call to solar_zenith_analytical(latitude, hourangle, declination)
please replace code for azimuth calc, etc, with a call to solar_azimuth_analytical(latitude, hourangle, declination, zenith)
can you please provide a warning, similar to the one in the existing analytical functions, that the Spencer calculation does not include atmospheric refraction, and therefore it only returns extraterrestrial zenith, which can be up to 2-degrees different from apparent zenith.

thanks! still think this is a valuable contribution, would like to see it v0.7.0 !

mikofski · 2019-07-10T23:29:38Z

pvlib/solarposition.py

+    zenith = np.pi / 2 - elevation
+
+    # Compute the sun's azimuth angle.
+    y = -(np.sin(latitude_radians) * np.sin(elevation) - np.sin(declination)) \


@meteoDaniel can you please replace these lines

pvlib-python/pvlib/solarposition.py

Lines 1517 to 1524 in 03bb143

# Compute the sun's azimuth angle.

y = -(np.sin(latitude_radians) * np.sin(elevation) - np.sin(declination)) \

/ (np.cos(latitude_radians) * np.cos(elevation))

azimuth = np.arccos(y)

# Convert azimuth angle from 0-pi to 0-2pi.

tlt_filter = 0 <= tlt

azimuth[tlt_filter] = 2 * np.pi - azimuth[tlt_filter]

with a call to solar_azimuth_analytical(latitude, hourangle, declination, zenith) ?

mikofski · 2019-07-10T23:30:51Z

pvlib/solarposition.py

+     Springer Science & Business Media, 2008.
+    """
+
+    julians, julian_2000 = datetime_to_julian(times)


@meteoDaniel can you please replace all of these lines

pvlib-python/pvlib/solarposition.py

Lines 1490 to 1511 in 03bb143

julians, julian_2000 = datetime_to_julian(times)

julians_2000 = np.asarray(julians, dtype=np.float) - julian_2000

latitude_radians = np.radians(latitude)

day_time = (julians_2000 % 1) * 24

declination = np.array(declination_spencer71(times.dayofyear))

# Equation of time (difference between standard time and solar time).

eot = np.array(equation_of_time_spencer71(times.dayofyear))

# True local time

tlt = (day_time + longitude / 15 + eot / 60) % 24 - 12

# Solar hour angle

ha = np.radians(tlt * 15)

# Calculate sun elevation.

sin_sun_elevation = (

np.sin(declination) * np.sin(latitude_radians) +

np.cos(declination) * np.cos(latitude_radians) * np.cos(ha)

)

with a call to solar_zenith_analytical(latitude, hourangle, declination) ?

adriesse · 2019-07-11T08:13:37Z

Maybe the title could be changed to something like: "Convenience function to calculate the complete solar position using the analytical/geometric methods."

Not sure whether the name Spencer should be used.

adriesse · 2019-07-11T08:40:46Z

By the way, refraction reduces the apparent zenith angle at the horizon by about 0.5 degrees (not 2 degrees as noted above).

wholmgren · 2020-02-13T16:29:29Z

Let's start fresh in a new PR if there's anything we want to continue here.

Daniel Lassahn and others added 5 commits May 21, 2019 00:28

Add a new way to calculate the solarposition

d97cf22

Tidy

07c9ea3

Tidy

98f6539

Merge branch 'master' into master

19531ad

Use existent declination function

5d46472

meteoDaniel changed the title ~~Add a new way to calculate the solarposition~~ Add a new and faster way to calculate the solarposition May 21, 2019

Daniel Lassahn added 2 commits May 21, 2019 10:07

Documentation stuff

6ae7dae

Tidy

a43cf78

cwhanse reviewed May 21, 2019

View reviewed changes

pvlib/test/test_solarposition.py Outdated Show resolved Hide resolved

cwhanse reviewed May 21, 2019

View reviewed changes

Daniel Lassahn added 4 commits May 21, 2019 14:02

Adapt tests and place datetime_to_julian in tools

1c1b2c6

Linter optimisation

0749c55

Linter optimisation

bf88e3c

Linter optimisation

b20ca53

cwhanse approved these changes May 22, 2019

View reviewed changes

wholmgren reviewed May 22, 2019

View reviewed changes

mikofski reviewed May 22, 2019

View reviewed changes

Daniel Lassahn added 2 commits May 22, 2019 18:34

Rename function, use existend equation of time function, whatsnew in …

fe5bf63

…7.0.0.

Avoid duplicate import of pvlib tools

71dc8c5

cwhanse added this to the 0.6.4 milestone Jun 11, 2019

Handle timezones in timeRange

03bb143

cwhanse modified the milestones: 0.6.4, 0.7.0 Jun 25, 2019

mikofski requested changes Jul 10, 2019

View reviewed changes

wholmgren removed this from the 0.7.0 milestone Nov 26, 2019

wholmgren closed this Feb 13, 2020

mikofski mentioned this pull request May 17, 2025

New default solar position algorithm #2455

Open

	def test_ephemeris_physical(expected_solpos, golden_mst):
	times = pd.date_range(datetime.datetime(2003, 10, 17, 12, 30, 30),
	periods=1, freq='D', tz=golden_mst.tz)
	ephem_data = solarposition.ephemeris(times, golden_mst.latitude,
	golden_mst.longitude,
	pressure=82000,
	temperature=11)
	expected_solpos.index = times
	expected_solpos = np.round(expected_solpos, 2)
	ephem_data = np.round(ephem_data, 2)
	assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])


	def test_ephemeris_physical_dst(expected_solpos, golden):
	times = pd.date_range(datetime.datetime(2003, 10, 17, 13, 30, 30),
	periods=1, freq='D', tz=golden.tz)
	ephem_data = solarposition.ephemeris(times, golden.latitude,
	golden.longitude, pressure=82000,
	temperature=11)
	expected_solpos.index = times
	expected_solpos = np.round(expected_solpos, 2)
	ephem_data = np.round(ephem_data, 2)
	assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])


	def test_ephemeris_physical_no_tz(expected_solpos, golden_mst):
	times = pd.date_range(datetime.datetime(2003, 10, 17, 19, 30, 30),
	periods=1, freq='D')
	ephem_data = solarposition.ephemeris(times, golden_mst.latitude,
	golden_mst.longitude,
	pressure=82000,
	temperature=11)
	expected_solpos.index = times
	expected_solpos = np.round(expected_solpos, 2)
	ephem_data = np.round(ephem_data, 2)
	assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])

	# Compute the sun's azimuth angle.
	y = -(np.sin(latitude_radians) * np.sin(elevation) - np.sin(declination)) \
	/ (np.cos(latitude_radians) * np.cos(elevation))
	azimuth = np.arccos(y)

	# Convert azimuth angle from 0-pi to 0-2pi.
	tlt_filter = 0 <= tlt
	azimuth[tlt_filter] = 2 * np.pi - azimuth[tlt_filter]

	julians, julian_2000 = datetime_to_julian(times)
	julians_2000 = np.asarray(julians, dtype=np.float) - julian_2000

	latitude_radians = np.radians(latitude)
	day_time = (julians_2000 % 1) * 24

	declination = np.array(declination_spencer71(times.dayofyear))

	# Equation of time (difference between standard time and solar time).
	eot = np.array(equation_of_time_spencer71(times.dayofyear))

	# True local time
	tlt = (day_time + longitude / 15 + eot / 60) % 24 - 12

	# Solar hour angle
	ha = np.radians(tlt * 15)

	# Calculate sun elevation.
	sin_sun_elevation = (
	np.sin(declination) * np.sin(latitude_radians) +
	np.cos(declination) * np.cos(latitude_radians) * np.cos(ha)
	)

Add a new and faster way to calculate the solarposition #730

Add a new and faster way to calculate the solarposition #730

Uh oh!

Conversation

meteoDaniel commented May 21, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

pvlib python pull request guidelines

Uh oh!

meteoDaniel commented May 21, 2019

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cwhanse left a comment

Choose a reason for hiding this comment

Uh oh!

meteoDaniel commented May 22, 2019

Uh oh!

mikofski commented May 22, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

cwhanse commented May 22, 2019

Uh oh!

cwhanse left a comment

Choose a reason for hiding this comment

Uh oh!

wholmgren left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mikofski left a comment

Choose a reason for hiding this comment

Uh oh!

mikofski May 22, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mikofski commented May 23, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

cwhanse commented May 30, 2019

Uh oh!

adriesse commented Jul 10, 2019

Uh oh!

mikofski commented Jul 10, 2019

azimuth

zenith

Uh oh!

meteoDaniel commented May 21, 2019 •

edited

Loading

mikofski commented May 22, 2019 •

edited

Loading

mikofski May 22, 2019 •

edited

Loading

mikofski commented May 23, 2019 •

edited

Loading