fix read_tmy3 with year coerced not monotonic, breaks soiling

docs/examples/plot_greensboro_kimber_soiling.py

@@ -42,7 +42,7 @@
 # calculate soiling with no wash dates
 THRESHOLD = 25.0
 soiling_no_wash = soiling_kimber(
-    greensboro_rain, cleaning_threshold=THRESHOLD, istmy=True)
+    greensboro_rain, cleaning_threshold=THRESHOLD, is_tmy=True)
 soiling_no_wash.name = 'soiling'
 # daily rain totals
 daily_rain = greensboro_rain.iloc[:-1].resample('D').sum()

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

@@ -39,6 +39,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/losses.py

@@ -9,13 +9,20 @@
 import pandas as pd
 from pvlib.tools import cosd
 
+_RAIN_ACC_PERIOD = pd.Timedelta('1h')
+_DEPO_VELOCITY = {'2_5': 0.004, '10': 0.0009}
+
 
 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=_RAIN_ACC_PERIOD,
+                is_tmy=False):
     """
     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
     ----------
@@ -46,6 +53,9 @@ def soiling_hsu(rainfall, cleaning_threshold, tilt, pm2_5, pm10,
         It is recommended that `rain_accum_period` be between 1 hour and
         24 hours.
 
+    is_tmy : bool, default False
+        Fix last timestep in TMY so that it is monotonically increasing.
+
     Returns
     -------
     soiling_ratio : Series
@@ -65,6 +75,22 @@ 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 = _DEPO_VELOCITY
+
+    # 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 is_tmy:
+        # get rainfall timezone, timestep as timedelta64
+        rain_tz = rainfall.index.tz
+        rain_index = rainfall.index.values
+        timestep_interval = (rain_index[1] - rain_index[0])
+        rain_vals = rainfall.values
+        rain_name = rainfall.name
+        rainfall = _fix_tmy_monotonicity(
+            rain_index, rain_vals, timestep_interval, rain_tz, rain_name)
+
     # 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
@@ -89,14 +115,26 @@ def soiling_hsu(rainfall, cleaning_threshold, tilt, pm2_5, pm10,
     return soiling_ratio
 
 
+def _fix_tmy_monotonicity(rain_index, rain_vals, timestep_interval, rain_tz,
+                          rain_name):
+    # fix TMY to be monotonically increasing by rolling index by 1 interval
+    # and then adding 1 interval, while the values stay the same
+    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
+    return pd.Series(rain_vals, index=rain_index, name=rain_name)
+
+
 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, is_tmy=False):
     """
-    Calculate soiling ratio with rainfall data and a daily soiling rate using
-    the Kimber soiling model [1]_.
+    Calculates fraction of energy lossed 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,7 +165,7 @@ 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
+    is_tmy : bool, default False
         Fix last timestep in TMY so that it is monotonically increasing.
 
     Returns
@@ -166,22 +204,19 @@ 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
+    # get rainfall timezone, timestep as timedelta64, and timestep as day-frac
     rain_tz = rainfall.index.tz
