From b9e57f02f433d32f6feea0340e17ca2ea0d19b19 Mon Sep 17 00:00:00 2001 From: Joe Ranalli Date: Thu, 22 Apr 2021 13:09:18 -0400 Subject: [PATCH 1/7] Add theoretical DWT test for scaling.py --- pvlib/tests/test_scaling.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pvlib/tests/test_scaling.py b/pvlib/tests/test_scaling.py index 09c81ae3c7..f00a1237dd 100644 --- a/pvlib/tests/test_scaling.py +++ b/pvlib/tests/test_scaling.py @@ -118,6 +118,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) From 73e7d08c78a748f525d13ad756d0b06f46a56f83 Mon Sep 17 00:00:00 2001 From: Joe Ranalli Date: Thu, 22 Apr 2021 20:26:06 -0400 Subject: [PATCH 2/7] Adjust test fixtures relative to outputs from corrected MATLAB version. --- pvlib/tests/test_scaling.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/pvlib/tests/test_scaling.py b/pvlib/tests/test_scaling.py index f00a1237dd..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): From a6efc7016b00812ae6cdd75264ba0ad1d751e335 Mon Sep 17 00:00:00 2001 From: Joe Ranalli Date: Thu, 22 Apr 2021 21:05:40 -0400 Subject: [PATCH 3/7] Implement fix to Python version of code. Some additional work needed to be done here to get the wavelet levels to match the MATLAB output. The only place python and MATLAB don't agree now is at the edges, where there are different assumptions about how the moving average should be computed for the beginning/end of the time series. In principle, this could actually matter for time series that aren't "much" longer than 4096 seconds. --- pvlib/scaling.py | 34 +++++++++++++++++++++++----------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/pvlib/scaling.py b/pvlib/scaling.py index b6a6caaf66..b42c980cdc 100644 --- a/pvlib/scaling.py +++ b/pvlib/scaling.py @@ -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,7 +179,8 @@ 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] @@ -185,7 +191,7 @@ def _compute_wavelet(clearsky_index, dt=None): 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: + [2] Wavelet Variability Model - Matlab Code: https://pvpmc.sandia.gov/applications/wavelet-variability-model/ """ @@ -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 From 979cbdbaa1cff7376fc83f48e29d8941c8a429ca Mon Sep 17 00:00:00 2001 From: Joe Ranalli Date: Thu, 22 Apr 2021 21:38:30 -0400 Subject: [PATCH 4/7] Update docs --- docs/sphinx/source/whatsnew/v0.9.0.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/sphinx/source/whatsnew/v0.9.0.rst b/docs/sphinx/source/whatsnew/v0.9.0.rst index 9a8f3878b7..570ca49024 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 scaling.py WVM. Tracks with fix in MATLAB ver. + (:issue:`1206`, :pull:`1213`) Testing ~~~~~~~ @@ -148,3 +150,4 @@ Contributors * Joshua Stein (:ghuser:`jsstein`) * Tony Lorenzo (:ghuser:`alorenzo175`) * Damjan Postolovski (:ghuser:`dpostolovski`) +* Joe Ranalli (:ghuser:`jranalli`) From 54928eb62b71a3c72c89aa3a3543f93e65b1bef3 Mon Sep 17 00:00:00 2001 From: Joe Ranalli Date: Mon, 26 Apr 2021 16:05:54 -0400 Subject: [PATCH 5/7] Update reference to WVM in scaling.py and formatting thereof. --- pvlib/scaling.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/pvlib/scaling.py b/pvlib/scaling.py index b42c980cdc..3d46f8e8bd 100644 --- a/pvlib/scaling.py +++ b/pvlib/scaling.py @@ -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 @@ -187,12 +187,12 @@ def _compute_wavelet(clearsky_index, 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] 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 From 56a892885d6a991c8d5757ff11e35661b337363a Mon Sep 17 00:00:00 2001 From: Joe Ranalli Date: Thu, 6 May 2021 11:33:51 -0400 Subject: [PATCH 6/7] Suggested update to the documentation language. Co-authored-by: Cliff Hansen --- docs/sphinx/source/whatsnew/v0.9.0.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.0.rst b/docs/sphinx/source/whatsnew/v0.9.0.rst index 570ca49024..151426ba65 100644 --- a/docs/sphinx/source/whatsnew/v0.9.0.rst +++ b/docs/sphinx/source/whatsnew/v0.9.0.rst @@ -117,7 +117,7 @@ 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 scaling.py WVM. Tracks with fix in MATLAB ver. +* Corrected methodology error in scaling.py WVM. Tracks with fix in PVLib for MATLAB. (:issue:`1206`, :pull:`1213`) Testing From 0053e59c39e5a1f66fbfb0bd2e6f2aeb4d1cd093 Mon Sep 17 00:00:00 2001 From: Cliff Hansen Date: Fri, 7 May 2021 09:57:25 -0600 Subject: [PATCH 7/7] docstring and whatsnew adjustments --- docs/sphinx/source/whatsnew/v0.9.0.rst | 5 ++--- pvlib/scaling.py | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/docs/sphinx/source/whatsnew/v0.9.0.rst b/docs/sphinx/source/whatsnew/v0.9.0.rst index 151426ba65..5a2fa4ed02 100644 --- a/docs/sphinx/source/whatsnew/v0.9.0.rst +++ b/docs/sphinx/source/whatsnew/v0.9.0.rst @@ -117,8 +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 scaling.py WVM. Tracks with fix in PVLib for MATLAB. - (:issue:`1206`, :pull:`1213`) +* Corrected methodology error in :py:func:`~pvlib.scaling.wvm`. Tracks with + fix in PVLib for MATLAB. (:issue:`1206`, :pull:`1213`) Testing ~~~~~~~ @@ -128,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` diff --git a/pvlib/scaling.py b/pvlib/scaling.py index 3d46f8e8bd..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 ----------