refactor ineichen. good riddance

pvlib · wholmgren · Jul 7, 2016 · Jun 24, 2016 · Jul 6, 2016 · Jul 6, 2016
commit b0bfc0d192b3d9bc3d3df9dba8528e1ff16d9931
diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py
@@ -20,69 +20,45 @@
 from pvlib import solarposition
 
 
-def ineichen(time, latitude, longitude, altitude=0, linke_turbidity=None,
-             solarposition_method='nrel_numpy', zenith_data=None,
-             airmass_model='young1994', airmass_data=None,
-             interp_turbidity=True):
+def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
+             dni_extra=1364., altitude=0):
     '''
-    Determine clear sky GHI, DNI, and DHI from Ineichen/Perez model
+    Determine clear sky GHI, DNI, and DHI from Ineichen/Perez model.
 
-    Implements the Ineichen and Perez clear sky model for global horizontal
-    irradiance (GHI), direct normal irradiance (DNI), and calculates
-    the clear-sky diffuse horizontal (DHI) component as the difference
-    between GHI and DNI*cos(zenith) as presented in [1, 2]. A report on clear
-    sky models found the Ineichen/Perez model to have excellent performance
-    with a minimal input data set [3].
+    Implements the Ineichen and Perez clear sky model for global
+    horizontal irradiance (GHI), direct normal irradiance (DNI), and
+    calculates the clear-sky diffuse horizontal (DHI) component as the
+    difference between GHI and DNI*cos(zenith) as presented in [1, 2]. A
+    report on clear sky models found the Ineichen/Perez model to have
+    excellent performance with a minimal input data set [3].
 
     Default values for montly Linke turbidity provided by SoDa [4, 5].
 
     Parameters
     -----------
-    time : pandas.DatetimeIndex
-
-    latitude : float
-
-    longitude : float
+    apparent_zenith: numeric
 
-    altitude : float
+    airmass_absolute: numeric
 
-    linke_turbidity : None or float
-        If None, uses ``LinkeTurbidities.mat`` lookup table.
+    linke_turbidity: numeric
 
-    solarposition_method : string
-        Sets the solar position algorithm.
-        See solarposition.get_solarposition()
-
-    zenith_data : None or Series
-        If None, ephemeris data will be calculated using ``solarposition_method``.
-
-    airmass_model : string
-        See pvlib.airmass.relativeairmass().
+    dni_extra: numeric
+        Extraterrestrial irradiance. The units of ``dni_extra``
+        determine the units of the output.
 
-    airmass_data : None or Series
-        If None, absolute air mass data will be calculated using
-        ``airmass_model`` and location.alitude.
-
-    interp_turbidity : bool
-        If ``True``, interpolates the monthly Linke turbidity values
-        found in ``LinkeTurbidities.mat`` to daily values.
+    altitude: numeric
 
     Returns
    --------
-    DataFrame with the following columns: ``ghi, dni, dhi``.
-
-    Notes
-    -----
-    If you are using this function
-    in a loop, it may be faster to load LinkeTurbidities.mat outside of
-    the loop and feed it in as a keyword argument, rather than
-    having the function open and process the file each time it is called.
+    clearsky : DataFrame (if Series input) or OrderedDict of arrays
+        DataFrame/OrderedDict contains the columns/keys
+        ``'dhi', 'dni', 'ghi'``.
 
     References
     ----------
-
     [1] P. Ineichen and R. Perez, "A New airmass independent formulation for
-        the Linke turbidity coefficient", Solar Energy, vol 73, pp. 151-157, 2002.
+        the Linke turbidity coefficient", Solar Energy, vol 73, pp. 151-157,
+        2002.
 
     [2] R. Perez et. al., "A New Operational Model for Satellite-Derived
         Irradiances: Description and Validation", Solar Energy, vol 73, pp.
@@ -98,97 +74,66 @@ def ineichen(time, latitude, longitude, altitude=0, linke_turbidity=None,
     [5] J. Remund, et. al., "Worldwide Linke Turbidity Information", Proc.
         ISES Solar World Congress, June 2003. Goteborg, Sweden.
     '''
