Refine the repairs and add test

adriesse · adriesse · commit bcdf6263a363 · 2018-02-23T10:25:10.000+01:00
diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py
@@ -1096,17 +1096,27 @@ def solar_azimuth_analytical(latitude, hour_angle, declination, zenith):
     hour_angle
     solar_zenith_analytical
     """
+
+    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)
 
-        cos_azi = ((np.cos(zenith) * np.sin(latitude) - np.sin(declination)) /
-                   (np.sin(zenith) * np.cos(latitude)))
+    # 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)
 
-    cos_azi = np.where(np.isclose(zenith,   0.0),  1.0, cos_azi)
-    cos_azi = np.where(np.isclose(latitude, 0.0),  1.0, cos_azi)
-    cos_azi = np.where(np.isclose(cos_azi,  1.0),  1.0, cos_azi)
-    cos_azi = np.where(np.isclose(cos_azi, -1.0), -1.0, cos_azi)
+    # when NaN values occur in input, ignore and pass to output
+    with np.errstate(invalid='ignore'):
+        sign_ha = np.sign(hour_angle)
 
-    return (np.sign(hour_angle) * np.arccos(cos_azi) + np.pi)
+    return (sign_ha * np.arccos(cos_azi) + np.pi)
 
 
 def solar_zenith_analytical(latitude, hour_angle, declination):
diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py
@@ -341,7 +341,7 @@ def test_get_solarposition_deltat(delta_t, method, expected):
     this_expected = np.round(this_expected, 5)
     ephem_data = np.round(ephem_data, 5)
     assert_frame_equal(this_expected, ephem_data[this_expected.columns])
-   
+
 
 def test_get_solarposition_no_kwargs(expected_solpos):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
@@ -463,3 +463,27 @@ def test_analytical_azimuth():
                        atol=0.01)
     assert np.allclose(azimuth_2[idx], solar_azimuth.as_matrix()[idx],
                        atol=0.017)
+
+    # test for NaN values at boundary conditions (PR #431)
+    # use degrees rather than radians for clarity
+    test_angles = np.radians(np.array(
+                   [[   0., -180.,  -20.],
+                    [   0.,    0.,   -5.],
+                    [   0.,    0.,    0.],
+                    [   0.,    0.,   15.],
+                    [   0.,  180.,   20.],
+                    [  30.,    0.,  -20.],
+                    [  30.,    0.,   -5.],
+                    [  30.,    0.,    0.],
+                    [  30.,  180.,    5.],
+                    [  30.,    0.,   10.],
+                    [ -30.,    0.,  -20.],
+                    [ -30.,    0.,  -15.],
+                    [ -30.,    0.,    0.],
+                    [ -30., -180.,    5.],
+                    [ -30.,  180.,   10.]]))
+
+    zeniths  = solarposition.solar_zenith_analytical(*test_angles.T)
+    azimuths = solarposition.solar_azimuth_analytical(*test_angles.T, zeniths)
+
+    assert not np.isnan(azimuths).any()
