avoid runtime warnings from nan comparisons and divide by zero (#429)

* replace manual perez binning with np.digitize

* update whatsnew

* add invalid=ignore

* refactor solis masking to avoid nan comparisons

* avoid iam and adrinverter warnings. added tests to make sure it works

* avoid sapm log warnings

* avoid runtime warning in sapm aoi

* add another adrinverter test, fix night bug

* document adrinverter warnings

* minimum numpy version to 1.10.1

* ignore over warnings in v_from_i

* use specific names for mask arrays

* comments, use errstate for clarity

* ensure adrinverter handles nan inputs

Author: wholmgren

.travis.yml

@@ -53,7 +53,7 @@ install:
     - if [[ "$CONDA_ENV" == "py27-min" ]]; then
         pip uninstall numpy --yes;
         pip uninstall pandas --yes;
-        pip install --no-cache-dir numpy==1.9.0;
+        pip install --no-cache-dir numpy==1.10.1;
         pip install --no-cache-dir pandas==0.14.0;
       fi
     - conda list

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

@@ -17,6 +17,10 @@ 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')
+* 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
 ~~~~~~~~~~~~~

pvlib/clearsky.py

@@ -433,11 +433,7 @@ def simplified_solis(apparent_elevation, aod700=0.1, precipitable_water=1.,
     w = precipitable_water
 
     # algorithm fails for pw < 0.2
-    if np.isscalar(w):
-        w = 0.2 if w < 0.2 else w
-    else:
-        w = w.copy()
-        w[w < 0.2] = 0.2
+    w = np.maximum(w, 0.2)
 
     # this algorithm is reasonably fast already, but it could be made
     # faster by precalculating the powers of aod700, the log(p/p0), and
@@ -542,8 +538,10 @@ def _calc_taud(w, aod700, p):
     elif np.isscalar(aod700):
         aod700 = np.full_like(w, aod700)
 
-    aod700_mask = aod700 < 0.05
-    aod700_mask = np.array([aod700_mask, ~aod700_mask], dtype=np.int)
+    # set up nan-tolerant masks
+    aod700_lt_0p05 = np.full_like(aod700, False, dtype='bool')
+    np.less(aod700, 0.05, where=~np.isnan(aod700), out=aod700_lt_0p05)
+    aod700_mask = np.array([aod700_lt_0p05, ~aod700_lt_0p05], dtype=np.int)
 
     # create tuples of coefficients for
     # aod700 < 0.05, aod700 >= 0.05

pvlib/irradiance.py

@@ -995,7 +995,8 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     delta = dhi * airmass / dni_extra
 
     # epsilon is the sky's "clearness"
-    eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
+    with np.errstate(invalid='ignore'):
+        eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
 
     # numpy indexing below will not work with a Series
     if isinstance(eps, pd.Series):
@@ -1005,15 +1006,8 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     # rules. 1 = overcast ... 8 = clear (these names really only make
     # sense for small zenith angles, but...) these values will
     # eventually be used as indicies for coeffecient look ups
-    ebin = np.zeros_like(eps, dtype=np.int8)
-    ebin[eps < 1.065] = 1
-    ebin[(eps >= 1.065) & (eps < 1.23)] = 2
-    ebin[(eps >= 1.23) & (eps < 1.5)] = 3
-    ebin[(eps >= 1.5) & (eps < 1.95)] = 4
-    ebin[(eps >= 1.95) & (eps < 2.8)] = 5
-    ebin[(eps >= 2.8) & (eps < 4.5)] = 6
-    ebin[(eps >= 4.5) & (eps < 6.2)] = 7
-    ebin[eps >= 6.2] = 8
+    ebin = np.digitize(eps, (0., 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2))
+    ebin[np.isnan(eps)] = 0
 
     # correct for 0 indexing in coeffecient lookup
     # later, ebin = -1 will yield nan coefficients
@@ -1333,7 +1327,7 @@ def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
     Determine DNI from GHI using the DIRINDEX model, which is a modification of
     the DIRINT model with information from a clear sky model.
 
-    DIRINDEX [1] improves upon the DIRINT model by taking into account 
+    DIRINDEX [1] improves upon the DIRINT model by taking into account
     turbidity when used with the Ineichen clear sky model results.
 
     Parameters
@@ -1396,7 +1390,7 @@ def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
                         use_delta_kt_prime=use_delta_kt_prime,
                         temp_dew=temp_dew)
 
-    dni_dirint_clearsky = dirint(ghi_clearsky, zenith, times, 
+    dni_dirint_clearsky = dirint(ghi_clearsky, zenith, times,
                                  pressure=pressure,
                                  use_delta_kt_prime=use_delta_kt_prime,
                                  temp_dew=temp_dew)

pvlib/pvsystem.py

@@ -732,9 +732,10 @@ def ashraeiam(aoi, b=0.05):
     physicaliam
     '''
 
-    iam = 1 - b*((1/np.cos(np.radians(aoi)) - 1))
-
-    iam = np.where(np.abs(aoi) >= 90, 0, iam)
+    iam = 1 - b * ((1 / np.cos(np.radians(aoi)) - 1))
+    aoi_gte_90 = np.full_like(aoi, False, dtype='bool')
+    np.greater_equal(np.abs(aoi), 90, where=~np.isnan(aoi), out=aoi_gte_90)
+    iam = np.where(aoi_gte_90, 0, iam)
     iam = np.maximum(0, iam)
 
     if isinstance(iam, pd.Series):
@@ -836,16 +837,18 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
     # after deducting the reflected portion of each
     iam = ((1 - (rho_para + rho_perp) / 2) / (1 - rho_zero) * tau / tau_zero)
 
-    # angles near zero produce nan, but iam is defined as one
-    small_angle = 1e-06
-    iam = np.where(np.abs(aoi) < small_angle, 1.0, iam)
+    with np.errstate(invalid='ignore'):
+        # angles near zero produce nan, but iam is defined as one
+        small_angle = 1e-06
+        iam = np.where(np.abs(aoi) < small_angle, 1.0, iam)
 
-    # angles at 90 degrees can produce tiny negative values, which should be zero
-    # this is a result of calculation precision rather than the physical model
-    iam = np.where(iam < 0, 0, iam)
+        # angles at 90 degrees can produce tiny negative values,
+        # which should be zero. this is a result of calculation precision
+        # rather than the physical model
+        iam = np.where(iam < 0, 0, iam)
 
-    # for light coming from behind the plane, none can enter the module
-    iam = np.where(aoi > 90, 0, iam)
+        # for light coming from behind the plane, none can enter the module
+        iam = np.where(aoi > 90, 0, iam)
 
     if isinstance(aoi_input, pd.Series):
         iam = pd.Series(iam, index=aoi_input.index)
@@ -1296,12 +1299,27 @@ def sapm(effective_irradiance, temp_cell, module):
     q = 1.60218e-19  # Elementary charge in units of coulombs
     kb = 1.38066e-23  # Boltzmann's constant in units of J/K
 
-    Ee = effective_irradiance
+    # avoid problem with integer input
+    Ee = np.array(effective_irradiance, dtype='float64')
+
+    # set up masking for 0, positive, and nan inputs
+    Ee_gt_0 = np.full_like(Ee, False, dtype='bool')
+    Ee_eq_0 = np.full_like(Ee, False, dtype='bool')
+    notnan = ~np.isnan(Ee)
+    np.greater(Ee, 0, where=notnan, out=Ee_gt_0)
+    np.equal(Ee, 0, where=notnan, out=Ee_eq_0)
 
     Bvmpo = module['Bvmpo'] + module['Mbvmp']*(1 - Ee)
     Bvoco = module['Bvoco'] + module['Mbvoc']*(1 - Ee)
     delta = module['N'] * kb * (temp_cell + 273.15) / q
 
+    # avoid repeated computation
+    logEe = np.full_like(Ee, np.nan)
+    np.log(Ee, where=Ee_gt_0, out=logEe)
+    logEe = np.where(Ee_eq_0, -np.inf, logEe)
+    # avoid repeated __getitem__
+    cells_in_series = module['Cells_in_Series']
+
     out = OrderedDict()
 
     out['i_sc'] = (
@@ -1312,13 +1330,13 @@ def sapm(effective_irradiance, temp_cell, module):
         (1 + module['Aimp']*(temp_cell - T0)))
 
     out['v_oc'] = np.maximum(0, (
-        module['Voco'] + module['Cells_in_Series']*delta*np.log(Ee) +
+        module['Voco'] + cells_in_series * delta * logEe +
         Bvoco*(temp_cell - T0)))
 
     out['v_mp'] = np.maximum(0, (
         module['Vmpo'] +
-        module['C2']*module['Cells_in_Series']*delta*np.log(Ee) +
-        module['C3']*module['Cells_in_Series']*((delta*np.log(Ee)) ** 2) +
+        module['C2'] * cells_in_series * delta * logEe +
+        module['C3'] * cells_in_series * ((delta * logEe) ** 2) +
         Bvmpo*(temp_cell - T0)))
 
     out['p_mp'] = out['i_mp'] * out['v_mp']
@@ -1518,7 +1536,10 @@ def sapm_aoi_loss(aoi, module, upper=None):
 
     aoi_loss = np.polyval(aoi_coeff, aoi)
     aoi_loss = np.clip(aoi_loss, 0, upper)
-    aoi_loss = np.where(aoi < 0, 0, aoi_loss)
+    # nan tolerant masking
+    aoi_lt_0 = np.full_like(aoi, False, dtype='bool')
+    np.less(aoi, 0, where=~np.isnan(aoi), out=aoi_lt_0)
+    aoi_loss = np.where(aoi_lt_0, 0, aoi_loss)
 
     if isinstance(aoi, pd.Series):
         aoi_loss = pd.Series(aoi_loss, aoi.index)
@@ -1912,9 +1933,11 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     # Only compute using LambertW if there are cases with Gsh>0
     if np.any(idx_p):
         # LambertW argument, cannot be float128, may overflow to np.inf
-        argW = I0[idx_p] / (Gsh[idx_p]*a[idx_p]) * \
-            np.exp((-I[idx_p] + IL[idx_p] + I0[idx_p]) /
-                   (Gsh[idx_p]*a[idx_p]))
+        # overflow is explicitly handled below, so ignore warnings here
+        with np.errstate(over='ignore'):
+            argW = (I0[idx_p] / (Gsh[idx_p]*a[idx_p]) *
+                    np.exp((-I[idx_p] + IL[idx_p] + I0[idx_p]) /
+                           (Gsh[idx_p]*a[idx_p])))
 
         # lambertw typically returns complex value with zero imaginary part
         # may overflow to np.inf
@@ -2274,22 +2297,39 @@ def adrinverter(v_dc, p_dc, inverter, vtol=0.10):
     mppt_hi = inverter['MPPTHi']
     mppt_low = inverter['MPPTLow']
 
-    v_lim_upper = np.nanmax([v_max, vdc_max, mppt_hi])*(1+vtol)
-    v_lim_lower = np.nanmax([v_min, mppt_low])*(1-vtol)
-
-    pdc = p_dc/p_nom
-    vdc = v_dc/v_nom
-    poly = np.array([pdc**0, pdc, pdc**2, vdc-1, pdc*(vdc-1),
-                     pdc**2*(vdc-1), 1/vdc-1, pdc*(1./vdc-1),
-                     pdc**2*(1./vdc-1)])
+    v_lim_upper = np.nanmax([v_max, vdc_max, mppt_hi]) * (1 + vtol)
+    v_lim_lower = np.nanmax([v_min, mppt_low]) * (1 - vtol)
+
+    pdc = p_dc / p_nom
+    vdc = v_dc / v_nom
+    # zero voltage will lead to division by zero, but since power is
+    # set to night time value later, these errors can be safely ignored
+    with np.errstate(invalid='ignore', divide='ignore'):
+        poly = np.array([pdc**0,  # replace with np.ones_like?
+                         pdc,
+                         pdc**2,
+                         vdc - 1,
+                         pdc * (vdc - 1),
+                         pdc**2 * (vdc - 1),
+                         1. / vdc - 1,  # divide by 0
+                         pdc * (1. / vdc - 1),  # invalid 0./0. --> nan
+                         pdc**2 * (1. / vdc - 1)])  # divide by 0
     p_loss = np.dot(np.array(ce_list), poly)
     ac_power = p_nom * (pdc-p_loss)
-    p_nt = -1*np.absolute(p_nt)
+    p_nt = -1 * np.absolute(p_nt)
+
+    # set output to nan where input is outside of limits
+    # errstate silences case where input is nan
+    with np.errstate(invalid='ignore'):
+        invalid = (v_lim_upper < v_dc) | (v_dc < v_lim_lower)
+    ac_power = np.where(invalid, np.nan, ac_power)
+
+    # set night values
+    ac_power = np.where(vdc == 0, p_nt, ac_power)
+    ac_power = np.maximum(ac_power, p_nt)
 
-    ac_power = np.where((v_lim_upper < v_dc) | (v_dc < v_lim_lower),
-                        np.nan, ac_power)
-    ac_power = np.where((ac_power < p_nt) | (vdc == 0), p_nt, ac_power)
-    ac_power = np.where(ac_power > pac_max, pac_max, ac_power)
+    # set max ac output
+    ac_power = np.minimum(ac_power, pac_max)
 
     if isinstance(p_dc, pd.Series):
         ac_power = pd.Series(ac_power, index=pdc.index)

pvlib/test/test_pvsystem.py

@@ -101,6 +101,18 @@ def test_ashraeiam():
     assert_allclose(iam, expected, equal_nan=True)
 
 
+@needs_numpy_1_10
+def test_ashraeiam_scalar():
+    thetas = -45.
+    iam = pvsystem.ashraeiam(thetas, .05)
+    expected = 0.97928932
+    assert_allclose(iam, expected, equal_nan=True)
+    thetas = np.nan
+    iam = pvsystem.ashraeiam(thetas, .05)
+    expected = np.nan
+    assert_allclose(iam, expected, equal_nan=True)
+
+
 @needs_numpy_1_10
 def test_PVSystem_ashraeiam():
     module_parameters = pd.Series({'b': 0.05})
@@ -127,6 +139,18 @@ def test_physicaliam():
     assert_series_equal(iam, expected)
 
 
+@needs_numpy_1_10
+def test_physicaliam_scalar():
+    aoi = -45.
+    iam = pvsystem.physicaliam(aoi, 1.526, 0.002, 4)
+    expected = 0.98797788
+    assert_allclose(iam, expected, equal_nan=True)
+    aoi = np.nan
+    iam = pvsystem.physicaliam(aoi, 1.526, 0.002, 4)
+    expected = np.nan
+    assert_allclose(iam, expected, equal_nan=True)
+
+
 @needs_numpy_1_10
 def test_PVSystem_physicaliam():
     module_parameters = pd.Series({'K': 4, 'L': 0.002, 'n': 1.526})
@@ -815,6 +839,16 @@ def test_adrinverter_float(sam_data):
     assert_allclose(pacs, 1161.5745)
 
 
+def test_adrinverter_invalid_and_night(sam_data):
+    inverters = sam_data['adrinverter']
+    testinv = 'Zigor__Sunzet_3_TL_US_240V__CEC_2011_'
+    vdcs = np.array([39.873036, 0., np.nan, 420])
+    pdcs = np.array([188.09182, 0., 420, np.nan])
+
+    pacs = pvsystem.adrinverter(vdcs, pdcs, inverters[testinv])
+    assert_allclose(pacs, np.array([np.nan, -0.25, np.nan, np.nan]))
+
+
 def test_snlinverter(sam_data):
     inverters = sam_data['cecinverter']
     testinv = 'ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_'

setup.py

@@ -40,7 +40,7 @@
 MAINTAINER_EMAIL = '[email protected]'
 URL = 'https://github.com/pvlib/pvlib-python'
 
-INSTALL_REQUIRES = ['numpy >= 1.9.0',
+INSTALL_REQUIRES = ['numpy >= 1.10.1',
                     'pandas >= 0.14.0',
                     'pytz',
                     'six',