fix dropped nighttime values in perez, propagate airmass to perez in ModelChain (#192)

* send am through mc to perez function. fix perez nans

* really fix the nan/dropped data problem in perez

* update whatsnew

* rebase on master

Author: wholmgren

docs/sphinx/source/whatsnew/v0.3.3.txt

@@ -37,6 +37,12 @@ Bug fixes
   dict-like object. (:issue:`157`)
 * Fix numpy 1.11 deprecation warnings caused by some functions using
   non-integer indices.
+* Propagate airmass data through ModelChain's get_irradiance call so that
+  the perez model can use it, if necessary. (:issue:`172`)
+* Fix problem in which the perez function dropped nighttime values.
+  Nighttime values are now set to 0.
+  (:issue:`191`)
+
 
 Documentation
 ~~~~~~~~~~~~~

docs/tutorials/irradiance.ipynb

@@ -967,32 +967,30 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
           solar_zenith, solar_azimuth, airmass,
           modelt='allsitescomposite1990'):
     '''
-    Determine diffuse irradiance from the sky on a tilted surface using one of
-    the Perez models.
+    Determine diffuse irradiance from the sky on a tilted surface using
+    one of the Perez models.
 
     Perez models determine the diffuse irradiance from the sky (ground
     reflected irradiance is not included in this algorithm) on a tilted
     surface using the surface tilt angle, surface azimuth angle, diffuse
     horizontal irradiance, direct normal irradiance, extraterrestrial
     irradiance, sun zenith angle, sun azimuth angle, and relative (not
-    pressure-corrected) airmass. Optionally a selector may be used to use
-    any of Perez's model coefficient sets.
+    pressure-corrected) airmass. Optionally a selector may be used to
+    use any of Perez's model coefficient sets.
 
 
     Parameters
     ----------
 
     surface_tilt : float or Series
-        Surface tilt angles in decimal degrees.
-        surface_tilt must be >=0 and <=180. The tilt angle is defined as
-        degrees from horizontal (e.g. surface facing up = 0, surface facing
-        horizon = 90)
+        Surface tilt angles in decimal degrees. surface_tilt must be >=0
+        and <=180. The tilt angle is defined as degrees from horizontal
+        (e.g. surface facing up = 0, surface facing horizon = 90)
 
     surface_azimuth : float or Series
-        Surface azimuth angles in decimal degrees.
-        surface_azimuth must be >=0 and <=360. The Azimuth convention is defined
-        as degrees east of north (e.g. North = 0, South=180 East = 90,
-        West = 270).
+        Surface azimuth angles in decimal degrees. surface_azimuth must
+        be >=0 and <=360. The Azimuth convention is defined as degrees
+        east of north (e.g. North = 0, South=180 East = 90, West = 270).
 
     dhi : float or Series
         Diffuse horizontal irradiance in W/m^2.
@@ -1011,22 +1009,20 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
         solar_zenith must be >=0 and <=180.
 
     solar_azimuth : float or Series
-        Sun azimuth angles in decimal degrees.
-        solar_azimuth must be >=0 and <=360. The Azimuth convention is defined
-        as degrees east of north (e.g. North = 0, East = 90, West = 270).
+        Sun azimuth angles in decimal degrees. solar_azimuth must be >=0
+        and <=360. The Azimuth convention is defined as degrees east of
+        north (e.g. North = 0, East = 90, West = 270).
 
     airmass : float or Series
-        relative (not pressure-corrected) airmass
-        values. If AM is a DataFrame it must be of the same size as all other
-        DataFrame inputs. AM must be >=0 (careful using the 1/sec(z) model of
-        AM generation)
+        Relative (not pressure-corrected) airmass values. If AM is a
+        DataFrame it must be of the same size as all other DataFrame
+        inputs. AM must be >=0 (careful using the 1/sec(z) model of AM
+        generation)
 
     model : string (optional, default='allsitescomposite1990')
-
-        A string which selects the desired set of Perez
-        coefficients. If model is not provided as an input, the default,
-        '1990' will be used.
-        All possible model selections are:
+        A string which selects the desired set of Perez coefficients. If
+        model is not provided as an input, the default, '1990' will be
+        used. All possible model selections are:
 
         * '1990'
         * 'allsitescomposite1990' (same as '1990')
@@ -1046,44 +1042,42 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
 
     float or Series
 
-        The diffuse component of the solar radiation  on an
-        arbitrarily tilted surface defined by the Perez model as given in
-        reference [3].
-        SkyDiffuse is the diffuse component ONLY and does not include the
-        ground reflected irradiance or the irradiance due to the beam.
+        The diffuse component of the solar radiation on an arbitrarily
+        tilted surface defined by the Perez model as given in reference
+        [3]. SkyDiffuse is the diffuse component ONLY and does not
+        include the ground reflected irradiance or the irradiance due to
+        the beam.
 
 
     References
     ----------
 
-    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute
-    solar irradiance on inclined surfaces for building energy simulation"
-    2007, Solar Energy vol. 81. pp. 254-267
+    [1] Loutzenhiser P.G. et. al. "Empirical validation of models to
+    compute solar irradiance on inclined surfaces for building energy
+    simulation" 2007, Solar Energy vol. 81. pp. 254-267
 
-    [2] Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D., 1987.
-    A new simplified version of the Perez diffuse irradiance model for tilted
-    surfaces. Solar Energy 39(3), 221-232.
+    [2] Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D.,
+    1987. A new simplified version of the Perez diffuse irradiance model
+    for tilted surfaces. Solar Energy 39(3), 221-232.
 
-    [3] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R., 1990.
-    Modeling daylight availability and irradiance components from direct
-    and global irradiance. Solar Energy 44 (5), 271-289.
+    [3] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R.,
+    1990. Modeling daylight availability and irradiance components from
+    direct and global irradiance. Solar Energy 44 (5), 271-289.
 
     [4] Perez, R. et. al 1988. "The Development and Verification of the
     Perez Diffuse Radiation Model". SAND88-7030
     '''
 
