Merge pull request #1 from wholmgren/int_indexing

make dirint work with nans again

Author: alorenzo175

pvlib/irradiance.py

@@ -1313,12 +1313,12 @@ def _get_perez_coefficients(perezmodelt):
 
 def disc(ghi, zenith, times, pressure=101325):
     '''
-    Estimate Direct Normal Irradiance from Global Horizontal Irradiance 
+    Estimate Direct Normal Irradiance from Global Horizontal Irradiance
     using the DISC model.
 
     The DISC algorithm converts global horizontal irradiance to direct
     normal irradiance through empirical relationships between the global
-    and direct clearness indices. 
+    and direct clearness indices.
 
     Parameters
     ----------
@@ -1327,62 +1327,62 @@ def disc(ghi, zenith, times, pressure=101325):
         Global horizontal irradiance in W/m^2.
 
     solar_zenith : Series
-        True (not refraction - corrected) solar zenith 
-        angles in decimal degrees. 
+        True (not refraction - corrected) solar zenith
+        angles in decimal degrees.
 
     times : DatetimeIndex
 
     pressure : float or Series
         Site pressure in Pascal.
 
-    Returns   
+    Returns
     -------
     DataFrame with the following keys:
-        * ``dni``: The modeled direct normal irradiance 
+        * ``dni``: The modeled direct normal irradiance
           in W/m^2 provided by the
-          Direct Insolation Simulation Code (DISC) model. 
-        * ``kt``: Ratio of global to extraterrestrial 
+          Direct Insolation Simulation Code (DISC) model.
+        * ``kt``: Ratio of global to extraterrestrial
           irradiance on a horizontal plane.
         * ``airmass``: Airmass
 
     References
     ----------
 
-    [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly 
-    Global Horizontal to Direct Normal Insolation", Technical 
-    Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research 
+    [1] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
+    Global Horizontal to Direct Normal Insolation", Technical
+    Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
     Institute, 1987.
 
-    [2] J.W. "Fourier series representation of the position of the sun". 
+    [2] J.W. "Fourier series representation of the position of the sun".
     Found at:
     http://www.mail-archive.com/[email protected]/msg01050.html on
     January 12, 2012
 
-    See Also 
-    -------- 
-    atmosphere.alt2pres 
+    See Also
+    --------
+    atmosphere.alt2pres
     dirint
     '''
 
     pvl_logger.debug('clearsky.disc')
-    
+
     temp = pd.DataFrame(index=times, columns=['A','B','C'], dtype=float)
 
     doy = times.dayofyear
-    
+
     DayAngle = 2. * np.pi*(doy - 1) / 365
-    
+
     re = (1.00011 + 0.034221*np.cos(DayAngle) + 0.00128*np.sin(DayAngle)
           + 0.000719*np.cos(2.*DayAngle) + 7.7e-05*np.sin(2.*DayAngle) )
-          
+
     I0 = re * 1370.
     I0h = I0 * np.cos(np.radians(zenith))
-    
+
     Ztemp = zenith.copy()
     Ztemp[zenith > 87] = np.NaN
-    
+
     AM = 1.0 / ( np.cos(np.radians(Ztemp)) + 0.15*( (93.885 - Ztemp)**(-1.253) ) ) * (pressure / 101325)
-    
+
     Kt = ghi / I0h
     Kt[Kt < 0] = 0
     Kt[Kt > 2] = np.NaN
@@ -1395,10 +1395,10 @@ def disc(ghi, zenith, times, pressure=101325):
     temp.C[Kt <= 0.6] = -0.28 + 0.932*(Kt[Kt <= 0.6]) - 2.048*(Kt[Kt <= 0.6] ** 2)
 
     delKn = temp.A + temp.B * np.exp(temp.C*AM)
-   
+
     Knc = 0.866 - 0.122*(AM) + 0.0121*(AM ** 2) - 0.000653*(AM ** 3) + 1.4e-05*(AM ** 4)
     Kn = Knc - delKn
-    
+
     dni = Kn * I0
 
     dni[zenith > 87] = np.NaN
