From 8da3850eb4a8879be0edae4af6e1df5132199cdc Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Fri, 1 Nov 2019 09:57:35 +0100
Subject: [PATCH 01/13] create initial files for scaling package

---
 pvlib/scaling.py           | 4 ++++
 pvlib/test/test_scaling.py | 0
 2 files changed, 4 insertions(+)
 create mode 100644 pvlib/scaling.py
 create mode 100644 pvlib/test/test_scaling.py

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
new file mode 100644
index 0000000000..54db8e40dc
--- /dev/null
+++ b/pvlib/scaling.py
@@ -0,0 +1,4 @@
+"""
+The ``scaling`` module contains functions for manipulating irradiance
+or other variables to account for temporal or spatial characteristics.
+"""
\ No newline at end of file
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
new file mode 100644
index 0000000000..e69de29bb2

From 0bde78adfbbea06a5081c9efe44143cba0205e98 Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Fri, 1 Nov 2019 15:47:04 +0100
Subject: [PATCH 02/13] Completed draft of WVM using only discrete distance
 algorithm.

---
 pvlib/scaling.py           | 311 ++++++++++++++++++++++++++++++++++++-
 pvlib/test/test_scaling.py | 112 +++++++++++++
 2 files changed, 422 insertions(+), 1 deletion(-)

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index 54db8e40dc..b5b32a0f33 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -1,4 +1,313 @@
 """
 The ``scaling`` module contains functions for manipulating irradiance
 or other variables to account for temporal or spatial characteristics.
-"""
\ No newline at end of file
+"""
+
+import numpy as np
+import pandas as pd
+from scipy.spatial.distance import pdist
+import scipy.optimize
+
+
+def wvm(cs_series, latitude, longitude, cloud_speed,
+        method, capacity=None, density=41, dt=None):
+    """
+    Compute spatial aggregation time series smoothing on clear sky index based
+    on the Wavelet Variability model of Lave et al [1-2]. Implementation is
+    basically a port of the Matlab version of the code [3].
+
+    Parameters
+    ----------
+    cs_series : numeric or pandas.Series
+        Clear Sky Index time series that will be smoothed.
+
+    latitude : numeric
+        Latitudes making up the plant in degrees. Type depends on method used.
+        'discrete' uses a list of lat/lon points representing discrete sites.
+        'square' uses a single valued lat/lon pair for the plant center.
+        'polygon' uses a list of lat/lon points representing the corners of
+        a polygon making up the site.
+
+    longitude : numeric
+        Longitudes making up the plant in degrees. Type depends on method used.
+        'discrete' uses a list of lat/lon points representing discrete sites.
+        'square' uses a single valued lat/lon pair for the plant center.
+        'polygon' uses a list of lat/lon points representing the corners of
+        a polygon making up the site.
+
+    cloud_speed : numeric
+        Speed of cloud movement in meters per second
+
+    method : string
+        The type of plant distribution to model.
+        Options are ``'discrete', 'square', 'polygon'``.
+
+    capacity : numeric
+        The plant capacity in MW. Must be specified for 'square' method,
+        ignored otherwise.
+
+    density : numeric, default 41
+        The density of installed PV in W/m^2. Must be specified for 'square'
+        method, ignored otherwise. Default value of 41 W/m^2 is 1MW per 6 acres.
+
+    dt : numeric
+        The time series time delta. By default, is inferred from the cs_series.
+        Must be specified for a time series that doesn't include an index.
+
+    Returns
+    -------
+    smoothed : pandas.Series or numeric, depending on inlet type
+        The clear sky index time series smoothed for the described plant.
+
+    wavelet: numeric
+        The individual wavelets for the time series before smoothing
+
+    tmscales: numeric
+        The timescales (in sec) associated with the wavelets
+
+    References
+    ----------
+    [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
+    Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
+    Energy, vol. 4, no. 2, pp. 501-509, 2013.
+
+    [2] M. Lave and J. Kleissl. Cloud speed impact on solar variability
+    scaling - Application to the wavelet variability model. Solar Energy,
+    vol. 91, pp. 11-21, 2013.
+
+    [3] Wavelet Variability Model - Matlab Code:
+    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    """
+
+    # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
+
+    dist = _compute_distances(longitude, latitude, method, capacity, density)
+    wavelet, tmscales = _compute_wavelet(cs_series, dt)
+
+    n_pairs = len(dist)
+    # Find eff length of position vector, 'dist' is full pairwise
+    fn = lambda x: np.abs((x ** 2 - x) / 2 - n_pairs)
+    n_dist = np.round(scipy.optimize.fmin(fn, np.sqrt(n_pairs), disp=False))
+
+    # Compute VR
+    A = cloud_speed / 2  # Resultant fit for A from [2]
+    vr = np.zeros(tmscales.shape)
+    for i, tmscale in enumerate(tmscales):
+        rho = np.exp(-1 / A * dist / tmscale)  # Eq 5 from [1]
+
+        # 2*rho is because rho_ij = rho_ji. +n_dist accounts for sum(rho_ii=1)
+        denominator = 2 * np.sum(rho) + n_dist
+        vr[i] = n_dist ** 2 / denominator  # Eq 6 of [1]
+
+    # Scale each wavelet by VR (Eq 7 in [1])
+    wavelet_smooth = np.zeros_like(wavelet)
+    for i in np.arange(len(tmscales)):
+        if i < len(tmscales) - 1:  # Treat the lowest freq differently
+            wavelet_smooth[i, :] = wavelet[i, :] / np.sqrt(vr[i])
+        else:
+            wavelet_smooth[i, :] = wavelet[i, :]
+
+    outsignal = np.sum(wavelet_smooth, 0)
+
+    try:  # See if there's an index already, if so, return as a pandas Series
+        smoothed = pd.Series(outsignal, index=cs_series.index)
+    except AttributeError:
+        smoothed = outsignal  # just output the numpy signal
+
+    return smoothed, wavelet, tmscales
+
+
+def _latlon_to_dist(latitude, longitude):
+    """
+    Convert latitude and longitude in degrees to a coordinate system measured
+    in meters from zero deg latitude, zero deg longitude.
+
+    Parameters
+    ----------
+    latitude : numeric
+        Latitude in degrees
+
+    longitude : numeric
+        Longitude in degrees
+
+    Returns
+    -------
+    ypos : numeric
+        the northward distance in meters
+
+    xpos: numeric
+        the eastward distance in meters.
+
+    References
+    ----------
+    [1] H. Moritz. Geodetic Reference System 1980, Journal of Geodesy, vol. 74,
+    no. 1, pp 128–133, 2000.
+
+    [2] Wavelet Variability Model - Matlab Code:
+    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    """
+
+    # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
+
+    r_earth = 6371008.7714  # mean radius of Earth, in meters
+    m_per_deg_lat = r_earth * np.pi / 180
+    m_per_deg_lon = r_earth * np.cos(np.pi/180*np.mean(latitude)) * np.pi/180
+
+    try:
+        xpos = m_per_deg_lon * longitude
+        ypos = m_per_deg_lat * latitude
+    except TypeError:  # When longitude and latitude are a list
+        xpos = m_per_deg_lon * np.array(longitude)
+        ypos = m_per_deg_lat * np.array(latitude)
+
+    return ypos, xpos
+
+
+def _compute_distances(longitude, latitude, method, capacity=None, density=41):
+    """
+    Compute points representing a plant and the pairwise distances between them
+
+    Parameters
+    ----------
+    latitude : numeric
+        Latitudes making up the plant in degrees. Type depends on method used.
+        'discrete' uses a list of lat/lon points representing discrete sites.
+        'square' uses a single valued lat/lon pair for the plant center.
+        'polygon' uses a list of lat/lon points representing the corners of
+        a polygon making up the site.
+
+    longitude : numeric
+        Longitudes making up the plant in degrees. Type depends on method used.
+        'discrete' uses a list of lat/lon points representing discrete sites.
+        'square' uses a single valued lat/lon pair for the plant center.
+        'polygon' uses a list of lat/lon points representing the corners of
+        a polygon making up the site.
+
+    cloud_speed : numeric
+        Speed of cloud movement in meters per second
+
+    method : string
+        The type of plant distribution to model.
+        Options are ``'discrete', 'square', 'polygon'``.
+
+    capacity : numeric
+        The plant capacity in MW. Must be specified for 'square' method,
+        ignored otherwise.
+
+    density : numeric, default 41
+        The density of installed PV in W/m^2. Must be specified for 'square'
+        method, ignored otherwise. Default value of 41 W/m^2 is 1MW per 6 acres.
+
+    Returns
+    -------
+    dist : numeric
+        The complete set of pairwise distances of all points representing the
+        plant being modelled.
+
+
+    References
+    ----------
+    [1] Wavelet Variability Model - Matlab Code:
+    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    """
+
+    # Convert latitude and longitude points to distances in meters
+    ypos, xpos = _latlon_to_dist(latitude, longitude)
+
+    if method == "discrete":
+        # Positions are individual plant centers.
+        # Treat as existing subscale points.
+        pos = np.array([xpos, ypos]).transpose()
+
+    elif method == "square":
+        raise NotImplementedError("To be implemented")
+
+    elif method == "polygon":
+        raise NotImplementedError("To be implemented")
+
+    else:
+        raise ValueError("Plant distance calculation method must be one of: "
+                         "discrete, square, polygon")
+
+    # Compute the full list of point-to-point distances
+    return pdist(pos, 'euclidean')
+
+
+def _compute_wavelet(cs_series, dt=None):
+    """
+    Compute the wavelet transform on the input clear_sky time series.
+
+    Parameters
+    ----------
+    cs_series : numeric or pandas.Series
+        Clear Sky Index time series that will be smoothed.
+
+    dt : numeric
+        The time series time delta. By default, is inferred from the cs_series.
+        Must be specified for a time series that doesn't include an index.
+
+    Returns
+    -------
+    wavelet: numeric
+        The individual wavelets for the time series before smoothing
+
+    tmscales: numeric
+        The timescales (in sec) associated with the wavelets
+
+    References
+    ----------
+    [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
+    Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
+    Energy, vol. 4, no. 2, pp. 501-509, 2013.
+
+    [3] Wavelet Variability Model - Matlab Code:
+    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    """
+
+    # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
+
+    try:  # Assume it's a pandas type
+        vals = cs_series.values.flatten()
+        try:  # Assume it's a time series type index
+            dt = (cs_series.index[1] - cs_series.index[0]).seconds
+        except AttributeError:  # It must just be a numeric index
+            dt = (cs_series.index[1] - cs_series.index[0])
+    except AttributeError:  # Assume it's a numpy type
+        vals = cs_series.flatten()
+        if dt is None:
+            raise ValueError("dt must be specified for numpy type inputs.")
+
+    # Pad the series on both ends in time and place in a dataframe
+    cs_long = np.pad(vals, (len(vals), len(vals)), 'symmetric')
+    cs_long = pd.DataFrame(cs_long)
+
+    # Compute wavelet time scales
+    mindt = np.ceil(np.log(dt)/np.log(2))  # Minimum wavelet dt
+    maxdt = int(12 - mindt)  # maximum wavelet dt
+
+    tmscales = np.zeros(maxdt)
+    csi_mean = np.zeros([maxdt, len(cs_long)])
+    # Loop for all time scales we will consider
+    for i in np.arange(0, maxdt):
+        j = i+1
+        tmscales[i] = 2**j * dt  # Wavelet integration time scale
+        intvllen = 2**j  # Wavelet integration time series interval
+        # Rolling average, retains only lower frequencies than interval
+        df = cs_long.rolling(window=intvllen, center=True, min_periods=1).mean()
+        # Fill nan's in both directions
+        df = df.fillna(method='bfill').fillna(method='ffill')
+        # Pop values back out of the dataframe and store
+        csi_mean[i, :] = df.values.flatten()
+
+    # Calculate the wavelets by isolating the rolling mean frequency ranges
+    wavelet_long = np.zeros(csi_mean.shape)
+    for i in np.arange(0, maxdt-1):
+        wavelet_long[i, :] = csi_mean[i, :] - csi_mean[i+1, :]
+    wavelet_long[maxdt-1, :] = csi_mean[maxdt-1, :]  # Lowest frequency
+
+    # Clip off the padding and just return the original time window
+    wavelet = np.zeros([maxdt, len(vals)])
+    for i in np.arange(0, maxdt):
+        wavelet[i, :] = wavelet_long[i, len(vals)+1: 2*len(vals)+1]
+
+    return wavelet, tmscales
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index e69de29bb2..227bba3b01 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -0,0 +1,112 @@
+import numpy as np
+import pandas as pd
+
+import pytest
+from numpy.testing import assert_almost_equal
+
+from pvlib import scaling
+
+# All expected_xxxxxx variable results computed in Matlab code
+
+# Sample positions
+lat = np.array((9.99, 10, 10.01))
+lon = np.array((4.99, 5, 5.01))
+# Sample cloud speed
+cloud_speed = 5
+# Generate a sample clear_sky_index and time vector.
+clear_sky_index = np.ones(10000)
+clear_sky_index[5000:5005] = np.array([1, 1, 1.1, 0.9, 1])
+time = np.arange(0, len(clear_sky_index))
+# Sample dt
+dt = 1
+
+# Expected distance and positions for sample lat/lon given above
+expect_dist = np.array((1560.6, 3121.3, 1560.6))
+expect_xpos = np.array([554863.4, 555975.4, 557087.3])
+expect_ypos = np.array([1106611.8, 1107719.5, 1108827.2])
+
+# Expected timescales for dt = 1
+expect_tmscale = [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+
+# Expected wavelet for indices 5000:5004 using clear_sky_index above
+expect_wavelet = np.array([[-0.025, 0.05, 0. ,-0.05, 0.025],
+                           [ 0.025, 0.  , 0. , 0.  ,-0.025],
+                           [ 0.,    0.  , 0. , 0.  , 0.]])
+
+# Expected smoothed clear sky index for indices 5000:5004 using inputs above
+expect_cs_smooth = np.array([1., 1.0289, 1., 0.9711, 1.])
+
+
+def test_latlon_to_dist_zero():
+    lat = 0
+    lon = 0
+    xpos_e = 0
+    ypos_e = 0
+    xpos, ypos = scaling._latlon_to_dist(lon, lat)
+    assert_almost_equal(xpos, xpos_e, decimal=1)
+    assert_almost_equal(ypos, ypos_e, decimal=1)
+
+
+def test_latlon_to_dist_single():
+    # Must test against central value, because latlon_to_dist uses the mean
+    xpos, ypos = scaling._latlon_to_dist(lon[1], lat[1])
+    assert_almost_equal(xpos, expect_xpos[1], decimal=1)
+    assert_almost_equal(ypos, expect_ypos[1], decimal=1)
+
+
+def test_latlon_to_dist_array():
+    xpos, ypos = scaling._latlon_to_dist(lon, lat)
+    assert_almost_equal(xpos, expect_xpos, decimal=1)
+    assert_almost_equal(ypos, expect_ypos, decimal=1)
+
+
+def test_compute_distances_invalid():
+    with pytest.raises(ValueError):
+        scaling._compute_distances(0, 0, method='invalid')
+
+
+def test_compute_distances_discrete_zero():
+    lat = np.array((0, 0))
+    lon = np.array((0, 0))
+    assert_almost_equal(scaling._compute_distances(lon, lat, 'discrete'), 0)
+
+
+def test_compute_distances_discrete_array():
+
+    dist = scaling._compute_distances(lon, lat, 'discrete')
+    assert_almost_equal(dist, expect_dist, decimal=1)
+
+
+def test_compute_wavelet_series():
+    csi_series = pd.Series(clear_sky_index, index=time)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(tmscale, expect_tmscale)
+    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+
+
+def test_compute_wavelet_array():
+    wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
+    assert_almost_equal(tmscale, expect_tmscale)
+    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+
+
+def test_compute_wavelet_array_invalid():
+    with pytest.raises(ValueError):
+        scaling._compute_wavelet(clear_sky_index)
+
+
+def test_wvm_series():
+    csi_series = pd.Series(clear_sky_index, index=time)
+    cs_sm, _, _ = scaling.wvm(csi_series, lat, lon, cloud_speed, "discrete")
+    assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
+
+
+def test_compute_wvm_array():
+    cs_sm, _, _ = scaling.wvm(clear_sky_index, lat, lon, cloud_speed,
+                              "discrete", dt=dt)
+    assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
+
+
+def test_compute_wvm_array_invalid():
+    with pytest.raises(ValueError):
+        scaling.wvm(clear_sky_index, lat, lon, cloud_speed, "discrete")

From 67fc83a56a23777cba220ae509c01cdf9d96d4bb Mon Sep 17 00:00:00 2001
From: Joe Ranalli <jar339@psu.edu>
Date: Fri, 1 Nov 2019 21:22:22 +0100
Subject: [PATCH 03/13] fixed lint errors, moved scipy imports and marked tests
 @requires_scipy

---
 pvlib/scaling.py           | 27 +++++++++++++++++----------
 pvlib/test/test_scaling.py | 17 ++++++++++++-----
 2 files changed, 29 insertions(+), 15 deletions(-)

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index b5b32a0f33..5ab5015f2e 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -5,8 +5,6 @@
 
 import numpy as np
 import pandas as pd
-from scipy.spatial.distance import pdist
-import scipy.optimize
 
 
 def wvm(cs_series, latitude, longitude, cloud_speed,
@@ -48,7 +46,7 @@ def wvm(cs_series, latitude, longitude, cloud_speed,
 
     density : numeric, default 41
         The density of installed PV in W/m^2. Must be specified for 'square'
-        method, ignored otherwise. Default value of 41 W/m^2 is 1MW per 6 acres.
+        method, ignored otherwise. Default of 41 W/m^2 is 1MW per 6 acres.
 
     dt : numeric
         The time series time delta. By default, is inferred from the cs_series.
@@ -81,12 +79,19 @@ def wvm(cs_series, latitude, longitude, cloud_speed,
 
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
 
+    try:
+        import scipy.optimize
+    except ImportError:
+        raise ImportError("The WVM function requires scipy.")
+
     dist = _compute_distances(longitude, latitude, method, capacity, density)
     wavelet, tmscales = _compute_wavelet(cs_series, dt)
 
     n_pairs = len(dist)
     # Find eff length of position vector, 'dist' is full pairwise
-    fn = lambda x: np.abs((x ** 2 - x) / 2 - n_pairs)
+
+    def fn(x):
+        return np.abs((x ** 2 - x) / 2 - n_pairs)
     n_dist = np.round(scipy.optimize.fmin(fn, np.sqrt(n_pairs), disp=False))
 
     # Compute VR
@@ -183,9 +188,6 @@ def _compute_distances(longitude, latitude, method, capacity=None, density=41):
         'polygon' uses a list of lat/lon points representing the corners of
         a polygon making up the site.
 
-    cloud_speed : numeric
-        Speed of cloud movement in meters per second
-
     method : string
         The type of plant distribution to model.
         Options are ``'discrete', 'square', 'polygon'``.
@@ -196,7 +198,7 @@ def _compute_distances(longitude, latitude, method, capacity=None, density=41):
 
     density : numeric, default 41
         The density of installed PV in W/m^2. Must be specified for 'square'
-        method, ignored otherwise. Default value of 41 W/m^2 is 1MW per 6 acres.
+        method, ignored otherwise. Default of 41 W/m^2 is 1MW per 6 acres.
 
     Returns
     -------
@@ -211,6 +213,11 @@ def _compute_distances(longitude, latitude, method, capacity=None, density=41):
     https://pvpmc.sandia.gov/applications/wavelet-variability-model/
     """
 
+    try:
+        from scipy.spatial.distance import pdist
+    except ImportError:
+        raise ImportError("The WVM function requires scipy.")
+
     # Convert latitude and longitude points to distances in meters
     ypos, xpos = _latlon_to_dist(latitude, longitude)
 
@@ -291,9 +298,9 @@ def _compute_wavelet(cs_series, dt=None):
     for i in np.arange(0, maxdt):
         j = i+1
         tmscales[i] = 2**j * dt  # Wavelet integration time scale
-        intvllen = 2**j  # Wavelet integration time series interval
+        intvlen = 2**j  # Wavelet integration time series interval
         # Rolling average, retains only lower frequencies than interval
-        df = cs_long.rolling(window=intvllen, center=True, min_periods=1).mean()
+        df = cs_long.rolling(window=intvlen, center=True, min_periods=1).mean()
         # Fill nan's in both directions
         df = df.fillna(method='bfill').fillna(method='ffill')
         # Pop values back out of the dataframe and store
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index 227bba3b01..bc96a86dbc 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -5,6 +5,7 @@
 from numpy.testing import assert_almost_equal
 
 from pvlib import scaling
+from conftest import requires_scipy
 
 # All expected_xxxxxx variable results computed in Matlab code
 
@@ -29,9 +30,9 @@
 expect_tmscale = [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
 
 # Expected wavelet for indices 5000:5004 using clear_sky_index above
-expect_wavelet = np.array([[-0.025, 0.05, 0. ,-0.05, 0.025],
-                           [ 0.025, 0.  , 0. , 0.  ,-0.025],
-                           [ 0.,    0.  , 0. , 0.  , 0.]])
+expect_wavelet = np.array([[-0.025, 0.05, 0., -0.05, 0.025],
+                           [0.025, 0., 0., 0., -0.025],
+                           [0., 0., 0., 0., 0.]])
 
 # Expected smoothed clear sky index for indices 5000:5004 using inputs above
 expect_cs_smooth = np.array([1., 1.0289, 1., 0.9711, 1.])
@@ -60,17 +61,20 @@ def test_latlon_to_dist_array():
     assert_almost_equal(ypos, expect_ypos, decimal=1)
 
 
+@requires_scipy
 def test_compute_distances_invalid():
     with pytest.raises(ValueError):
         scaling._compute_distances(0, 0, method='invalid')
 
 
+@requires_scipy
 def test_compute_distances_discrete_zero():
     lat = np.array((0, 0))
     lon = np.array((0, 0))
     assert_almost_equal(scaling._compute_distances(lon, lat, 'discrete'), 0)
 
 
+@requires_scipy
 def test_compute_distances_discrete_array():
 
     dist = scaling._compute_distances(lon, lat, 'discrete')
@@ -95,18 +99,21 @@ def test_compute_wavelet_array_invalid():
         scaling._compute_wavelet(clear_sky_index)
 
 
+@requires_scipy
 def test_wvm_series():
     csi_series = pd.Series(clear_sky_index, index=time)
     cs_sm, _, _ = scaling.wvm(csi_series, lat, lon, cloud_speed, "discrete")
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
-def test_compute_wvm_array():
+@requires_scipy
+def test_wvm_array():
     cs_sm, _, _ = scaling.wvm(clear_sky_index, lat, lon, cloud_speed,
                               "discrete", dt=dt)
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
-def test_compute_wvm_array_invalid():
+@requires_scipy
+def test_wvm_array_invalid():
     with pytest.raises(ValueError):
         scaling.wvm(clear_sky_index, lat, lon, cloud_speed, "discrete")

From 78b8c1ab5e0522e83185557868f3f7c5b9c0c0cc Mon Sep 17 00:00:00 2001
From: Joe Ranalli <jar339@psu.edu>
Date: Fri, 1 Nov 2019 22:15:20 +0100
Subject: [PATCH 04/13] added tests to complete coverage as far as I know how

---
 pvlib/test/test_scaling.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index bc96a86dbc..16afd4e6e3 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -61,12 +61,30 @@ def test_latlon_to_dist_array():
     assert_almost_equal(ypos, expect_ypos, decimal=1)
 
 
+def test_latlon_to_dist_list():
+    xpos, ypos = scaling._latlon_to_dist(lon.tolist(), lat.tolist())
+    assert_almost_equal(xpos, expect_xpos, decimal=1)
+    assert_almost_equal(ypos, expect_ypos, decimal=1)
+
+
 @requires_scipy
 def test_compute_distances_invalid():
     with pytest.raises(ValueError):
         scaling._compute_distances(0, 0, method='invalid')
 
 
+@requires_scipy
+def test_compute_distances_square():
+    with pytest.raises(NotImplementedError):
+        scaling._compute_distances(0, 0, method='square')
+
+
+@requires_scipy
+def test_compute_distances_polygon():
+    with pytest.raises(NotImplementedError):
+        scaling._compute_distances(0, 0, method='polygon')
+
+
 @requires_scipy
 def test_compute_distances_discrete_zero():
     lat = np.array((0, 0))

From 69bd4173e261d7d66f434ab9b6aac15632e33fdf Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Mon, 4 Nov 2019 08:45:06 +0100
Subject: [PATCH 05/13] added lines to the sphinx documentation

---
 docs/sphinx/source/api.rst             | 12 ++++++++++++
 docs/sphinx/source/whatsnew/v0.7.0.rst |  2 ++
 2 files changed, 14 insertions(+)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 26786d4f88..e2e019af4f 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -565,3 +565,15 @@ Methods for calculating back surface irradiance
    :toctree: generated/
 
    bifacial.pvfactors_timeseries
+
+
+Scaling
+========
+
+Methods for manipulating irradiance for temporal or spatial considerations
+--------------------------------------------------------------------------
+
+.. autosummary::
+   :toctree: generated/
+
+   scaling.wvm
\ No newline at end of file
diff --git a/docs/sphinx/source/whatsnew/v0.7.0.rst b/docs/sphinx/source/whatsnew/v0.7.0.rst
index fdc419ef6e..7668f11fcc 100644
--- a/docs/sphinx/source/whatsnew/v0.7.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.0.rst
@@ -114,6 +114,8 @@ Enhancements
 * Add :py:func:`~pvlib.ivtools.fit_sdm_desoto`, a method to fit the De Soto single
   diode model to the typical specifications given in manufacturers datasheets.
 * Add `timeout` to :py:func:`pvlib.iotools.get_psm3`.
+* Add :py:func:`~pvlib.scaling.wvm`, a port of the wavelet variability model for
+  computing reductions in variability due to a spatially distributed plant.
 
 Bug fixes
 ~~~~~~~~~

From ce33db97aa9c3477ac6ef3902d78cbf10e7fe5a0 Mon Sep 17 00:00:00 2001
From: Joe Ranalli <jar339@psu.edu>
Date: Wed, 13 Nov 2019 10:11:18 +0100
Subject: [PATCH 06/13] Changes following review. WVM accepts list of xy tuples
 as inputs. Cleanup of docs.

---
 pvlib/scaling.py           | 201 +++++++++++--------------------------
 pvlib/test/test_scaling.py |  99 +++++++-----------
 2 files changed, 97 insertions(+), 203 deletions(-)

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index 5ab5015f2e..23e2dbacbc 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -7,8 +7,7 @@
 import pandas as pd
 
 
-def wvm(cs_series, latitude, longitude, cloud_speed,
-        method, capacity=None, density=41, dt=None):
+def wvm(clearsky_index, positions, cloud_speed, dt=None):
     """
     Compute spatial aggregation time series smoothing on clear sky index based
     on the Wavelet Variability model of Lave et al [1-2]. Implementation is
@@ -16,52 +15,33 @@ def wvm(cs_series, latitude, longitude, cloud_speed,
 
     Parameters
     ----------
-    cs_series : numeric or pandas.Series
+    clearsky_index : numeric or pandas.Series
         Clear Sky Index time series that will be smoothed.
 
-    latitude : numeric
-        Latitudes making up the plant in degrees. Type depends on method used.
-        'discrete' uses a list of lat/lon points representing discrete sites.
-        'square' uses a single valued lat/lon pair for the plant center.
-        'polygon' uses a list of lat/lon points representing the corners of
-        a polygon making up the site.
-
-    longitude : numeric
-        Longitudes making up the plant in degrees. Type depends on method used.
-        'discrete' uses a list of lat/lon points representing discrete sites.
-        'square' uses a single valued lat/lon pair for the plant center.
-        'polygon' uses a list of lat/lon points representing the corners of
-        a polygon making up the site.
+    positions : numeric
+        Array of coordinate distances as (x,y) pairs representing the
+        easting, northing of the site positions in meters [m]. Distributed
+        plants could be simulated by gridded points throughout the plant
+        footprint.
 
     cloud_speed : numeric
-        Speed of cloud movement in meters per second
-
-    method : string
-        The type of plant distribution to model.
-        Options are ``'discrete', 'square', 'polygon'``.
-
-    capacity : numeric
-        The plant capacity in MW. Must be specified for 'square' method,
-        ignored otherwise.
+        Speed of cloud movement in meters per second [m/s].
 
-    density : numeric, default 41
-        The density of installed PV in W/m^2. Must be specified for 'square'
-        method, ignored otherwise. Default of 41 W/m^2 is 1MW per 6 acres.
-
-    dt : numeric
-        The time series time delta. By default, is inferred from the cs_series.
-        Must be specified for a time series that doesn't include an index.
+    dt : float, default None
+        The time series time delta. By default, is inferred from the
+        clearsky_index. Must be specified for a time series that doesn't
+        include an index. Units of seconds [s].
 
     Returns
     -------
-    smoothed : pandas.Series or numeric, depending on inlet type
-        The clear sky index time series smoothed for the described plant.
+    smoothed : numeric or pandas.Series
+        The Clear Sky Index time series smoothed for the described plant.
 
     wavelet: numeric
-        The individual wavelets for the time series before smoothing
+        The individual wavelets for the time series before smoothing.
 
     tmscales: numeric
-        The timescales (in sec) associated with the wavelets
+        The timescales associated with the wavelets in seconds [s].
 
     References
     ----------
@@ -81,14 +61,16 @@ def wvm(cs_series, latitude, longitude, cloud_speed,
 
     try:
         import scipy.optimize
+        from scipy.spatial.distance import pdist
     except ImportError:
         raise ImportError("The WVM function requires scipy.")
 
-    dist = _compute_distances(longitude, latitude, method, capacity, density)
-    wavelet, tmscales = _compute_wavelet(cs_series, dt)
+    pos = np.array(positions)
+    dist = pdist(pos, 'euclidean')
+    wavelet, tmscales = _compute_wavelet(clearsky_index, dt)
 
+    # Find effective length of position vector, 'dist' is full pairwise
     n_pairs = len(dist)
-    # Find eff length of position vector, 'dist' is full pairwise
 
     def fn(x):
         return np.abs((x ** 2 - x) / 2 - n_pairs)
@@ -115,40 +97,46 @@ def fn(x):
     outsignal = np.sum(wavelet_smooth, 0)
 
     try:  # See if there's an index already, if so, return as a pandas Series
-        smoothed = pd.Series(outsignal, index=cs_series.index)
+        smoothed = pd.Series(outsignal, index=clearsky_index.index)
     except AttributeError:
         smoothed = outsignal  # just output the numpy signal
 
     return smoothed, wavelet, tmscales
 
 
-def _latlon_to_dist(latitude, longitude):
+def latlon_to_xy(coordinates):
     """
     Convert latitude and longitude in degrees to a coordinate system measured
     in meters from zero deg latitude, zero deg longitude.
 
+    This is a convenience method to support inputs to wvm. Note that the
+    methodology used is only suitable for short distances. For conversions of
+    longer distances, users should consider use of Universal Transverse
+    Mercator (UTM) or other suitable cartographic projection. Consider
+    packages built for cartographic projection such as pyproj (e.g.
+    pyproj.transform()) [2].
+
     Parameters
     ----------
-    latitude : numeric
-        Latitude in degrees
 
-    longitude : numeric
-        Longitude in degrees
+    coordinates : numeric
+        Array or list of (latitude, longitude) coordinate pairs. Use decimal
+        degrees notation.
 
     Returns
     -------
-    ypos : numeric
-        the northward distance in meters
-
-    xpos: numeric
-        the eastward distance in meters.
+    xypos : numeric
+        Array of coordinate distances as (x,y) pairs representing the
+        easting, northing of the position in meters [m].
 
     References
     ----------
     [1] H. Moritz. Geodetic Reference System 1980, Journal of Geodesy, vol. 74,
     no. 1, pp 128–133, 2000.
 
-    [2] Wavelet Variability Model - Matlab Code:
+    [2] https://pypi.org/project/pyproj/
+
+    [3] Wavelet Variability Model - Matlab Code:
     https://pvpmc.sandia.gov/applications/wavelet-variability-model/
     """
 
@@ -156,110 +144,43 @@ def _latlon_to_dist(latitude, longitude):
 
     r_earth = 6371008.7714  # mean radius of Earth, in meters
     m_per_deg_lat = r_earth * np.pi / 180
-    m_per_deg_lon = r_earth * np.cos(np.pi/180*np.mean(latitude)) * np.pi/180
-
     try:
-        xpos = m_per_deg_lon * longitude
-        ypos = m_per_deg_lat * latitude
-    except TypeError:  # When longitude and latitude are a list
-        xpos = m_per_deg_lon * np.array(longitude)
-        ypos = m_per_deg_lat * np.array(latitude)
-
-    return ypos, xpos
-
-
-def _compute_distances(longitude, latitude, method, capacity=None, density=41):
-    """
-    Compute points representing a plant and the pairwise distances between them
-
-    Parameters
-    ----------
-    latitude : numeric
-        Latitudes making up the plant in degrees. Type depends on method used.
-        'discrete' uses a list of lat/lon points representing discrete sites.
-        'square' uses a single valued lat/lon pair for the plant center.
-        'polygon' uses a list of lat/lon points representing the corners of
-        a polygon making up the site.
-
-    longitude : numeric
-        Longitudes making up the plant in degrees. Type depends on method used.
-        'discrete' uses a list of lat/lon points representing discrete sites.
-        'square' uses a single valued lat/lon pair for the plant center.
-        'polygon' uses a list of lat/lon points representing the corners of
-        a polygon making up the site.
-
-    method : string
-        The type of plant distribution to model.
-        Options are ``'discrete', 'square', 'polygon'``.
-
-    capacity : numeric
-        The plant capacity in MW. Must be specified for 'square' method,
-        ignored otherwise.
-
-    density : numeric, default 41
-        The density of installed PV in W/m^2. Must be specified for 'square'
-        method, ignored otherwise. Default of 41 W/m^2 is 1MW per 6 acres.
-
-    Returns
-    -------
-    dist : numeric
-        The complete set of pairwise distances of all points representing the
-        plant being modelled.
+        meanlat = np.mean([lat for (lat, lon) in coordinates])  # Mean latitude
+    except TypeError:  # Assume it's a single value?
+        meanlat = coordinates[0]
+    m_per_deg_lon = r_earth * np.cos(np.pi/180 * meanlat) * np.pi/180
 
+    # Conversion
+    pos = coordinates * np.array(m_per_deg_lat, m_per_deg_lon)
 
-    References
-    ----------
-    [1] Wavelet Variability Model - Matlab Code:
-    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
-    """
-
+    # reshape as (x,y) pairs to return
     try:
-        from scipy.spatial.distance import pdist
-    except ImportError:
-        raise ImportError("The WVM function requires scipy.")
-
-    # Convert latitude and longitude points to distances in meters
-    ypos, xpos = _latlon_to_dist(latitude, longitude)
-
-    if method == "discrete":
-        # Positions are individual plant centers.
-        # Treat as existing subscale points.
-        pos = np.array([xpos, ypos]).transpose()
-
-    elif method == "square":
-        raise NotImplementedError("To be implemented")
-
-    elif method == "polygon":
-        raise NotImplementedError("To be implemented")
-
-    else:
-        raise ValueError("Plant distance calculation method must be one of: "
-                         "discrete, square, polygon")
-
-    # Compute the full list of point-to-point distances
-    return pdist(pos, 'euclidean')
+        return np.column_stack([pos[:, 1], pos[:, 0]])
+    except IndexError:  # Assume it's a single value, which has a 1D shape
+        return np.array((pos[1], pos[0]))
 
 
-def _compute_wavelet(cs_series, dt=None):
+def _compute_wavelet(clearsky_index, dt=None):
     """
     Compute the wavelet transform on the input clear_sky time series.
 
     Parameters
     ----------
-    cs_series : numeric or pandas.Series
+    clearsky_index : numeric or pandas.Series
         Clear Sky Index time series that will be smoothed.
 
-    dt : numeric
-        The time series time delta. By default, is inferred from the cs_series.
-        Must be specified for a time series that doesn't include an index.
+    dt : float, default None
+        The time series time delta. By default, is inferred from the
+        clearsky_index. Must be specified for a time series that doesn't
+        include an index. Units of seconds [s].
 
     Returns
     -------
     wavelet: numeric
-        The individual wavelets for the time series before smoothing
+        The individual wavelets for the time series
 
     tmscales: numeric
-        The timescales (in sec) associated with the wavelets
+        The timescales associated with the wavelets in seconds [s]
 
     References
     ----------
@@ -274,13 +195,13 @@ def _compute_wavelet(cs_series, dt=None):
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
 
     try:  # Assume it's a pandas type
-        vals = cs_series.values.flatten()
+        vals = clearsky_index.values.flatten()
         try:  # Assume it's a time series type index
-            dt = (cs_series.index[1] - cs_series.index[0]).seconds
+            dt = (clearsky_index.index[1] - clearsky_index.index[0]).seconds
         except AttributeError:  # It must just be a numeric index
-            dt = (cs_series.index[1] - cs_series.index[0])
+            dt = (clearsky_index.index[1] - clearsky_index.index[0])
     except AttributeError:  # Assume it's a numpy type
-        vals = cs_series.flatten()
+        vals = clearsky_index.flatten()
         if dt is None:
             raise ValueError("dt must be specified for numpy type inputs.")
 
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index 16afd4e6e3..3f856da3fe 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -7,96 +7,63 @@
 from pvlib import scaling
 from conftest import requires_scipy
 
-# All expected_xxxxxx variable results computed in Matlab code
 
 # Sample positions
 lat = np.array((9.99, 10, 10.01))
 lon = np.array((4.99, 5, 5.01))
+coordinates = np.array([(lati, loni) for (lati, loni) in zip(lat, lon)])
+
 # Sample cloud speed
 cloud_speed = 5
 # Generate a sample clear_sky_index and time vector.
 clear_sky_index = np.ones(10000)
 clear_sky_index[5000:5005] = np.array([1, 1, 1.1, 0.9, 1])
 time = np.arange(0, len(clear_sky_index))
+
 # Sample dt
 dt = 1
 
-# Expected distance and positions for sample lat/lon given above
-expect_dist = np.array((1560.6, 3121.3, 1560.6))
+# Expected positions for sample lat/lon given above (calculated manually)
 expect_xpos = np.array([554863.4, 555975.4, 557087.3])
-expect_ypos = np.array([1106611.8, 1107719.5, 1108827.2])
+expect_ypos = np.array([1110838.8, 1111950.8, 1113062.7])
+
+# Sample positions based on the previous lat/lon
+positions = np.array([pt for pt in zip(expect_xpos, expect_ypos)])
 
 # Expected timescales for dt = 1
 expect_tmscale = [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
 
-# Expected wavelet for indices 5000:5004 using clear_sky_index above
+# Expected wavelet for indices 5000:5004 using clear_sky_index above (Matlab)
 expect_wavelet = np.array([[-0.025, 0.05, 0., -0.05, 0.025],
                            [0.025, 0., 0., 0., -0.025],
                            [0., 0., 0., 0., 0.]])
 
-# Expected smoothed clear sky index for indices 5000:5004 using inputs above
+# Expected smoothed clear sky index for indices 5000:5004 (Matlab)
 expect_cs_smooth = np.array([1., 1.0289, 1., 0.9711, 1.])
 
 
-def test_latlon_to_dist_zero():
-    lat = 0
-    lon = 0
-    xpos_e = 0
-    ypos_e = 0
-    xpos, ypos = scaling._latlon_to_dist(lon, lat)
-    assert_almost_equal(xpos, xpos_e, decimal=1)
-    assert_almost_equal(ypos, ypos_e, decimal=1)
+def test_latlon_to_xy_zero():
+    coord = [0, 0]
+    pos_e = [0, 0]
+    pos = scaling.latlon_to_xy(coord)
+    assert_almost_equal(pos, pos_e, decimal=1)
 
 
-def test_latlon_to_dist_single():
+def test_latlon_to_xy_single():
     # Must test against central value, because latlon_to_dist uses the mean
-    xpos, ypos = scaling._latlon_to_dist(lon[1], lat[1])
-    assert_almost_equal(xpos, expect_xpos[1], decimal=1)
-    assert_almost_equal(ypos, expect_ypos[1], decimal=1)
-
+    coord = (lat[1], lon[1])
+    pos = scaling.latlon_to_xy(coord)
+    assert_almost_equal(pos, (expect_xpos[1], expect_ypos[1]), decimal=1)
 
-def test_latlon_to_dist_array():
-    xpos, ypos = scaling._latlon_to_dist(lon, lat)
-    assert_almost_equal(xpos, expect_xpos, decimal=1)
-    assert_almost_equal(ypos, expect_ypos, decimal=1)
 
+def test_latlon_to_xy_array():
+    pos = scaling.latlon_to_xy(coordinates)
+    assert_almost_equal(pos, positions, decimal=1)
 
-def test_latlon_to_dist_list():
-    xpos, ypos = scaling._latlon_to_dist(lon.tolist(), lat.tolist())
-    assert_almost_equal(xpos, expect_xpos, decimal=1)
-    assert_almost_equal(ypos, expect_ypos, decimal=1)
 
-
-@requires_scipy
-def test_compute_distances_invalid():
-    with pytest.raises(ValueError):
-        scaling._compute_distances(0, 0, method='invalid')
-
-
-@requires_scipy
-def test_compute_distances_square():
-    with pytest.raises(NotImplementedError):
-        scaling._compute_distances(0, 0, method='square')
-
-
-@requires_scipy
-def test_compute_distances_polygon():
-    with pytest.raises(NotImplementedError):
-        scaling._compute_distances(0, 0, method='polygon')
-
-
-@requires_scipy
-def test_compute_distances_discrete_zero():
-    lat = np.array((0, 0))
-    lon = np.array((0, 0))
-    assert_almost_equal(scaling._compute_distances(lon, lat, 'discrete'), 0)
-
-
-@requires_scipy
-def test_compute_distances_discrete_array():
-
-    dist = scaling._compute_distances(lon, lat, 'discrete')
-    assert_almost_equal(dist, expect_dist, decimal=1)
+def test_latlon_to_xy_list():
+    pos = scaling.latlon_to_xy(coordinates.tolist())
+    assert_almost_equal(pos, positions, decimal=1)
 
 
 def test_compute_wavelet_series():
@@ -120,18 +87,24 @@ def test_compute_wavelet_array_invalid():
 @requires_scipy
 def test_wvm_series():
     csi_series = pd.Series(clear_sky_index, index=time)
-    cs_sm, _, _ = scaling.wvm(csi_series, lat, lon, cloud_speed, "discrete")
+    cs_sm, _, _ = scaling.wvm(csi_series, positions, cloud_speed)
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
 @requires_scipy
 def test_wvm_array():
-    cs_sm, _, _ = scaling.wvm(clear_sky_index, lat, lon, cloud_speed,
-                              "discrete", dt=dt)
+    cs_sm, _, _ = scaling.wvm(clear_sky_index, positions, cloud_speed, dt=dt)
+    assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
+
+
+@requires_scipy
+def test_wvm_series_xyaslist():
+    csi_series = pd.Series(clear_sky_index, index=time)
+    cs_sm, _, _ = scaling.wvm(csi_series, positions.tolist(), cloud_speed)
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
 @requires_scipy
-def test_wvm_array_invalid():
+def test_wvm_invalid():
     with pytest.raises(ValueError):
-        scaling.wvm(clear_sky_index, lat, lon, cloud_speed, "discrete")
+        scaling.wvm(clear_sky_index, positions, cloud_speed)

From 20cc3c9e2dfee79fc6ea219f68a37a1e56d6887c Mon Sep 17 00:00:00 2001
From: Joe Ranalli <jar339@psu.edu>
Date: Wed, 13 Nov 2019 10:17:32 +0100
Subject: [PATCH 07/13] Added my name to the list in the what's new docs. Hope
 this is the intent.

---
 docs/sphinx/source/whatsnew/v0.7.0.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/sphinx/source/whatsnew/v0.7.0.rst b/docs/sphinx/source/whatsnew/v0.7.0.rst
index 74a1699249..0bb7c478a0 100644
--- a/docs/sphinx/source/whatsnew/v0.7.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.0.rst
@@ -170,3 +170,4 @@ Contributors
 * Miguel Sánchez de León Peque (:ghuser:`Peque`)
 * Tanguy Lunel (:ghuser:`tylunel`)
 * Veronica Guo (:ghuser:`veronicaguo`)
+* Joseph Ranalli (:ghuser:`jranalli`)

From 4135a505f39e13ef68334c68c219e2bd6ecddc5a Mon Sep 17 00:00:00 2001
From: Joe Ranalli <jar339@psu.edu>
Date: Wed, 27 Nov 2019 19:51:19 +0100
Subject: [PATCH 08/13] rename internal variables, truncate comment for lint

---
 pvlib/scaling.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index 23e2dbacbc..cc09df64fb 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -210,13 +210,13 @@ def _compute_wavelet(clearsky_index, dt=None):
     cs_long = pd.DataFrame(cs_long)
 
     # Compute wavelet time scales
-    mindt = np.ceil(np.log(dt)/np.log(2))  # Minimum wavelet dt
-    maxdt = int(12 - mindt)  # maximum wavelet dt
+    min_tmscale = np.ceil(np.log(dt)/np.log(2))  # Minimum wavelet timescale
+    max_tmscale = int(12 - min_tmscale)  # maximum wavelet timescale
 
-    tmscales = np.zeros(maxdt)
-    csi_mean = np.zeros([maxdt, len(cs_long)])
+    tmscales = np.zeros(max_tmscale)
+    csi_mean = np.zeros([max_tmscale, len(cs_long)])
     # Loop for all time scales we will consider
-    for i in np.arange(0, maxdt):
+    for i in np.arange(0, max_tmscale):
         j = i+1
         tmscales[i] = 2**j * dt  # Wavelet integration time scale
         intvlen = 2**j  # Wavelet integration time series interval
@@ -229,13 +229,13 @@ def _compute_wavelet(clearsky_index, dt=None):
 
     # Calculate the wavelets by isolating the rolling mean frequency ranges
     wavelet_long = np.zeros(csi_mean.shape)
-    for i in np.arange(0, maxdt-1):
+    for i in np.arange(0, max_tmscale-1):
         wavelet_long[i, :] = csi_mean[i, :] - csi_mean[i+1, :]
-    wavelet_long[maxdt-1, :] = csi_mean[maxdt-1, :]  # Lowest frequency
+    wavelet_long[max_tmscale-1, :] = csi_mean[max_tmscale-1, :]  # Lowest freq
 
     # Clip off the padding and just return the original time window
-    wavelet = np.zeros([maxdt, len(vals)])
-    for i in np.arange(0, maxdt):
+    wavelet = np.zeros([max_tmscale, len(vals)])
+    for i in np.arange(0, max_tmscale):
         wavelet[i, :] = wavelet_long[i, len(vals)+1: 2*len(vals)+1]
 
     return wavelet, tmscales

From 6bd1c09eda830569a9681ebe9a7aaffafcd7c5da Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Thu, 28 Nov 2019 09:02:27 +0100
Subject: [PATCH 09/13] Removed section headings from api.rst to shorten

---
 docs/sphinx/source/api.rst | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index b511cc8a6b..b9c7a82651 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -560,7 +560,6 @@ Bifacial
 ========
 
 Methods for calculating back surface irradiance
------------------------------------------------
 
 .. autosummary::
    :toctree: generated/
@@ -569,10 +568,9 @@ Methods for calculating back surface irradiance
 
 
 Scaling
-========
+=======
 
 Methods for manipulating irradiance for temporal or spatial considerations
---------------------------------------------------------------------------
 
 .. autosummary::
    :toctree: generated/

From 341814ee9c204156a200a2194b71cc1e1bc3baf7 Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Thu, 28 Nov 2019 09:27:15 +0100
Subject: [PATCH 10/13] Change to try catch block. Added additional test to
 match.

---
 pvlib/scaling.py           | 9 +++++----
 pvlib/test/test_scaling.py | 8 ++++++++
 2 files changed, 13 insertions(+), 4 deletions(-)

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index cc09df64fb..453d1c80f3 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -196,14 +196,15 @@ def _compute_wavelet(clearsky_index, dt=None):
 
     try:  # Assume it's a pandas type
         vals = clearsky_index.values.flatten()
-        try:  # Assume it's a time series type index
-            dt = (clearsky_index.index[1] - clearsky_index.index[0]).seconds
-        except AttributeError:  # It must just be a numeric index
-            dt = (clearsky_index.index[1] - clearsky_index.index[0])
     except AttributeError:  # Assume it's a numpy type
         vals = clearsky_index.flatten()
         if dt is None:
             raise ValueError("dt must be specified for numpy type inputs.")
+    else: # flatten() succeeded, thus it's a pandas type, so get its dt
+        try:  # Assume it's a time series type index
+            dt = (clearsky_index.index[1] - clearsky_index.index[0]).seconds
+        except AttributeError:  # It must just be a numeric index
+            dt = (clearsky_index.index[1] - clearsky_index.index[0])
 
     # Pad the series on both ends in time and place in a dataframe
     cs_long = np.pad(vals, (len(vals), len(vals)), 'symmetric')
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index 3f856da3fe..c8103e8002 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -73,6 +73,14 @@ def test_compute_wavelet_series():
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
 
 
+def test_compute_wavelet_series_numindex():
+    dtindex = pd.to_datetime(time, unit='s')
+    csi_series = pd.Series(clear_sky_index, index=dtindex)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(tmscale, expect_tmscale)
+    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+
+
 def test_compute_wavelet_array():
     wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
     assert_almost_equal(tmscale, expect_tmscale)

From 84cd0fb28282e11e0bf1ed6f31c3883b364d9ced Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Thu, 28 Nov 2019 09:55:55 +0100
Subject: [PATCH 11/13] convert tests to use fixtures

---
 pvlib/test/test_scaling.py | 99 ++++++++++++++++++++++++--------------
 1 file changed, 63 insertions(+), 36 deletions(-)

diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index c8103e8002..ac0f0ad8be 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -8,38 +8,62 @@
 from conftest import requires_scipy
 
 
-# Sample positions
-lat = np.array((9.99, 10, 10.01))
-lon = np.array((4.99, 5, 5.01))
-coordinates = np.array([(lati, loni) for (lati, loni) in zip(lat, lon)])
-
 # Sample cloud speed
 cloud_speed = 5
-# Generate a sample clear_sky_index and time vector.
-clear_sky_index = np.ones(10000)
-clear_sky_index[5000:5005] = np.array([1, 1, 1.1, 0.9, 1])
-time = np.arange(0, len(clear_sky_index))
 
 # Sample dt
 dt = 1
 
-# Expected positions for sample lat/lon given above (calculated manually)
-expect_xpos = np.array([554863.4, 555975.4, 557087.3])
-expect_ypos = np.array([1110838.8, 1111950.8, 1113062.7])
 
-# Sample positions based on the previous lat/lon
-positions = np.array([pt for pt in zip(expect_xpos, expect_ypos)])
+@pytest.fixture
+def coordinates():
+    # Sample positions in lat/lon
+    lat = np.array((9.99, 10, 10.01))
+    lon = np.array((4.99, 5, 5.01))
+    coordinates = np.array([(lati, loni) for (lati, loni) in zip(lat, lon)])
+    return coordinates
+
+
+@pytest.fixture
+def clear_sky_index():
+    # Generate a sample clear_sky_index
+    clear_sky_index = np.ones(10000)
+    clear_sky_index[5000:5005] = np.array([1, 1, 1.1, 0.9, 1])
+    return clear_sky_index
+
+
+@pytest.fixture
+def time(clear_sky_index):
+    # Sample time vector
+    return np.arange(0, len(clear_sky_index))
+
+
+@pytest.fixture
+def positions():
+    # Sample positions based on the previous lat/lon (calculated manually)
+    expect_xpos = np.array([554863.4, 555975.4, 557087.3])
+    expect_ypos = np.array([1110838.8, 1111950.8, 1113062.7])
+    return np.array([pt for pt in zip(expect_xpos, expect_ypos)])
+
+
+@pytest.fixture
+def expect_tmscale():
+    # Expected timescales for dt = 1
+    return [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+
 
-# Expected timescales for dt = 1
-expect_tmscale = [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+@pytest.fixture
+def expect_wavelet():
+    # Expected wavelet for indices 5000:5004 with clear_sky_index above (Matlab)
+    return np.array([[-0.025, 0.05, 0., -0.05, 0.025],
+                     [0.025, 0., 0., 0., -0.025],
+                     [0., 0., 0., 0., 0.]])
 
-# Expected wavelet for indices 5000:5004 using clear_sky_index above (Matlab)
-expect_wavelet = np.array([[-0.025, 0.05, 0., -0.05, 0.025],
-                           [0.025, 0., 0., 0., -0.025],
-                           [0., 0., 0., 0., 0.]])
 
-# Expected smoothed clear sky index for indices 5000:5004 (Matlab)
-expect_cs_smooth = np.array([1., 1.0289, 1., 0.9711, 1.])
+@pytest.fixture
+def expect_cs_smooth():
+    # Expected smoothed clear sky index for indices 5000:5004 (Matlab)
+    return np.array([1., 1.0289, 1., 0.9711, 1.])
 
 
 def test_latlon_to_xy_zero():
@@ -49,31 +73,33 @@ def test_latlon_to_xy_zero():
     assert_almost_equal(pos, pos_e, decimal=1)
 
 
-def test_latlon_to_xy_single():
-    # Must test against central value, because latlon_to_dist uses the mean
-    coord = (lat[1], lon[1])
+def test_latlon_to_xy_single(coordinates, positions):
+    # Must test against central value, because latlon_to_xy uses the mean
+    coord = coordinates[1]
     pos = scaling.latlon_to_xy(coord)
-    assert_almost_equal(pos, (expect_xpos[1], expect_ypos[1]), decimal=1)
+    assert_almost_equal(pos, positions[1], decimal=1)
 
 
-def test_latlon_to_xy_array():
+def test_latlon_to_xy_array(coordinates, positions):
     pos = scaling.latlon_to_xy(coordinates)
     assert_almost_equal(pos, positions, decimal=1)
 
 
-def test_latlon_to_xy_list():
+def test_latlon_to_xy_list(coordinates, positions):
     pos = scaling.latlon_to_xy(coordinates.tolist())
     assert_almost_equal(pos, positions, decimal=1)
 
 
-def test_compute_wavelet_series():
+def test_compute_wavelet_series(clear_sky_index, time,
+                                expect_tmscale, expect_wavelet):
     csi_series = pd.Series(clear_sky_index, index=time)
     wavelet, tmscale = scaling._compute_wavelet(csi_series)
     assert_almost_equal(tmscale, expect_tmscale)
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
 
 
-def test_compute_wavelet_series_numindex():
+def test_compute_wavelet_series_numindex(clear_sky_index, time,
+                                         expect_tmscale, expect_wavelet):
     dtindex = pd.to_datetime(time, unit='s')
     csi_series = pd.Series(clear_sky_index, index=dtindex)
     wavelet, tmscale = scaling._compute_wavelet(csi_series)
@@ -81,38 +107,39 @@ def test_compute_wavelet_series_numindex():
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
 
 
-def test_compute_wavelet_array():
+def test_compute_wavelet_array(clear_sky_index, expect_tmscale, expect_wavelet):
     wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
     assert_almost_equal(tmscale, expect_tmscale)
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
 
 
-def test_compute_wavelet_array_invalid():
+def test_compute_wavelet_array_invalid(clear_sky_index):
     with pytest.raises(ValueError):
         scaling._compute_wavelet(clear_sky_index)
 
 
 @requires_scipy
-def test_wvm_series():
+def test_wvm_series(clear_sky_index, time, positions, expect_cs_smooth):
     csi_series = pd.Series(clear_sky_index, index=time)
     cs_sm, _, _ = scaling.wvm(csi_series, positions, cloud_speed)
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
 @requires_scipy
-def test_wvm_array():
+def test_wvm_array(clear_sky_index, positions, expect_cs_smooth):
     cs_sm, _, _ = scaling.wvm(clear_sky_index, positions, cloud_speed, dt=dt)
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
 @requires_scipy
-def test_wvm_series_xyaslist():
+def test_wvm_series_xyaslist(clear_sky_index, time, positions,
+                             expect_cs_smooth):
     csi_series = pd.Series(clear_sky_index, index=time)
     cs_sm, _, _ = scaling.wvm(csi_series, positions.tolist(), cloud_speed)
     assert_almost_equal(cs_sm[5000:5005], expect_cs_smooth, decimal=4)
 
 
 @requires_scipy
-def test_wvm_invalid():
+def test_wvm_invalid(clear_sky_index, positions):
     with pytest.raises(ValueError):
         scaling.wvm(clear_sky_index, positions, cloud_speed)

From b62f16f5a1d41f0575b62a1a52549dcc53a5aa94 Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Thu, 28 Nov 2019 10:03:12 +0100
Subject: [PATCH 12/13] fix lint errors

---
 pvlib/scaling.py           | 2 +-
 pvlib/test/test_scaling.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index 453d1c80f3..c7ec11b7c3 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -200,7 +200,7 @@ def _compute_wavelet(clearsky_index, dt=None):
         vals = clearsky_index.flatten()
         if dt is None:
             raise ValueError("dt must be specified for numpy type inputs.")
-    else: # flatten() succeeded, thus it's a pandas type, so get its dt
+    else:  # flatten() succeeded, thus it's a pandas type, so get its dt
         try:  # Assume it's a time series type index
             dt = (clearsky_index.index[1] - clearsky_index.index[0]).seconds
         except AttributeError:  # It must just be a numeric index
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index ac0f0ad8be..964aad947c 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -54,7 +54,7 @@ def expect_tmscale():
 
 @pytest.fixture
 def expect_wavelet():
-    # Expected wavelet for indices 5000:5004 with clear_sky_index above (Matlab)
+    # Expected wavelet for indices 5000:5004 for clear_sky_index above (Matlab)
     return np.array([[-0.025, 0.05, 0., -0.05, 0.025],
                      [0.025, 0., 0., 0., -0.025],
                      [0., 0., 0., 0., 0.]])

From a310f1a76082600c88b030baa3685c944393256e Mon Sep 17 00:00:00 2001
From: "Ranalli, Joe" <jar339@psu.edu>
Date: Thu, 28 Nov 2019 10:10:09 +0100
Subject: [PATCH 13/13] one more lint error

---
 pvlib/test/test_scaling.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
index 964aad947c..62f37794d2 100644
--- a/pvlib/test/test_scaling.py
+++ b/pvlib/test/test_scaling.py
@@ -107,7 +107,8 @@ def test_compute_wavelet_series_numindex(clear_sky_index, time,
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
 
 
-def test_compute_wavelet_array(clear_sky_index, expect_tmscale, expect_wavelet):
+def test_compute_wavelet_array(clear_sky_index,
+                               expect_tmscale, expect_wavelet):
     wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
     assert_almost_equal(tmscale, expect_tmscale)
     assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)