-    pvl_logger.debug('diffuse_sky.perez()')
-
     kappa = 1.041  # for solar_zenith in radians
     z = np.radians(solar_zenith)  # convert to radians
 
     # epsilon is the sky's "clearness"
     eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
 
-    # Perez et al define clearness bins according to the following 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
+    # Perez et al define clearness bins according to the following
+    # 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 = eps.copy()
     ebin[(eps < 1.065)] = 1
     ebin[(eps >= 1.065) & (eps < 1.23)] = 2
@@ -1100,31 +1094,29 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     ebin = ebin.dropna().astype(int)
 
     # This is added because in cases where the sun is below the horizon
-    # (var.solar_zenith > 90) but there is still diffuse horizontal light
-    # (var.DHI>0), it is possible that the airmass (var.AM) could be NaN, which
-    # messes up later calculations. Instead, if the sun is down, and there is
-    # still var.DHI, we set the airmass to the airmass value on the horizon
-    # (approximately 37-38).
+    # (var.solar_zenith > 90) but there is still diffuse horizontal
+    # light (var.DHI>0), it is possible that the airmass (var.AM) could
+    # be NaN, which messes up later calculations. Instead, if the sun is
+    # down, and there is still var.DHI, we set the airmass to the
+    # airmass value on the horizon (approximately 37-38).
     # var.AM(var.solar_zenith >=90 & var.DHI >0) = 37;
 
     # var.DNI_ET[var.DNI_ET==0] = .00000001 #very hacky, fix this
 
     # delta is the sky's "brightness"
     delta = dhi * airmass / dni_extra
 
-    # keep only valid times
-    delta = delta[ebin.index]
-    z = z[ebin.index]
-
     # The various possible sets of Perez coefficients are contained
     # in a subfunction to clean up the code.
     F1c, F2c = _get_perez_coefficients(modelt)
 
-    F1 = F1c[ebin, 0] + F1c[ebin, 1] * delta + F1c[ebin, 2] * z
+    F1 = (F1c[ebin, 0] + F1c[ebin, 1] * delta[ebin.index] +
+          F1c[ebin, 2] * z[ebin.index])
     F1[F1 < 0] = 0
     F1 = F1.astype(float)
 
-    F2 = F2c[ebin, 0] + F2c[ebin, 1] * delta + F2c[ebin, 2] * z
+    F2 = (F2c[ebin, 0] + F2c[ebin, 1] * delta[ebin.index] +
+          F2c[ebin, 2] * z[ebin.index])
     F2[F2 < 0] = 0
     F2 = F2.astype(float)
 
@@ -1141,8 +1133,9 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     term2 = F1 * A[ebin.index] / B[ebin.index]
     term3 = F2 * tools.sind(surface_tilt)
 
-    sky_diffuse = dhi[ebin.index] * (term1 + term2 + term3)
+    sky_diffuse = dhi * (term1 + term2 + term3)
     sky_diffuse[sky_diffuse < 0] = 0
+    sky_diffuse[airmass.isnull()] = 0
 
     return sky_diffuse
 

pvlib/irradiance.py

