fix read_tmy3 with year coerced not monotonic, breaks soiling (#910)

* fix soiling_hsu for tmy not monotonic

* if tmy is coerced to year, then the last hour is in the same year,
 before the other timestamps, so timestamps are no monotonic
* closes #889
* add underscore to is_tmy
* add is_tmy to soiling_hsu
* add private _fix_tmy_monotonicity
* also fix reference footnote in first line of docstring,
 since used in api.rst, causes doc warning
* don't use mutable input args
* don't create pd.timedelta object every single time soiling_hsu() called
* always create soiling pd.Series for consistency and robust to future edits
* fix rain_index is not a tz-aware pd.DatetimeIndex series,
 so rename as rain_index_vals and use rainfall.index in Soiling
* add new test for gh889, check soiling ratio for 24 hours,
 and check timestamp of first and last hours, and last day

* update what's new with fix for TMY3 in soiling_hsu

* also change 1st line of kimber soiling docstring, to say calculates
 energy loss, not soiling ratio, since kimber is 1-hsu,
  they're opposite!

* move fix_tmy_monotonicity to tmy in iotools

* and remove it from losses, remove is_tmy args, and remove calls to
 fix tmy monotonicity
* test to make sure if coerced year, last record is Jan 1 of next year
* add fix_tmy3_coerce_year_monotonicity to iotools api
* revert changes to soiling_hsu, don't use private constants for
 rain accumulation period or deposition velocity
* update kimber soiling example in gallery
  - to use fix_tmy3_coerce_year_monotonicity
  - and not use conftest, instead build path
  - also better comments
* remove test for gh889 in test_losses
* apply fix to greensboro rain data fixture and remove is_tmy args from kimber tests
* add needs_pandas_0_22 to soiling_hsu tests to show that it is required,
 since requies_scipy masks the requirement

* fix path to data in kimber example

* add test for HSU soiling defaults for coverage
* fix iotools api stickler warning
* remove unused pytz import in test_losses

* update what's new and api.rst with fix_tmy3_coerce_year_monotonicity

* parameters spelling typo in fix_tmy3_monotonicity docstring

* also move reference in kimber soiling example to comment in first code cell

* add monotonic_index kwarg to read_tmy3

* change name to tmy3_monotonic_index
* set monotonic_index default to false but warn that
 this will change to True in v0.8
* improve tmy3_monotonic_index docstring with reviewer comments, and add
 note to call only after read_tmy3 with coerce year,
  to use monotonic_index=True to combine calls, and that there's no validation
  that input TMY3 has been coerced or if it has already been fixed
* also if only used for TMY3 no need to calculate timestep interval,
 it's always 1-hour
* that's all folks

* update what's new with monotonic_index kwarg

* Update pvlib/tests/iotools/test_tmy.py

* change test name to `test_tmy3_monotonic_index`

Co-Authored-By: Cliff Hansen <cwhanse@sandia.gov>

* fix bug in read_tmy3 when year coerced

- so the indices are monotonically increasing
- add note to API changes and update bug fix in what's new
- update docstring in read_tmy3
  * explain in coerce_year parmaeter that all of the indices except the
  last one will have this, value, but the last index will be the value+1
- update the warnings to explain that if the year is coerced that the last index
 will be the next year
- do coerce year on date_ymd which is a Series of Timestamps, not an
 index so mutable (can assign by integer index) which makes it trivial
 to change the last index to be the next year
- revert most of the changes from the previous commits in this PR
  * remove `tmy3_monotonic_index` from `api.rst`, `iotools/__init__.py`, and from `tmy.py`
  * remove `monotonic_index` kwarg from read_tmy3
  * remove use of `tmy3_monotonic_index` from sphinx gallery example,
   test_losses, and test_tmy
  * remove test for tmy3_monotonic_index and the monotonic_index kwarg
  * fix test_tmy test to treat last index differently, except for the
  leap year test which now has a constant diff, yay!

* stickler: rm unused imports from tmy, test

* suggested changes in read_tmy3 for monotonic index

* sed s/data/index/g
* wordsmithing docstring: coerce_year kwarg
* wordsmithing what's new bug fixes
* also consistently use periods for bullets in what's new

* Apply suggestions from code review to test last hour in read_tmy3

* and minor wordsmithing in what's new

Co-Authored-By: Will Holmgren <william.holmgren@gmail.com>

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
Co-authored-by: Will Holmgren <william.holmgren@gmail.com>

Author: 3 people

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,22 @@
 # step.
 
 from datetime import datetime
+import pathlib
 from matplotlib import pyplot as plt
 from pvlib.iotools import read_tmy3
 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)
+# 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/whatsnew/v0.7.2.rst

