pvlib · cwhanse · Aug 27, 2018 · Aug 31, 2017 · Aug 31, 2017 · Aug 31, 2017
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
@@ -166,6 +166,7 @@ DNI estimation models
    irradiance.dirindex
    irradiance.erbs
    irradiance.liujordan
+   irradiance.gti_dirint
 
 
 PV Modeling

diff --git a/docs/sphinx/source/whatsnew/v0.5.2.rst b/docs/sphinx/source/whatsnew/v0.5.2.rst
@@ -6,10 +6,13 @@ v0.5.2 (__, 2017)
 API Changes
 ~~~~~~~~~~~
 * removed unused 'pressure' arg from irradiance.liujordan function (:issue:`386`)
+* Add kt and kt_prime keyword arguments to dirint function. (:issue:`396`)
+* Add kt and min_cos_zenith keyword arguments to disc function. (:issue:`396`)
+* Change the disc dni filter thresholds. (:issue:`396`)
 
 Enhancements
 ~~~~~~~~~~~~
-*
+* Add gti_dirint function. (:issue:`396`)
 
 Bug fixes
 ~~~~~~~~~

diff --git a/docs/tutorials/gti_dirint.ipynb b/docs/tutorials/gti_dirint.ipynb
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
@@ -1070,7 +1070,8 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
         return sky_diffuse
 
 
-def disc(ghi, zenith, datetime_or_doy, pressure=101325):
+def disc(ghi, zenith, datetime_or_doy, pressure=101325, min_cos_zenith=0.065,
+         kt=None):
     """
     Estimate Direct Normal Irradiance from Global Horizontal Irradiance
     using the DISC model.
@@ -1127,13 +1128,17 @@ def disc(ghi, zenith, datetime_or_doy, pressure=101325):
 
     # this is the I0 calculation from the reference
     I0 = extraradiation(datetime_or_doy, 1370, 'spencer')
-    I0h = I0 * np.cos(np.radians(zenith))
+    cos_zenith = np.maximum(min_cos_zenith, np.cos(np.radians(zenith)))
+    I0h = I0 * cos_zenith
 
     am = atmosphere.relativeairmass(zenith, model='kasten1966')
     am = atmosphere.absoluteairmass(am, pressure)
 
-    kt = ghi / I0h
+    if kt is None:
+        kt = ghi / I0h
+
     kt = np.maximum(kt, 0)
+    kt = np.minimum(kt, 0.82)
     # powers of kt will be used repeatedly, so compute only once
     kt2 = kt * kt  # about the same as kt ** 2
     kt3 = kt2 * kt  # 5-10x faster than kt ** 3
@@ -1156,7 +1161,8 @@ def disc(ghi, zenith, datetime_or_doy, pressure=101325):
 
     dni = Kn * I0
 
-    dni = np.where((zenith > 87) | (ghi < 0) | (dni < 0), 0, dni)
+    dni = np.where((zenith > 90) | (ghi < 0) | (dni < 0) | (dni > I0*0.82),
+                   0, dni)
 
     output = OrderedDict()
     output['dni'] = dni
@@ -1170,7 +1176,7 @@ def disc(ghi, zenith, datetime_or_doy, pressure=101325):
 
 
 def dirint(ghi, zenith, times, pressure=101325., use_delta_kt_prime=True,
-           temp_dew=None):
+           temp_dew=None, kt=None, kt_prime=None, return_kt_prime=False):
     """
     Determine DNI from GHI using the DIRINT modification of the DISC
     model.
