From 677de4c190eca946ac868abd523afca6aeb7fe84 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Thu, 11 Feb 2021 14:04:06 -0700
Subject: [PATCH 01/31] First cut porting pvlib_detect_shadows from matlab

This is a rought implementation that returns almost the correct shadow
image. There is still some work to be done. The output needs to be
returned as a series instead of an image and steps should be taken to
minimize memory use, especially in the _clean_wires() function, by
passing ndarrays rather than allocating a new array for every step.

Currently this implementation appears to be more sensitive than the
matlab version, and includes more remnants of clouds in the output.
---
 pvanalytics/features/shading.py | 306 ++++++++++++++++++++++++++++++++
 1 file changed, 306 insertions(+)
 create mode 100644 pvanalytics/features/shading.py

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
new file mode 100644
index 00000000..c6583fd3
--- /dev/null
+++ b/pvanalytics/features/shading.py
@@ -0,0 +1,306 @@
+"""Functions for labeling shading and shodows."""
+import numpy as np
+import pandas as pd
+from scipy import ndimage
+from skimage import morphology, measure
+import pvlib
+from pvanalytics import util
+
+
+def _to_image(data, width):
+    """Convert data to an image.
+
+    Parameters
+    ----------
+    data : Series
+        Values of the pixels
+    width : int
+        Width of the image in pixels
+
+    Returns
+    -------
+    ndarray
+        Data reshaped into an ndarray with width `width`
+    """
+    return data.reshape(len(data) // width, width)
+
+
+def _smooth(image):
+    """Smooth `image` by local mean filtering.
+
+    Local means are computed in a 3x3 square surrounding each pixel.
+    """
+    return ndimage.uniform_filter(image, size=15)
+
+
+def _prepare_images(ghi, clearsky, daytime, interval):
+    """Prepare data as images.
+
+    Performs pre-processing steps on `ghi` and `clearsky` before
+    returning images for use in the shadow detection algorithm.
+
+    Returns
+    -------
+    ghi_image : np.ndarray
+        Image form of `ghi`
+    clearsky_image : np.ndarray
+        Image form of `clearsky`
+    clouds_image : np.ndarray
+        Image of the cloudy periods in `ghi`
+    image_times : pandas.DatetimeIndex
+        Index for the data included in the returned images. Leading
+        and trailing days with incomplete data are not included in the
+        image, these times are needed to build a Series from the image
+        later on.
+
+    """
+    # Fill missing days by interpolation. Interpolation is limited to
+    # times with valid data before and after; missing data at the
+    # beginning or end of the series is not filled in.
+    image_width = 1440 // interval
+    ghi = ghi.interpolate(limit_area='inside')
+    # drop incomplete days.
+    ghi = ghi[ghi.resample('D').transform('count') == image_width]
+    image_times = ghi.index
+    ghi_image = _to_image(ghi.to_numpy(), image_width)
+    scaled_ghi = (ghi * 1000) / _smooth(ghi_image).max()
+    scaled_clearsky = (clearsky * 1000) / clearsky.max()
+    scaled_clearsky = scaled_clearsky.reindex_like(scaled_ghi)
+    daytime = daytime.reindex_like(scaled_ghi)
+    # Detect clouds.
+    clouds = _detect_clouds(scaled_ghi, scaled_clearsky)
+    # Interpolate across days (i.e. along columns) to remove clouds
+    # replace clouds with nans
+    #
+    # This could probably be done directly with scipy.interpolate.inter1d,
+    # but the easiest approach is to turn the image into a dataframe and
+    # interpolate along the columns.
+    cloudless_image = ghi_image.copy()
+    cloudless_image[_to_image(clouds.to_numpy(), image_width)] = np.nan
+    clouds_image = ghi_image.copy()
+    clouds_image[~_to_image(clouds.to_numpy(), image_width)] = np.nan
+    ghi_image = pd.DataFrame(cloudless_image).interpolate(
+        axis=0,
+        limit_direction='both'
+    ).to_numpy()
+    # set night to nan
+    ghi_image[~_to_image(daytime.to_numpy(), image_width)] = np.nan
+    return (
+        ghi_image,
+        _to_image(scaled_clearsky.to_numpy(), image_width),
+        clouds_image,
+        image_times
+    )
+
+
+def _detect_clouds(ghi, clearsky_ghi):
+    """Use :py:func:`pvanalytics.clearsky.reno` to detect clouds"""
+    clouds = pvlib.clearsky.detect_clearsky(
+        ghi,
+        clearsky_ghi,
+        ghi.index,
+        window_length=10
+    )
+    return clouds.rolling('60T').sum() == 0
+
+
+def _remove_pillars(wires):
+    j = np.array([[0, 0, 0, 0, 0],
+                  [0, 0, 1, 0, 0],
+                  [0, 0, 1, 0, 0]])
+    k = np.array([[1, 0, 0, 0, 1],
+                  [1, 0, 0, 0, 1],
+                  [1, 0, 0, 0, 1]])
+    a1 = np.logical_and(
+        wires,
+        np.logical_not(ndimage.binary_hit_or_miss(wires, k, j))
+    )
+    a2 = np.logical_and(
+        wires,
+        np.logical_not(
+            ndimage.binary_hit_or_miss(wires, np.flipud(k), np.flipud(j))
+        )
+    )
+    return np.logical_and(a1, a2)
+
+
+def _fill_gaps(wires):
+    """Fill 1-minute gaps."""
+    mask = np.array([[1, 0, 1]])
+    return np.logical_or(wires, ndimage.binary_hit_or_miss(wires, mask))
+
+
+def _remove_spikes(wires):
+    """Remove 1-day spikes."""
+    mask = np.array([[0, 1, 0]]).T
+    return np.logical_and(
+        wires,
+        np.logical_not(ndimage.binary_hit_or_miss(wires, mask))
+    )
+
+
+def _restore_gaps(wires):
+    """Restore 2-minute and 1-minute gaps"""
+    wires = np.logical_or(
+        wires,
+        ndimage.binary_hit_or_miss(
+            wires,
+            np.array([[0, 1, 0, 0, 1]]), np.array([[0, 0, 1, 1, 0]])
+        )
+    )
+    return np.logical_or(
+        wires,
+        ndimage.binary_hit_or_miss(wires, np.array([[1, 0, 1]]))
+    )
+
+
+def _width(bbox):
+    """Return the width of a bounding box.
+
+    Parameters
+    ----------
+    bbox : tuple
+        4-tuple specifying ``(min_row, min_col, max_row, max_col)``.
+
+    Returns
+    -------
+    int
+        The width of the bounding box in pixels.
+    """
+    return bbox[3] - bbox[1] + 1
+
+
+def _filter_blobs(wires, blob_length, connectivity):
+    """Remove blobs shorter than blob_length.
+
+    Examines all connected components in the image and removes any
+    that span fewer than `blob_length` columns.
+    """
+    components = measure.label(wires, connectivity=connectivity, background=0)
+    keep = [props.label for props in measure.regionprops(components)
+            if _width(props.bbox) > blob_length]
+    return np.isin(components, keep)
+
+
+def _tand(degrees):
+    """Tangent of an angle in degrees."""
+    return np.tan(np.deg2rad(degrees))
+
+
+def _filter_bars(wires):
+    """Only keep shadows that are 'bar-shaped' and span multiple days."""
+    for angle in range(0,90,5):
+        if angle < 75:
+            dim = 21
+            mid = 10
+        else:
+            dim = round(10 + 2 * angle // 5)
+            mid = round(6 + angle // 5)
+        cc = np.arange(0, dim)
+        se = np.zeros((dim, dim))
+        if angle <= 45:
+            rr = mid + np.round((cc - mid) * _tand(angle))
+            se[cc, rr.astype('int')] = 1
+        else:
+            rr = mid + np.round((cc - mid) * _tand(90 - angle))
+            se[cc, rr.astype('int')] = 1
+        wires = np.logical_or(
+            wires,
+            morphology.binary_opening(wires, se.T)
+        )
+        wires = np.logical_or(
+            wires,
+            morphology.binary_opening(wires, np.flipud(se.T))
+        )
+    return wires
+
+
+def _clean_wires(wires):
+    """Clean up clouds that are connected to wires."""
+    wires = _remove_pillars(wires)
+    wires = _fill_gaps(wires)
+    wires = _remove_spikes(wires)
+    wires = _fill_gaps(wires)
+    wires = _restore_gaps(wires)
+    wires = _filter_blobs(wires, 20, connectivity=1)
+    wires = _filter_bars(wires)
+    wires = _filter_blobs(wires, 15, connectivity=2)
+    return wires
+
+
+def fixed(ghi, daytime, clearsky, interval=None):
+    """Detects shadows form fixed structures such as wires and poles.
+
+    Uses morphological image processing methods to identify shadows
+    from fixed local objects in GHI data. GHI data are are assumed to
+    be reasonably complete with relatively few missing values and at a
+    fixed time interval nominally of 1 minute over the course of
+    several months. Shadows are assumed to be relatively short
+    duration. The algorithm forms a 2D image of the GHI data by
+    arranging time of day along the x-axis and day of year along the
+    y-axis. Rapid change in GHI in the x-direction is used to identify
+    edges of shadows; continuity in the y-direction is used to
+    separate local object shading from cloud shadows.
+
+    Parameters
+    ----------
+    ghi : Series
+        Time series of GHI measurements. Should be at 1-minute frequency
+        for best results. Must be 1 year of data.
+    daytime : Series
+        Boolean series with True for times when the sun is up.
+    clearsky : Series
+        Clearsky GHI.
+    interval : int, optional
+        Interval between data points in minutes. If not specified the
+        interval is inferred from the frequency of the index of `ghi`.
+
+    Returns
+    -------
+    Series
+        Boolean series with true for times that are impacted by shadows.
+
+    References
+    ----------
+    References
+    .. [1] Martin, C. E., Hansen, C. W., An Image Processing Algorithm to
+       Identify Near-Field Shading in Irradiance Measurements, preprint 2016
+    .. [2] Reno, M.J. and C.W. Hansen, "Identification of periods of clear sky
+       irradiance in time series of GHI measurements" Renewable Energy, v90,
+       p. 520-531, 2016.
+    """
+    if interval is None:
+        interval = util.freq_to_timedelta(
+            pd.infer_freq(ghi.index)
+        ).seconds // 60
+    # TODO Handle invalid intervals
+    ghi_image, clearsky_image, clouds_image, index = _prepare_images(
+        ghi,
+        clearsky,
+        daytime,
+        interval
+    )
+    # normalize the GHI and dampen the dynamic range where the clear
+    # sky model may have large errors (e.g. at very low sun elevation)
+    alpha = 2000
+    ghi_boosted = 1000 * (ghi_image + alpha) / (clearsky_image + alpha)
+
+    # We must use scipy.ndimage here because skimage does not support
+    # floating point the data outside the range [-1, 1]
+    gradient = ndimage.morphological_gradient(ghi_boosted, size=(1, 3))
+    min_gradient = 2  # TODO consider making this a parameter
+    threshold = gradient > min_gradient  # binary image of wire candidates
+
+    # From here we CAN use skimage we are working with binary
+    # images. We use the grayscale version of these funcitons rather
+    # than the binary version since the binary version do not support
+    # the 'area_threshold' parameter.
+    three_minute_mask = morphology.rectangle(1, 3)
+    wires = morphology.area_opening(
+        morphology.closing(threshold, three_minute_mask),
+        area_threshold=200,
+        connectivity=2  # all neighbors (including diagonals)
+    )
+    wires = _clean_wires(wires)
+    # TODO convert wires back to a time-series indexed by ghi.index
+    return wires, ghi_image, gradient

From 2fcfb9e3b8d6d3fd7650dc30a4e372541d0fd61c Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Mon, 1 Mar 2021 11:03:53 -0700
Subject: [PATCH 02/31] Add test module for features.shading

Initially adds basic tests for data without any fixed shadows.
---
 pvanalytics/features/__init__.py           |  1 +
 pvanalytics/tests/features/test_shading.py | 51 ++++++++++++++++++++++
 2 files changed, 52 insertions(+)
 create mode 100644 pvanalytics/tests/features/test_shading.py

diff --git a/pvanalytics/features/__init__.py b/pvanalytics/features/__init__.py
index dfb7a12d..9a8ff7c5 100644
--- a/pvanalytics/features/__init__.py
+++ b/pvanalytics/features/__init__.py
@@ -1,2 +1,3 @@
 from pvanalytics.features import clearsky  # noqa: F401
 from pvanalytics.features import clipping  # noqa: F401
+from pvanalytics.features import shading   # noqa: F401
diff --git a/pvanalytics/tests/features/test_shading.py b/pvanalytics/tests/features/test_shading.py
new file mode 100644
index 00000000..421adfd1
--- /dev/null
+++ b/pvanalytics/tests/features/test_shading.py
@@ -0,0 +1,51 @@
+import pandas as pd
+from pandas.testing import assert_series_equal
+import pytest
+from pvanalytics.features import shading
+
+
+# TODO testing plan
+#      - [ ] index of output same as input
+#      - [ ] no shadows -> all False
+#      - [ ] 5-minute timestamps
+#      - [ ] unaligned timestamps (?)
+#      - [ ] partial days at start and end of series
+
+
+@pytest.fixture(scope='module')
+def times():
+    return pd.date_range(
+        start='1/1/2020',
+        end='12/31/2020',
+        closed='left',
+        freq='T',
+        tz='MST'
+    )
+
+
+@pytest.fixture(scope='module')
+def daytime(times, albuquerque):
+    solar_position = albuquerque.get_solarposition(times)
+    return solar_position['zenith'] < 87
+
+
+@pytest.fixture(scope='module')
+def clearsky_ghi(times, albuquerque):
+    clearsky = albuquerque.get_clearsky(times, model='simplified_solis')
+    return clearsky['ghi']
+
+
+def test_fixed_no_shadows(daytime, clearsky_ghi):
+    shadows = shading.fixed(clearsky_ghi, daytime, clearsky_ghi)
+    assert not shadows.any()
+
+
+def test_fixed_same_index(daytime, clearsky_ghi):
+    daytime = daytime['1/2/2020 11:00':'11/2/2020 13:00']
+    clearsky_ghi = clearsky_ghi['1/2/2020 11:00':'11/2/2020 13:00']
+    shadows = shading.fixed(daytime, clearsky_ghi, clearsky_ghi)
+    assert_series_equal(
+        pd.series(False, index=daytime.index),
+        shadows,
+        check_names=False
+    )

From 5facb7732f2ce81e788ce772876b27d77875eaba Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Mon, 1 Mar 2021 11:08:40 -0700
Subject: [PATCH 03/31] Correct order of masks passed to binary hit/miss filter

In the remove_pillars step the masks were passed in the wrong
order.
---
 pvanalytics/features/shading.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index c6583fd3..80b5c7fc 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -113,12 +113,12 @@ def _remove_pillars(wires):
                   [1, 0, 0, 0, 1]])
     a1 = np.logical_and(
         wires,
-        np.logical_not(ndimage.binary_hit_or_miss(wires, k, j))
+        np.logical_not(ndimage.binary_hit_or_miss(wires, j, k))
     )
     a2 = np.logical_and(
         wires,
         np.logical_not(
-            ndimage.binary_hit_or_miss(wires, np.flipud(k), np.flipud(j))
+            ndimage.binary_hit_or_miss(wires, np.flipud(j), np.flipud(k))
         )
     )
     return np.logical_and(a1, a2)

From df75c49657b80412295763da1bfa1e6ae87591c2 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Mon, 1 Mar 2021 11:15:08 -0700
Subject: [PATCH 04/31] Filter blobs by height (e.g. days) rather than width

This is consistent with the comments in the original matlab
code, although it is not 100% clear to me that it matches the
implementation. Since it produces results that are closer to
the output (in addition to matching the comments) from the MATLAB
code, I will assume that my reading of the original
implementation was not correct.
---
 pvanalytics/features/shading.py | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 80b5c7fc..3324971e 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -170,6 +170,22 @@ def _width(bbox):
     return bbox[3] - bbox[1] + 1
 
 
+def _height(bbox):
+    """Return the height of a bounding box.
+
+    Parameters
+    ----------
+    bbox : tuple
+        4-tuple specifying ``(min_row, min_col, max_row, max_col)``.
+
+    Returns
+    -------
+    int
+        The height of the bounding box in pixels.
+    """
+    return bbox[2] - bbox[0] + 1
+
+
 def _filter_blobs(wires, blob_length, connectivity):
     """Remove blobs shorter than blob_length.
 
@@ -178,7 +194,7 @@ def _filter_blobs(wires, blob_length, connectivity):
     """
     components = measure.label(wires, connectivity=connectivity, background=0)
     keep = [props.label for props in measure.regionprops(components)
-            if _width(props.bbox) > blob_length]
+            if _height(props.bbox) > blob_length]
     return np.isin(components, keep)
 
 

