Repair solar_azimuth_analytical

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

@@ -17,11 +17,15 @@ Bug fixes
 * physicaliam now returns a Series if called with a Series as an
   argument. (:issue:`397`)
 * corrected docstring for irradiance.total_irrad (:issue: '423')
+* modified solar_azimuth_analytical to handle some borderline cases, thereby
+  avoiding the NaN values and/or warnings that were previously produced
+  (:issue: '420')
 * removed RuntimeWarnings due to divide by 0 or nans in input data within
   irradiance.perez, clearsky.simplified_solis, pvsystem.adrinverter,
   pvsystem.ashraeiam, pvsystem.physicaliam, pvsystem.sapm_aoi_loss,
   pvsystem.v_from_i. (:issue:`428`)
 
+
 Documentation
 ~~~~~~~~~~~~~
 * Improve physicaliam doc string. (:issue:`397`)
@@ -37,5 +41,6 @@ Contributors
 * Cliff Hansen
 * Will Holmgren
 * KonstantinTr
+* Anton Driesse
 
 

pvlib/solarposition.py

@@ -1096,9 +1096,27 @@ def solar_azimuth_analytical(latitude, hour_angle, declination, zenith):
     hour_angle
     solar_zenith_analytical
     """
-    return np.sign(hour_angle) * np.abs(np.arccos((np.cos(zenith) * np.sin(
-        latitude) - np.sin(declination)) / (np.sin(zenith) * np.cos(
-        latitude)))) + np.pi
+
+    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(hour_angle)
+
+    return (sign_ha * np.arccos(cos_azi) + np.pi)
 
 
 def solar_zenith_analytical(latitude, hour_angle, declination):

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,26 @@ 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)
+    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, zenith=zeniths)
+
+    assert not np.isnan(azimuths).any()