diff --git a/.travis.yml b/.travis.yml index eee4301467..66dcd4c1cd 100644 --- a/.travis.yml +++ b/.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 diff --git a/docs/sphinx/source/whatsnew/v0.5.2.rst b/docs/sphinx/source/whatsnew/v0.5.2.rst index b02691448f..eab450fcba 100644 --- a/docs/sphinx/source/whatsnew/v0.5.2.rst +++ b/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 ~~~~~~~~~~~~~ diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py index 3ca7c181bc..b43849e311 100644 --- a/pvlib/clearsky.py +++ b/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 diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 49b51a9088..b8ce4e5cb2 100755 --- a/pvlib/irradiance.py +++ b/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) diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py index 40ef7c2f0e..4517d5806b 100644 --- a/pvlib/pvsystem.py +++ b/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) diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py index 7d8dbf8b32..f5cb78f651 100644 --- a/pvlib/test/test_pvsystem.py +++ b/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_' diff --git a/setup.py b/setup.py index 62789514d8..26c7612aac 100755 --- a/setup.py +++ b/setup.py @@ -40,7 +40,7 @@ MAINTAINER_EMAIL = 'holmgren@email.arizona.edu' 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',