diff --git a/docs/sphinx/source/whatsnew/v0.9.0.rst b/docs/sphinx/source/whatsnew/v0.9.0.rst index 9a8f3878b7..5a2fa4ed02 100644 --- a/docs/sphinx/source/whatsnew/v0.9.0.rst +++ b/docs/sphinx/source/whatsnew/v0.9.0.rst @@ -117,6 +117,8 @@ Bug fixes * Reindl model fixed to generate sky_diffuse=0 when GHI=0. (:issue:`1153`, :pull:`1154`) * Update GFS product names for GFS v16. (:issue:`1202`, :pull:`1203`) +* Corrected methodology error in :py:func:`~pvlib.scaling.wvm`. Tracks with + fix in PVLib for MATLAB. (:issue:`1206`, :pull:`1213`) Testing ~~~~~~~ @@ -126,7 +128,6 @@ Testing Documentation ~~~~~~~~~~~~~ - * Update intro tutorial to highlight the use of historical meteorological data and to make the procedural and object oriented results match exactly. * Update documentation links in :py:func:`pvlib.iotools.get_psm3` @@ -148,3 +149,4 @@ Contributors * Joshua Stein (:ghuser:`jsstein`) * Tony Lorenzo (:ghuser:`alorenzo175`) * Damjan Postolovski (:ghuser:`dpostolovski`) +* Joe Ranalli (:ghuser:`jranalli`) diff --git a/pvlib/scaling.py b/pvlib/scaling.py index b6a6caaf66..a288912e8d 100644 --- a/pvlib/scaling.py +++ b/pvlib/scaling.py @@ -13,8 +13,8 @@ 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 - basically a port of the Matlab version of the code [3]. + 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 ---------- @@ -48,16 +48,16 @@ def wvm(clearsky_index, positions, cloud_speed, dt=None): 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. + .. [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. + .. [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/ + .. [3] Wavelet Variability Model - Matlab Code: + https://github.com/sandialabs/wvm """ # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019 @@ -128,13 +128,13 @@ def latlon_to_xy(coordinates): References ---------- - [1] H. Moritz. Geodetic Reference System 1980, Journal of Geodesy, vol. 74, - no. 1, pp 128–133, 2000. + .. [1] H. Moritz. Geodetic Reference System 1980, Journal of Geodesy, vol. + 74, no. 1, pp 128–133, 2000. - [2] https://pypi.org/project/pyproj/ + .. [2] https://pypi.org/project/pyproj/ - [3] Wavelet Variability Model - Matlab Code: - https://pvpmc.sandia.gov/applications/wavelet-variability-model/ + .. [3] Wavelet Variability Model - Matlab Code: + https://github.com/sandialabs/wvm """ # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019 @@ -159,7 +159,12 @@ def latlon_to_xy(coordinates): def _compute_wavelet(clearsky_index, dt=None): """ - Compute the wavelet transform on the input clear_sky time series. + Compute the wavelet transform on the input clear_sky time series. Uses a + top hat wavelet [-1,1,1,-1] shape, based on the difference of successive + centered moving averages. Smallest scale (filter size of 2) is a degenerate + case that resembles a Haar wavelet. Returns one level of approximation + coefficient (CAn) and n levels of detail coefficients (CD1, CD2, ..., + CDn-1, CDn). Parameters ---------- @@ -174,19 +179,20 @@ def _compute_wavelet(clearsky_index, dt=None): Returns ------- wavelet: numeric - The individual wavelets for the time series + The individual wavelets for the time series. Format follows increasing + scale (decreasing frequency): [CD1, CD2, ..., CDn, CAn] tmscales: numeric The timescales associated with the wavelets in seconds [s] 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. + .. [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/ + .. [2] Wavelet Variability Model - Matlab Code: + https://github.com/sandialabs/wvm """ # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019 @@ -209,31 +215,37 @@ def _compute_wavelet(clearsky_index, dt=None): # Compute wavelet time scales min_tmscale = np.ceil(np.log(dt)/np.log(2)) # Minimum wavelet timescale - max_tmscale = int(12 - min_tmscale) # maximum wavelet timescale + max_tmscale = int(13 - min_tmscale) # maximum wavelet timescale tmscales = np.zeros(max_tmscale) 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 # Loop for all time scales we will consider - 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 + for i in np.arange(1, max_tmscale): + tmscales[i] = 2**i * dt # Wavelet integration time scale + intvlen = 2**i # Wavelet integration time series interval # Rolling average, retains only lower frequencies than interval + # Produces slightly different end effects than the MATLAB version 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 csi_mean[i, :] = df.values.flatten() + # Shift to account for different indexing in MATLAB moving average + csi_mean[i, :] = np.roll(csi_mean[i, :], -1) + csi_mean[i, -1] = csi_mean[i, -2] - # Calculate the wavelets by isolating the rolling mean frequency ranges + # Calculate detail coefficients by difference between successive averages wavelet_long = np.zeros(csi_mean.shape) for i in np.arange(0, max_tmscale-1): wavelet_long[i, :] = csi_mean[i, :] - csi_mean[i+1, :] - wavelet_long[max_tmscale-1, :] = csi_mean[max_tmscale-1, :] # Lowest freq + wavelet_long[-1, :] = csi_mean[-1, :] # Lowest freq (CAn) # Clip off the padding and just return the original time window 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] + wavelet[i, :] = wavelet_long[i, len(vals): 2*len(vals)] return wavelet, tmscales diff --git a/pvlib/tests/test_scaling.py b/pvlib/tests/test_scaling.py index 09c81ae3c7..ba2c00ae3f 100644 --- a/pvlib/tests/test_scaling.py +++ b/pvlib/tests/test_scaling.py @@ -48,21 +48,24 @@ def positions(): @pytest.fixture def expect_tmscale(): # Expected timescales for dt = 1 - return [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096] + return [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) - return np.array([[-0.025, 0.05, 0., -0.05, 0.025], - [0.025, 0., 0., 0., -0.025], - [0., 0., 0., 0., 0.]]) + e = np.zeros([13, 5]) + e[0, :] = np.array([0, -0.05, 0.1, -0.05, 0]) + e[1, :] = np.array([-0.025, 0.05, 0., -0.05, 0.025]) + e[2, :] = np.array([0.025, 0., 0., 0., -0.025]) + e[-1, :] = np.array([1, 1, 1, 1, 1]) + return e @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.]) + return np.array([1., 1., 1.05774, 0.94226, 1.]) def test_latlon_to_xy_zero(): @@ -94,7 +97,7 @@ def test_compute_wavelet_series(clear_sky_index, time, 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) + assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet) def test_compute_wavelet_series_numindex(clear_sky_index, time, @@ -103,14 +106,14 @@ def test_compute_wavelet_series_numindex(clear_sky_index, time, 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) + assert_almost_equal(wavelet[:, 5000:5005], 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) + assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet) def test_compute_wavelet_array_invalid(clear_sky_index): @@ -118,6 +121,14 @@ def test_compute_wavelet_array_invalid(clear_sky_index): scaling._compute_wavelet(clear_sky_index) +def test_compute_wavelet_dwttheory(clear_sky_index, time, + expect_tmscale, expect_wavelet): + # Confirm detail coeffs sum to original signal + csi_series = pd.Series(clear_sky_index, index=time) + wavelet, tmscale = scaling._compute_wavelet(csi_series) + assert_almost_equal(np.sum(wavelet, 0), csi_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)