add test for hour_angle, vectorize (#599)

* add test for hour_angle, vectorize

* closes #598 vectorize to make it more efficient
* closes #597 add test

Signed-off-by: Mark Mikofski <[email protected]>

* BUG: use np.int64 works better for older numpy version than python int

* also converting times to int before subtracting works better for older
pandas versions which were not calculating the timedeltas correctly
* remove comment with missing space after hash, add FIXME that explains
why the expected values are slightly different than the SPA calculator
output

* fix hanging indent

* BUG: replace utcoffset() with reliable, efficient approach ...

* ... suggested by @wholmgren (thx!)
* utcoffset() is unpredictable when used with pandas datetime indices
* only predictable with Python datetime objects or pandas Timestamps
* instead replace tzinfo with None to get naive local times, and
calculate difference from tz-aware times to get timezones

* BUG: combine arithmetic to make calculation more efficient

* also use asarray wrapper at return to ensure consistency
* add comments to explain to future maintainers

* stickler

* remove try-except, doesn't work anyway, wait for pandas >=0.15.0

Author: mikofski

pvlib/solarposition.py

@@ -1327,8 +1327,11 @@ def hour_angle(times, longitude, equation_of_time):
     equation_of_time_Spencer71
     equation_of_time_pvcdrom
     """
-    hours = np.array([(t - t.tz.localize(
-        dt.datetime(t.year, t.month, t.day)
-    )).total_seconds() / 3600. for t in times])
-    timezone = times.tz.utcoffset(times).total_seconds() / 3600.
-    return 15. * (hours - 12. - timezone) + longitude + equation_of_time / 4.
+    naive_times = times.tz_localize(None)  # naive but still localized
+    # hours - timezone = (times - normalized_times) - (naive_times - times)
+    hrs_minus_tzs = 1 / (3600. * 1.e9) * (
+        2 * times.astype(np.int64) - times.normalize().astype(np.int64) -
+        naive_times.astype(np.int64))
+    # ensure array return instead of a version-dependent pandas <T>Index
+    return np.asarray(
+        15. * (hrs_minus_tzs - 12.) + longitude + equation_of_time / 4.)

pvlib/test/test_solarposition.py

@@ -694,3 +694,27 @@ def test_analytical_azimuth():
     azimuths = solarposition.solar_azimuth_analytical(*test_angles.T, zenith=zeniths)
 
     assert not np.isnan(azimuths).any()
+
+
+def test_hour_angle():
+    """
+    Test conversion from hours to hour angles in degrees given the following
+    inputs from NREL SPA calculator at Golden, CO
+    date,times,eot,sunrise,sunset
+    1/2/2015,7:21:55,-3.935172,-70.699400,70.512721
+    1/2/2015,16:47:43,-4.117227,-70.699400,70.512721
+    1/2/2015,12:04:45,-4.026295,-70.699400,70.512721
+    """
+    longitude = -105.1786  # degrees
+    times = pd.DatetimeIndex([
+        '2015-01-02 07:21:55.2132',
+        '2015-01-02 16:47:42.9828',
+        '2015-01-02 12:04:44.6340'
+    ]).tz_localize('Etc/GMT+7')
+    eot = np.array([-3.935172, -4.117227, -4.026295])
+    hours = solarposition.hour_angle(times, longitude, eot)
+    expected = (-70.682338, 70.72118825000001, 0.000801250)
+    # FIXME: there are differences from expected NREL SPA calculator values
+    # sunrise: 4 seconds, sunset: 48 seconds, transit: 0.2 seconds
+    # but the differences may be due to other SPA input parameters
+    assert np.allclose(hours, expected)