-    # Initial implementation of this algorithm by Matthew Reno.
-    # Ported to python by Rob Andrews
-    # Added functionality by Will Holmgren (@wholmgren)
-
-    I0 = irradiance.extraradiation(time.dayofyear)
-
-    if zenith_data is None:
-        ephem_data = solarposition.get_solarposition(time,
-                                                     latitude=latitude,
-                                                     longitude=longitude,
-                                                     altitude=altitude,
-                                                     method=solarposition_method)
-        time = ephem_data.index # fixes issue with time possibly not being tz-aware
-        try:
-            ApparentZenith = ephem_data['apparent_zenith']
-        except KeyError:
-            ApparentZenith = ephem_data['zenith']
-            logger.warning('could not find apparent_zenith. using zenith')
-    else:
-        ApparentZenith = zenith_data
-    #ApparentZenith[ApparentZenith >= 90] = 90 # can cause problems in edge cases
-
 
-    if linke_turbidity is None:
-        TL = lookup_linke_turbidity(time, latitude, longitude,
-                                    interp_turbidity=interp_turbidity)
-    else:
-        TL = linke_turbidity
+    # Dan's note on the TL correction: By my reading of the publication
+    # on pages 151-157, Ineichen and Perez introduce (among other
+    # things) three things. 1) Beam model in eqn. 8, 2) new turbidity
+    # factor in eqn 9 and appendix A, and 3) Global horizontal model in
+    # eqn. 11. They do NOT appear to use the new turbidity factor (item
+    # 2 above) in either the beam or GHI models. The phrasing of
+    # appendix A seems as if there are two separate corrections, the
+    # first correction is used to correct the beam/GHI models, and the
+    # second correction is used to correct the revised turibidity
+    # factor. In my estimation, there is no need to correct the
+    # turbidity factor used in the beam/GHI models.
+
+    # Create the corrected TL for TL < 2
+    # TLcorr = TL;
+    # TLcorr(TL < 2) = TLcorr(TL < 2) - 0.25 .* (2-TLcorr(TL < 2)) .^ (0.5);
+
+    # This equation is found in Solar Energy 73, pg 311. Full ref: Perez
+    # et. al., Vol. 73, pp. 307-317 (2002). It is slightly different
+    # than the equation given in Solar Energy 73, pg 156. We used the
+    # equation from pg 311 because of the existence of known typos in
+    # the pg 156 publication (notably the fh2-(TL-1) should be fh2 *
+    # (TL-1)).
 
-    # Get the absolute airmass assuming standard local pressure (per
-    # alt2pres) using Kasten and Young's 1989 formula for airmass.
+    cos_zenith = tools.cosd(apparent_zenith)
 
-    if airmass_data is None:
-        AMabsolute = atmosphere.absoluteairmass(airmass_relative=atmosphere.relativeairmass(ApparentZenith, airmass_model),
-                                                pressure=atmosphere.alt2pres(altitude))
-    else:
-        AMabsolute = airmass_data
+    tl = linke_turbidity
 
     fh1 = np.exp(-altitude/8000.)
     fh2 = np.exp(-altitude/1250.)
     cg1 = 5.09e-05 * altitude + 0.868
     cg2 = 3.92e-05 * altitude + 0.0387
-    logger.debug('fh1=%s, fh2=%s, cg1=%s, cg2=%s', fh1, fh2, cg1, cg2)
-
-    #  Dan's note on the TL correction: By my reading of the publication on
-    #  pages 151-157, Ineichen and Perez introduce (among other things) three
-    #  things. 1) Beam model in eqn. 8, 2) new turbidity factor in eqn 9 and
-    #  appendix A, and 3) Global horizontal model in eqn. 11. They do NOT appear
-    #  to use the new turbidity factor (item 2 above) in either the beam or GHI
-    #  models. The phrasing of appendix A seems as if there are two separate
-    #  corrections, the first correction is used to correct the beam/GHI models,
-    #  and the second correction is used to correct the revised turibidity
-    #  factor. In my estimation, there is no need to correct the turbidity
-    #  factor used in the beam/GHI models.
-
-    #  Create the corrected TL for TL < 2
-    #  TLcorr = TL;
-    #  TLcorr(TL < 2) = TLcorr(TL < 2) - 0.25 .* (2-TLcorr(TL < 2)) .^ (0.5);
-
-    #  This equation is found in Solar Energy 73, pg 311.
-    #  Full ref: Perez et. al., Vol. 73, pp. 307-317 (2002).
-    #  It is slightly different than the equation given in Solar Energy 73, pg 156.
-    #  We used the equation from pg 311 because of the existence of known typos
-    #  in the pg 156 publication (notably the fh2-(TL-1) should be fh2 * (TL-1)).
-
-    cos_zenith = tools.cosd(ApparentZenith)
-
-    clearsky_GHI = ( cg1 * I0 * cos_zenith *
-                     np.exp(-cg2*AMabsolute*(fh1 + fh2*(TL - 1))) *
-                     np.exp(0.01*AMabsolute**1.8) )
-    clearsky_GHI[clearsky_GHI < 0] = 0
 