@@ -268,12 +268,12 @@ def __init__(self, system, location,
 
         # calls setter
         self.orientation_strategy = orientation_strategy
-        
+
     def __repr__(self):
-        return ('ModelChain for: '+ str(self.system) + 
+        return ('ModelChain for: '+ str(self.system) +
                 ' orientation_startegy: ' + str(self.orientation_strategy) +
                 ' clearsky_model: ' + str(self.clearsky_model) +
-                'transposition_model: ' + str(self.transposition_model) + 
+                'transposition_model: ' + str(self.transposition_model) +
                 ' solar_position_method: ' + str(self.solar_position_method) +
                 'airmass_model: ' + str(self.airmass_model))
 
@@ -360,7 +360,8 @@ def run_model(self, times, irradiance=None, weather=None):
             self.irradiance['dni'],
             self.irradiance['ghi'],
             self.irradiance['dhi'],
-            model=self.transposition_model)
+            model=self.transposition_model,
+            airmass=self.airmass['airmass_relative'])
 
         if weather is None:
             weather = {'wind_speed': 0, 'temp_air': 20}

pvlib/modelchain.py

@@ -1,8 +1,6 @@
 import logging
 pvl_logger = logging.getLogger('pvlib')
 
-import datetime
-
 import numpy as np
 import pandas as pd
 
@@ -20,15 +18,14 @@
 from . import requires_ephem
 
 # setup times and location to be tested.
-times = pd.date_range(start=datetime.datetime(2014, 6, 24),
-                      end=datetime.datetime(2014, 6, 26), freq='1Min')
-
 tus = Location(32.2, -111, 'US/Arizona', 700)
 
-times_localized = times.tz_localize(tus.tz)
+# must include night values
+times = pd.date_range(start='20140624', end='20140626', freq='1Min',
+                      tz=tus.tz)
 
-ephem_data = solarposition.get_solarposition(times, tus.latitude,
-                                             tus.longitude, method='nrel_numpy')
+ephem_data = solarposition.get_solarposition(
+    times, tus.latitude, tus.longitude, method='nrel_numpy')
 
 irrad_data = clearsky.ineichen(times, tus.latitude, tus.longitude,
                                altitude=tus.altitude, linke_turbidity=3,
@@ -141,10 +138,11 @@ def test_king():
 
 
 def test_perez():
-    AM = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
-    irradiance.perez(40, 180, irrad_data['dhi'], irrad_data['dni'],
+    am = atmosphere.relativeairmass(ephem_data['apparent_zenith'])
+    out = irradiance.perez(40, 180, irrad_data['dhi'], irrad_data['dni'],
                      dni_et, ephem_data['apparent_zenith'],
-                     ephem_data['azimuth'], AM)
+                     ephem_data['azimuth'], am)
+    assert not out.isnull().any()
 
 # klutcher (misspelling) will be removed in 0.3
 def test_total_irrad():

pvlib/test/test_irradiance.py

@@ -90,6 +90,33 @@ def test_run_model_with_irradiance():
     assert_series_equal(ac, expected)
 
 
+def test_run_model_perez():
+    system, location = mc_setup()
+    mc = ModelChain(system, location, transposition_model='perez')
+    times = pd.date_range('20160101 1200-0700', periods=2, freq='6H')
+    irradiance = pd.DataFrame({'dni':900, 'ghi':600, 'dhi':150},
+                              index=times)
+    ac = mc.run_model(times, irradiance=irradiance).ac
+
+    expected = pd.Series(np.array([  190.194545796,  -2.00000000e-02]),
+                         index=times)
+    assert_series_equal(ac, expected)
+
+
+def test_run_model_gueymard_perez():
+    system, location = mc_setup()
+    mc = ModelChain(system, location, airmass_model='gueymard1993',
+                    transposition_model='perez')
+    times = pd.date_range('20160101 1200-0700', periods=2, freq='6H')
+    irradiance = pd.DataFrame({'dni':900, 'ghi':600, 'dhi':150},
+                              index=times)
+    ac = mc.run_model(times, irradiance=irradiance).ac
+
+    expected = pd.Series(np.array([  190.194760203,  -2.00000000e-02]),
+                         index=times)
+    assert_series_equal(ac, expected)
+
+
 def test_run_model_with_weather():
     system, location = mc_setup()
     mc = ModelChain(system, location)
@@ -220,8 +247,8 @@ def test_basic_chain_altitude_pressure():
     expected = pd.Series(np.array([  1.15771428788e+02,  -2.00000000e-02]),
                          index=times)
     assert_series_equal(ac, expected)
-    
-    
+
+
 def test_ModelChain___repr__():
     system = PVSystem()
     location = Location(32.2, -111, altitude=700)