From a532ea65981f8ce4f2307ded419d7849b4c6e2f7 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Mon, 1 Mar 2021 11:16:37 -0700
Subject: [PATCH 05/31] Clean up variable naming and comments in
 shading.fixed()

---
 pvanalytics/features/shading.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 3324971e..bef6c135 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -308,12 +308,12 @@ def fixed(ghi, daytime, clearsky, interval=None):
     threshold = gradient > min_gradient  # binary image of wire candidates
 
     # From here we CAN use skimage we are working with binary
-    # images. We use the grayscale version of these funcitons rather
+    # images. We use the grayscale version of these functions rather
     # than the binary version since the binary version do not support
     # the 'area_threshold' parameter.
-    three_minute_mask = morphology.rectangle(1, 3)
+    three_day_mask = morphology.rectangle(1, 3)
     wires = morphology.area_opening(
-        morphology.closing(threshold, three_minute_mask),
+        morphology.closing(threshold, three_day_mask),
         area_threshold=200,
         connectivity=2  # all neighbors (including diagonals)
     )

From 370fc06a9422cb038ebc6195e651f53143b2d5ec Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Mon, 1 Mar 2021 11:18:31 -0700
Subject: [PATCH 06/31] Tweak shading._filter_bars()

Changes to better match the original implementation. Still not
completely sure that this correct though.
---
 pvanalytics/features/shading.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index bef6c135..83725bd4 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -101,7 +101,7 @@ def _detect_clouds(ghi, clearsky_ghi):
         ghi.index,
         window_length=10
     )
