bugfix for high frequency time series in scaling.py, regarding issue 1257. Includes refactor to separate out VR computation within pvlib.scaling.wvm to a new function pvlib.scaling._compute_vr.

Author: jranalli

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

@@ -175,6 +175,10 @@ Bug fixes
 * Corrected an error affecting :py:func:`~pvlib.clearsky.detect_clearsky`
   when data time step is not one minute. Error was introduced in v0.8.1.
   (:issue:`1241`, :pull:`1242`)
+* Corrected error affecting :py:func:`~pvlib.scaling._compute_wavelet` when
+  passing a pandas time series with a sampling rate faster than 1 second.
+  (:issue:`1257`, :pull:`1258`)
+  (:issue:`1257`, :pull:``)
 
 Testing
 ~~~~~~~

pvlib/scaling.py

@@ -62,26 +62,9 @@ def wvm(clearsky_index, positions, cloud_speed, dt=None):
 
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
 
-    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)
-
-    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
-    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]
+    vr = _compute_vr(positions, cloud_speed, tmscales)
 
     # Scale each wavelet by VR (Eq 7 in [1])
     wavelet_smooth = np.zeros_like(wavelet)
@@ -101,6 +84,68 @@ def fn(x):
     return smoothed, wavelet, tmscales
 
 
+def _compute_vr(positions, cloud_speed, tmscales):
+    """
+    Compute the variability reduction factors for each wavelet mode for the
+    Wavelet Variability Model [1-3].
+
+    Parameters
+    ----------
+    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 [m/s].
+
+    tmscales: numeric
+        The timescales associated with the wavelets in seconds [s].
+
+    Returns
+    -------
+    vr : numeric
+        an array of variability reduction factors for each tmscale.
+
+    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://github.com/sandialabs/wvm
+    """
+
+    # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2021
+
+    pos = np.array(positions)
+    dist = pdist(pos, 'euclidean')
+
+    # Find effective length of position vector, 'dist' is full pairwise
+    n_pairs = len(dist)
+
+    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
+    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]
+    return vr
+
+
 def latlon_to_xy(coordinates):
     """
     Convert latitude and longitude in degrees to a coordinate system measured
@@ -205,7 +250,8 @@ def _compute_wavelet(clearsky_index, dt=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
+            dt = clearsky_index.index[1] - clearsky_index.index[0]
+            dt = dt.seconds + dt.microseconds/1e6
         except AttributeError:  # It must just be a numeric index
             dt = (clearsky_index.index[1] - clearsky_index.index[0])
 
@@ -221,7 +267,7 @@ def _compute_wavelet(clearsky_index, dt=None):
     csi_mean = np.zeros([max_tmscale, len(cs_long)])
     # Skip averaging for the 0th scale
     csi_mean[0, :] = cs_long.values.flatten()
-    tmscales[0] = 1
+    tmscales[0] = dt
     # Loop for all time scales we will consider
     for i in np.arange(1, max_tmscale):
         tmscales[i] = 2**i * dt  # Wavelet integration time scale

pvlib/tests/test_scaling.py

@@ -37,6 +37,18 @@ def time(clear_sky_index):
     return np.arange(0, len(clear_sky_index))
 
 
+@pytest.fixture
+def time_60s(clear_sky_index):
+    # Sample time vector 60s resolution
+    return np.arange(0, len(clear_sky_index))*60
+
+
+@pytest.fixture
+def time_500ms(clear_sky_index):
+    # Sample time vector 0.5s resolution
+    return np.arange(0, len(clear_sky_index))*0.5
+
+
 @pytest.fixture
 def positions():
     # Sample positions based on the previous lat/lon (calculated manually)
@@ -51,6 +63,18 @@ def expect_tmscale():
     return [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
 
 
+@pytest.fixture
+def expect_tmscale_1min():
+    # Expected timescales for dt = 60
+    return [60, 120, 240, 480, 960, 1920, 3840]
+
+
+@pytest.fixture
+def expect_tmscale_500ms():
+    # Expected timescales for dt = 0.5
+    return [0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+
+
 @pytest.fixture
 def expect_wavelet():
     # Expected wavelet for indices 5000:5004 for clear_sky_index above (Matlab)
@@ -68,6 +92,13 @@ def expect_cs_smooth():
     return np.array([1., 1., 1.05774, 0.94226, 1.])
 
 
+@pytest.fixture
+def expect_vr():
+    # Expected VR for expecttmscale
+    return np.array([3., 3., 3., 3., 3., 3., 2.9997844, 2.9708118, 2.6806291,
+                     2.0726611, 1.5653324, 1.2812714, 1.1389995])
+
+
 def test_latlon_to_xy_zero():
     coord = [0, 0]
     pos_e = [0, 0]
@@ -109,6 +140,25 @@ def test_compute_wavelet_series_numindex(clear_sky_index, time,
     assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
+def test_compute_wavelet_series_highres(clear_sky_index, time_500ms,
+                                        expect_tmscale_500ms, expect_wavelet):
+    dtindex = pd.to_datetime(time_500ms, unit='s')
+    csi_series = pd.Series(clear_sky_index, index=dtindex)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(tmscale, expect_tmscale_500ms)
+    assert_almost_equal(wavelet[:, 5000:5005].shape, (14, 5))
+
+
+def test_compute_wavelet_series_minuteres(clear_sky_index, time_60s,
+                                          expect_tmscale_1min, expect_wavelet):
+    dtindex = pd.to_datetime(time_60s, unit='s')
+    csi_series = pd.Series(clear_sky_index, index=dtindex)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(tmscale, expect_tmscale_1min)
+    assert_almost_equal(wavelet[:, 5000:5005].shape,
+                        expect_wavelet[0:len(tmscale), :].shape)
+
+
 def test_compute_wavelet_array(clear_sky_index,
                                expect_tmscale, expect_wavelet):
     wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
@@ -129,6 +179,11 @@ def test_compute_wavelet_dwttheory(clear_sky_index, time,
     assert_almost_equal(np.sum(wavelet, 0), csi_series)
 
 
+def test_compute_vr(positions, expect_tmscale, expect_vr):
+    vr = scaling._compute_vr(positions, cloud_speed, np.array(expect_tmscale))
+    assert_almost_equal(vr, expect_vr)
+
+
 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)