Repository: pvlib/pvlib-python
Pull request #807 (2019-11-01 14:54:52)

Create scaling.py and implement WVM model
 - [X] Closes  #806 
 - [X] I am familiar with the [contributing guidelines](http://pvlib-python.readthedocs.io/en/latest/contributing.html).
 - [x] Fully tested. Added and/or modified tests to ensure correct behavior for all reasonable inputs. Tests (usually) must pass on the TravisCI and Appveyor testing services.
 - [x] Updates entries to `docs/sphinx/source/api.rst` for API changes.
 - [x] Adds description and name entries in the appropriate `docs/sphinx/source/whatsnew` file for all changes.
 - [x] Code quality and style is sufficient. Passes LGTM and SticklerCI checks.
 - [x] New code is fully documented. Includes sphinx/numpydoc compliant docstrings and comments in the code where necessary.
 - [x] Pull request is nearly complete and ready for detailed review.

Brief description of the problem and proposed solution (if not already fully described in the issue linked to above):

I put in just the "discrete" version of distance finding for now, because I'm already looking at quite a bit of code, and the guidelines suggest keeping it short as is reasonable. I have the other distance finding modes complete, but they're an additional 75 lines of code within the main algorithm, and they're likely to prompt a bit of discussion as to best approach.

cwhanse: @jranalli thanks! If you are OK with it, I'd like to use this PR as a way for us to explore improving our review process - something we've been discussing but haven't yet tried out.

We should review the API (first priority) and numerical performance. I'd like to communicate who plans to review, so please reply here.

I will ask the originator of the algorithm to check the porting.

jranalli: @cwhanse I'm completely ok with a new process. Any process would be new to me, so it's all the same.

@kevinsa5 oops you're right. I'll straighten that out.
wholmgren: Thanks @jranalli for the submission and thanks @cwhanse for coordinating. I'd use this code and I am interested in reviewing this PR. I skimmed the code and it looks like a great start. I agree that this scope is appropriate and other methods would be best implemented in future PRs.
adriesse: @cwhanse Thanks for the invitation, but I need a breather; will try to resist temptations and sit this one out.
cwhanse: @jranalli I pushed a merge to your branch to fix the conflict in the whatsnew file (the conflict was nothing you did, it happens sometimes). I should have notified you first, my apologies. Just pull your branch back from GitHub to update your local files and all should be fine.

The test failures are also nothing you did - we have some functions in pvlib that fetch data from servers, and when the servers don't respond during tests, the test fail. Often, these issues are resolved by restarting the build from the test report page. 
jranalli: Being new to this, I had been assuming I should wait to make any changes until discussion is finished, but I also saw your comment about "moving review ahead." Is there anything I should start addressing (e.g. the obvious docstring corrections), or should I just continue to hang tight?
cwhanse: > Being new to this, I had been assuming I should wait to make any changes until discussion is finished, but I also saw your comment about "moving review ahead." Is there anything I should start addressing (e.g. the obvious docstring corrections), or should I just continue to hang tight?

The main issue is changing from lat, long to x, y. I'd like others to weigh in.
mikofski: I haven't looked over the code in great detail, but I agree with @cwhanse API suggestions:
- [ ] change inputs from separate latitude and longitude lists to a single list of point tuples using `(x, y)`
- [ ] remove the `method` argument and explain and/or provide helper to generate grid for "square" plant type
- [ ] remove lat/lon(wgs84) to easting/northing (bng) private function, and just input points as `(x, y)` pairs, and/or possibly make it public helper in `tools` with caveats and links to other projection conversion methods
cwhanse: > Being new to this, I had been assuming I should wait to make any changes until discussion is finished, but I also saw your comment about "moving review ahead." Is there anything I should start addressing (e.g. the obvious docstring corrections), or should I just continue to hang tight?
> 
> The main issue is changing from lat, long to x, y. I'd like others to weigh in.

@jranalli I think we have consensus on x,y as input rather than lat/long, and on removing the `method` kwarg. I would use meters as the unit on the coordinates.

I'm still waiting to hear from Matt Lave about a technical review, I'll let you know soon.
jranalli: These are the API changes that were requested and I think everything is back to working correctly.
wholmgren: What do we need to get this PR across the finish line? A technical review? 
cwhanse: > What do we need to get this PR across the finish line? A technical review?

Yes. I'll take care of it today, I don't know when Matt will get to it.
jranalli: @cwhanse My pleasure, thanks to you all for the opportunity to contribute. 

The renaming is done. If there's anything else needed from me on this please just let me know.
cwhanse: 


> @cwhanse My pleasure, thanks to you all for the opportunity to contribute.
> The renaming is done. If there's anything else needed from me on this please just let me know.

Can you merge pvlib-python/master into your branch, and resolve the merge conflict in the whatsnew file?  Or you can do it from this page.
jranalli: I hate to say it, but this kind of thing is where I start pushing the boundaries of my git expertise. I think I did it right, but if not, I might need some guidance.
cwhanse: > I hate to say it, but this kind of thing is where I start pushing the boundaries of my git expertise. I think I did it right, but if not, I might need some guidance.

Looks good now.

Patch:
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 7241f78..b9c7a82 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -560,9 +560,19 @@ Bifacial
 ========
 
 Methods for calculating back surface irradiance
------------------------------------------------
 
 .. autosummary::
    :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 3b83ac9..6de109c 100644
--- a/docs/sphinx/source/whatsnew/v0.7.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.0.rst
@@ -136,6 +136,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.
 * Created one new incidence angle modifier (IAM) function for diffuse irradiance:
   :py:func:`pvlib.iam.martin_ruiz_diffuse`. (:issue:`751`)
 
@@ -197,5 +199,6 @@ Contributors
 * Miguel Sánchez de León Peque (:ghuser:`Peque`)
 * Tanguy Lunel (:ghuser:`tylunel`)
 * Veronica Guo (:ghuser:`veronicaguo`)
+* Joseph Ranalli (:ghuser:`jranalli`)
 * Tony Lorenzo (:ghuser:`alorenzo175`)
 * Todd Karin (:ghuser:`toddkarin`)
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 93491a5..850339b 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1537,26 +1537,41 @@ def retrieve_sam(name=None, path=None):
     return _parse_raw_sam_df(csvdata)
 
 
+def _normalize_sam_product_names(names):
+    '''
+    Replace special characters within the product names to make them more
+    suitable for use as Dataframe column names.
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
+
+    import warnings
+
+    BAD_CHARS  = ' -.()[]:+/",'
+    GOOD_CHARS = '____________'
+
+    mapping = str.maketrans(BAD_CHARS, GOOD_CHARS)
+    names = pd.Series(data=names)
+    norm_names = names.str.translate(mapping)
+
+    n_duplicates = names.duplicated().sum()
+    if n_duplicates > 0:
+        warnings.warn('Original names contain %d duplicate(s).' % n_duplicates)
+
+    n_duplicates = norm_names.duplicated().sum()
+    if n_duplicates > 0:
+        warnings.warn('Normalized names contain %d duplicate(s).' % n_duplicates)
+
+    return norm_names.values
+
+
 def _parse_raw_sam_df(csvdata):
+
     df = pd.read_csv(csvdata, index_col=0, skiprows=[1, 2])
-    colnames = df.columns.values.tolist()
-    parsedcolnames = []
-    for cn in colnames:
-        parsedcolnames.append(cn.replace(' ', '_'))
-
-    df.columns = parsedcolnames
-
-    parsedindex = []
-    for index in df.index:
-        parsedindex.append(index.replace(' ', '_').replace('-', '_')
-                                .replace('.', '_').replace('(', '_')
-                                .replace(')', '_').replace('[', '_')
-                                .replace(']', '_').replace(':', '_')
-                                .replace('+', '_').replace('/', '_')
-                                .replace('"', '_').replace(',', '_'))
-
-    df.index = parsedindex
+
+    df.columns = df.columns.str.replace(' ', '_')
+    df.index = _normalize_sam_product_names(df.index)
     df = df.transpose()
+
     if 'ADRCoefficients' in df.index:
         ad_ce = 'ADRCoefficients'
         # for each inverter, parses a string of coefficients like
diff --git a/pvlib/scaling.py b/pvlib/scaling.py
new file mode 100644
index 0000000..c7ec11b
--- /dev/null
+++ b/pvlib/scaling.py
@@ -0,0 +1,242 @@
+"""
+The ``scaling`` module contains functions for manipulating irradiance
+or other variables to account for temporal or spatial characteristics.
+"""
+
+import numpy as np
+import pandas as pd
+
+
+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].
+
+    Parameters
+    ----------
+    clearsky_index : numeric or pandas.Series
+        Clear Sky Index time series that will be smoothed.
+
+    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].
+
+    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 : 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.
+
+    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.
+
+    [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
+
+    try:
+        import scipy.optimize
+        from scipy.spatial.distance import pdist
+    except ImportError:
+        raise ImportError("The WVM function requires scipy.")
+
+    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]
+
+    # 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=clearsky_index.index)
+    except AttributeError:
+        smoothed = outsignal  # just output the numpy signal
+
+    return smoothed, wavelet, tmscales
+
+
+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
+    ----------
+
+    coordinates : numeric
+        Array or list of (latitude, longitude) coordinate pairs. Use decimal
+        degrees notation.
+
+    Returns
+    -------
+    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] https://pypi.org/project/pyproj/
+
+    [3] 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
+    try:
+        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)
+
+    # reshape as (x,y) pairs to return
+    try:
+        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(clearsky_index, dt=None):
+    """
+    Compute the wavelet transform on the input clear_sky time series.
+
+    Parameters
+    ----------
+    clearsky_index : numeric or pandas.Series
+        Clear Sky Index time series that will be smoothed.
+
+    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
+
+    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.
+
+    [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 = clearsky_index.values.flatten()
+    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')
+    cs_long = pd.DataFrame(cs_long)
+
+    # 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
+
+    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, max_tmscale):
+        j = i+1
+        tmscales[i] = 2**j * dt  # Wavelet integration time scale
+        intvlen = 2**j  # Wavelet integration time series interval
+        # Rolling average, retains only lower frequencies than interval
+        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()
+
+    # Calculate the wavelets by isolating the rolling mean frequency ranges
+    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
+
+    # 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]
+
+    return wavelet, tmscales


Test patch:
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index c9bec95..ecf68cd 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -107,6 +107,28 @@ def test_PVSystem_get_iam_interp(sapm_module_params, mocker):
     with pytest.raises(ValueError):
         system.get_iam(45, iam_model='interp')
 
+def test__normalize_sam_product_names():
+
+    BAD_NAMES  = [' -.()[]:+/",', 'Module[1]']
+    NORM_NAMES = ['____________', 'Module_1_']
+
+    norm_names = pvsystem._normalize_sam_product_names(BAD_NAMES)
+    assert list(norm_names) == NORM_NAMES
+
+    BAD_NAMES  = ['Module[1]', 'Module(1)']
+    NORM_NAMES = ['Module_1_', 'Module_1_']
+
+    with pytest.warns(UserWarning):
+        norm_names = pvsystem._normalize_sam_product_names(BAD_NAMES)
+    assert list(norm_names) == NORM_NAMES
+
+    BAD_NAMES  = ['Module[1]', 'Module[1]']
+    NORM_NAMES = ['Module_1_', 'Module_1_']
+
+    with pytest.warns(UserWarning):
+        norm_names = pvsystem._normalize_sam_product_names(BAD_NAMES)
+    assert list(norm_names) == NORM_NAMES
+
 
 def test_PVSystem_get_iam_invalid(sapm_module_params, mocker):
     system = pvsystem.PVSystem(module_parameters=sapm_module_params)
diff --git a/pvlib/test/test_scaling.py b/pvlib/test/test_scaling.py
new file mode 100644
index 0000000..62f3779
--- /dev/null
+++ b/pvlib/test/test_scaling.py
@@ -0,0 +1,146 @@
+import numpy as np
+import pandas as pd
+
+import pytest
+from numpy.testing import assert_almost_equal
+
+from pvlib import scaling
+from conftest import requires_scipy
+
+
+# Sample cloud speed
+cloud_speed = 5
+
+# Sample dt
+dt = 1
+
+
+@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]
+
+
+@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.]])
+
+
+@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():
+    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_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, positions[1], decimal=1)
+
+
+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(coordinates, positions):
+    pos = scaling.latlon_to_xy(coordinates.tolist())
+    assert_almost_equal(pos, positions, decimal=1)
+
+
+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(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)
+    assert_almost_equal(tmscale, expect_tmscale)
+    assert_almost_equal(wavelet[0:3, 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)
+
+
+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(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(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(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(clear_sky_index, positions):
+    with pytest.raises(ValueError):
+        scaling.wvm(clear_sky_index, positions, cloud_speed)