@@ -7,42 +7,47 @@ API Changes
 ~~~~~~~~~~~
 * :py:class:`pvlib.forecast.ForecastModel` now requires ``start`` and ``end``
   arguments to be tz-localized. (:issue:`877`, :pull:`879`)
+* :py:func:`pvlib.iotools.read_tmy3` when coerced to a single year now returns
+  indices that are monotonically increasing. Therefore, the last index will be
+  January 1, 00:00 of the *next* year. (:pull:`910`)
 
 Enhancements
 ~~~~~~~~~~~~
 * TMY3 dataframe returned by :py:func:`~pvlib.iotools.read_tmy3` now contains
   the original ``Date (MM/DD/YYYY)`` and ``Time (HH:MM)`` columns that the
-  indices were parsed from (:pull:`866`)
+  indices were parsed from. (:pull:`866`)
 * Add :py:func:`~pvlib.pvsystem.PVSystem.faiman` and added 
   ``temperature_model='faiman'`` option to :py:class:`~pvlib.modelchain.ModelChain`
   (:pull:`897`) (:issue:`836`).
-* Add Kimber soiling model :py:func:`pvlib.losses.soiling_kimber` (:pull:`860`)
+* Add Kimber soiling model :py:func:`pvlib.losses.soiling_kimber`. (:pull:`860`)
 
 Bug fixes
 ~~~~~~~~~
 * Fix :py:func:`~pvlib.iotools.read_tmy3` parsing when February contains
-  a leap year (:pull:`866`)
+  a leap year. (:pull:`866`)
 * Implement NREL Developer Network API key for consistent success with API
-  calls in :py:mod:`pvlib.tests.iotools.test_psm3` (:pull:`873`)
+  calls in :py:mod:`pvlib.tests.iotools.test_psm3`. (:pull:`873`)
 * Fix issue with :py:class:`pvlib.location.Location` creation when
-  passing ``tz=datetime.timezone.utc`` (:pull:`879`)
+  passing ``tz=datetime.timezone.utc``. (:pull:`879`)
 * Fix documentation homepage title to "pvlib python" based on first heading on
   the page. (:pull:`890`) (:issue:`888`)
 * Implement `pytest-remotedata <https://github.com/astropy/pytest-remotedata>`_
   to increase test suite speed. Requires ``--remote-data`` pytest flag to
-  execute data retrieval tests over a network.(:issue:`882`)(:pull:`896`)
+  execute data retrieval tests over a network. (:issue:`882`)(:pull:`896`)
 * Fix missing
   `0.7.0 what's new <https://pvlib-python.readthedocs.io/en/stable/whatsnew.html#v0-7-0-december-18-2019>`_
   entries about changes to ``PVSystem.pvwatts_ac``. Delete unreleased
   0.6.4 what's new file. (:issue:`898`)
 * Compatibility with cftime 1.1. (:issue:`895`)
-* Add Python3.8 to Azure Pipelines CI (:issue:`903`)(:pull:`904`)
-* Add documentation build test to Azure Pipelines CI (:pull:`909`)
+* Add Python3.8 to Azure Pipelines CI. (:issue:`903`)(:pull:`904`)
+* Add documentation build test to Azure Pipelines CI. (:pull:`909`)
 * Minor implemention changes to avoid runtime and deprecation warnings in
   :py:func:`~pvlib.clearsky.detect_clearsky`,
   :py:func:`~pvlib.iam.martin_ruiz_diffuse`,
   :py:func:`~pvlib.losses.soiling_hsu`,
   and various test functions.
+* Fix :py:func:`~pvlib.iotools.read_tmy3` so that when coerced to a single year
+  the TMY3 index will be monotonically increasing. (:pull:`910`)
 
 Testing
 ~~~~~~~

pvlib/iotools/__init__.py

@@ -1,5 +1,4 @@
-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.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

@@ -28,10 +28,12 @@ def read_tmy3(filename=None, coerce_year=None, recolumn=True):
         a relative file path, absolute file path, or url.
 
     coerce_year : None or int, default None
-        If supplied, the year of the data will be set to this value.
+        If supplied, the year of the index will be set to `coerce_year`, except
+        for the last index value which will be set to the *next* year so that
+        the index increases monotonically.
 
     recolumn : bool, default True
-        If True, apply standard names to TMY3 columns. Typically this
+        If ``True``, apply standard names to TMY3 columns. Typically this
         results in stripping the units from the column name.
 
     Returns
@@ -49,7 +51,6 @@ def read_tmy3(filename=None, coerce_year=None, recolumn=True):
 
     Notes
     -----
-
     The returned structures have the following fields.
 
     ===============   ======  ===================
@@ -139,15 +140,16 @@ def read_tmy3(filename=None, coerce_year=None, recolumn=True):
     TMYData.PresWthUncertainty          Present weather code uncertainty, see [2]_.
     =============================       ======================================================================================================================================================
 