-    rain_index = rainfall.index.values
-    timestep_interval = (rain_index[1] - rain_index[0])
+    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)
+    if is_tmy:
+        rain_vals = rainfall.values
+        rain_name = rainfall.name
+        rainfall = _fix_tmy_monotonicity(
+            rain_index_vals, rain_vals, timestep_interval, rain_tz, rain_name)
 
     # accumulate rainfall
     accumulated_rainfall = rainfall.rolling(
@@ -191,6 +226,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
@@ -206,7 +242,6 @@ def soiling_kimber(rainfall, cleaning_threshold=6, soiling_loss_rate=0.0015,
     # manual wash dates
     if manual_wash_dates is not None:
         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/test_losses.py

@@ -10,6 +10,7 @@
 from conftest import (
     requires_scipy, needs_pandas_0_22, DATA_DIR)
 import pytest
+import pytz
 
 
 @pytest.fixture
@@ -121,7 +122,7 @@ def test_kimber_soiling_nowash(greensboro_rain,
     # Greensboro typical expected annual rainfall is 8345mm
     assert greensboro_rain.sum() == 8345
     # calculate soiling with no wash dates
-    soiling_nowash = soiling_kimber(greensboro_rain, istmy=True)
+    soiling_nowash = soiling_kimber(greensboro_rain, is_tmy=True)
     # test no washes
     assert np.allclose(
         soiling_nowash.values,
@@ -143,7 +144,7 @@ def test_kimber_soiling_manwash(greensboro_rain,
     manwash = [datetime.date(1990, 2, 15), ]
     # calculate soiling with manual wash
     soiling_manwash = soiling_kimber(
-        greensboro_rain, manual_wash_dates=manwash, istmy=True)
+        greensboro_rain, manual_wash_dates=manwash, is_tmy=True)
     # test manual wash
     assert np.allclose(
         soiling_manwash.values,
@@ -168,7 +169,7 @@ def test_kimber_soiling_norain(greensboro_rain,
     # a year with no rain
     norain = pd.Series(0, index=greensboro_rain.index)
     # calculate soiling with no rain
-    soiling_norain = soiling_kimber(norain, istmy=True)
+    soiling_norain = soiling_kimber(norain, is_tmy=True)
     # test no rain, soiling reaches maximum
     assert np.allclose(soiling_norain.values, expected_kimber_soiling_norain)
 
@@ -191,7 +192,41 @@ def test_kimber_soiling_initial_soil(greensboro_rain,
     # a year with no rain
     norain = pd.Series(0, index=greensboro_rain.index)
     # calculate soiling with no rain
-    soiling_norain = soiling_kimber(norain, initial_soiling=0.1, istmy=True)
+    soiling_norain = soiling_kimber(norain, initial_soiling=0.1, is_tmy=True)
     # test no rain, soiling reaches maximum
     assert np.allclose(
         soiling_norain.values, expected_kimber_soiling_initial_soil)
+
+
+@pytest.fixture
+def expected_greensboro_hsu_soil():
+    return np.array([
+        0.99927224, 0.99869067, 0.99815393, 0.99764437, 0.99715412,
+        0.99667873, 0.99621536, 0.99576203, 0.99531731, 0.99488010,
+        0.99444954, 0.99402494, 0.99360572, 0.99319142, 1.00000000,
+        1.00000000, 0.99927224, 0.99869067, 0.99815393, 0.99764437,
+        0.99715412, 1.00000000, 0.99927224, 0.99869067])
+
+
+@requires_scipy
+def test_gh889_soiing_hsu_tmy_not_monotonic(expected_greensboro_hsu_soil):
+    """doesn't raise value error"""
+    greensboro = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
+    greensboro_rain = greensboro[0].Lprecipdepth
+    soiling_ratio = soiling_hsu(
+        greensboro_rain, cleaning_threshold=10.0, tilt=0.0, pm2_5=1.0,
+        pm10=2.0, is_tmy=True)
+    # check first day of soiling ratio, actually (1 - transmission loss)
+    # greensboro rains hours 3pm, 4pm, and 10pm, so expect soiling ratio of one
+    assert np.allclose(expected_greensboro_hsu_soil, soiling_ratio.values[:24])
+    # greensboro timezone is UTC-5 or Eastern time
+    gmt_5 = pytz.timezone('Etc/GMT+5')
+    # check last day, should be 1991 now
+    lastday = datetime.datetime(1991, 1, 1, 0, 0, 0)
+    assert gmt_5.localize(lastday) == soiling_ratio.index[-1]
+    # check last hour is still 1990
+    lasthour = datetime.datetime(1990, 12, 31, 23, 0, 0)
+    assert gmt_5.localize(lasthour) == soiling_ratio.index[-2]
+    # check first hour is still 1990
+    firsthour = datetime.datetime(1990, 1, 1, 1, 0, 0)
+    assert gmt_5.localize(firsthour) == soiling_ratio.index[0]