From 80c4d4b0d5b15d0ea8f4242a01a75a0723d39b45 Mon Sep 17 00:00:00 2001 From: wdwettin Date: Tue, 20 Jul 2021 10:40:00 +0200 Subject: [PATCH 1/7] Implemented linear blending Added linear blending function, example file and test --- doc/source/pysteps_reference/nowcasts.rst | 1 + examples/plot_linear_blending.py | 141 ++++++++++++++++++ pysteps/nowcasts/interface.py | 6 +- pysteps/nowcasts/lagrangian_probability.py | 2 +- pysteps/nowcasts/linear_blending.py | 139 +++++++++++++++++ .../tests/test_nowcasts_linear_blending.py | 107 +++++++++++++ 6 files changed, 394 insertions(+), 2 deletions(-) create mode 100644 examples/plot_linear_blending.py create mode 100644 pysteps/nowcasts/linear_blending.py create mode 100644 pysteps/tests/test_nowcasts_linear_blending.py diff --git a/doc/source/pysteps_reference/nowcasts.rst b/doc/source/pysteps_reference/nowcasts.rst index 9747d1658..12bc406c5 100644 --- a/doc/source/pysteps_reference/nowcasts.rst +++ b/doc/source/pysteps_reference/nowcasts.rst @@ -13,3 +13,4 @@ Implementation of deterministic and ensemble nowcasting methods. .. automodule:: pysteps.nowcasts.sseps .. automodule:: pysteps.nowcasts.steps .. automodule:: pysteps.nowcasts.utils +.. automodule:: pysteps.nowcasts.linear_blending diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py new file mode 100644 index 000000000..ccaa9a872 --- /dev/null +++ b/examples/plot_linear_blending.py @@ -0,0 +1,141 @@ +# -*- coding: utf-8 -*- + +""" +Plot linear blending +==================== +""" + +from matplotlib import cm, pyplot as plt +import numpy as np +from pprint import pprint +from pysteps.motion.lucaskanade import dense_lucaskanade +from pysteps import io, rcparams, nowcasts +from pysteps.utils import conversion, transformation +from pysteps.visualization import plot_precip_field +from datetime import datetime +from pysteps.utils import dimension +from pysteps.nowcasts.linear_blending import forecast + + +def gaussian(x, max, mean, sigma): + return max * np.exp(-(x - mean) * (x - mean) / sigma / sigma / 2) + + +def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): + """Generates dummy NWP data with the same dimension as the input + precipitation field R. The NWP data is a vertical line with a Gaussian + profile moving to the left""" + + # R is original radar image + rows = R.shape[0] + cols = R.shape[1] + + # Initialise the dummy NWP data + R_nwp = np.zeros((n_leadtimes, rows, cols)) + x = np.linspace(-5, 5, cols) + + for n in range(n_leadtimes): + for i in range(rows): + R_nwp[n, i, :] = gaussian(x, max, mean, sigma) + mean -= speed / rows + + return R_nwp + + +############################################################################### +# Set nowcast parameters +n_ens_members = 20 +n_leadtimes = 12 +start_blending = 20 # in minutes +end_blending = 40 # in minutes +seed = 24 + +# Read precipitation field +# ------------------------ +# +# First thing, the sequence of Swiss radar composites is imported, converted and +# transformed into units of dBR. + +date = datetime.strptime("201701311200", "%Y%m%d%H%M") +data_source = "mch" + +# Load data source config +root_path = rcparams.data_sources[data_source]["root_path"] +path_fmt = rcparams.data_sources[data_source]["path_fmt"] +fn_pattern = rcparams.data_sources[data_source]["fn_pattern"] +fn_ext = rcparams.data_sources[data_source]["fn_ext"] +importer_name = rcparams.data_sources[data_source]["importer"] +importer_kwargs = rcparams.data_sources[data_source]["importer_kwargs"] +timestep = rcparams.data_sources[data_source]["timestep"] + +# Find the radar files in the archive +fns = io.find_by_date( + date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2 +) + +# Read the data from the archive +importer = io.get_method(importer_name, "importer") +R, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) + +# Convert to rain rate +R, metadata = conversion.to_rainrate(R, metadata) + +# Upscale data to 2 km to limit memory usage +R, metadata = dimension.aggregate_fields_space(R, metadata, 2000) + +# Import the dummy NWP data (vertical line moving to the left) +R_nwp = dummy_nwp(R[-1, :, :], n_leadtimes + 1, max=7, mean=4, speed=0.2 * 350) + +# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h, +# set the fill value to -15 dBR +R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0) + +# Set missing values with the fill value +R[~np.isfinite(R)] = -15.0 + +# Nicely print the metadata +pprint(metadata) + +############################################################################### +# Linear blending of nowcast and NWP data + +# Estimate the motion field +V = dense_lucaskanade(R) + +# Define nowcast keyword arguments +nowcast_kwargs = { + "n_ens_members": n_ens_members, + "n_cascade_levels": 6, + "R_thr": -10.0, + "kmperpixel": 2.0, + "timestep": timestep, + "noise_method": "nonparametric", + "vel_pert_method": "bps", + "mask_method": "incremental", +} + +# Calculate the blended precipitation field +R_blended = forecast( + R[-3:, :, :], + V, + n_leadtimes, + timestep, + "steps", + R_nwp=R_nwp[1:, :, :], + start_blending=start_blending, + end_blending=end_blending, + nowcast_kwargs=nowcast_kwargs, +) + +# Calculate the ensemble average +R_blended_mean = np.mean(R_blended[:, :, :, :], axis=0) + +# Plot the blended field +for i in range(n_leadtimes): + plot_precip_field( + R_blended_mean[i, :, :], + geodata=metadata, + title="Blended field (+ %i min)" % ((i + 1) * timestep), + ) + plt.show() + plt.close() diff --git a/pysteps/nowcasts/interface.py b/pysteps/nowcasts/interface.py index 7daa4e7ad..0ca68ae62 100644 --- a/pysteps/nowcasts/interface.py +++ b/pysteps/nowcasts/interface.py @@ -32,7 +32,7 @@ """ from pysteps.extrapolation.interface import eulerian_persistence -from pysteps.nowcasts import anvil, sprog, steps, sseps, extrapolation +from pysteps.nowcasts import anvil, sprog, steps, sseps, extrapolation, linear_blending from pysteps.nowcasts import lagrangian_probability _nowcast_methods = dict() @@ -45,6 +45,7 @@ _nowcast_methods["sprog"] = sprog.forecast _nowcast_methods["sseps"] = sseps.forecast _nowcast_methods["steps"] = steps.forecast +_nowcast_methods["linear blending"] = linear_blending.forecast def get_method(name): @@ -80,6 +81,9 @@ def get_method(name): | sseps | short-space ensemble prediction system (SSEPS). | | | Essentially, this is a localization of STEPS. | +-----------------+-------------------------------------------------------+ + | linear | the linear blending of a nowcast method with other | + | blending | data (e.g. NWP data). | + +-----------------+-------------------------------------------------------+ """ if isinstance(name, str): name = name.lower() diff --git a/pysteps/nowcasts/lagrangian_probability.py b/pysteps/nowcasts/lagrangian_probability.py index 48e04bb98..90c93586b 100644 --- a/pysteps/nowcasts/lagrangian_probability.py +++ b/pysteps/nowcasts/lagrangian_probability.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """ pysteps.nowcasts.lagrangian_probability -====================================== +======================================= Implementation of the local Lagrangian probability nowcasting technique described in :cite:`GZ2004`. diff --git a/pysteps/nowcasts/linear_blending.py b/pysteps/nowcasts/linear_blending.py new file mode 100644 index 000000000..119e4287b --- /dev/null +++ b/pysteps/nowcasts/linear_blending.py @@ -0,0 +1,139 @@ +# -*- coding: utf-8 -*- +""" +pysteps.nowcasts.linear_blending +================================ + +Implementation of the linear blending between nowcast and NWP data. + +.. autosummary:: + :toctree: ../generated/ + + forecast +""" + +import numpy as np +from pysteps import nowcasts +from pysteps.utils import transformation + + +def forecast( + precip, + velocity, + timesteps, + timestep, + nowcast_method, + R_nwp=None, + start_blending=120, + end_blending=240, + nowcast_kwargs=dict(), +): + + """Generate a forecast by linearly blending nowcasts with NWP data + + Parameters + ---------- + precip: array_like + Array containing the input precipitation field(s) ordered by timestamp + from oldest to newest. The time steps between the inputs are assumed + to be regular. + velocity; array_like + Array of shape (2, m, n) containing the x- and y-components of the advection + field. The velocities are assumed to represent one time step between the + inputs. All values are required to be finite. + timesteps: int + Number of time steps to forecast. + timestep: int or float + The time difference (in minutes) between consecutive forecast fields. + nowcast_method: {'anvil', 'eulerian', 'lagrangian', 'extrapolation', 'lagrangian_probability', 'sprog', 'steps', 'sseps'} + Name of the nowcast method used to forecast the precipitation. + R_nwp: array_like or NoneType, optional + Array of shape (timesteps, m, n) containing the NWP precipitation fields + ordered by timestamp from oldest to newest. The time steps between the + inputs are assumed to be regular (and identical to the time step between + the nowcasts). If no NWP data is given the value of R_nwp is None and no + blending will be performed. + start_blending: int, optional + Time stamp (in minutes) after which the blending should start. Before this + only the nowcast data is used. + end_blending: int, optional + Time stamp (in minutes) after which the blending should end. Between + start_blending and end_blending the nowcasts and NWP data are linearly + merged with each other. After end_blending only the NWP data is used. + nowcast_kwargs: dict, optional + Dictionary containing keyword arguments for the nowcast method. + + Returns + ------- + R_blended: ndarray + Array of shape (timesteps, m, n) in the case of no ensemble method or + of shape (n_ens_members, timesteps, m, n) in the case of an ensemble + method containing the precipation forecast generated by linearly blending + the nowcasts and the NWP data. + """ + + # Calculate the nowcasts + nowcast_method_func = nowcasts.get_method(nowcast_method) + R_nowcast = nowcast_method_func( + precip, + velocity, + timesteps, + **nowcast_kwargs, + ) + + # Transform the precipitation back to mm/h + R_nowcast = transformation.dB_transform(R_nowcast, threshold=-10.0, inverse=True)[0] + + # Check if NWP data is given as input + if R_nwp is not None: + + # Initialise output + R_blended = np.zeros_like(R_nowcast) + + # Check if the nowcast method is an ensemble method, in order to + # determine the shape of the NWP data + if nowcast_method in ["steps", "sseps"]: + n_ens_members = R_nowcast.shape[0] + + R_nwp = np.repeat(R_nwp[np.newaxis, :, :], n_ens_members, axis=0) + + # Check if dimensions are correct + assert ( + R_nwp.shape == R_nowcast.shape + ), "The dimensions of R_nowcast and R_nwp need to be identical: dimension of R_nwp = {} and dimension of R_nowcast = {}".format( + R_nwp.shape, R_nowcast.shape + ) + + # Calculate the weights + for i in range(timesteps): + # Calculate what time we are at + t = (i + 1) * timestep + + # Calculate the weight with a linear relation (weight_nwp at start_blending = 0.0) + # and (weight_nwp at end_blending = 1.0) + weight_nwp = (t - start_blending) / (end_blending - start_blending) + + # Set weights at times before start_blending and after end_blending + if weight_nwp < 0.0: + weight_nwp = 0.0 + elif weight_nwp > 1.0: + weight_nwp = 1.0 + + # Calculate weight_nowcast + weight_nowcast = 1.0 - weight_nwp + + # Calculate output by combining R_nwp and R_nowcast, + # while distinguishing between ensemble and non-ensemble methods + if nowcast_method in ["steps", "sseps"]: + R_blended[:, i, :, :] = ( + weight_nwp * R_nwp[:, i, :, :] + + weight_nowcast * R_nowcast[:, i, :, :] + ) + else: + R_blended[i, :, :] = ( + weight_nwp * R_nwp[i, :, :] + weight_nowcast * R_nowcast[i, :, :] + ) + else: + # If no NWP data is given, the blended field is simply equal to the nowcast field + R_blended = R_nowcast + + return R_blended diff --git a/pysteps/tests/test_nowcasts_linear_blending.py b/pysteps/tests/test_nowcasts_linear_blending.py new file mode 100644 index 000000000..781c03d7c --- /dev/null +++ b/pysteps/tests/test_nowcasts_linear_blending.py @@ -0,0 +1,107 @@ +# -*- coding: utf-8 -*- + +import numpy as np +import pytest +from pysteps.nowcasts.linear_blending import forecast +from numpy.testing import assert_array_almost_equal +from pysteps.utils import conversion, transformation + +# Test function arguments +linear_arg_values = [ + (5, 30, 60, 20, 45, "eulerian", None), + (4, 23, 33, 9, 28, "eulerian", None), + (3, 18, 36, 13, 27, "eulerian", None), + (7, 30, 68, 11, 49, "eulerian", None), + (10, 100, 160, 25, 130, "eulerian", None), + (6, 60, 180, 22, 120, "eulerian", None), + (5, 100, 200, 40, 150, "eulerian", None), + (5, 30, 60, 20, 45, "extrapolation", np.zeros((2, 200, 200))), + (4, 23, 33, 9, 28, "extrapolation", np.zeros((2, 200, 200))), + (3, 18, 36, 13, 27, "extrapolation", np.zeros((2, 200, 200))), + (7, 30, 68, 11, 49, "extrapolation", np.zeros((2, 200, 200))), + (10, 100, 160, 25, 130, "extrapolation", np.zeros((2, 200, 200))), + (6, 60, 180, 22, 120, "extrapolation", np.zeros((2, 200, 200))), + (5, 100, 200, 40, 150, "extrapolation", np.zeros((2, 200, 200))), +] + + +@pytest.mark.parametrize( + "timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V", + linear_arg_values, +) +def test_linear_blending( + timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V +): + """Tests if the linear blending function is correct. For the nowcast data a precipitation field + which is constant over time is taken where one half of the field has no rain and the other half + has a set value. For the NWP data a similar field is taken with the only difference + being that now the other half of the field is zero. After blending the blended field should + have a constant value over the entire field at the timestep right in the middle between the start + of the blending and the end of the blending. This assertion is checked to see if the + linear blending function works well.""" + + # The argument controltime gives the timestep at which the field is assumed to be + # entirely constant + + # Assert that the control time step is in the range of the forecasted time steps + assert controltime <= ( + n_timesteps * timestep + ), "Control time needs to be within reach of forecasts, controltime = {} and n_timesteps = {}".format( + controltime, n_timesteps + ) + + # Assert that the start time of the blending comes before the end time of the blending + assert ( + start_blending < end_blending + ), "Start time of blending needs to be smaller than end time of blending" + + # Assert that the control time is a multiple of the time step + assert ( + not controltime % timestep + ), "Control time needs to be a multiple of the time step" + + # Initialise dummy NWP data + R_NWP = np.zeros((n_timesteps, 200, 200)) + + for i in range(100): + R_NWP[:, i, :] = 11.0 + + # Define nowcast input data + R_input = np.zeros((200, 200)) + + for i in range(100, 200): + R_input[i, :] = 11.0 + + # Transform from mm/h to dB + R_input, _ = transformation.dB_transform( + R_input, None, threshold=0.1, zerovalue=-15.0 + ) + + # Calculate the blended field + R_blended = forecast( + R_input, + V, + n_timesteps, + timestep, + nowcast_method, + R_NWP, + start_blending=start_blending, + end_blending=end_blending, + ) + + # Assert that the blended field has the expected dimension + assert R_blended.shape == ( + n_timesteps, + 200, + 200, + ), "The shape of the blended array does not have the expected value. The shape is {}".format( + R_blended.shape + ) + + # Assert that the blended field at the control time step is equal to + # a constant field with the expected value. + assert_array_almost_equal( + R_blended[controltime // timestep - 1], + np.ones((200, 200)) * 5.5, + err_msg="The blended array does not have the expected value", + ) From f01f4f2c97b4ffcb323a9f52cd0dffe134bdd9ee Mon Sep 17 00:00:00 2001 From: Wout Dewettinck Date: Thu, 5 Aug 2021 15:21:03 +0200 Subject: [PATCH 2/7] Replaced missing data in nowcasting with NWP data The replacement is done at the end --- examples/plot_linear_blending.py | 10 +++++----- pysteps/nowcasts/linear_blending.py | 4 ++++ 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py index ccaa9a872..6a7ed095e 100644 --- a/examples/plot_linear_blending.py +++ b/examples/plot_linear_blending.py @@ -45,9 +45,9 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): ############################################################################### # Set nowcast parameters n_ens_members = 20 -n_leadtimes = 12 -start_blending = 20 # in minutes -end_blending = 40 # in minutes +n_leadtimes = 18 +start_blending = 30 # in minutes +end_blending = 60 # in minutes seed = 24 # Read precipitation field @@ -75,7 +75,7 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): # Read the data from the archive importer = io.get_method(importer_name, "importer") -R, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to rain rate R, metadata = conversion.to_rainrate(R, metadata) @@ -91,7 +91,7 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0) # Set missing values with the fill value -R[~np.isfinite(R)] = -15.0 +# R[~np.isfinite(R)] = -15.0 # Nicely print the metadata pprint(metadata) diff --git a/pysteps/nowcasts/linear_blending.py b/pysteps/nowcasts/linear_blending.py index 119e4287b..6f6567b70 100644 --- a/pysteps/nowcasts/linear_blending.py +++ b/pysteps/nowcasts/linear_blending.py @@ -132,6 +132,10 @@ def forecast( R_blended[i, :, :] = ( weight_nwp * R_nwp[i, :, :] + weight_nowcast * R_nowcast[i, :, :] ) + + # Find where the NaN values are and replace them with NWP data + nan_indices = np.isnan(R_blended) + R_blended[nan_indices] = R_nwp[nan_indices] else: # If no NWP data is given, the blended field is simply equal to the nowcast field R_blended = R_nowcast From 7d6008b222e2582438f0398746024f7d54a7a57a Mon Sep 17 00:00:00 2001 From: Wout Dewettinck Date: Thu, 5 Aug 2021 18:12:58 +0200 Subject: [PATCH 3/7] Added compatibility for NWP ensembles --- examples/plot_linear_blending.py | 25 +++++- pysteps/nowcasts/linear_blending.py | 81 +++++++++++++------ .../tests/test_nowcasts_linear_blending.py | 8 +- 3 files changed, 82 insertions(+), 32 deletions(-) diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py index 6a7ed095e..40986d234 100644 --- a/examples/plot_linear_blending.py +++ b/examples/plot_linear_blending.py @@ -44,7 +44,7 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): ############################################################################### # Set nowcast parameters -n_ens_members = 20 +n_ens_members = 10 n_leadtimes = 18 start_blending = 30 # in minutes end_blending = 60 # in minutes @@ -127,8 +127,29 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): nowcast_kwargs=nowcast_kwargs, ) +""" +extrap_kwargs = {"allow_nonfinite_values": True} +nowcast_kwargs = {"extrap_kwargs": extrap_kwargs} + +# Calculate the blended precipitation field +R_blended = forecast( + R[-1, :, :], + V, + n_leadtimes, + timestep, + "extrapolation", + R_nwp=R_nwp[1: :, :], + start_blending=start_blending, + end_blending=end_blending, + nowcast_kwargs=nowcast_kwargs, +) +""" + # Calculate the ensemble average -R_blended_mean = np.mean(R_blended[:, :, :, :], axis=0) +if len(R_blended.shape) == 4: + R_blended_mean = np.mean(R_blended[:, :, :, :], axis=0) +else: + R_blended_mean = np.copy(R_blended) # Plot the blended field for i in range(n_leadtimes): diff --git a/pysteps/nowcasts/linear_blending.py b/pysteps/nowcasts/linear_blending.py index 6f6567b70..35ab0ce19 100644 --- a/pysteps/nowcasts/linear_blending.py +++ b/pysteps/nowcasts/linear_blending.py @@ -47,11 +47,12 @@ def forecast( nowcast_method: {'anvil', 'eulerian', 'lagrangian', 'extrapolation', 'lagrangian_probability', 'sprog', 'steps', 'sseps'} Name of the nowcast method used to forecast the precipitation. R_nwp: array_like or NoneType, optional - Array of shape (timesteps, m, n) containing the NWP precipitation fields - ordered by timestamp from oldest to newest. The time steps between the - inputs are assumed to be regular (and identical to the time step between - the nowcasts). If no NWP data is given the value of R_nwp is None and no - blending will be performed. + Array of shape (timesteps, m, n) in the case of no ensemble or + of shape (n_ens_members, timesteps, m, n) in the case of an ensemble + containing the NWP precipitation fields ordered by timestamp from oldest + to newest. The time steps between the inputs are assumed to be regular + (and identical to the time step between the nowcasts). If no NWP + data is given the value of R_nwp is None and no blending will be performed. start_blending: int, optional Time stamp (in minutes) after which the blending should start. Before this only the nowcast data is used. @@ -65,9 +66,9 @@ def forecast( Returns ------- R_blended: ndarray - Array of shape (timesteps, m, n) in the case of no ensemble method or + Array of shape (timesteps, m, n) in the case of no ensemble or of shape (n_ens_members, timesteps, m, n) in the case of an ensemble - method containing the precipation forecast generated by linearly blending + containing the precipation forecast generated by linearly blending the nowcasts and the NWP data. """ @@ -86,23 +87,51 @@ def forecast( # Check if NWP data is given as input if R_nwp is not None: + if len(R_nowcast.shape) == 4: + n_ens_members_nowcast = R_nowcast.shape[0] + if n_ens_members_nowcast == 1: + R_nowcast = np.squeeze(R_nowcast) + else: + n_ens_members_nowcast = 1 + + if len(R_nwp.shape) == 4: + n_ens_members_nwp = R_nwp.shape[0] + if n_ens_members_nwp == 1: + R_nwp = np.squeeze(R_nwp) + else: + n_ens_members_nwp = 1 + + n_ens_members_max = max(n_ens_members_nowcast, n_ens_members_nwp) + n_ens_members_min = min(n_ens_members_nowcast, n_ens_members_nwp) + + if n_ens_members_min != n_ens_members_max: + if n_ens_members_nwp == 1: + R_nwp = np.repeat(R_nwp[np.newaxis, :, :], n_ens_members_max, axis=0) + elif n_ens_members_nowcast == 1: + R_nowcast = np.repeat( + R_nowcast[np.newaxis, :, :], n_ens_members_max, axis=0 + ) + else: + repeats = [ + (n_ens_members_max + i) // n_ens_members_min + for i in range(n_ens_members_min) + ] + + if n_ens_members_nwp == n_ens_members_min: + R_nwp = np.repeat(R_nwp, repeats, axis=0) + elif n_ens_members_nowcast == n_ens_members_min: + R_nowcast = np.repeat(R_nowcast, repeats, axis=0) + + # Check if dimensions are correct + assert ( + R_nwp.shape == R_nowcast.shape + ), "The dimensions of R_nowcast and R_nwp need to be identical: dimension of R_nwp = {} and dimension of R_nowcast = {}".format( + R_nwp.shape, R_nowcast.shape + ) + # Initialise output R_blended = np.zeros_like(R_nowcast) - # Check if the nowcast method is an ensemble method, in order to - # determine the shape of the NWP data - if nowcast_method in ["steps", "sseps"]: - n_ens_members = R_nowcast.shape[0] - - R_nwp = np.repeat(R_nwp[np.newaxis, :, :], n_ens_members, axis=0) - - # Check if dimensions are correct - assert ( - R_nwp.shape == R_nowcast.shape - ), "The dimensions of R_nowcast and R_nwp need to be identical: dimension of R_nwp = {} and dimension of R_nowcast = {}".format( - R_nwp.shape, R_nowcast.shape - ) - # Calculate the weights for i in range(timesteps): # Calculate what time we are at @@ -123,15 +152,15 @@ def forecast( # Calculate output by combining R_nwp and R_nowcast, # while distinguishing between ensemble and non-ensemble methods - if nowcast_method in ["steps", "sseps"]: + if n_ens_members_max == 1: + R_blended[i, :, :] = ( + weight_nwp * R_nwp[i, :, :] + weight_nowcast * R_nowcast[i, :, :] + ) + else: R_blended[:, i, :, :] = ( weight_nwp * R_nwp[:, i, :, :] + weight_nowcast * R_nowcast[:, i, :, :] ) - else: - R_blended[i, :, :] = ( - weight_nwp * R_nwp[i, :, :] + weight_nowcast * R_nowcast[i, :, :] - ) # Find where the NaN values are and replace them with NWP data nan_indices = np.isnan(R_blended) diff --git a/pysteps/tests/test_nowcasts_linear_blending.py b/pysteps/tests/test_nowcasts_linear_blending.py index 781c03d7c..ca6053f30 100644 --- a/pysteps/tests/test_nowcasts_linear_blending.py +++ b/pysteps/tests/test_nowcasts_linear_blending.py @@ -33,10 +33,10 @@ def test_linear_blending( timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V ): """Tests if the linear blending function is correct. For the nowcast data a precipitation field - which is constant over time is taken where one half of the field has no rain and the other half - has a set value. For the NWP data a similar field is taken with the only difference - being that now the other half of the field is zero. After blending the blended field should - have a constant value over the entire field at the timestep right in the middle between the start + which is constant over time is taken. One half of the field has no rain and the other half + has a set value. For the NWP data a similar field is taken, the only difference + being that now the other half of the field is zero. The blended field should have a + constant value over the entire field at the timestep right in the middle between the start of the blending and the end of the blending. This assertion is checked to see if the linear blending function works well.""" From 086c41dfec1371498f503b54a2c2ccfd9429cbb9 Mon Sep 17 00:00:00 2001 From: Wout Dewettinck Date: Tue, 24 Aug 2021 20:25:06 +0200 Subject: [PATCH 4/7] Add fill_nwp as parameter Add use_nwp as parameter Add commentary about use_nwp Add fill_nwp as parameter --- pysteps/nowcasts/linear_blending.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pysteps/nowcasts/linear_blending.py b/pysteps/nowcasts/linear_blending.py index 35ab0ce19..8057baacc 100644 --- a/pysteps/nowcasts/linear_blending.py +++ b/pysteps/nowcasts/linear_blending.py @@ -25,6 +25,7 @@ def forecast( R_nwp=None, start_blending=120, end_blending=240, + fill_nwp=True, nowcast_kwargs=dict(), ): @@ -60,6 +61,9 @@ def forecast( Time stamp (in minutes) after which the blending should end. Between start_blending and end_blending the nowcasts and NWP data are linearly merged with each other. After end_blending only the NWP data is used. + fill_nwp: bool, optional + Standard value is True. If True, the NWP data will be used to fill in the + no data mask of the nowcast. nowcast_kwargs: dict, optional Dictionary containing keyword arguments for the nowcast method. @@ -163,8 +167,9 @@ def forecast( ) # Find where the NaN values are and replace them with NWP data - nan_indices = np.isnan(R_blended) - R_blended[nan_indices] = R_nwp[nan_indices] + if fill_nwp: + nan_indices = np.isnan(R_blended) + R_blended[nan_indices] = R_nwp[nan_indices] else: # If no NWP data is given, the blended field is simply equal to the nowcast field R_blended = R_nowcast From d9d07772911f2295df4a4097042ebe96c6ce72e5 Mon Sep 17 00:00:00 2001 From: RubenImhoff Date: Tue, 19 Oct 2021 13:30:47 +0200 Subject: [PATCH 5/7] snake_case implementation and move linear blending to blending folder --- .../{nowcasts => blending}/linear_blending.py | 61 ++++++++++--------- 1 file changed, 33 insertions(+), 28 deletions(-) rename pysteps/{nowcasts => blending}/linear_blending.py (73%) diff --git a/pysteps/nowcasts/linear_blending.py b/pysteps/blending/linear_blending.py similarity index 73% rename from pysteps/nowcasts/linear_blending.py rename to pysteps/blending/linear_blending.py index 8057baacc..e7297f1b6 100644 --- a/pysteps/nowcasts/linear_blending.py +++ b/pysteps/blending/linear_blending.py @@ -22,7 +22,7 @@ def forecast( timesteps, timestep, nowcast_method, - R_nwp=None, + precip_nwp=None, start_blending=120, end_blending=240, fill_nwp=True, @@ -47,13 +47,13 @@ def forecast( The time difference (in minutes) between consecutive forecast fields. nowcast_method: {'anvil', 'eulerian', 'lagrangian', 'extrapolation', 'lagrangian_probability', 'sprog', 'steps', 'sseps'} Name of the nowcast method used to forecast the precipitation. - R_nwp: array_like or NoneType, optional + precip_nwp: array_like or NoneType, optional Array of shape (timesteps, m, n) in the case of no ensemble or of shape (n_ens_members, timesteps, m, n) in the case of an ensemble containing the NWP precipitation fields ordered by timestamp from oldest to newest. The time steps between the inputs are assumed to be regular (and identical to the time step between the nowcasts). If no NWP - data is given the value of R_nwp is None and no blending will be performed. + data is given the value of precip_nwp is None and no blending will be performed. start_blending: int, optional Time stamp (in minutes) after which the blending should start. Before this only the nowcast data is used. @@ -78,30 +78,32 @@ def forecast( # Calculate the nowcasts nowcast_method_func = nowcasts.get_method(nowcast_method) - R_nowcast = nowcast_method_func( + precip_nowcast = nowcast_method_func( precip, velocity, timesteps, **nowcast_kwargs, ) - # Transform the precipitation back to mm/h - R_nowcast = transformation.dB_transform(R_nowcast, threshold=-10.0, inverse=True)[0] + # Make sure that precip_nowcast and precip_nwp are in mm/h + precip_nowcast = transformation.dB_transform( + precip_nowcast, threshold=-10.0, inverse=True + )[0] # Check if NWP data is given as input - if R_nwp is not None: + if precip_nwp is not None: - if len(R_nowcast.shape) == 4: - n_ens_members_nowcast = R_nowcast.shape[0] + if len(precip_nowcast.shape) == 4: + n_ens_members_nowcast = precip_nowcast.shape[0] if n_ens_members_nowcast == 1: - R_nowcast = np.squeeze(R_nowcast) + precip_nowcast = np.squeeze(precip_nowcast) else: n_ens_members_nowcast = 1 - if len(R_nwp.shape) == 4: - n_ens_members_nwp = R_nwp.shape[0] + if len(precip_nwp.shape) == 4: + n_ens_members_nwp = precip_nwp.shape[0] if n_ens_members_nwp == 1: - R_nwp = np.squeeze(R_nwp) + precip_nwp = np.squeeze(precip_nwp) else: n_ens_members_nwp = 1 @@ -110,10 +112,12 @@ def forecast( if n_ens_members_min != n_ens_members_max: if n_ens_members_nwp == 1: - R_nwp = np.repeat(R_nwp[np.newaxis, :, :], n_ens_members_max, axis=0) + precip_nwp = np.repeat( + precip_nwp[np.newaxis, :, :], n_ens_members_max, axis=0 + ) elif n_ens_members_nowcast == 1: - R_nowcast = np.repeat( - R_nowcast[np.newaxis, :, :], n_ens_members_max, axis=0 + precip_nowcast = np.repeat( + precip_nowcast[np.newaxis, :, :], n_ens_members_max, axis=0 ) else: repeats = [ @@ -122,19 +126,19 @@ def forecast( ] if n_ens_members_nwp == n_ens_members_min: - R_nwp = np.repeat(R_nwp, repeats, axis=0) + precip_nwp = np.repeat(precip_nwp, repeats, axis=0) elif n_ens_members_nowcast == n_ens_members_min: - R_nowcast = np.repeat(R_nowcast, repeats, axis=0) + precip_nowcast = np.repeat(precip_nowcast, repeats, axis=0) # Check if dimensions are correct assert ( - R_nwp.shape == R_nowcast.shape - ), "The dimensions of R_nowcast and R_nwp need to be identical: dimension of R_nwp = {} and dimension of R_nowcast = {}".format( - R_nwp.shape, R_nowcast.shape + precip_nwp.shape == precip_nowcast.shape + ), "The dimensions of precip_nowcast and precip_nwp need to be identical: dimension of precip_nwp = {} and dimension of precip_nowcast = {}".format( + precip_nwp.shape, precip_nowcast.shape ) # Initialise output - R_blended = np.zeros_like(R_nowcast) + R_blended = np.zeros_like(precip_nowcast) # Calculate the weights for i in range(timesteps): @@ -154,24 +158,25 @@ def forecast( # Calculate weight_nowcast weight_nowcast = 1.0 - weight_nwp - # Calculate output by combining R_nwp and R_nowcast, + # Calculate output by combining precip_nwp and precip_nowcast, # while distinguishing between ensemble and non-ensemble methods if n_ens_members_max == 1: R_blended[i, :, :] = ( - weight_nwp * R_nwp[i, :, :] + weight_nowcast * R_nowcast[i, :, :] + weight_nwp * precip_nwp[i, :, :] + + weight_nowcast * precip_nowcast[i, :, :] ) else: R_blended[:, i, :, :] = ( - weight_nwp * R_nwp[:, i, :, :] - + weight_nowcast * R_nowcast[:, i, :, :] + weight_nwp * precip_nwp[:, i, :, :] + + weight_nowcast * precip_nowcast[:, i, :, :] ) # Find where the NaN values are and replace them with NWP data if fill_nwp: nan_indices = np.isnan(R_blended) - R_blended[nan_indices] = R_nwp[nan_indices] + R_blended[nan_indices] = precip_nwp[nan_indices] else: # If no NWP data is given, the blended field is simply equal to the nowcast field - R_blended = R_nowcast + R_blended = precip_nowcast return R_blended From 6d4a37517d0acd84fd4f5232a1a318947e21689c Mon Sep 17 00:00:00 2001 From: RubenImhoff Date: Tue, 19 Oct 2021 14:55:00 +0200 Subject: [PATCH 6/7] Changes after review by Carlos, Daniele and Ruben --- doc/source/pysteps_reference/blending.rst | 8 ++ doc/source/pysteps_reference/nowcasts.rst | 1 - examples/plot_linear_blending.py | 100 +++++++++++------- pysteps/blending/__init__.py | 4 + pysteps/blending/interface.py | 48 +++++++++ pysteps/blending/linear_blending.py | 28 ++++- pysteps/nowcasts/interface.py | 5 - ...ng.py => test_blending_linear_blending.py} | 30 +++--- 8 files changed, 163 insertions(+), 61 deletions(-) create mode 100644 doc/source/pysteps_reference/blending.rst create mode 100644 pysteps/blending/__init__.py create mode 100644 pysteps/blending/interface.py rename pysteps/tests/{test_nowcasts_linear_blending.py => test_blending_linear_blending.py} (85%) diff --git a/doc/source/pysteps_reference/blending.rst b/doc/source/pysteps_reference/blending.rst new file mode 100644 index 000000000..b827e1aa7 --- /dev/null +++ b/doc/source/pysteps_reference/blending.rst @@ -0,0 +1,8 @@ +================ +pysteps.blending +================ + +Implementation of blending methods for blending (ensemble) nowcasts with Numerical Weather Prediction (NWP) models. + + +.. automodule:: pysteps.blending.linear_blending diff --git a/doc/source/pysteps_reference/nowcasts.rst b/doc/source/pysteps_reference/nowcasts.rst index 1487c4524..35d95cbbe 100644 --- a/doc/source/pysteps_reference/nowcasts.rst +++ b/doc/source/pysteps_reference/nowcasts.rst @@ -14,4 +14,3 @@ Implementation of deterministic and ensemble nowcasting methods. .. automodule:: pysteps.nowcasts.sseps .. automodule:: pysteps.nowcasts.steps .. automodule:: pysteps.nowcasts.utils -.. automodule:: pysteps.nowcasts.linear_blending diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py index 40986d234..dada60256 100644 --- a/examples/plot_linear_blending.py +++ b/examples/plot_linear_blending.py @@ -3,43 +3,48 @@ """ Plot linear blending ==================== + +This tutorial shows how to construct a simple linear blending between an ensemble +nowcast and a dummy Numerical Weather Prediction (NWP) rainfall forecast. """ -from matplotlib import cm, pyplot as plt +# TODO: Make xarray ready for pysteps-v2 + +from matplotlib import pyplot as plt import numpy as np from pprint import pprint +from datetime import datetime from pysteps.motion.lucaskanade import dense_lucaskanade -from pysteps import io, rcparams, nowcasts +from pysteps import io, rcparams from pysteps.utils import conversion, transformation from pysteps.visualization import plot_precip_field -from datetime import datetime from pysteps.utils import dimension -from pysteps.nowcasts.linear_blending import forecast +from pysteps.blending.linear_blending import forecast def gaussian(x, max, mean, sigma): return max * np.exp(-(x - mean) * (x - mean) / sigma / sigma / 2) -def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): +def dummy_nwp(precip, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): """Generates dummy NWP data with the same dimension as the input - precipitation field R. The NWP data is a vertical line with a Gaussian + precipitation field precip. The NWP data is a vertical line with a Gaussian profile moving to the left""" - # R is original radar image - rows = R.shape[0] - cols = R.shape[1] + # precip is original radar image + rows = precip.shape[0] + cols = precip.shape[1] # Initialise the dummy NWP data - R_nwp = np.zeros((n_leadtimes, rows, cols)) + precip_nwp = np.zeros((n_leadtimes, rows, cols)) x = np.linspace(-5, 5, cols) for n in range(n_leadtimes): for i in range(rows): - R_nwp[n, i, :] = gaussian(x, max, mean, sigma) + precip_nwp[n, i, :] = gaussian(x, max, mean, sigma) mean -= speed / rows - return R_nwp + return precip_nwp ############################################################################### @@ -75,23 +80,44 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): # Read the data from the archive importer = io.get_method(importer_name, "importer") -R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) +precip = io.read_timeseries(fns, importer, legacy=False, **importer_kwargs) + +# Get the metadata +metadata = precip.x.attrs.copy() +metadata.update(**precip.y.attrs) +metadata.update(**precip.attrs) # Convert to rain rate -R, metadata = conversion.to_rainrate(R, metadata) +precip, metadata = conversion.to_rainrate(precip, metadata) # Upscale data to 2 km to limit memory usage -R, metadata = dimension.aggregate_fields_space(R, metadata, 2000) +precip, metadata = dimension.aggregate_fields_space(precip, metadata, 2000) # Import the dummy NWP data (vertical line moving to the left) -R_nwp = dummy_nwp(R[-1, :, :], n_leadtimes + 1, max=7, mean=4, speed=0.2 * 350) +precip_nwp = dummy_nwp( + precip[-1, :, :], n_leadtimes + 1, max=7, mean=4, speed=0.2 * 350 +) +metadata_nwp = metadata.copy() + +# Plot the radar rainfall field and the first time step of the dummy NWP forecast. +date_str = datetime.strftime(date, "%Y-%m-%d %H:%M") +plt.figure(figsize=(10, 5)) +plt.subplot(121) +plot_precip_field( + precip[-1, :, :], geodata=metadata, title=f"Radar observation at {date_str}" +) +plt.subplot(122) +plot_precip_field( + precip_nwp[0, :, :], geodata=metadata_nwp, title=f"Dummy NWP forecast at {date_str}" +) +plt.tight_layout() +plt.show() # Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h, # set the fill value to -15 dBR -R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0) - -# Set missing values with the fill value -# R[~np.isfinite(R)] = -15.0 +precip, metadata = transformation.dB_transform( + precip, metadata, threshold=0.1, zerovalue=-15.0 +) # Nicely print the metadata pprint(metadata) @@ -100,7 +126,7 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): # Linear blending of nowcast and NWP data # Estimate the motion field -V = dense_lucaskanade(R) +velocity = dense_lucaskanade(precip) # Define nowcast keyword arguments nowcast_kwargs = { @@ -115,13 +141,15 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): } # Calculate the blended precipitation field -R_blended = forecast( - R[-3:, :, :], - V, - n_leadtimes, - timestep, - "steps", - R_nwp=R_nwp[1:, :, :], +precip_blended = forecast( + precip=precip[-3:, :, :], + precip_metadata=metadata, + velocity=velocity, + timesteps=n_leadtimes, + timestep=timestep, + nowcast_method="steps", + precip_nwp=precip_nwp[1:, :, :], + precip_nwp_metadata=metadata_nwp, start_blending=start_blending, end_blending=end_blending, nowcast_kwargs=nowcast_kwargs, @@ -132,13 +160,13 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): nowcast_kwargs = {"extrap_kwargs": extrap_kwargs} # Calculate the blended precipitation field -R_blended = forecast( - R[-1, :, :], +precip_blended = forecast( + precip[-1, :, :], V, n_leadtimes, timestep, "extrapolation", - R_nwp=R_nwp[1: :, :], + precip_nwp=precip_nwp[1: :, :], start_blending=start_blending, end_blending=end_blending, nowcast_kwargs=nowcast_kwargs, @@ -146,15 +174,15 @@ def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100): """ # Calculate the ensemble average -if len(R_blended.shape) == 4: - R_blended_mean = np.mean(R_blended[:, :, :, :], axis=0) +if len(precip_blended.shape) == 4: + precip_blended_mean = np.mean(precip_blended[:, :, :, :], axis=0) else: - R_blended_mean = np.copy(R_blended) + precip_blended_mean = np.copy(precip_blended) # Plot the blended field -for i in range(n_leadtimes): +for i in range(0, n_leadtimes, 3): plot_precip_field( - R_blended_mean[i, :, :], + precip_blended_mean[i, :, :], geodata=metadata, title="Blended field (+ %i min)" % ((i + 1) * timestep), ) diff --git a/pysteps/blending/__init__.py b/pysteps/blending/__init__.py new file mode 100644 index 000000000..420a11bd0 --- /dev/null +++ b/pysteps/blending/__init__.py @@ -0,0 +1,4 @@ +# -*- coding: utf-8 -*- +"""Methods for blending NWP model(s) with nowcasts.""" + +from pysteps.blending.interface import get_method diff --git a/pysteps/blending/interface.py b/pysteps/blending/interface.py new file mode 100644 index 000000000..6764e95ba --- /dev/null +++ b/pysteps/blending/interface.py @@ -0,0 +1,48 @@ +# -*- coding: utf-8 -*- +""" +pysteps.blending.interface +========================== +Interface for the blending module. It returns a callable function for computing +blended nowcasts with NWP models. +.. autosummary:: + :toctree: ../generated/ + get_method +""" + +from pysteps.blending import linear_blending + +_blending_methods = dict() +_blending_methods["linear_blending"] = linear_blending.forecast + + +def get_method(name): + """Return a callable function for computing nowcasts. + Description: + Return a callable function for computing deterministic or ensemble + precipitation nowcasts. + Implemented methods: + +-----------------+-------------------------------------------------------+ + | Name | Description | + +=================+=======================================================+ + +-----------------+-------------------------------------------------------+ + | linear | the linear blending of a nowcast method with other | + | blending | data (e.g. NWP data). | + +-----------------+-------------------------------------------------------+ + """ + if isinstance(name, str): + name = name.lower() + else: + raise TypeError( + "Only strings supported for the method's names.\n" + + "Available names:" + + str(list(_blending_methods.keys())) + ) from None + + try: + return _blending_methods[name] + except KeyError: + raise ValueError( + "Unknown blending method {}\n".format(name) + + "The available methods are:" + + str(list(_blending_methods.keys())) + ) from None diff --git a/pysteps/blending/linear_blending.py b/pysteps/blending/linear_blending.py index e7297f1b6..5b8ea5a21 100644 --- a/pysteps/blending/linear_blending.py +++ b/pysteps/blending/linear_blending.py @@ -13,16 +13,20 @@ import numpy as np from pysteps import nowcasts -from pysteps.utils import transformation +from pysteps.utils import conversion + +# TODO: Make xarray ready for pysteps-v2 def forecast( precip, + precip_metadata, velocity, timesteps, timestep, nowcast_method, precip_nwp=None, + precip_nwp_metadata=None, start_blending=120, end_blending=240, fill_nwp=True, @@ -37,6 +41,9 @@ def forecast( Array containing the input precipitation field(s) ordered by timestamp from oldest to newest. The time steps between the inputs are assumed to be regular. + precip_metadata: dict + Metadata dictionary containing (at least) the transform, unit and threshold + attributes as described in the documentation of :py:mod:`pysteps.io.importers`. velocity; array_like Array of shape (2, m, n) containing the x- and y-components of the advection field. The velocities are assumed to represent one time step between the @@ -54,6 +61,9 @@ def forecast( to newest. The time steps between the inputs are assumed to be regular (and identical to the time step between the nowcasts). If no NWP data is given the value of precip_nwp is None and no blending will be performed. + precip_nwp_metadata: dict or NoneType, optional + NWP metadata dictionary containing (at least) the transform, unit and threshold + attributes as described in the documentation of :py:mod:`pysteps.io.importers`. start_blending: int, optional Time stamp (in minutes) after which the blending should start. Before this only the nowcast data is used. @@ -73,7 +83,8 @@ def forecast( Array of shape (timesteps, m, n) in the case of no ensemble or of shape (n_ens_members, timesteps, m, n) in the case of an ensemble containing the precipation forecast generated by linearly blending - the nowcasts and the NWP data. + the nowcasts and the NWP data. n_ens_members equals the maximum no. of + ensemble members in either the nowcast or nwp model(s). """ # Calculate the nowcasts @@ -86,9 +97,8 @@ def forecast( ) # Make sure that precip_nowcast and precip_nwp are in mm/h - precip_nowcast = transformation.dB_transform( - precip_nowcast, threshold=-10.0, inverse=True - )[0] + precip_nowcast, _ = conversion.to_rainrate(precip_nowcast, metadata=precip_metadata) + precip_nwp, _ = conversion.to_rainrate(precip_nwp, metadata=precip_nwp_metadata) # Check if NWP data is given as input if precip_nwp is not None: @@ -107,6 +117,14 @@ def forecast( else: n_ens_members_nwp = 1 + # Now, repeat the nowcast ensemble members or the nwp models/members until + # it has the same amount of members as n_ens_members_max. For instance, if + # you have 10 ensemble nowcasts members and 3 NWP members, the output will + # be an ensemble of 10 members. Hence, the three NWP members are blended + # with the first three members of the nowcast (member one with member one, + # two with two, etc.), subsequently, the same NWP members are blended with + # the next three members (NWP member one with member 4, NWP member 2 with + # member 5, etc.), until 10 is reached. n_ens_members_max = max(n_ens_members_nowcast, n_ens_members_nwp) n_ens_members_min = min(n_ens_members_nowcast, n_ens_members_nwp) diff --git a/pysteps/nowcasts/interface.py b/pysteps/nowcasts/interface.py index 2f5908fd0..0b7b94bce 100644 --- a/pysteps/nowcasts/interface.py +++ b/pysteps/nowcasts/interface.py @@ -39,7 +39,6 @@ sprog, steps, sseps, - linear_blending, ) from pysteps.nowcasts import lagrangian_probability @@ -54,7 +53,6 @@ _nowcast_methods["sprog"] = sprog.forecast _nowcast_methods["sseps"] = sseps.forecast _nowcast_methods["steps"] = steps.forecast -_nowcast_methods["linear blending"] = linear_blending.forecast def get_method(name): @@ -92,9 +90,6 @@ def get_method(name): | sseps | short-space ensemble prediction system (SSEPS). | | | Essentially, this is a localization of STEPS | +-----------------+-------------------------------------------------------+ - | linear | the linear blending of a nowcast method with other | - | blending | data (e.g. NWP data). | - +-----------------+-------------------------------------------------------+ """ if isinstance(name, str): name = name.lower() diff --git a/pysteps/tests/test_nowcasts_linear_blending.py b/pysteps/tests/test_blending_linear_blending.py similarity index 85% rename from pysteps/tests/test_nowcasts_linear_blending.py rename to pysteps/tests/test_blending_linear_blending.py index ca6053f30..56a8b5031 100644 --- a/pysteps/tests/test_nowcasts_linear_blending.py +++ b/pysteps/tests/test_blending_linear_blending.py @@ -2,9 +2,9 @@ import numpy as np import pytest -from pysteps.nowcasts.linear_blending import forecast +from pysteps.blending.linear_blending import forecast from numpy.testing import assert_array_almost_equal -from pysteps.utils import conversion, transformation +from pysteps.utils import transformation # Test function arguments linear_arg_values = [ @@ -61,47 +61,49 @@ def test_linear_blending( ), "Control time needs to be a multiple of the time step" # Initialise dummy NWP data - R_NWP = np.zeros((n_timesteps, 200, 200)) + r_nwp = np.zeros((n_timesteps, 200, 200)) for i in range(100): - R_NWP[:, i, :] = 11.0 + r_nwp[:, i, :] = 11.0 # Define nowcast input data - R_input = np.zeros((200, 200)) + r_input = np.zeros((200, 200)) for i in range(100, 200): - R_input[i, :] = 11.0 + r_input[i, :] = 11.0 # Transform from mm/h to dB - R_input, _ = transformation.dB_transform( - R_input, None, threshold=0.1, zerovalue=-15.0 + r_input, _ = transformation.dB_transform( + r_input, None, threshold=0.1, zerovalue=-15.0 ) # Calculate the blended field - R_blended = forecast( - R_input, + r_blended = forecast( + r_input, + dict({"unit": "mm/h", "transform": "dB"}), V, n_timesteps, timestep, nowcast_method, - R_NWP, + r_nwp, + dict({"unit": "mm/h"}), start_blending=start_blending, end_blending=end_blending, ) # Assert that the blended field has the expected dimension - assert R_blended.shape == ( + assert r_blended.shape == ( n_timesteps, 200, 200, ), "The shape of the blended array does not have the expected value. The shape is {}".format( - R_blended.shape + r_blended.shape ) # Assert that the blended field at the control time step is equal to # a constant field with the expected value. assert_array_almost_equal( - R_blended[controltime // timestep - 1], + r_blended[controltime // timestep - 1], np.ones((200, 200)) * 5.5, err_msg="The blended array does not have the expected value", ) From 758f88b3b70cecb3774d5d17ff95e942ecfa9fba Mon Sep 17 00:00:00 2001 From: RubenImhoff Date: Tue, 19 Oct 2021 16:26:08 +0200 Subject: [PATCH 7/7] Fixes for tests --- pysteps/tests/test_blending_linear_blending.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/tests/test_blending_linear_blending.py b/pysteps/tests/test_blending_linear_blending.py index 56a8b5031..712963eba 100644 --- a/pysteps/tests/test_blending_linear_blending.py +++ b/pysteps/tests/test_blending_linear_blending.py @@ -86,7 +86,7 @@ def test_linear_blending( timestep, nowcast_method, r_nwp, - dict({"unit": "mm/h"}), + dict({"unit": "mm/h", "transform": None}), start_blending=start_blending, end_blending=end_blending, )