diff --git a/docs/sphinx/source/whatsnew/v0.9.0.rst b/docs/sphinx/source/whatsnew/v0.9.0.rst index fa7ed689b5..745903fcce 100644 --- a/docs/sphinx/source/whatsnew/v0.9.0.rst +++ b/docs/sphinx/source/whatsnew/v0.9.0.rst @@ -187,6 +187,9 @@ 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`) * Changed deprecated use of ``.astype()`` to ``.view()`` in :py:mod:`~pvlib.solarposition`. (:pull:`1256`, :issue:`1261`, :pull:`1262`) diff --git a/pvlib/scaling.py b/pvlib/scaling.py index a288912e8d..dca2ca4935 100644 --- a/pvlib/scaling.py +++ b/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 diff --git a/pvlib/tests/test_scaling.py b/pvlib/tests/test_scaling.py index ba2c00ae3f..344e2209b5 100644 --- a/pvlib/tests/test_scaling.py +++ b/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)