fix read_tmy3 with year coerced not monotonic, breaks soiling

docs/examples/plot_greensboro_kimber_soiling.py

@@ -2,20 +2,20 @@
 Kimber Soiling Model
 ====================
 
-Examples of soiling using the Kimber model [1]_.
-
-References
-----------
-.. [1] "The Effect of Soiling on Large Grid-Connected Photovoltaic Systems
-   in California and the Southwest Region of the United States," Adrianne
-   Kimber, et al., IEEE 4th World Conference on Photovoltaic Energy
-   Conference, 2006, :doi:`10.1109/WCPEC.2006.279690`
+Examples of soiling using the Kimber model.
 """
 
 # %%
-# This example shows basic usage of pvlib's Kimber Soiling model with
+# This example shows basic usage of pvlib's Kimber Soiling model [1]_ with
 # :py:meth:`pvlib.losses.soiling_kimber`.
 #
+# References
+# ----------
+# .. [1] "The Effect of Soiling on Large Grid-Connected Photovoltaic Systems
+#    in California and the Southwest Region of the United States," Adrianne
+#    Kimber, et al., IEEE 4th World Conference on Photovoltaic Energy
+#    Conference, 2006, :doi:`10.1109/WCPEC.2006.279690`
+#
 # The Kimber Soiling model assumes that soiling builds up at a constant rate
 # until cleaned either manually or by rain. The rain must reach a threshold to
 # clean the panels. When rains exceeds the threshold, it's assumed the earth is
@@ -30,19 +30,24 @@
 # step.
 
 from datetime import datetime
+import pathlib
 from matplotlib import pyplot as plt
-from pvlib.iotools import read_tmy3
+from pvlib.iotools import read_tmy3, fix_tmy3_coerce_year_monotonicity
 from pvlib.losses import soiling_kimber
-from pvlib.tests.conftest import DATA_DIR
+import pvlib
+
+# get full path to the data directory
+DATA_DIR = pathlib.Path(pvlib.__file__).parent / 'data'
 
 # get TMY3 data with rain
-greensboro = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
-# NOTE: can't use Sand Point, AK b/c Lprecipdepth is -9900, ie: missing
-greensboro_rain = greensboro[0].Lprecipdepth
-# calculate soiling with no wash dates
+greensboro, _ = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
+# fix TMY3 index to be monotonically increasing
+greensboro = fix_tmy3_coerce_year_monotonicity(greensboro)
+# get the rain data
+greensboro_rain = greensboro.Lprecipdepth
+# calculate soiling with no wash dates and cleaning threshold of 25-mm of rain
 THRESHOLD = 25.0
-soiling_no_wash = soiling_kimber(
-    greensboro_rain, cleaning_threshold=THRESHOLD, istmy=True)
+soiling_no_wash = soiling_kimber(greensboro_rain, cleaning_threshold=THRESHOLD)
 soiling_no_wash.name = 'soiling'
 # daily rain totals
 daily_rain = greensboro_rain.iloc[:-1].resample('D').sum()

docs/sphinx/source/api.rst

@@ -369,6 +369,7 @@ relevant to solar energy modeling.
 
    iotools.read_tmy2
    iotools.read_tmy3
+   iotools.fix_tmy3_coerce_year_monotonicity
    iotools.read_epw
    iotools.parse_epw
    iotools.read_srml

docs/sphinx/source/whatsnew/v0.7.2.rst

@@ -14,6 +14,8 @@ Enhancements
   the original ``Date (MM/DD/YYYY)`` and ``Time (HH:MM)`` columns that the
   indices were parsed from (:pull:`866`)
 * Add Kimber soiling model :py:func:`pvlib.losses.soiling_kimber` (:pull:`860`)
+* Add :py:func:`~pvlib.iotools.fix_tmy3_coerce_year_monotonicity` to fix TMY3
+  data coerced to a single year to be monotonically increasing (:pull:`910`)
 
 Bug fixes
 ~~~~~~~~~
@@ -39,6 +41,8 @@ Bug fixes
   :py:func:`~pvlib.iam.martin_ruiz_diffuse`,
   :py:func:`~pvlib.losses.soiling_hsu`,
   and various test functions.
+* Fix :py:func:`~pvlib.losses.soiling_hsu` TMY3 data not monotonic exception
+  (:pull:`910`)
 
 Documentation
 ~~~~~~~~~~~~~

pvlib/iotools/__init__.py

@@ -1,5 +1,5 @@
-from pvlib.iotools.tmy import read_tmy2  # noqa: F401
-from pvlib.iotools.tmy import read_tmy3  # noqa: F401
+from pvlib.iotools.tmy import read_tmy2, read_tmy3  # noqa: F401
+from pvlib.iotools.tmy import fix_tmy3_coerce_year_monotonicity  # noqa: F401
 from pvlib.iotools.epw import read_epw, parse_epw  # noqa: F401
 from pvlib.iotools.srml import read_srml  # noqa: F401
 from pvlib.iotools.srml import read_srml_month_from_solardat  # noqa: F401

pvlib/iotools/tmy.py

@@ -7,6 +7,7 @@
 import re
 from urllib.request import urlopen, Request
 import pandas as pd
+import numpy as np
 
 
 def read_tmy3(filename=None, coerce_year=None, recolumn=True):
@@ -526,3 +527,44 @@ def _read_tmy2(string, columns, hdr_columns, fname):
         columns=columns.split(',')).tz_localize(int(meta['TZ'] * 3600))
 
     return data, meta
+
+
+def fix_tmy3_coerce_year_monotonicity(tmy3):
+    """
+    Fix TMY3 coerced to a single year to be monotonically increasing by
+    changing the last record to be in January 1st of the next year.
+
+    Parameters
+    ----------
+    tmy3 : pandas.DataFrame
+        TMY3 data frame from :func:`pvlib.iotools.read_tmy3` with year coearced
+
+    Returns
+    -------
+    pandas.DataFrame
+        Copy of the TMY3 data frame with monotonically increasing index
+    """
+    # NOTE: pandas index is immutable, therefore it's not possible to change a
+    # single item in the index, so the entire index must be replaced
+
+    # get tmy3 timezone, index values as np.datetime64[ns], and index frequency
+    # as np.timedelta64[ns]
+    # NOTE: numpy converts index values to UTC
+    index_tz = tmy3.index.tz
+    index_values = tmy3.index.values
+    timestep_interval = index_values[1] - index_values[0]
+
+    # fix index to be monotonically increasing by rolling indices 1 interval,
+    # then adding 1 interval to all indices
+    index_values = np.roll(index_values, 1) + timestep_interval
+
+    # create new datetime index and convert it to the original timezone
+    new_index = pd.DatetimeIndex(index_values, tz='UTC').tz_convert(index_tz)
+
+    # copy the original TMY3 so it doesn't change, then replace the index of
+    # the copy with the new fixed monotonically increasing index
+    new_tmy3 = tmy3.copy()
+    new_tmy3.index = new_index
+
+    # fixed TMY3 with monotonically increasing index
+    return new_tmy3

pvlib/losses.py

@@ -11,11 +11,14 @@
 
 
 def soiling_hsu(rainfall, cleaning_threshold, tilt, pm2_5, pm10,
-                depo_veloc={'2_5': 0.004, '10': 0.0009},
-                rain_accum_period=pd.Timedelta('1h')):
+                depo_veloc=None, rain_accum_period=pd.Timedelta('1h')):
     """
     Calculates soiling ratio given particulate and rain data using the model