@@ -1409,14 +1409,14 @@ def disc(ghi, zenith, times, pressure=101325):
     dfout['airmass'] = AM
 
     return dfout
-    
 
-def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True, 
+
+def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
            temp_dew=None):
     """
-    Determine DNI from GHI using the DIRINT modification 
+    Determine DNI from GHI using the DIRINT modification
     of the DISC model.
-    
+
     Implements the modified DISC model known as "DIRINT" introduced in [1].
     DIRINT predicts direct normal irradiance (DNI) from measured global
     horizontal irradiance (GHI). DIRINT improves upon the DISC model by
@@ -1425,22 +1425,22 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
     information provided.
 
     Parameters
-    ----------  
+    ----------
     ghi : pd.Series
-        Global horizontal irradiance in W/m^2. 
-    
+        Global horizontal irradiance in W/m^2.
+
     zenith : pd.Series
         True (not refraction-corrected) zenith
         angles in decimal degrees. If Z is a vector it must be of the
         same size as all other vector inputs. Z must be >=0 and <=180.
-    
+
     times : DatetimeIndex
-        
+
     pressure : float or pd.Series
-        The site pressure in Pascal. 
-        Pressure may be measured or an average pressure may be 
+        The site pressure in Pascal.
+        Pressure may be measured or an average pressure may be
         calculated from site altitude.
-    
+
     use_delta_kt_prime : bool
         Indicates if the user would like to
         utilize the time-series nature of the GHI measurements. A value of ``False``
@@ -1450,13 +1450,13 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
         than 1.5 hours. If none of the input arguments are
         vectors, then time-series improvements are not used (because it's not
         a time-series).
-    
-    temp_dew : None, float, or pd.Series 
-        Surface dew point temperatures, in degrees C. 
+
+    temp_dew : None, float, or pd.Series
+        Surface dew point temperatures, in degrees C.
         Values of temp_dew may be numeric or NaN. Any
         single time period point with a DewPtTemp=NaN does not have dew point
-        improvements applied. If DewPtTemp is not provided, then dew point 
-        improvements are not applied.  
+        improvements applied. If DewPtTemp is not provided, then dew point
+        improvements are not applied.
 
     Returns
     -------
@@ -1467,57 +1467,58 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
     References
     ----------
     [1] Perez, R., P. Ineichen, E. Maxwell, R. Seals and A. Zelenka, (1992).
-    "Dynamic Global-to-Direct Irradiance Conversion Models".  ASHRAE 
+    "Dynamic Global-to-Direct Irradiance Conversion Models".  ASHRAE
     Transactions-Research Series, pp. 354-369
 
-    [2] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly 
-    Global Horizontal to Direct Normal Insolation", Technical 
-    Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research 
+    [2] Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
+    Global Horizontal to Direct Normal Insolation", Technical
+    Report No. SERI/TR-215-3087, Golden, CO: Solar Energy Research
     Institute, 1987.
 
     DIRINT model requires time series data (ie. one of the inputs must be a
     vector of length >2.
     """
-    
+
     pvl_logger.debug('clearsky.dirint')
-    
+
     disc_out = disc(ghi, zenith, times)
     kt = disc_out['kt']
-    
+
     # Absolute Airmass, per the DISC model
     # Note that we calculate the AM pressure correction slightly differently
     # than Perez. He uses altitude, we use pressure (which we calculate
     # slightly differently)