-    return clouds.rolling('60T').sum() == 0
+    return clouds.rolling('50T').sum() == 0
 
 
 def _remove_pillars(wires):
@@ -205,12 +205,12 @@ def _tand(degrees):
 
 def _filter_bars(wires):
     """Only keep shadows that are 'bar-shaped' and span multiple days."""
-    for angle in range(0,90,5):
+    for angle in range(0, 90, 5):
         if angle < 75:
             dim = 21
-            mid = 10
+            mid = 11
         else:
-            dim = round(10 + 2 * angle // 5)
+            dim = round(11 + 2 * angle // 5)
             mid = round(6 + angle // 5)
         cc = np.arange(0, dim)
         se = np.zeros((dim, dim))
@@ -220,13 +220,14 @@ def _filter_bars(wires):
         else:
             rr = mid + np.round((cc - mid) * _tand(90 - angle))
             se[cc, rr.astype('int')] = 1
+            se = se.T
         wires = np.logical_or(
             wires,
-            morphology.binary_opening(wires, se.T)
+            morphology.binary_opening(wires, se)
         )
         wires = np.logical_or(
             wires,
-            morphology.binary_opening(wires, np.flipud(se.T))
+            morphology.binary_opening(wires, np.flipud(se))
         )
     return wires
 

From 25144fb2a25b559fb6091d44c419295c30ae768d Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 10 Mar 2021 10:58:36 -0700
Subject: [PATCH 07/31] Apply _filter_bars() to the raw shadow image

Rather than applying it on successively more cleaned up images.
---
 pvanalytics/features/shading.py | 36 ++++++++++++++++-----------------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 83725bd4..0db7b2c4 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -203,14 +203,14 @@ def _tand(degrees):
     return np.tan(np.deg2rad(degrees))
 
 
-def _filter_bars(wires):
+def _filter_bars(wires, out):
     """Only keep shadows that are 'bar-shaped' and span multiple days."""
     for angle in range(0, 90, 5):
         if angle < 75:
-            dim = 21
-            mid = 11
+            dim = 20
+            mid = 10
         else:
-            dim = round(11 + 2 * angle // 5)
+            dim = round(10 + 2 * angle // 5)
             mid = round(6 + angle // 5)
         cc = np.arange(0, dim)
         se = np.zeros((dim, dim))
@@ -221,28 +221,28 @@ def _filter_bars(wires):
             rr = mid + np.round((cc - mid) * _tand(90 - angle))
             se[cc, rr.astype('int')] = 1
             se = se.T
-        wires = np.logical_or(
-            wires,
+        out = np.logical_or(
+            out,
             morphology.binary_opening(wires, se)
         )
-        wires = np.logical_or(
-            wires,
+        out = np.logical_or(
+            out,
             morphology.binary_opening(wires, np.flipud(se))
         )
-    return wires
+    return out
 
 
 def _clean_wires(wires):
     """Clean up clouds that are connected to wires."""
-    wires = _remove_pillars(wires)
-    wires = _fill_gaps(wires)
-    wires = _remove_spikes(wires)
-    wires = _fill_gaps(wires)
-    wires = _restore_gaps(wires)
-    wires = _filter_blobs(wires, 20, connectivity=1)
-    wires = _filter_bars(wires)
-    wires = _filter_blobs(wires, 15, connectivity=2)
-    return wires
+    out = _remove_pillars(wires)
+    out = _fill_gaps(out)
+    out = _remove_spikes(out)
+    out = _fill_gaps(out)
+    out = _restore_gaps(out)
+    out = _filter_blobs(out, 20, connectivity=1)
+    out = _filter_bars(wires, out)
+    out = _filter_blobs(out, 15, connectivity=2)
+    return out
 
 
 def fixed(ghi, daytime, clearsky, interval=None):

From e94b2f8e363b1a8edbc00a20f6230e9f40b48faf Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 10 Mar 2021 11:00:31 -0700
Subject: [PATCH 08/31] Make min_gradient a parameter to shading.fixed()

---
 pvanalytics/features/shading.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 0db7b2c4..b3ccee3a 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -245,7 +245,7 @@ def _clean_wires(wires):
     return out
 
 
-def fixed(ghi, daytime, clearsky, interval=None):
+def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     """Detects shadows form fixed structures such as wires and poles.
 
     Uses morphological image processing methods to identify shadows
@@ -305,7 +305,6 @@ def fixed(ghi, daytime, clearsky, interval=None):
     # We must use scipy.ndimage here because skimage does not support
     # floating point the data outside the range [-1, 1]
     gradient = ndimage.morphological_gradient(ghi_boosted, size=(1, 3))
-    min_gradient = 2  # TODO consider making this a parameter
     threshold = gradient > min_gradient  # binary image of wire candidates
 
     # From here we CAN use skimage we are working with binary

From 54ac08e10ee20b22ac334cb5f145b081be42f9da Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 10 Mar 2021 11:02:20 -0700
Subject: [PATCH 09/31] Use morphology.remove_small_objects in initial shadow
 detection

This should be slightly more efficient, plus it matched the intent
of the code better.

Also updates name of structuring element to match its shape
---
 pvanalytics/features/shading.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index b3ccee3a..873a97bd 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -307,14 +307,11 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     gradient = ndimage.morphological_gradient(ghi_boosted, size=(1, 3))
     threshold = gradient > min_gradient  # binary image of wire candidates
 
-    # From here we CAN use skimage we are working with binary
-    # images. We use the grayscale version of these functions rather
-    # than the binary version since the binary version do not support
-    # the 'area_threshold' parameter.
-    three_day_mask = morphology.rectangle(1, 3)
-    wires = morphology.area_opening(
-        morphology.closing(threshold, three_day_mask),
-        area_threshold=200,
+    # From here we CAN use skimage we are working with binary images.
+    three_minute_mask = morphology.rectangle(1, 3)
+    wires = morphology.remove_small_objects(
+        morphology.binary_closing(threshold, three_minute_mask),
+        min_size=200,
         connectivity=2  # all neighbors (including diagonals)
     )
     wires = _clean_wires(wires)

From 36f3f60ff6be9f58e7d2ea6209972b2bca9c24bb Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 10 Mar 2021 11:04:41 -0700
Subject: [PATCH 10/31] Make the minimum duration of cloudy times a parameter

---
 pvanalytics/features/shading.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 873a97bd..ee9c9e52 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -68,7 +68,7 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     scaled_clearsky = scaled_clearsky.reindex_like(scaled_ghi)
     daytime = daytime.reindex_like(scaled_ghi)
     # Detect clouds.
-    clouds = _detect_clouds(scaled_ghi, scaled_clearsky)
+    clouds = _detect_clouds(scaled_ghi, scaled_clearsky, '50T')
     # Interpolate across days (i.e. along columns) to remove clouds
     # replace clouds with nans
     #
@@ -93,7 +93,7 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     )
 
 
-def _detect_clouds(ghi, clearsky_ghi):
+def _detect_clouds(ghi, clearsky_ghi, window='50T'):
     """Use :py:func:`pvanalytics.clearsky.reno` to detect clouds"""
     clouds = pvlib.clearsky.detect_clearsky(
         ghi,
@@ -101,7 +101,7 @@ def _detect_clouds(ghi, clearsky_ghi):
         ghi.index,
         window_length=10
     )
-    return clouds.rolling('50T').sum() == 0
+    return clouds.rolling(window).sum() == 0
 
 
 def _remove_pillars(wires):

From f4828275ba1ac6bc065afb2a59fd5dfc9a1202c8 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 24 Mar 2021 09:55:21 -0600
Subject: [PATCH 11/31] Return shadows as a series and an image

Converts the image back to a series and re-indexes it like
the input series.
---
 pvanalytics/features/shading.py            | 24 ++++++++++++++++++++--
 pvanalytics/tests/features/test_shading.py |  4 ++--
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index ee9c9e52..a7dc9674 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -25,6 +25,24 @@ def _to_image(data, width):
     return data.reshape(len(data) // width, width)
 
 
+def _to_series(image, image_index):
+    """Convert a binary image to a boolean series.
+
+    Parameters
+    ----------
+    image : ndarray
+        Binary image.
+    image_index : pd.Index
+        Index for the returned series.
+
+    Returns
+    -------
+    Series
+        A pandas Series with index `index` and values from the image.
+    """
+    return pd.Series(image.flatten(), index=image_index)
+
+
 def _smooth(image):
     """Smooth `image` by local mean filtering.
 
@@ -314,6 +332,8 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
         min_size=200,
         connectivity=2  # all neighbors (including diagonals)
     )
-    wires = _clean_wires(wires)
+    wires_image = _clean_wires(wires)
     # TODO convert wires back to a time-series indexed by ghi.index
-    return wires, ghi_image, gradient
+    wires_series = _to_series(wires, index)
+    wires_series = wires_series.reindex(ghi.index, fill_value=False)
+    return wires_series, wires_image
diff --git a/pvanalytics/tests/features/test_shading.py b/pvanalytics/tests/features/test_shading.py
index 421adfd1..abf9dda9 100644
--- a/pvanalytics/tests/features/test_shading.py
+++ b/pvanalytics/tests/features/test_shading.py
@@ -36,14 +36,14 @@ def clearsky_ghi(times, albuquerque):
 
 
 def test_fixed_no_shadows(daytime, clearsky_ghi):
-    shadows = shading.fixed(clearsky_ghi, daytime, clearsky_ghi)
+    shadows, image = shading.fixed(clearsky_ghi, daytime, clearsky_ghi)
     assert not shadows.any()
 
 
 def test_fixed_same_index(daytime, clearsky_ghi):
     daytime = daytime['1/2/2020 11:00':'11/2/2020 13:00']
     clearsky_ghi = clearsky_ghi['1/2/2020 11:00':'11/2/2020 13:00']
-    shadows = shading.fixed(daytime, clearsky_ghi, clearsky_ghi)
+    shadows, image = shading.fixed(daytime, clearsky_ghi, clearsky_ghi)
     assert_series_equal(
         pd.series(False, index=daytime.index),
         shadows,

From 7ea77fb8a645918177e5a14532428426b47a4083 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 24 Mar 2021 10:10:57 -0600
Subject: [PATCH 12/31] Fix errors in test_shading.test_fixed_same_index

---
 pvanalytics/tests/features/test_shading.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvanalytics/tests/features/test_shading.py b/pvanalytics/tests/features/test_shading.py
index abf9dda9..9973acf5 100644
--- a/pvanalytics/tests/features/test_shading.py
+++ b/pvanalytics/tests/features/test_shading.py
@@ -43,9 +43,9 @@ def test_fixed_no_shadows(daytime, clearsky_ghi):
 def test_fixed_same_index(daytime, clearsky_ghi):
     daytime = daytime['1/2/2020 11:00':'11/2/2020 13:00']
     clearsky_ghi = clearsky_ghi['1/2/2020 11:00':'11/2/2020 13:00']
-    shadows, image = shading.fixed(daytime, clearsky_ghi, clearsky_ghi)
+    shadows, image = shading.fixed(clearsky_ghi, daytime, clearsky_ghi)
     assert_series_equal(
-        pd.series(False, index=daytime.index),
+        pd.Series(False, index=daytime.index),
         shadows,
         check_names=False
     )

From 574afa008a9ca8c5517103aadcfdc509d95297f1 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 24 Mar 2021 10:34:53 -0600
Subject: [PATCH 13/31] Test that a fake shadow at the same time every day is
 detected

---
 pvanalytics/tests/features/test_shading.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/pvanalytics/tests/features/test_shading.py b/pvanalytics/tests/features/test_shading.py
index 9973acf5..37e2db21 100644
--- a/pvanalytics/tests/features/test_shading.py
+++ b/pvanalytics/tests/features/test_shading.py
@@ -49,3 +49,10 @@ def test_fixed_same_index(daytime, clearsky_ghi):
         shadows,
         check_names=False
     )
+
+
+def test_simple_shadow(daytime, clearsky_ghi):
+    shadow_ghi = clearsky_ghi.copy()
+    shadow_ghi[shadow_ghi.between_time('11:00', '11:03').index] *= 0.5
+    shadows, image = shading.fixed(shadow_ghi, daytime, clearsky_ghi)
+    assert shadows.between_time('11:00', '11:03').all()

From 1040f8550042c207dcfd37f87a879a2440988556 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 24 Mar 2021 10:54:21 -0600
Subject: [PATCH 14/31] Fix grammar in comment

---
 pvanalytics/features/shading.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index a7dc9674..6a0904be 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -325,7 +325,7 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     gradient = ndimage.morphological_gradient(ghi_boosted, size=(1, 3))
     threshold = gradient > min_gradient  # binary image of wire candidates
 
-    # From here we CAN use skimage we are working with binary images.
+    # From here we CAN use skimage because we are working with binary images.
     three_minute_mask = morphology.rectangle(1, 3)
     wires = morphology.remove_small_objects(
         morphology.binary_closing(threshold, three_minute_mask),
@@ -333,7 +333,6 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
         connectivity=2  # all neighbors (including diagonals)
     )
     wires_image = _clean_wires(wires)
-    # TODO convert wires back to a time-series indexed by ghi.index
     wires_series = _to_series(wires, index)
     wires_series = wires_series.reindex(ghi.index, fill_value=False)
     return wires_series, wires_image

From a130e783fe8e095669dec58db3e2e71944f2694d Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Wed, 24 Mar 2021 11:03:16 -0600
Subject: [PATCH 15/31] Add scikit-image to dependencies

---
 requirements-min.txt | 1 +
 requirements.txt     | 1 +
 setup.py             | 3 ++-
 3 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/requirements-min.txt b/requirements-min.txt
index c27b656e..e7d61919 100644
--- a/requirements-min.txt
+++ b/requirements-min.txt
@@ -3,3 +3,4 @@ pandas~=0.23.0
 pvlib~=0.8.0
 scipy~=1.2.0
 statsmodels~=0.9.0
+scikit-image~=0.16.0
diff --git a/requirements.txt b/requirements.txt
index 07e9cf86..71c90084 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,3 +3,4 @@ pandas>=0.23.0,<1.1.0
 pvlib>=0.8.0
 scipy>=1.2.0
 statsmodels>=0.9.0
+scikit-image>=0.16.0
diff --git a/setup.py b/setup.py
index da684ae4..fe070ba7 100644
--- a/setup.py
+++ b/setup.py
@@ -35,7 +35,8 @@
     'pandas >= 0.23.0',
     'pvlib >= 0.8.0',
     'scipy >= 1.2.0',
-    'statsmodels >= 0.9.0'
+    'statsmodels >= 0.9.0',
+    'scikit-image >= 0.16.0'
 ]
 
 DOCS_REQUIRE = [

From d88e19c9644a4815e10cc41862dd76dc7f3d3198 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 30 Mar 2021 09:28:27 -0600
Subject: [PATCH 16/31] Use a 3x3 mask in shading._smooth()

---
 pvanalytics/features/shading.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 6a0904be..2252d1f3 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -48,7 +48,7 @@ def _smooth(image):
 
     Local means are computed in a 3x3 square surrounding each pixel.
     """
-    return ndimage.uniform_filter(image, size=15)
+    return ndimage.uniform_filter(image, size=(3, 3))
 
 
 def _prepare_images(ghi, clearsky, daytime, interval):

From ef32bece0077414df8ce3973ad6e77e33cd99999 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 30 Mar 2021 09:34:44 -0600
Subject: [PATCH 17/31] Remove extra calls to _to_image() in
 shading._prepare_images()

---
 pvanalytics/features/shading.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 2252d1f3..454c0fad 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -81,12 +81,13 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     ghi = ghi[ghi.resample('D').transform('count') == image_width]
     image_times = ghi.index
     ghi_image = _to_image(ghi.to_numpy(), image_width)
-    scaled_ghi = (ghi * 1000) / _smooth(ghi_image).max()
+    scaled_ghi = (ghi * 1000) / np.max(_smooth(ghi_image))
     scaled_clearsky = (clearsky * 1000) / clearsky.max()
     scaled_clearsky = scaled_clearsky.reindex_like(scaled_ghi)
     daytime = daytime.reindex_like(scaled_ghi)
     # Detect clouds.
     clouds = _detect_clouds(scaled_ghi, scaled_clearsky, '50T')
+    cloud_mask = _to_image(clouds.to_numpy(), image_width)
     # Interpolate across days (i.e. along columns) to remove clouds
     # replace clouds with nans
     #
@@ -94,9 +95,9 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     # but the easiest approach is to turn the image into a dataframe and
     # interpolate along the columns.
     cloudless_image = ghi_image.copy()
-    cloudless_image[_to_image(clouds.to_numpy(), image_width)] = np.nan
+    cloudless_image[cloud_mask] = np.nan
     clouds_image = ghi_image.copy()
-    clouds_image[~_to_image(clouds.to_numpy(), image_width)] = np.nan
+    clouds_image[~cloud_mask] = np.nan
     ghi_image = pd.DataFrame(cloudless_image).interpolate(
         axis=0,
         limit_direction='both'

From 1dd87da7c15bdb6957180a54dd500bbfea42c831 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 30 Mar 2021 10:09:29 -0600
Subject: [PATCH 18/31] Docstring and comment updates

---
 pvanalytics/features/shading.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 454c0fad..4d257601 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -72,9 +72,9 @@ def _prepare_images(ghi, clearsky, daytime, interval):
         later on.
 
     """
-    # Fill missing days by interpolation. Interpolation is limited to
-    # times with valid data before and after; missing data at the
-    # beginning or end of the series is not filled in.
+    # Fill missing times by interpolation. Missing data at the
+    # beginning or end of the series is not filled in, and will be
+    # excluded from the images used for shadow detection.
     image_width = 1440 // interval
     ghi = ghi.interpolate(limit_area='inside')
     # drop incomplete days.
@@ -281,12 +281,12 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     Parameters
     ----------
     ghi : Series
-        Time series of GHI measurements. Should be at 1-minute frequency
-        for best results. Must be 1 year of data.
+        Time series of GHI measurements. Data must be at 1-minute frequency
+        and cover at least several months.
     daytime : Series
         Boolean series with True for times when the sun is up.
     clearsky : Series
-        Clearsky GHI.
+        Clearsky GHI with same index as `ghi`.
     interval : int, optional
         Interval between data points in minutes. If not specified the
         interval is inferred from the frequency of the index of `ghi`.
@@ -322,7 +322,7 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     ghi_boosted = 1000 * (ghi_image + alpha) / (clearsky_image + alpha)
 
     # We must use scipy.ndimage here because skimage does not support
-    # floating point the data outside the range [-1, 1]
+    # floating point data outside the range [-1, 1].
     gradient = ndimage.morphological_gradient(ghi_boosted, size=(1, 3))
     threshold = gradient > min_gradient  # binary image of wire candidates
 

From b35d37b79d2299447cf871d95bb1a6f3ef2b19d0 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 30 Mar 2021 10:19:27 -0600
Subject: [PATCH 19/31] Raise a ValueError for data with invalid intervals

---
 pvanalytics/features/shading.py            |  3 ++-
 pvanalytics/tests/features/test_shading.py | 17 +++++++++--------
 2 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 4d257601..65049ff1 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -309,7 +309,8 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
         interval = util.freq_to_timedelta(
             pd.infer_freq(ghi.index)
         ).seconds // 60
-    # TODO Handle invalid intervals
+    if interval != 1:
+        raise ValueError("Data must be at 1-minute intervals")
     ghi_image, clearsky_image, clouds_image, index = _prepare_images(
         ghi,
         clearsky,
diff --git a/pvanalytics/tests/features/test_shading.py b/pvanalytics/tests/features/test_shading.py
index 37e2db21..311594b5 100644
--- a/pvanalytics/tests/features/test_shading.py
+++ b/pvanalytics/tests/features/test_shading.py
@@ -4,14 +4,6 @@
 from pvanalytics.features import shading
 
 
-# TODO testing plan
-#      - [ ] index of output same as input
-#      - [ ] no shadows -> all False
-#      - [ ] 5-minute timestamps
-#      - [ ] unaligned timestamps (?)
-#      - [ ] partial days at start and end of series
-
-
 @pytest.fixture(scope='module')
 def times():
     return pd.date_range(
@@ -56,3 +48,12 @@ def test_simple_shadow(daytime, clearsky_ghi):
     shadow_ghi[shadow_ghi.between_time('11:00', '11:03').index] *= 0.5
     shadows, image = shading.fixed(shadow_ghi, daytime, clearsky_ghi)
     assert shadows.between_time('11:00', '11:03').all()
+
+
+def test_invalid_interval(daytime, clearsky_ghi):
+    ghi = clearsky_ghi.resample('5T').first()
+    daytime_resampled = daytime.resample('5T').first()
+    with pytest.raises(ValueError, match="Data must be at 1-minute intervals"):
+        shading.fixed(ghi, daytime_resampled, ghi)
+    with pytest.raises(ValueError, match="Data must be at 1-minute intervals"):
+        shading.fixed(daytime, clearsky_ghi, clearsky_ghi, interval=2)

From 8522144cc2b43932662fe080a3716a2589a21ecf Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 30 Mar 2021 10:27:42 -0600
Subject: [PATCH 20/31] Apply docstring edits from code review

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 pvanalytics/features/shading.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 65049ff1..517fe8ca 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -265,13 +265,13 @@ def _clean_wires(wires):
 
 
 def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
-    """Detects shadows form fixed structures such as wires and poles.
+    """Detects shadows from fixed structures such as wires and poles.
 
     Uses morphological image processing methods to identify shadows
-    from fixed local objects in GHI data. GHI data are are assumed to
+    from fixed local objects in GHI data. GHI data are assumed to
     be reasonably complete with relatively few missing values and at a
     fixed time interval nominally of 1 minute over the course of
-    several months. Shadows are assumed to be relatively short
+    several months. Detection focuses on shadows with relatively short
     duration. The algorithm forms a 2D image of the GHI data by
     arranging time of day along the x-axis and day of year along the
     y-axis. Rapid change in GHI in the x-direction is used to identify

From 15c448aa1b1160b7567dc2c1e5fc72a54a563de0 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 30 Mar 2021 10:59:00 -0600
Subject: [PATCH 21/31] Rename `window` to `duration` in
 shading._detect_clouds()

Adds a better docstring for the returned series.
---
 pvanalytics/features/shading.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 517fe8ca..f65fe986 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -112,15 +112,22 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     )
 
 
-def _detect_clouds(ghi, clearsky_ghi, window='50T'):
-    """Use :py:func:`pvanalytics.clearsky.reno` to detect clouds"""
+def _detect_clouds(ghi, clearsky_ghi, duration='50T'):
+    """Use :py:func:`pvanalytics.clearsky.reno` to detect clouds.
+
+    Returns
+    -------
+    Series
+        Boolean series with true for cloudy periods that are at least as long
+        as `duration`.
+    """
     clouds = pvlib.clearsky.detect_clearsky(
         ghi,
         clearsky_ghi,
         ghi.index,
         window_length=10
     )
-    return clouds.rolling(window).sum() == 0
+    return clouds.rolling(duration).sum() == 0
 
 
 def _remove_pillars(wires):
@@ -282,7 +289,7 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     ----------
     ghi : Series
         Time series of GHI measurements. Data must be at 1-minute frequency
-        and cover at least several months.
+        and should cover at least 60 days.
     daytime : Series
         Boolean series with True for times when the sun is up.
     clearsky : Series

From f1d1bc431a6c23a436b10845b4bbcd20a12c5e24 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 6 Apr 2021 07:36:17 -0600
Subject: [PATCH 22/31] Bump minimum pandas version to 0.24.0

Following NEP 29, we can now require at least 0.24.0 which was
released 24 months ago.
---
 requirements-min.txt | 2 +-
 requirements.txt     | 2 +-
 setup.py             | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/requirements-min.txt b/requirements-min.txt
index e7d61919..46f4b391 100644
--- a/requirements-min.txt
+++ b/requirements-min.txt
@@ -1,5 +1,5 @@
 numpy~=1.15.0
-pandas~=0.23.0
+pandas~=0.24.0
 pvlib~=0.8.0
 scipy~=1.2.0
 statsmodels~=0.9.0
diff --git a/requirements.txt b/requirements.txt
index 71c90084..b195b452 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 numpy>=1.15.0
-pandas>=0.23.0,<1.1.0
+pandas>=0.24.0,<1.1.0
 pvlib>=0.8.0
 scipy>=1.2.0
 statsmodels>=0.9.0
diff --git a/setup.py b/setup.py
index fe070ba7..e6917256 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,7 @@
 
 INSTALL_REQUIRES = [
     'numpy >= 1.15.0',
-    'pandas >= 0.23.0',
+    'pandas >= 0.24.0',
     'pvlib >= 0.8.0',
     'scipy >= 1.2.0',
     'statsmodels >= 0.9.0',

From 10e6954f7bb5a939bd23b68dec535d3ffc17a89d Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Tue, 6 Apr 2021 08:14:13 -0600
Subject: [PATCH 23/31] Add parameter descriptions to shading._prepare_images()

---
 pvanalytics/features/shading.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index f65fe986..35acaf18 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -57,6 +57,17 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     Performs pre-processing steps on `ghi` and `clearsky` before
     returning images for use in the shadow detection algorithm.
 
+    Parameters
+    ----------
+    ghi : Series
+        Measured GHI.
+    clearsky : Series
+        Expected clearsky GHI.
+    daytime : Series
+        Boolean series with True for daytime and False for night.
+    interval : int
+        Time between data points in `ghi`. [minutes]
+
     Returns
     -------
     ghi_image : np.ndarray

From 556443f077e2519e6d4ddfd9ead926a4159ccdfe Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Thu, 15 Apr 2021 12:20:51 -0600
Subject: [PATCH 24/31] Specify output array for image operations in
 _filter_bars()

Without specifying an output array the image processing and numpy
functions will allocate a new array for every operation. This can
be a very expensive operation; fortunately we can avoid it by using
a temporary array for the morphological operations and specifying
an output array for the resulting image with horizontal bars
removed.
---
 pvanalytics/features/shading.py | 29 ++++++++++++++++++++++++-----
 1 file changed, 24 insertions(+), 5 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 35acaf18..76065490 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -241,7 +241,22 @@ def _tand(degrees):
 
 
 def _filter_bars(wires, out):
-    """Only keep shadows that are 'bar-shaped' and span multiple days."""
+    """Only keep shadows that are 'bar-shaped' and span multiple days.
+
+    Parameters
+    ----------
+    wires : ndarray
+        Binary image of candidate shadows.
+    out : ndarray
+        An ndarray that will hold the output image resulting from this
+        filtering operation. If None then a new array will be allocated.
+
+    Returns
+    -------
+    ndarray
+        The output array. If `out` is not none, will return `out`.
+    """
+    temp_image = np.ndarray(wires.shape)
     for angle in range(0, 90, 5):
         if angle < 75:
             dim = 20
@@ -258,13 +273,17 @@ def _filter_bars(wires, out):
             rr = mid + np.round((cc - mid) * _tand(90 - angle))
             se[cc, rr.astype('int')] = 1
             se = se.T
-        out = np.logical_or(
+        morphology.binary_opening(wires, se, out=temp_image)
+        np.logical_or(
             out,
-            morphology.binary_opening(wires, se)
+            temp_image,
+            out=out
         )
-        out = np.logical_or(
+        morphology.binary_opening(wires, np.flipud(se), out=temp_image)
+        np.logical_or(
             out,
-            morphology.binary_opening(wires, np.flipud(se))
+            temp_image,
+            out=out
         )
     return out
 

From 6485079d3707dddf737d0e4723d49750e8edf564 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Thu, 15 Apr 2021 14:26:21 -0600
Subject: [PATCH 25/31] Remove unused shading._width() function

---
 pvanalytics/features/shading.py | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 76065490..e1fb2146 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -191,22 +191,6 @@ def _restore_gaps(wires):
     )
 
 
-def _width(bbox):
-    """Return the width of a bounding box.
-
-    Parameters
-    ----------
-    bbox : tuple
-        4-tuple specifying ``(min_row, min_col, max_row, max_col)``.
-
-    Returns
-    -------
-    int
-        The width of the bounding box in pixels.
-    """
-    return bbox[3] - bbox[1] + 1
-
-
 def _height(bbox):
     """Return the height of a bounding box.
 

From 4b0d2e3794ba68f8ab34219a255d7dd9bd826f08 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Thu, 22 Apr 2021 11:40:08 -0600
Subject: [PATCH 26/31] Pass output array to _fill_gaps() and _remove_spikes()

Improves performance by reducing the number of new arrays that must
be allocated.
---
 pvanalytics/features/shading.py | 21 ++++++++++++++-------
 1 file changed, 14 insertions(+), 7 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index e1fb2146..e323e403 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -161,18 +161,25 @@ def _remove_pillars(wires):
     return np.logical_and(a1, a2)
 
 
-def _fill_gaps(wires):
+def _fill_gaps(wires, out):
     """Fill 1-minute gaps."""
     mask = np.array([[1, 0, 1]])
-    return np.logical_or(wires, ndimage.binary_hit_or_miss(wires, mask))
+    return np.logical_or(
+        wires, ndimage.binary_hit_or_miss(wires, mask), out=out)
 
 
-def _remove_spikes(wires):
+def _remove_spikes(wires, out):
     """Remove 1-day spikes."""
     mask = np.array([[0, 1, 0]]).T
+    temp_image = np.ndarray(wires.shape)
+    hit_miss = ndimage.binary_hit_or_miss(wires, mask)
     return np.logical_and(
         wires,
-        np.logical_not(ndimage.binary_hit_or_miss(wires, mask))
+        np.logical_not(
+            hit_miss,
+            out=temp_image
+        ),
+        out=out
     )
 
 
@@ -275,9 +282,9 @@ def _filter_bars(wires, out):
 def _clean_wires(wires):
     """Clean up clouds that are connected to wires."""
     out = _remove_pillars(wires)
-    out = _fill_gaps(out)
-    out = _remove_spikes(out)
-    out = _fill_gaps(out)
+    out = _fill_gaps(out, out)
+    out = _remove_spikes(out, out)
+    out = _fill_gaps(out, out)
     out = _restore_gaps(out)
     out = _filter_blobs(out, 20, connectivity=1)
     out = _filter_bars(wires, out)

From 469a4a108fe0ce5de7525e8c963b7aa080100798 Mon Sep 17 00:00:00 2001
From: Cliff Hansen <cwhanse@sandia.gov>
Date: Fri, 28 May 2021 10:37:35 -0600
Subject: [PATCH 27/31] debugging

---
 pvanalytics/features/shading.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index e323e403..1e2d8275 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -97,7 +97,8 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     scaled_clearsky = scaled_clearsky.reindex_like(scaled_ghi)
     daytime = daytime.reindex_like(scaled_ghi)
     # Detect clouds.
-    clouds = _detect_clouds(scaled_ghi, scaled_clearsky, '50T')
+    window_size = 50 // interval
+    clouds = _detect_clouds(scaled_ghi, scaled_clearsky, window_size)
     cloud_mask = _to_image(clouds.to_numpy(), image_width)
     # Interpolate across days (i.e. along columns) to remove clouds
     # replace clouds with nans
@@ -123,7 +124,7 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     )
 
 
-def _detect_clouds(ghi, clearsky_ghi, duration='50T'):
+def _detect_clouds(ghi, clearsky_ghi, window_size):
     """Use :py:func:`pvanalytics.clearsky.reno` to detect clouds.
 
     Returns
@@ -132,22 +133,24 @@ def _detect_clouds(ghi, clearsky_ghi, duration='50T'):
         Boolean series with true for cloudy periods that are at least as long
         as `duration`.
     """
-    clouds = pvlib.clearsky.detect_clearsky(
+    cleartimes = pvlib.clearsky.detect_clearsky(
         ghi,
         clearsky_ghi,
         ghi.index,
-        window_length=10
+        window_length=10,
+        max_iterations=1
     )
-    return clouds.rolling(duration).sum() == 0
+
+    return cleartimes.rolling(window_size, center=True).sum() == 0
 
 
 def _remove_pillars(wires):
     j = np.array([[0, 0, 0, 0, 0],
                   [0, 0, 1, 0, 0],
-                  [0, 0, 1, 0, 0]])
+                  [0, 0, 1, 0, 0]]).T
     k = np.array([[1, 0, 0, 0, 1],
                   [1, 0, 0, 0, 1],
-                  [1, 0, 0, 0, 1]])
+                  [1, 0, 0, 0, 1]]).T
     a1 = np.logical_and(
         wires,
         np.logical_not(ndimage.binary_hit_or_miss(wires, j, k))

From 1788e43cf2fd35161a2d9b2588e37c202be1dc50 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Fri, 28 May 2021 11:25:54 -0600
Subject: [PATCH 28/31] Tidy up doc string

Add image to return values
---
 pvanalytics/features/shading.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 1e2d8275..6bb61cf6 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -312,8 +312,8 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     Parameters
     ----------
     ghi : Series
-        Time series of GHI measurements. Data must be at 1-minute frequency
-        and should cover at least 60 days.
+        Time series of GHI measurements. Data must be in local time at
+        1-minute frequency and should cover at least 60 days.
     daytime : Series
         Boolean series with True for times when the sun is up.
     clearsky : Series
@@ -326,10 +326,13 @@ def fixed(ghi, daytime, clearsky, interval=None, min_gradient=2):
     -------
     Series
         Boolean series with true for times that are impacted by shadows.
+    ndarray
+        A boolean image (black and white) showing the shadows that were
+        detected.
 
     References
     ----------
-    References
+
     .. [1] Martin, C. E., Hansen, C. W., An Image Processing Algorithm to
        Identify Near-Field Shading in Irradiance Measurements, preprint 2016
     .. [2] Reno, M.J. and C.W. Hansen, "Identification of periods of clear sky

From ac18e807f7a736d93e6ec1151a4e58985a3daa88 Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Fri, 28 May 2021 11:26:59 -0600
Subject: [PATCH 29/31] Add shading to api.rst and whatsnew

---
 docs/api.rst            | 10 ++++++++++
 docs/whatsnew/0.1.1.rst |  3 +++
 2 files changed, 13 insertions(+)

diff --git a/docs/api.rst b/docs/api.rst
index 39433e9c..c608a53b 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -231,6 +231,16 @@ Functions that return a Boolean mask indicating day and night.
 
    features.daytime.power_or_irradiance
 
+Shading
+-------
+
+Functions for labeling shadows.
+
+.. autosummary::
+   :toctree: generated/
+
+   features.shading.fixed
+
 System
 ======
 
diff --git a/docs/whatsnew/0.1.1.rst b/docs/whatsnew/0.1.1.rst
index 632306f0..3be4681b 100644
--- a/docs/whatsnew/0.1.1.rst
+++ b/docs/whatsnew/0.1.1.rst
@@ -13,6 +13,8 @@ Enhancements
   :py:func:`pvanalytics.metrics.variability_index`. (:issue:`60`, :pull:`106`)
 * Internal refactor of `pvanalytics.metrics.performance_ratio_nrel` to support
   other performance ratio formulas. (:pull:`109`)
+* Detect shadows from fixed objects in GHI data using
+  :py:func:`pvanalytics.features.shading.fixed`. (:issue:`24`, :pull:`101`)
 
 Bug Fixes
 ~~~~~~~~~
@@ -22,3 +24,4 @@ Contributors
 
 * Kevin Anderson (:ghuser:`kanderso-nrel`)
 * Cliff Hansen (:ghuser:`cwhanse`)
+* Will Vining (:ghuser:`wfvining`)

From e97a3bf606860fb203f9658887de8b53e554ca6d Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Thu, 3 Jun 2021 13:14:21 -0600
Subject: [PATCH 30/31] Docstring edits

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 pvanalytics/features/shading.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pvanalytics/features/shading.py b/pvanalytics/features/shading.py
index 6bb61cf6..90c9a419 100644
--- a/pvanalytics/features/shading.py
+++ b/pvanalytics/features/shading.py
@@ -12,7 +12,7 @@ def _to_image(data, width):
 
     Parameters
     ----------
-    data : Series
+    data : array_like
         Values of the pixels
     width : int
         Width of the image in pixels
@@ -60,9 +60,9 @@ def _prepare_images(ghi, clearsky, daytime, interval):
     Parameters
     ----------
     ghi : Series
-        Measured GHI.
+        Measured GHI. [W/m^2]
     clearsky : Series
-        Expected clearsky GHI.
+        Expected clearsky GHI. [W/m^2]
     daytime : Series
         Boolean series with True for daytime and False for night.
     interval : int
@@ -131,7 +131,7 @@ def _detect_clouds(ghi, clearsky_ghi, window_size):
     -------
     Series
         Boolean series with true for cloudy periods that are at least as long
-        as `duration`.
+        as `window_size`.
     """
     cleartimes = pvlib.clearsky.detect_clearsky(
         ghi,

From 8b2b17d0eb4348ea29fd503ab6fad5db56c6c1be Mon Sep 17 00:00:00 2001
From: Will Vining <wfvining@gmail.com>
Date: Thu, 3 Jun 2021 13:23:50 -0600
Subject: [PATCH 31/31] Add features.shading to package overview in index.rst

---
 docs/index.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/index.rst b/docs/index.rst
index 898fe943..10261198 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -53,6 +53,7 @@ library status.
     functions in the :py:mod:`system` module in that we are identifying
     features of data rather than properties of the system that produced the
     data.
+  - :py:mod:`features.shading` functions for identifying shadows.
 
 - :py:mod:`system` identification of PV system characteristics from data
   (e.g. nameplate power, tilt, azimuth)