-    # BncI == "normal beam clear sky radiation"
+    ghi = (cg1 * dni_extra * cos_zenith *
+           np.exp(-cg2*airmass_absolute*(fh1 + fh2*(tl - 1))) *
+           np.exp(0.01*airmass_absolute**1.8))
+    ghi = np.maximum(ghi, 0)
+
+    # BncI = "normal beam clear sky radiation"
     b = 0.664 + 0.163/fh1
-    BncI = b * I0 * np.exp( -0.09 * AMabsolute * (TL - 1) )
-    logger.debug('b=%s', b)
+    BncI = b * dni_extra * np.exp(-0.09 * airmass_absolute * (tl - 1))
 
     # "empirical correction" SE 73, 157 & SE 73, 312.
-    BncI_2 = ( clearsky_GHI *
-               ( 1 - (0.1 - 0.2*np.exp(-TL))/(0.1 + 0.882/fh1) ) /
-               cos_zenith )
+    BncI_2 = (ghi *
+              (1 - (0.1 - 0.2*np.exp(-tl))/(0.1 + 0.882/fh1)) /
+              cos_zenith)
 
-    clearsky_DNI = np.minimum(BncI, BncI_2)
+    dni = np.minimum(BncI, BncI_2)
 
-    clearsky_DHI = clearsky_GHI - clearsky_DNI*cos_zenith
+    dhi = ghi - dni*cos_zenith
 
-    df_out = pd.DataFrame({'ghi':clearsky_GHI, 'dni':clearsky_DNI,
-                           'dhi':clearsky_DHI})
-    df_out.fillna(0, inplace=True)
+    irrads = OrderedDict()
+    irrads['ghi'] = ghi
+    irrads['dni'] = dni
+    irrads['dhi'] = dhi
 
-    return df_out
+    if isinstance(dni, pd.Series):
+        irrads = pd.DataFrame.from_dict(irrads)
+
+    return irrads
 
 
 def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