-    airmass = (1./(tools.cosd(zenith) + 0.15*((93.885-zenith)**(-1.253))) * 
+    airmass = (1./(tools.cosd(zenith) + 0.15*((93.885-zenith)**(-1.253))) *
                pressure/101325)
-    
+
     coeffs = _get_dirint_coeffs()
-    
+
     kt_prime = kt / (1.031 * np.exp(-1.4/(0.9+9.4/airmass)) + 0.1)
     kt_prime[kt_prime > 0.82] = 0.82 # From SRRL code. consider np.NaN
-    kt_prime.fillna(0, inplace=True)
     pvl_logger.debug('kt_prime:\n%s', kt_prime)
-    
-    # wholmgren: 
+
+    # wholmgren:
     # the use_delta_kt_prime statement is a port of the MATLAB code.
     # I am confused by the abs() in the delta_kt_prime calculation.
     # It is not the absolute value of the central difference.
     if use_delta_kt_prime:
         delta_kt_prime = 0.5*( (kt_prime - kt_prime.shift(1)).abs()
                               .add(
-                               (kt_prime - kt_prime.shift(-1)).abs(), 
+                               (kt_prime - kt_prime.shift(-1)).abs(),
                                    fill_value=0))
     else:
         delta_kt_prime = pd.Series(-1, index=times)
-    
+
     if temp_dew is not None:
         w = pd.Series(np.exp(0.07 * temp_dew - 0.075), index=times)
     else:
         w = pd.Series(-1, index=times)
-    
+
     # @wholmgren: the following bin assignments use MATLAB's 1-indexing.
     # Later, we'll subtract 1 to conform to Python's 0-indexing.
-    
+    # The initial filled 0s will result in nans.
+    # see discussion at https://github.com/pvlib/pvlib-python/pull/164
+
     # Create kt_prime bins
     kt_prime_bin = pd.Series(0, index=times, dtype=np.int64)
     kt_prime_bin[(kt_prime>=0) & (kt_prime<0.24)] = 1
@@ -1527,7 +1528,7 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
     kt_prime_bin[(kt_prime>=0.7) & (kt_prime<0.8)] = 5
     kt_prime_bin[(kt_prime>=0.8) & (kt_prime<=1)] = 6
     pvl_logger.debug('kt_prime_bin:\n%s', kt_prime_bin)
-    
+
     # Create zenith angle bins
     zenith_bin = pd.Series(0, index=times, dtype=np.int64)
     zenith_bin[(zenith>=0) & (zenith<25)] = 1
@@ -1537,7 +1538,7 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
     zenith_bin[(zenith>=70) & (zenith<80)] = 5
     zenith_bin[(zenith>=80)] = 6
     pvl_logger.debug('zenith_bin:\n%s', zenith_bin)
-    
+
     # Create the bins for w based on dew point temperature
     w_bin = pd.Series(0, index=times, dtype=np.int64)
     w_bin[(w>=0) & (w<1)] = 1
@@ -1557,11 +1558,16 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
     delta_kt_prime_bin[(delta_kt_prime>=0.3) & (delta_kt_prime<=1)] = 6
     delta_kt_prime_bin[delta_kt_prime == -1] = 7
     pvl_logger.debug('delta_kt_prime_bin:\n%s', delta_kt_prime_bin)
-    
+
     # subtract 1 to account for difference between MATLAB-style bin
     # assignment and Python-style array lookup.
     dirint_coeffs = coeffs[kt_prime_bin-1, zenith_bin-1,
                            delta_kt_prime_bin-1, w_bin-1]
+    # convert unassigned bins to nan
+    # use where to avoid boolean indexing deprecation
+    dirint_coeffs[np.where((kt_prime_bin == 0) | (zenith_bin == 0) |
+                           (w_bin == 0) | (delta_kt_prime_bin == 0))] = np.nan
+
     dni = disc_out['dni'] * dirint_coeffs
 
     dni.name = 'dni'
@@ -1572,19 +1578,19 @@ def dirint(ghi, zenith, times, pressure=101325, use_delta_kt_prime=True,
 def _get_dirint_coeffs():
     """
     A place to stash the dirint coefficients.
-    
+
     Returns
     -------
     np.array with shape ``(6, 6, 7, 5)``.
     Ordering is ``[kt_prime_bin, zenith_bin, delta_kt_prime_bin, w_bin]``
     """
 
-    # To allow for maximum copy/paste from the MATLAB 1-indexed code, 
+    # To allow for maximum copy/paste from the MATLAB 1-indexed code,
     # we create and assign values to an oversized array.
     # Then, we return the [1:, 1:, :, :] slice.
-    
+
     coeffs = np.zeros((7,7,7,5))
-    
+
     coeffs[1,1,:,:] = [
         [0.385230, 0.385230, 0.385230, 0.462880, 0.317440],
         [0.338390, 0.338390, 0.221270, 0.316730, 0.503650],
@@ -1908,6 +1914,6 @@ def _get_dirint_coeffs():
         [0.570000, 0.550000, 0.598800, 0.400000, 0.560150],
         [0.475230, 0.500000, 0.518640, 0.339970, 0.520230],
         [0.743440, 0.592190, 0.603060, 0.316930, 0.794390 ]]
-        
+
     return coeffs[1:,1:,:,:]
 

pvlib/test/test_irradiance.py

@@ -152,14 +152,14 @@ def test_total_irrad():
 
     for model in models:
         total = irradiance.total_irrad(
-            32, 180, 
+            32, 180,
             ephem_data['apparent_zenith'], ephem_data['azimuth'],
             dni=irrad_data['dni'], ghi=irrad_data['ghi'],
             dhi=irrad_data['dhi'],
             dni_extra=dni_et, airmass=AM,
             model=model,
             surface_type='urban')
-        
+
         assert total.columns.tolist() == ['poa_global', 'poa_direct',
                                           'poa_diffuse', 'poa_sky_diffuse',
                                           'poa_ground_diffuse']
@@ -181,7 +181,7 @@ def test_globalinplane():
 def test_disc_keys():
     clearsky_data = clearsky.ineichen(times, tus.latitude, tus.longitude,
                                       linke_turbidity=3)
-    disc_data = irradiance.disc(clearsky_data['ghi'], ephem_data['zenith'], 
+    disc_data = irradiance.disc(clearsky_data['ghi'], ephem_data['zenith'],
                                 ephem_data.index)
     assert 'dni' in disc_data.columns
     assert 'kt' in disc_data.columns
@@ -202,7 +202,7 @@ def test_dirint():
     clearsky_data = clearsky.ineichen(times, tus.latitude, tus.longitude,
                                       linke_turbidity=3)
     pressure = 93193.
-    dirint_data = irradiance.dirint(clearsky_data['ghi'], ephem_data['zenith'], 
+    dirint_data = irradiance.dirint(clearsky_data['ghi'], ephem_data['zenith'],
                                     ephem_data.index, pressure=pressure)
 
 def test_dirint_value():
@@ -214,6 +214,19 @@ def test_dirint_value():
     assert_almost_equal(dirint_data.values,
                         np.array([928.85, 688.26]), 1)
 
+
+def test_dirint_nans():
+    times = pd.DatetimeIndex(start='2014-06-24T12-0700', periods=5, freq='6H')
+    ghi = pd.Series([np.nan, 1038.62, 1038.62, 1038.62, 1038.62], index=times)
+    zenith = pd.Series([10.567, np.nan, 10.567, 10.567, 10.567,], index=times)
+    pressure = pd.Series([93193., 93193., np.nan, 93193., 93193.], index=times)
+    temp_dew = pd.Series([10, 10, 10, np.nan, 10], index=times)
+    dirint_data = irradiance.dirint(ghi, zenith, times, pressure=pressure,
+                                    temp_dew=temp_dew)
+    assert_almost_equal(dirint_data.values,
+                        np.array([np.nan, np.nan, np.nan, np.nan, 934.2]), 1)
+
+
 def test_dirint_tdew():
     times = pd.DatetimeIndex(['2014-06-24T12-0700','2014-06-24T18-0700'])
     ghi = pd.Series([1038.62, 254.53], index=times)
@@ -238,4 +251,4 @@ def test_dirint_coeffs():
     coeffs = irradiance._get_dirint_coeffs()
     assert coeffs[0,0,0,0] == 0.385230
     assert coeffs[0,1,2,1] == 0.229970
-    assert coeffs[3,2,6,3] == 1.032260
+    assert coeffs[3,2,6,3] == 1.032260