-    from Humboldt State University [1]_.
+    from Humboldt State University (HSU).
+
+    The HSU soiling model [1]_ returns the soiling ratio, a value between zero
+    and one which is equivalent to (1 - transmission loss). Therefore a soiling
+    ratio of 1.0 is equivalent to zero transmission loss.
 
     Parameters
     ----------
@@ -65,6 +68,10 @@ def soiling_hsu(rainfall, cleaning_threshold, tilt, pm2_5, pm10,
     except ImportError:
         raise ImportError("The soiling_hsu function requires scipy.")
 
+    # never use mutable input arguments
+    if depo_veloc is None:
+        depo_veloc = {'2_5': 0.004, '10': 0.0009}
+
     # accumulate rainfall into periods for comparison with threshold
     accum_rain = rainfall.rolling(rain_accum_period, closed='right').sum()
     # cleaning is True for intervals with rainfall greater than threshold
@@ -91,12 +98,12 @@ def soiling_hsu(rainfall, cleaning_threshold, tilt, pm2_5, pm10,
 
 def soiling_kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
                    grace_period=14, max_soiling=0.3, manual_wash_dates=None,
-                   initial_soiling=0, rain_accum_period=24, istmy=False):
+                   initial_soiling=0, rain_accum_period=24):
     """
-    Calculate soiling ratio with rainfall data and a daily soiling rate using
-    the Kimber soiling model [1]_.
+    Calculates fraction of energy lost due to soiling given rainfall data and
+    daily loss rate using the Kimber model.
 
-    Kimber soiling model assumes soiling builds up at a daily rate unless
+    Kimber soiling model [1]_ assumes soiling builds up at a daily rate unless
     the daily rainfall is greater than a threshold. The model also assumes that
     if daily rainfall has exceeded the threshold within a grace period, then
     the ground is too damp to cause soiling build-up. The model also assumes
@@ -127,8 +134,6 @@ def soiling_kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
     rain_accum_period : int, default 24
         Period for accumulating rainfall to check against `cleaning_threshold`.
         The Kimber model defines this period as one day. [hours]
-    istmy : bool, default False
-        Fix last timestep in TMY so that it is monotonically increasing.
 
     Returns
     -------
@@ -166,23 +171,12 @@ def soiling_kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
     # convert grace_period to timedelta
     grace_period = datetime.timedelta(days=grace_period)
 
-    # get rainfall timezone, timestep as timedelta64, and timestep in int days
-    rain_tz = rainfall.index.tz
-    rain_index = rainfall.index.values
-    timestep_interval = (rain_index[1] - rain_index[0])
+    # get indices as numpy datetime64, calculate timestep as numpy timedelta64,
+    # and convert timestep to fraction of days
+    rain_index_vals = rainfall.index.values
+    timestep_interval = (rain_index_vals[1] - rain_index_vals[0])
     day_fraction = timestep_interval / np.timedelta64(24, 'h')
 
-    # if TMY fix to be monotonically increasing by rolling index by 1 interval
-    # and then adding 1 interval, while the values stay the same
-    if istmy:
-        rain_index = np.roll(rain_index, 1) + timestep_interval
-        # NOTE: numpy datetim64[ns] has no timezone
-        # convert to datetimeindex at UTC and convert to original timezone
-        rain_index = pd.DatetimeIndex(rain_index, tz='UTC').tz_convert(rain_tz)
-        # fixed rainfall timeseries with monotonically increasing index
-        rainfall = pd.Series(
-            rainfall.values, index=rain_index, name=rainfall.name)
-
     # accumulate rainfall
     accumulated_rainfall = rainfall.rolling(
         rain_accum_period, closed='right').sum()
@@ -191,6 +185,7 @@ def soiling_kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
     soiling = np.ones_like(rainfall.values) * soiling_loss_rate * day_fraction
     soiling[0] = initial_soiling
     soiling = np.cumsum(soiling)
+    soiling = pd.Series(soiling, index=rainfall.index, name='soiling')
 
     # rainfall events that clean the panels
     rain_events = accumulated_rainfall > cleaning_threshold
@@ -205,8 +200,9 @@ def soiling_kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
 
     # manual wash dates
     if manual_wash_dates is not None:
+        rain_tz = rainfall.index.tz
+        # convert manual wash dates to datetime index in the timezone of rain
         manual_wash_dates = pd.DatetimeIndex(manual_wash_dates, tz=rain_tz)
-        soiling = pd.Series(soiling, index=rain_index, name='soiling')
         cleaning[manual_wash_dates] = soiling[manual_wash_dates]
 
     # remove soiling by foward filling cleaning where NaN

pvlib/tests/iotools/test_tmy.py

@@ -1,8 +1,11 @@
+import datetime
 from pandas.util.testing import network
 import numpy as np
 import pandas as pd
 import pytest
+import pytz
 from pvlib.iotools import tmy
+from pvlib.iotools import read_tmy3, fix_tmy3_coerce_year_monotonicity
 from conftest import DATA_DIR
 
 # test the API works
@@ -77,3 +80,34 @@ def test_gh865_read_tmy3_feb_leapyear_hr24():
     # hour so check that the 1st hour is 1AM and the last hour is midnite
     assert data.index[0].hour == 1
     assert data.index[-1].hour == 0
+
+
+def test_fix_tmy_coerce_year_monotonicity():
+    # greensboro timezone is UTC-5 or Eastern time
+    gmt_5 = pytz.timezone('Etc/GMT+5')
+
+    # tmy3 coerced to year is not monotonically increasing
+    greensboro, _ = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
+
+    # check first hour was coerced to 1990
+    firsthour = gmt_5.localize(datetime.datetime(1990, 1, 1, 1, 0, 0))
+    assert firsthour == greensboro.index[0]
+
+    # check last hour was coerced to 1990
+    lasthour = gmt_5.localize(datetime.datetime(1990, 12, 31, 23, 0, 0))
+    assert lasthour == greensboro.index[-2]
+
+    # check last day, was coerced to 1990
+    lastday1990 = gmt_5.localize(datetime.datetime(1990, 1, 1, 0, 0, 0))
+    assert lastday1990 == greensboro.index[-1]
+
+    # fix the index to be monotonically increasing
+    greensboro = fix_tmy3_coerce_year_monotonicity(greensboro)
+
+    # check first and last hours are still 1990
+    assert firsthour == greensboro.index[0]
+    assert lasthour == greensboro.index[-2]
+
+    # check last day, should be 1991 now
+    lastday1991 = lastday1990.replace(year=1991)
+    assert lastday1991 == greensboro.index[-1]