@@ -360,16 +305,15 @@ def simplified_solis(apparent_elevation, aod700=0.1, precipitable_water=1.,
         or 101325 and 41000 Pascals.
 
     dni_extra: numeric
-        Extraterrestrial irradiance.
+        Extraterrestrial irradiance. The units of ``dni_extra``
+        determine the units of the output.
 
     Returns
     --------
     clearsky : DataFrame (if Series input) or OrderedDict of arrays
         DataFrame/OrderedDict contains the columns/keys
         ``'dhi', 'dni', 'ghi'``.
 
-        The units of ``dni_extra`` determine the units of the output.
-
     References
     ----------
     .. [1] P. Ineichen, "A broadband simplified version of the

diff --git a/pvlib/location.py b/pvlib/location.py
@@ -167,62 +167,74 @@ def get_solarposition(self, times, pressure=None, temperature=12
6D40
,
                                                **kwargs)
 
 
-    def get_clearsky(self, times, model='ineichen', **kwargs):
+    def get_clearsky(self, times, model='ineichen', solar_position=None,
+                     dni_extra=None, **kwargs):
         """
         Calculate the clear sky estimates of GHI, DNI, and/or DHI
         at this location.
 
         Parameters
         ----------
-        times : DatetimeIndex
-
-        model : str
+        times: DatetimeIndex
+        model: str
             The clear sky model to use. Must be one of
             'ineichen', 'haurwitz', 'simplified_solis'.
+        solar_position : None or DataFrame
+            DataFrame with with columns 'apparent_zenith', 'zenith',
+            'apparent_elevation'.
+        dni_extra: None or numeric
+            If None, will be calculated from times.
 
         kwargs passed to the relevant functions. Climatological values
-        are assumed in many cases. See code for details.
+        are assumed in many cases. See source code for details!
 
         Returns
         -------
         clearsky : DataFrame
             Column names are: ``ghi, dni, dhi``.
         """
+        if dni_extra is None:
+            dni_extra = irradiance.extraradiation(times.dayofyear)
 
-        if model == 'ineichen':
-            cs = clearsky.ineichen(times, latitude=self.latitude,
-                                   longitude=self.longitude,
-                                   altitude=self.altitude,
-                                   **kwargs)
-        elif model == 'haurwitz':
-            solpos = self.get_solarposition(times, **kwargs)
-            cs = clearsky.haurwitz(solpos['apparent_zenith'])
-        elif model == 'simplified_solis':
+        try:
+            pressure = kwargs.pop('pressure')
+        except KeyError:
+            pressure = atmosphere.alt2pres(self.altitude)
 
-            # these try/excepts define default values that are only
-            # evaluated if necessary. ineichen does some of this internally
-            try:
-                dni_extra = kwargs.pop('dni_extra')
-            except KeyError:
-                dni_extra = irradiance.extraradiation(times.dayofyear)
+        if solar_position is None:
+            solar_position = self.get_solarposition(times, pressure=pressure,
+                                                    **kwargs)
+
+        apparent_zenith = solar_position['apparent_zenith']
+        apparent_elevation = solar_position['apparent_elevation']
 
+        if model == 'ineichen':
             try:
-                pressure = kwargs.pop('pressure')
+                linke_turbidity = kwargs.pop('linke_turbidity')
             except KeyError:
-                pressure = atmosphere.alt2pres(self.altitude)
+                interp_turbidity = kwargs.pop('interp_turbidity', True)
+                linke_turbidity = clearsky.lookup_linke_turbidity(
+                    times, self.latitude, self.longitude,
+                    interp_turbidity=interp_turbidity)
 
             try:
-                apparent_elevation = kwargs.pop('apparent_elevation')
+                airmass_absolute = kwargs.pop('airmass_absolute')
             except KeyError:
-                solpos = self.get_solarposition(
-                    times, pressure=pressure, **kwargs)
-                apparent_elevation = solpos['apparent_elevation']
+                airmass_absolute = self.get_airmass(
+                    times, solar_position=solar_position)['airmass_absolute']
 
+            cs = clearsky.ineichen(apparent_zenith, airmass_absolute,
+                                   linke_turbidity, dni_extra=dni_extra,
+                                   altitude=self.altitude)
+        elif model == 'haurwitz':
+            cs = clearsky.haurwitz(apparent_zenith)
+        elif model == 'simplified_solis':
             cs = clearsky.simplified_solis(
                 apparent_elevation, pressure=pressure, dni_extra=dni_extra,
                 **kwargs)
         else:
-            raise ValueError('{} is not a valid clear sky model'
+            raise ValueError(('{} is not a valid clear sky model. Must be ' +
+                              'one of ineichen, simplified_solis, haurwitz')
                              .format(model))
 
         return cs

diff --git a/pvlib/test/test_location.py b/pvlib/test/test_location.py
@@ -107,9 +107,10 @@ def test_get_clearsky_simplified_solis_apparent_elevation():
     times = pd.DatetimeIndex(start='20160101T0600-0700',
                              end='20160101T1800-0700',
                              freq='3H')
-    apparent_elevation = pd.Series(80, index=times)
+    solar_position = {'apparent_elevation': pd.Series(80, index=times),
+                      'apparent_zenith': pd.Series(10, index=times)}
     clearsky = tus.get_clearsky(times, model='simplified_solis',
-                                apparent_elevation=apparent_elevation)
+                                solar_position=solar_position)
     expected = pd.DataFrame(data=np.
         array([[  131.3124497 ,  1001.14754036,  1108.14147919],
                [  131.3124497 ,  1001.14754036,  1108.14147919],