-    .. warning:: TMY3 irradiance data corresponds to the previous hour, so the
-        first hour is 1AM, corresponding to the net irradiance from midnite to
-        1AM, and the last hour is midnite of the *next* year, unless the year
-        has been coerced. EG: if TMY3 was 1988-12-31 24:00:00 this becomes
-        1989-01-01 00:00:00
+    .. warning:: TMY3 irradiance data corresponds to the *previous* hour, so
+        the first index is 1AM, corresponding to the irradiance from midnight
+        to 1AM, and the last index is midnight of the *next* year. For example,
+        if the last index in the TMY3 file was 1988-12-31 24:00:00 this becomes
+        1989-01-01 00:00:00 after calling :func:`~pvlib.iotools.read_tmy3`.
 
     .. warning:: When coercing the year, the last index in the dataframe will
-        be the first hour of the same year, EG: if TMY3 was 1988-12-31 24:00:00
-        and year is coerced to 1990 this becomes 1990-01-01
+        become midnight of the *next* year. For example, if the last index in
+        the TMY3 was 1988-12-31 24:00:00, and year is coerced to 1990 then this
+        becomes 1991-01-01 00:00:00.
 
     References
     ----------
@@ -214,11 +216,12 @@ def read_tmy3(filename=None, coerce_year=None, recolumn=True):
     data_ymd[leapday] += datetime.timedelta(days=1)
     # shifted_hour is a pd.Series, so use pd.to_timedelta to get a pd.Series of
     # timedeltas
+    if coerce_year is not None:
+        data_ymd = data_ymd.map(lambda dt: dt.replace(year=coerce_year))
+        data_ymd.iloc[-1] = data_ymd.iloc[-1].replace(year=coerce_year+1)
     # NOTE: as of pvlib-0.6.3, min req is pandas-0.18.1, so pd.to_timedelta
     # unit must be in (D,h,m,s,ms,us,ns), but pandas>=0.24 allows unit='hour'
     data.index = data_ymd + pd.to_timedelta(shifted_hour, unit='h')
-    if coerce_year is not None:
-        data.index = data.index.map(lambda dt: dt.replace(year=coerce_year))
 
     if recolumn:
         data = _recolumn(data)  # rename to standard column names

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

@@ -3,6 +3,7 @@
 import pandas as pd
 import pytest
 from pvlib.iotools import tmy
+from pvlib.iotools import read_tmy3
 from conftest import DATA_DIR
 
 # test the API works
@@ -30,21 +31,23 @@ def test_read_tmy3_recolumn():
 
 
 def test_read_tmy3_norecolumn():
-    data, meta = tmy.read_tmy3(TMY3_TESTFILE, recolumn=False)
+    data, _ = tmy.read_tmy3(TMY3_TESTFILE, recolumn=False)
     assert 'GHI source' in data.columns
 
 
 def test_read_tmy3_coerce_year():
     coerce_year = 1987
-    data, meta = tmy.read_tmy3(TMY3_TESTFILE, coerce_year=coerce_year)
-    assert (data.index.year == 1987).all()
+    data, _ = tmy.read_tmy3(TMY3_TESTFILE, coerce_year=coerce_year)
+    assert (data.index[:-1].year == 1987).all()
+    assert data.index[-1].year == 1988
 
 
 def test_read_tmy3_no_coerce_year():
     coerce_year = None
-    data, meta = tmy.read_tmy3(TMY3_TESTFILE, coerce_year=coerce_year)
+    data, _ = tmy.read_tmy3(TMY3_TESTFILE, coerce_year=coerce_year)
     assert 1997 and 1999 in data.index.year
-
+    assert data.index[-2] == pd.Timestamp('1998-12-31 23:00:00-09:00')
+    assert data.index[-1] == pd.Timestamp('1999-01-01 00:00:00-09:00')
 
 def test_read_tmy2():
     tmy.read_tmy2(TMY2_TESTFILE)
@@ -70,9 +73,10 @@ def test_gh865_read_tmy3_feb_leapyear_hr24():
     # now check if it parses correctly when we try to coerce the year
     data, _ = read_tmy3(TMY3_FEB_LEAPYEAR, coerce_year=1990)
     # if get's here w/o an error, then gh865 is fixed, but let's check anyway
-    assert all(data.index.year == 1990)
+    assert all(data.index[:-1].year == 1990)
+    assert data.index[-1].year == 1991
     # let's do a quick sanity check, are the indices monotonically increasing?
-    assert all(np.diff(data.index[:-1].astype(int)) == 3600000000000)
+    assert all(np.diff(data.index.astype(int)) == 3600000000000)
     # according to the TMY3 manual, each record corresponds to the previous
     # hour so check that the 1st hour is 1AM and the last hour is midnite
     assert data.index[0].hour == 1