@@ -1237,19 +1243,15 @@ def dirint(ghi, zenith, times, pressure=101325., use_delta_kt_prime=True,
     SERI/TR-215-3087, Golden, CO: Solar Energy Research Institute, 1987.
     """
 
-    disc_out = disc(ghi, zenith, times, pressure=pressure)
+    disc_out = disc(ghi, zenith, times, pressure=pressure, kt=kt)
     dni = disc_out['dni']
     kt = disc_out['kt']
     am = disc_out['airmass']
 
-    kt_prime = kt / (1.031 * np.exp(-1.4 / (0.9 + 9.4 / am)) + 0.1)
-    kt_prime = np.minimum(kt_prime, 0.82)  # From SRRL code
+    if kt_prime is None:
+        kt_prime = kt / (1.031 * np.exp(-1.4 / (0.9 + 9.4 / am)) + 0.1)
+        kt_prime = np.minimum(kt_prime, 0.82)  # From SRRL code
 
-    # 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.
-    # current implementation requires that kt_prime is a Series
     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(),
@@ -1317,7 +1319,10 @@ def dirint(ghi, zenith, times, pressure=101325., use_delta_kt_prime=True,
 
     dni *= dirint_coeffs
 
-    return dni
+    if return_kt_prime:
+        return dni, kt_prime
+    else:
+        return dni
 
 
 def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
@@ -1400,6 +1405,268 @@ def dirindex(ghi, ghi_clearsky, dni_clearsky, zenith, times, pressure=101325.,
     return dni_dirindex
 
 
+def gti_dirint(poa_global, aoi, solar_zenith, solar_azimuth,
+               surface_tilt, surface_azimuth,
+               times, pressure=101325.,
+               use_delta_kt_prime=True, temp_dew=None, albedo=.25,
+               model='perez', model_perez='allsitescomposite1990',
+               calculate_gt_90=True, max_iterations=30):
+    """
+    Determine GHI, DNI, DHI from POA global using the GTI DIRINT model.
+
+    Parameters
+    ----------
+    poa_global : array-like
+        Plane of array global irradiance in W/m^2.
+
+    aoi : array-like
+        Angle of incidence of solar rays with respect to the module
+        surface normal.
+
+    zenith : array-like
+        True (not refraction-corrected) zenith angles in decimal
+        degrees.
+
+    surface_tilt : numeric
+        Surface tilt angles in decimal degrees. 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).
+
+    times : DatetimeIndex
+
+    pressure : numeric, default 101325.0
+        The site pressure in Pascal. Pressure may be measured or an
+        average pressure may be calculated from site altitude.
+
+    use_delta_kt_prime : bool, default True
+        Indicates if the user would like to utilize the time-series
+        nature of the GHI measurements. A value of ``False`` will not
+        use the time-series improvements, any other numeric value will
+        use time-series improvements. It is recommended that time-series
+        data only be used if the time between measured data points is
+        less 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). If True, input data must be Series.
+
+    temp_dew : numeric, default None
+        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.
+
+    calculate_gt_90 : bool, default True
+        Controls if the algorithm evaluates inputs with AOI >= 90 degrees.
+        If False, returns nan for AOI >= 90 degrees. Significant speed ups
+        can be achieved by setting this parameter to False.
+
+    Returns
+    -------
+    data : OrderedDict or DataFrame
+        Contains the following keys/columns:
+
+            * ``ghi``: the modeled global horizontal irradiance in W/m^2.
+            * ``dni``: the modeled direct normal irradiance in W/m^2.
+            * ``dhi``: the modeled diffuse horizontal irradiance in
+              W/m^2.
+
+    References
+    ----------
+    .. [1] B. Marion, A model for deriving the direct normal and
+           diffuse horizontal irradiance from the global tilted
+           irradiance, Solar Energy 122, 1037-1046.
+           http://dx.doi.org/10.1016/j.solener.2015.10.024
+    """
+
+    aoi_lt_90 = aoi < 90
+
+    # for AOI less than 90 degrees
+
+    ghi, dni, dhi, kt_prime = _gti_dirint_lt_90(
+        poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth,
+        surface_tilt, surface_azimuth, times, pressure=pressure,
+        use_delta_kt_prime=use_delta_kt_prime, temp_dew=temp_dew,
+        albedo=albedo, model=model, model_perez=model_perez,
+        max_iterations=max_iterations)
+
+    # for AOI greater than or equal to 90 degrees
+
+    if calculate_gt_90:
+        ghi_gte_90, dni_gte_90, dhi_gte_90 = _gti_dirint_gte_90(
+            poa_global, aoi, solar_zenith, solar_azimuth,
+            surface_tilt, times, kt_prime,
+            pressure=pressure, temp_dew=temp_dew, albedo=albedo)
+    else:
+        ghi_gte_90, dni_gte_90, dhi_gte_90 = np.nan, np.nan, np.nan
+
+    # put the AOI < 90 and AOI >= 90 conditions together
+
+    output = OrderedDict()
+    output['ghi'] = ghi.where(aoi_lt_90, ghi_gte_90)
+    output['dni'] = dni.where(aoi_lt_90, dni_gte_90)
+    output['dhi'] = dhi.where(aoi_lt_90, dhi_gte_90)
+
+    output = pd.DataFrame(output, index=times)
+
+    return output
+
+
+def _gti_dirint_lt_90(poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth,
+                      surface_tilt, surface_azimuth,
+                      times, pressure=101325.,
+                      use_delta_kt_prime=True, temp_dew=None, albedo=.25,
+                      model='perez', model_perez='allsitescomposite1990',
+                      max_iterations=30):
+
+    I0 = extraradiation(times, 1370, 'spencer')
+    # cos_zenith = np.maximum(0.065, tools.cosd(solar_zenith))
+    cos_zenith = tools.cosd(solar_zenith)
+    I0h = I0 * cos_zenith
+
+    # these coeffs and diff variables and the loop below
+    # implement figure 1 of Marion 2015
+
+    # make coeffs that is at least 30 elements long
+    # for loop below will limit iterations if necessary
+    coeffs = np.empty(max(30, max_iterations))
+    coeffs[0:3] = 1
+    coeffs[3:10] = 0.5
+    coeffs[10:20] = 0.25
+    coeffs[20:] = 0.125
+
+    # initialize diff
+    diff = pd.Series(9999, index=times)
+    best_diff = diff
+
+    # initialize poa_global_i
+    poa_global_i = poa_global
+
+    for coeff in coeffs[:max_iterations]:
+        # use slightly > 1 for float errors
+        best_diff_lte_1 = best_diff <= 1.000001
+
+        # only test for aoi less than 90 deg
+        best_diff_lte_1_lt_90 = best_diff_lte_1[aoi_lt_90]
+        if best_diff_lte_1_lt_90.all():
+            break
+
+        # calculate kt and DNI using GTI and pvlib implementations
+        # of Marion eqn 2 and DIRINT
+        disc_out = disc(poa_global_i, aoi, times, pressure=pressure)
+        dni, kt_prime = dirint(
+            poa_global_i, aoi, times, pressure=pressure,
+            use_delta_kt_prime=use_delta_kt_prime,
+            temp_dew=temp_dew, return_kt_prime=True)
+
+        # calculate DHI using Marion eqn 3 (identify 1st term as GHI)
+        #ghi = np.minimum(disc_out['kt'], 0.82) * I0h
+        ghi = disc_out['kt'] * I0h
+        dhi = ghi - dni * cos_zenith
+        bad_values = (dhi < 0) | (dni < 0) | (ghi < 0)
+        dni[bad_values] = np.nan
+        ghi[bad_values] = np.nan
+        dhi[bad_values] = np.nan
+
+        # use DNI and DHI to model GTI
+        all_irrad = total_irrad(
+            surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
+            dni, ghi, dhi, dni_extra=I0, airmass=disc_out['airmass'],
+            albedo=albedo, model=model, model_perez=model_perez)
+
+        gti_model = all_irrad['poa_global']
+
+        # calculate new diff
+        diff = gti_model - poa_global
+
+        # determine if the new diff is smaller in magnitude
+        # than the old diff
+        diff_abs = diff.abs()
+        smallest_diff = diff_abs < best_diff
+
+        # save the best differences
+        best_diff = diff_abs.where(smallest_diff, best_diff)
+
+        # save DNI, DHI, DHI if they provide the best consistency
+        # otherwise use the older values.
+        # try/except accounts for first iteration
+        try:
+            best_ghi = ghi.where(smallest_diff, best_ghi)
+            best_dni = dni.where(smallest_diff, best_dni)
+            best_dhi = dhi.where(smallest_diff, best_dhi)
+            best_kt_prime = kt_prime.where(smallest_diff, best_kt_prime)
+        except UnboundLocalError:
+            best_ghi = ghi
+            best_dni = dni
+            best_dhi = dhi
+            best_kt_prime = kt_prime
+
+        # calculate adjusted inputs for next iteration
+        poa_global_i = np.maximum(1.0, poa_global_i - coeff * diff)
+    else:
+        import warnings
+        failed_points = best_diff[aoi_lt_90][best_diff_lte_1_lt_90 == False]
+        warnings.warn(
+            '%s points failed to converge after %s iterations.'
+            'best_diff:\n%s' %
+            (len(failed_points), max_iterations, failed_points),
+            RuntimeWarning)
+
+    return best_ghi, best_dni, best_dhi, best_kt_prime
+
+
+def _gti_dirint_gte_90(poa_global, aoi, solar_zenith, solar_azimuth,
+                       surface_tilt, times, kt_prime,
+                       pressure=101325., temp_dew=None, albedo=.25,):
+
+    # set the kt_prime for sunrise to AOI=90 to be equal to
+    # the kt_prime for 65 < AOI < 80 during the morning.
+    # similar for the afternoon. repeat for every day.
+    aoi_gte_90 = aoi >= 90
+    aoi_lt_90 = aoi < 90
+    aoi_65_80 = (aoi > 65) & (aoi < 80)
+    zenith_lt_90 = solar_zenith < 90
+    morning = solar_azimuth < 180
+    afternoon = solar_azimuth > 180
+    aoi_65_80_morning = aoi_65_80 & morning
+    aoi_65_80_afternoon = aoi_65_80 & afternoon
+    zenith_lt_90_aoi_gte_90_morning = zenith_lt_90 & aoi_gte_90 & morning
+    zenith_lt_90_aoi_gte_90_afternoon = zenith_lt_90 & aoi_gte_90 & afternoon
+
+    kt_prime_90s = []
+    for date, data in kt_prime.groupby(times.date):
+        kt_prime_am_avg = data[aoi_65_80_morning].mean()
+        kt_prime_pm_avg = data[aoi_65_80_afternoon].mean()
+
+        kt_prime_90 = pd.Series(np.nan, index=data.index)
+        kt_prime_90[zenith_lt_90_aoi_gte_90_morning] = kt_prime_am_avg
+        kt_prime_90[zenith_lt_90_aoi_gte_90_afternoon] = kt_prime_pm_avg
+        kt_prime_90s.append(kt_prime_90)
+    kt_prime_90s = pd.concat(kt_prime_90s)
+
+    am = atmosphere.relativeairmass(solar_zenith, model='kasten1966')
+    am = atmosphere.absoluteairmass(am, pressure)
+    kt = kt_prime_90s * (1.031 * np.exp(-1.4 / (0.9 + 9.4 / am)) + 0.1)
+
+    dni_gte_90 = dirint(
+        poa_global, solar_zenith, times, pressure=pressure,
+        use_delta_kt_prime=False,
+        temp_dew=temp_dew, kt=kt, kt_prime=kt_prime_90s)
+
+    cos_zenith = tools.cosd(solar_zenith)
+
+    dni_gte_90_proj = dni_gte_90 * cos_zenith
+
+    cos_surface_tilt = tools.cosd(surface_tilt)
+    # isotropic sky plus ground diffuse
+    dhi_gte_90 = (
+        (2 * poa_global - dni_gte_90_proj * albedo * (1 - cos_surface_tilt)) /
+        (1 + cos_surface_tilt + albedo * (1 - cos_surface_tilt)))
+
+    ghi_gte_90 = dni_gte_90_proj + dhi_gte_90
+
+    return ghi_gte_90, dni_gte_90, dhi_gte_90
+
+
 def erbs(ghi, zenith, doy):
     r"""
     Estimate DNI and DHI from GHI using the Erbs model.