From d7f84fb4ef340f620b91c0384c542636bba16d1d Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Thu, 11 Sep 2025 13:45:36 +0200 Subject: [PATCH 1/8] fix blend utils and xarray dataformat --- pysteps/blending/utils.py | 276 ++++-------- pysteps/downscaling/rainfarm.py | 7 +- pysteps/io/importers.py | 13 +- pysteps/io/readers.py | 7 +- pysteps/tests/helpers.py | 14 +- .../tests/test_blending_linear_blending.py | 2 +- pysteps/tests/test_blending_utils.py | 410 +++++++----------- pysteps/tests/test_io_bom_rf3.py | 1 - pysteps/tests/test_io_dwd_hdf5.py | 1 - pysteps/tests/test_io_knmi_hdf5.py | 1 - pysteps/tests/test_io_mch_gif.py | 1 - pysteps/tests/test_io_mrms_grib.py | 1 - pysteps/tests/test_io_nowcast_importers.py | 4 - pysteps/tests/test_io_opera_hdf5.py | 3 - pysteps/tests/test_io_saf_crri.py | 1 - pysteps/tests/test_nowcasts_anvil.py | 7 +- pysteps/tests/test_nowcasts_sseps.py | 7 +- pysteps/tests/test_utils_conversion.py | 23 - pysteps/tests/test_utils_dimension.py | 1 - pysteps/tests/test_utils_reprojection.py | 2 +- pysteps/tests/test_utils_transformation.py | 9 - pysteps/utils/conversion.py | 8 +- pysteps/utils/transformation.py | 12 +- pysteps/xarray_helpers.py | 20 +- 24 files changed, 310 insertions(+), 521 deletions(-) diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py index aaed2cfa2..49e95619b 100644 --- a/pysteps/blending/utils.py +++ b/pysteps/blending/utils.py @@ -19,10 +19,12 @@ """ import datetime +from typing import Any import warnings from pathlib import Path import numpy as np +import xarray as xr from pysteps.cascade import get_method as cascade_get_method from pysteps.cascade.bandpass_filters import filter_gaussian @@ -241,12 +243,7 @@ def blend_optical_flows(flows, weights): def decompose_NWP( - R_NWP, - NWP_model, - analysis_time, - timestep, - valid_times, - output_path, + precip_nwp_dataset: xr.Dataset, num_cascade_levels=8, num_workers=1, decomp_method="fft", @@ -255,7 +252,7 @@ def decompose_NWP( normalize=True, compute_stats=True, compact_output=True, -): +) -> xr.Dataset: """Decomposes the NWP forecast data into cascades and saves it in a netCDF file @@ -269,11 +266,6 @@ def decompose_NWP( analysis_time: numpy.datetime64 The analysis time of the NWP forecast. The analysis time is assumed to be a numpy.datetime64 type as imported by the pysteps importer - timestep: int - Timestep in minutes between subsequent NWP forecast fields - valid_times: array_like - Array containing the valid times of the NWP forecast fields. The times are - assumed to be numpy.datetime64 types as imported by the pysteps importer. output_path: str The location where to save the file with the NWP cascade. Defaults to the path_workdir specified in the rcparams file. @@ -315,62 +307,41 @@ def decompose_NWP( Returns ------- - None + xarray.Dataset + The same dataset as was passed in but with the precip data replaced + with decomposed precip data and means and stds added """ - if not NETCDF4_IMPORTED: - raise MissingOptionalDependency( - "netCDF4 package is required to save the decomposed NWP data, " - "but it is not installed" - ) - - # Make a NetCDF file - output_date = f"{analysis_time.astype('datetime64[us]').astype(datetime.datetime):%Y%m%d%H%M%S}" - outfn = Path(output_path) / f"cascade_{NWP_model}_{output_date}.nc" - ncf = netCDF4.Dataset(outfn, "w", format="NETCDF4") - - # Express times relative to the zero time - zero_time = np.datetime64("1970-01-01T00:00:00", "ns") - valid_times = np.array(valid_times) - zero_time - analysis_time = analysis_time - zero_time - - # Set attributes of decomposition method - ncf.domain = domain - ncf.normalized = int(normalize) - ncf.compact_output = int(compact_output) - ncf.analysis_time = int(analysis_time) - ncf.timestep = int(timestep) - - # Create dimensions - ncf.createDimension("time", R_NWP.shape[0]) - ncf.createDimension("cascade_levels", num_cascade_levels) - ncf.createDimension("x", R_NWP.shape[2]) - ncf.createDimension("y", R_NWP.shape[1]) - - # Create variables (decomposed cascade, means and standard deviations) - R_d = ncf.createVariable( - "pr_decomposed", - np.float32, - ("time", "cascade_levels", "y", "x"), - zlib=True, - complevel=4, - ) - means = ncf.createVariable("means", np.float64, ("time", "cascade_levels")) - stds = ncf.createVariable("stds", np.float64, ("time", "cascade_levels")) - v_times = ncf.createVariable("valid_times", np.float64, ("time",)) - v_times.units = "nanoseconds since 1970-01-01 00:00:00" - - # The valid times are saved as an array of floats, because netCDF files can't handle datetime types - v_times[:] = np.array([np.float64(valid_times[i]) for i in range(len(valid_times))]) - + nwp_precip_var = precip_nwp_dataset.attrs["precip_var"] + precip_nwp = precip_nwp_dataset[nwp_precip_var].values # Decompose the NWP data - filter_g = filter_gaussian(R_NWP.shape[1:], num_cascade_levels) - fft = utils_get_method(fft_method, shape=R_NWP.shape[1:], n_threads=num_workers) + filter_g = filter_gaussian(precip_nwp.shape[1:], num_cascade_levels) + fft = utils_get_method( + fft_method, shape=precip_nwp.shape[1:], n_threads=num_workers + ) decomp_method, _ = cascade_get_method(decomp_method) - for i in range(R_NWP.shape[0]): - R_ = decomp_method( - field=R_NWP[i, :, :], + pr_decomposed = np.zeros( + ( + precip_nwp.shape[0], + num_cascade_levels, + precip_nwp.shape[1], + precip_nwp.shape[2], + ), + dtype=np.float32, + ) + means = np.zeros( + (precip_nwp.shape[0], num_cascade_levels), + dtype=np.float64, + ) + stds = np.zeros( + (precip_nwp.shape[0], num_cascade_levels), + dtype=np.float64, + ) + + for i in range(precip_nwp.shape[0]): + decomposed_precip_nwp = decomp_method( + field=precip_nwp[i, :, :], bp_filter=filter_g, fft_method=fft, input_domain=domain, @@ -380,157 +351,96 @@ def decompose_NWP( compact_output=compact_output, ) - # Save data to netCDF file - # print(R_["cascade_levels"]) - R_d[i, :, :, :] = R_["cascade_levels"] - means[i, :] = R_["means"] - stds[i, :] = R_["stds"] + pr_decomposed[i, :, :, :] = decomposed_precip_nwp["cascade_levels"] + means[i, :] = decomposed_precip_nwp["means"] + stds[i, :] = decomposed_precip_nwp["stds"] - # Close the file - ncf.close() - - -def compute_store_nwp_motion( - precip_nwp, - oflow_method, - analysis_time, - nwp_model, - output_path, + precip_nwp_dataset = precip_nwp_dataset.assign_coords( + cascade_level=( + "cascade_level", + np.arange(num_cascade_levels), + {"long_name": "cascade level", "units": ""}, + ) + ) + precip_nwp_dataset = precip_nwp_dataset.drop_vars(nwp_precip_var) + precip_nwp_dataset[nwp_precip_var] = ( + ["time", "cascade_level", "y", "x"], + pr_decomposed, + ) + precip_nwp_dataset["means"] = (["time", "cascade_level"], means) + precip_nwp_dataset["stds"] = (["time", "cascade_level"], stds) + return precip_nwp_dataset + + +def preprocess_and_store_nwp_data( + precip_nwp_dataset: xr.Dataset, + oflow_method: str, + nwp_model: str, + output_path: str | None, + decompose_nwp: bool, + decompose_kwargs: dict[str, Any] = {}, ): """Computes, per forecast lead time, the velocity field of an NWP model field. Parameters ---------- - precip_nwp: array-like - Array of dimension (n_timesteps, x, y) containing the precipitation forecast + precip_nwp_dataset: xarray.Dataset + xarray Dataset containing the precipitation forecast from some NWP model. oflow_method: {'constant', 'darts', 'lucaskanade', 'proesmans', 'vet'}, optional An optical flow method from pysteps.motion.get_method. - analysis_time: numpy.datetime64 - The analysis time of the NWP forecast. The analysis time is assumed to be a - numpy.datetime64 type as imported by the pysteps importer. nwp_model: str The name of the NWP model. output_path: str, optional - The location where to save the file with the NWP velocity fields. Defaults + The location where to save the netcdf file with the NWP velocity fields. Defaults to the path_workdir specified in the rcparams file. + decompose_nwp: bool + Defines wether or not the NWP needs to be decomposed before storing. This can + be beneficial for performance, because then the decomposition does not need + to happen during the blending anymore. It can however also be detrimental because + this increases the amount of storage and RAM required for the blending. + decompose_kwargs: dict + Keyword arguments passed to the decompose_NWP method. Returns ------- Nothing """ + if not NETCDF4_IMPORTED: + raise MissingOptionalDependency( + "netCDF4 package is required to save the NWP data, " + "but it is not installed" + ) + # Set the output file + analysis_time = precip_nwp_dataset.time.values[0] output_date = f"{analysis_time.astype('datetime64[us]').astype(datetime.datetime):%Y%m%d%H%M%S}" - outfn = Path(output_path) / f"motion_{nwp_model}_{output_date}.npy" + outfn = Path(output_path) / f"preprocessed_{nwp_model}_{output_date}.nc" + nwp_precip_var = precip_nwp_dataset.attrs["precip_var"] + precip_nwp = precip_nwp_dataset[nwp_precip_var].values # Get the velocity field per time step - v_nwp = np.zeros((precip_nwp.shape[0], 2, precip_nwp.shape[1], precip_nwp.shape[2])) + v_nwp_x = np.zeros((precip_nwp.shape[0], precip_nwp.shape[1], precip_nwp.shape[2])) + v_nwp_y = np.zeros((precip_nwp.shape[0], precip_nwp.shape[1], precip_nwp.shape[2])) # Loop through the timesteps. We need two images to construct a motion # field, so we can start from timestep 1. for t in range(1, precip_nwp.shape[0]): - v_nwp[t] = oflow_method(precip_nwp[t - 1 : t + 1, :, :]) + v_nwp_dataset = oflow_method(precip_nwp_dataset.isel(time=slice(t - 1, t + 1))) + v_nwp_x[t] = v_nwp_dataset.velocity_x + v_nwp_y[t] = v_nwp_dataset.velocity_y # Make timestep 0 the same as timestep 1. - v_nwp[0] = v_nwp[1] + v_nwp_x[0] = v_nwp_x[1] + v_nwp_y[0] = v_nwp_y[1] + precip_nwp_dataset["velocity_x"] = (["time", "y", "x"], v_nwp_x) + precip_nwp_dataset["velocity_y"] = (["time", "y", "x"], v_nwp_y) - assert v_nwp.ndim == 4, "v_nwp must be a four-dimensional array" + if decompose_nwp: + precip_nwp_dataset = decompose_NWP(precip_nwp_dataset, **decompose_kwargs) # Save it as a numpy array - np.save(outfn, v_nwp) - - -def load_NWP(input_nc_path_decomp, input_path_velocities, start_time, n_timesteps): - """Loads the decomposed NWP and velocity data from the netCDF files - - Parameters - ---------- - input_nc_path_decomp: str - Path to the saved netCDF file containing the decomposed NWP data. - input_path_velocities: str - Path to the saved numpy binary file containing the estimated velocity - fields from the NWP data. - start_time: numpy.datetime64 - The start time of the nowcasting. Assumed to be a numpy.datetime64 type - n_timesteps: int - Number of time steps to forecast - - Returns - ------- - R_d: list - A list of dictionaries with each element in the list corresponding to - a different time step. Each dictionary has the same structure as the - output of the decomposition function - uv: array-like - Array of shape (timestep,2,m,n) containing the x- and y-components - of the advection field for the (NWP) model field per forecast lead time. - """ - - if not NETCDF4_IMPORTED: - raise MissingOptionalDependency( - "netCDF4 package is required to load the decomposed NWP data, " - "but it is not installed" - ) - - # Open the file - ncf_decomp = netCDF4.Dataset(input_nc_path_decomp, "r", format="NETCDF4") - velocities = np.load(input_path_velocities) - - decomp_dict = { - "domain": ncf_decomp.domain, - "normalized": bool(ncf_decomp.normalized), - "compact_output": bool(ncf_decomp.compact_output), - } - - # Convert the start time and the timestep to datetime64 and timedelta64 type - zero_time = np.datetime64("1970-01-01T00:00:00", "ns") - analysis_time = np.timedelta64(int(ncf_decomp.analysis_time), "ns") + zero_time - - timestep = ncf_decomp.timestep - timestep = np.timedelta64(timestep, "m") - - valid_times = ncf_decomp.variables["valid_times"][:] - valid_times = np.array( - [np.timedelta64(int(valid_times[i]), "ns") for i in range(len(valid_times))] - ) - valid_times = valid_times + zero_time - - # Find the indices corresponding with the required start and end time - start_i = (start_time - analysis_time) // timestep - assert analysis_time + start_i * timestep == start_time - end_i = start_i + n_timesteps + 1 - - # Check if the requested end time (the forecast horizon) is in the stored data. - # If not, raise an error - if end_i > ncf_decomp.variables["pr_decomposed"].shape[0]: - raise IndexError( - "The requested forecast horizon is outside the stored NWP forecast horizon. Either request a shorter forecast horizon or store a longer NWP forecast horizon" - ) - - # Add the valid times to the output - decomp_dict["valid_times"] = valid_times[start_i:end_i] - - # Slice the velocity fields with the start and end indices - uv = velocities[start_i:end_i, :, :, :] - - # Initialise the list of dictionaries which will serve as the output (cf: the STEPS function) - R_d = list() - - pr_decomposed = ncf_decomp.variables["pr_decomposed"][start_i:end_i, :, :, :] - means = ncf_decomp.variables["means"][start_i:end_i, :] - stds = ncf_decomp.variables["stds"][start_i:end_i, :] - - for i in range(n_timesteps + 1): - decomp_dict["cascade_levels"] = np.ma.filled( - pr_decomposed[i], fill_value=np.nan - ) - decomp_dict["means"] = np.ma.filled(means[i], fill_value=np.nan) - decomp_dict["stds"] = np.ma.filled(stds[i], fill_value=np.nan) - - R_d.append(decomp_dict.copy()) - - ncf_decomp.close() - return R_d, uv + precip_nwp_dataset.to_netcdf(outfn) def check_norain(precip_arr, precip_thr=None, norain_thr=0.0): diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index b10e45c83..fddabac0d 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -318,7 +318,12 @@ def downscale( noise_dataset = xr.Dataset( data_vars={precip_var: (["time", "y", "x"], [noise_field])}, coords={ - "time": (["time"], precip_dataset.time.values, precip_dataset.time.attrs), + "time": ( + ["time"], + precip_dataset.time.values, + precip_dataset.time.attrs, + precip_dataset.time.encoding, + ), "y": ( ["y"], y_new, diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index 71a725842..cd45df034 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -66,6 +66,7 @@ +------------------+----------------------------------------------------------+ # XR: Move this to appropriate place +# XR: Add means, stds cascade level The data and metadata is then postprocessed into an xarray dataset. This dataset will always contain an x and y dimension, but can be extended with a time dimension and/or an ensemble member dimension over the course of the process. @@ -538,7 +539,6 @@ def _import_mrms_grib(filename, extent=None, window_size=4, **kwargs): ypixelsize=ysize, unit="mm/h", accutime=2.0, - transform=None, zerovalue=0, projection=proj_def.strip(), yorigin="upper", @@ -581,7 +581,6 @@ def import_bom_rf3(filename, **kwargs): precip, geodata = _import_bom_rf3_data(filename) metadata = geodata - metadata["transform"] = None metadata["zerovalue"] = np.nanmin(precip) metadata["threshold"] = _get_threshold_value(precip) @@ -1021,7 +1020,8 @@ def import_knmi_hdf5( metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute" metadata["accutime"] = accutime metadata["unit"] = unit - metadata["transform"] = transform + if transform is not None: + metadata["transform"] = transform metadata["zerovalue"] = 0.0 metadata["threshold"] = _get_threshold_value(precip) metadata["zr_a"] = 200.0 @@ -1152,7 +1152,6 @@ def import_mch_gif(filename, product, unit, accutime, **kwargs): metadata["accutime"] = accutime metadata["unit"] = unit - metadata["transform"] = None metadata["zerovalue"] = np.nanmin(precip) metadata["threshold"] = _get_threshold_value(precip) metadata["institution"] = "MeteoSwiss" @@ -1285,11 +1284,11 @@ def import_mch_hdf5(filename, qty="RATE", **kwargs): "institution": "MeteoSwiss", "accutime": 5.0, "unit": unit, - "transform": transform, "zerovalue": np.nanmin(precip), "threshold": thr, "zr_a": 316.0, "zr_b": 1.5, + **({"transform": transform} if transform is not None else {}), } ) @@ -1365,7 +1364,6 @@ def import_mch_metranet(filename, product, unit, accutime): metadata["institution"] = "MeteoSwiss" metadata["accutime"] = accutime metadata["unit"] = unit - metadata["transform"] = None metadata["zerovalue"] = np.nanmin(precip) metadata["threshold"] = _get_threshold_value(precip) metadata["zr_a"] = 316.0 @@ -1621,9 +1619,9 @@ def import_odim_hdf5(filename, qty="RATE", **kwargs): "institution": "Odyssey datacentre", "accutime": 15.0, "unit": unit, - "transform": transform, "zerovalue": np.nanmin(precip), "threshold": _get_threshold_value(precip), + **({"transform": transform} if transform is not None else {}), } metadata.update(kwargs) @@ -1723,7 +1721,6 @@ def import_saf_crri(filename, extent=None, **kwargs): precip, quality = _import_saf_crri_data(filename, idx_x, idx_y) - metadata["transform"] = None metadata["zerovalue"] = np.nanmin(precip) metadata["threshold"] = _get_threshold_value(precip) diff --git a/pysteps/io/readers.py b/pysteps/io/readers.py index 4fc29876d..1fd786042 100644 --- a/pysteps/io/readers.py +++ b/pysteps/io/readers.py @@ -94,11 +94,8 @@ def read_timeseries(inputfns, importer, timestep=None, **kwargs) -> xr.Dataset | time=( "time", [inputfns[1][i]], - { - "long_name": "forecast time", - "units": f"seconds since {startdate:%Y-%m-%d %H:%M:%S}", - "stepsize": timestep, - }, + {"long_name": "forecast time", "stepsize": timestep}, + {"units": f"seconds since {startdate:%Y-%m-%d %H:%M:%S}"}, ) ) datasets.append(dataset_) diff --git a/pysteps/tests/helpers.py b/pysteps/tests/helpers.py index e2dca3df8..a70c5a5af 100644 --- a/pysteps/tests/helpers.py +++ b/pysteps/tests/helpers.py @@ -42,12 +42,14 @@ def assert_dataset_equivalent(dataset1: xr.Dataset, dataset2: xr.Dataset) -> Non dataset2[precip_var].attrs["zerovalue"], ) assert dataset1[precip_var].attrs["units"] == dataset2[precip_var].attrs["units"] - assert ( - dataset1[precip_var].attrs["transform"] - == dataset2[precip_var].attrs["transform"] - or dataset1[precip_var].attrs["transform"] is None - and dataset2[precip_var].attrs["transform"] is None - ) + if ( + "transform" in dataset1[precip_var].attrs + or "transform" in dataset2[precip_var].attrs + ): + assert ( + dataset1[precip_var].attrs["transform"] + == dataset2[precip_var].attrs["transform"] + ) if ( "accutime" in dataset1[precip_var].attrs or "accutime" in dataset2[precip_var].attrs diff --git a/pysteps/tests/test_blending_linear_blending.py b/pysteps/tests/test_blending_linear_blending.py index 4dbbaf856..8dd139174 100644 --- a/pysteps/tests/test_blending_linear_blending.py +++ b/pysteps/tests/test_blending_linear_blending.py @@ -236,7 +236,7 @@ def test_linear_blending( timestep, nowcast_method, r_nwp, - dict({"unit": "mm/h", "transform": None}), + dict({"unit": "mm/h"}), start_blending=start_blending, end_blending=end_blending, fill_nwp=fill_nwp, diff --git a/pysteps/tests/test_blending_utils.py b/pysteps/tests/test_blending_utils.py index b423460a8..fccbe9067 100644 --- a/pysteps/tests/test_blending_utils.py +++ b/pysteps/tests/test_blending_utils.py @@ -6,13 +6,13 @@ import numpy as np import pytest from numpy.testing import assert_array_almost_equal - +import xarray as xr import pysteps from pysteps.blending.utils import ( blend_cascades, blend_optical_flows, compute_smooth_dilated_mask, - compute_store_nwp_motion, + preprocess_and_store_nwp_data, decompose_NWP, load_NWP, recompose_cascade, @@ -47,7 +47,6 @@ nwp_metadata = dict( projection=nwp_proj, institution="Royal Meteorological Institute of Belgium", - transform=None, zerovalue=0.0, threshold=0, unit="mm", @@ -62,38 +61,11 @@ y2=0.0, ) precip_nwp_dataset = convert_input_to_xarray_dataset( - precip_nwp, None, nwp_metadata, datetime.now(tz=timezone.utc) -) - -# Get the analysis time and valid time -times_nwp = np.array( - [ - "2021-07-04T16:05:00.000000000", - "2021-07-04T16:10:00.000000000", - "2021-07-04T16:15:00.000000000", - "2021-07-04T16:20:00.000000000", - "2021-07-04T16:25:00.000000000", - "2021-07-04T16:30:00.000000000", - "2021-07-04T16:35:00.000000000", - "2021-07-04T16:40:00.000000000", - "2021-07-04T16:45:00.000000000", - "2021-07-04T16:50:00.000000000", - "2021-07-04T16:55:00.000000000", - "2021-07-04T17:00:00.000000000", - "2021-07-04T17:05:00.000000000", - "2021-07-04T17:10:00.000000000", - "2021-07-04T17:15:00.000000000", - "2021-07-04T17:20:00.000000000", - "2021-07-04T17:25:00.000000000", - "2021-07-04T17:30:00.000000000", - "2021-07-04T17:35:00.000000000", - "2021-07-04T17:40:00.000000000", - "2021-07-04T17:45:00.000000000", - "2021-07-04T17:50:00.000000000", - "2021-07-04T17:55:00.000000000", - "2021-07-04T18:00:00.000000000", - ], - dtype="datetime64[ns]", + precip_nwp, + None, + nwp_metadata, + datetime.fromisoformat("2021-07-04T16:05:00.000000000"), + 300, ) @@ -112,12 +84,16 @@ # Transform the data transformer = pysteps.utils.get_method("dB") precip_nwp_dataset = transformer( - precip_nwp_dataset, threshold=nwp_metadata["threshold"] + precip_nwp_dataset, threshold=precip_nwp_dataset[nwp_precip_var].attrs["threshold"] ) # Set two issue times for testing -issue_time_first = times_nwp[0] -issue_time_second = times_nwp[3] +issue_time_first = np.datetime64( + datetime.fromisoformat("2021-07-04T16:05:00.000000000") +) +issue_time_second = np.datetime64( + datetime.fromisoformat("2021-07-04T16:20:00.000000000") +) # Set the blending weights (we'll blend with a 50-50 weight) weights = np.full((2, 8), fill_value=0.5) @@ -130,7 +106,6 @@ "issue_times", "timestep", "n_timesteps", - "valid_times", "shape", "weights", ) @@ -143,7 +118,6 @@ [issue_time_first, issue_time_second], 5.0, 3, - times_nwp, precip_nwp_dataset[nwp_precip_var].values.shape[1:], weights, ) @@ -177,7 +151,6 @@ def test_blending_utils( issue_times, timestep, n_timesteps, - valid_times, shape, weights, ): @@ -193,55 +166,29 @@ def test_blending_utils( ### # Compute and store the motion ### - compute_store_nwp_motion( - precip_nwp=precip_nwp_dataset, + preprocess_and_store_nwp_data( + precip_nwp_dataset=precip_nwp_dataset, oflow_method=oflow_method, - analysis_time=valid_times[0], nwp_model=nwp_model, output_path=tmpdir, + decompose_nwp=True, + decompose_kwargs=dict( + num_cascade_levels=8, + num_workers=1, + decomp_method="fft", + fft_method="numpy", + domain="spatial", + normalize=True, + compute_stats=True, + compact_output=False, + ), ) # Check if file exists - date_string = np.datetime_as_string(valid_times[0], "s") - motion_file = os.path.join( - tmpdir, - "motion_" - + nwp_model - + "_" - + date_string[:4] - + date_string[5:7] - + date_string[8:10] - + date_string[11:13] - + date_string[14:16] - + date_string[17:19] - + ".npy", - ) - assert os.path.exists(motion_file) - - ### - # Decompose and store NWP forecast - ### - decompose_NWP( - R_NWP=precip_nwp_dataset, - NWP_model=nwp_model, - analysis_time=valid_times[0], - timestep=timestep, - valid_times=valid_times, - num_cascade_levels=8, - num_workers=1, - output_path=tmpdir, - decomp_method="fft", - fft_method="numpy", - domain="spatial", - normalize=True, - compute_stats=True, - compact_output=False, - ) - - # Check if file exists - decomp_file = os.path.join( + date_string = np.datetime_as_string(precip_nwp_dataset.time.values[0], "s") + preprocessed_file = os.path.join( tmpdir, - "cascade_" + "preprocessed_" + nwp_model + "_" + date_string[:4] @@ -252,188 +199,159 @@ def test_blending_utils( + date_string[17:19] + ".nc", ) - assert os.path.exists(decomp_file) + assert os.path.exists(preprocessed_file) ### # Now check if files load correctly for two different issue times ### - precip_decomposed_nwp_first, v_nwp_first = load_NWP( - input_nc_path_decomp=os.path.join(decomp_file), - input_path_velocities=os.path.join(motion_file), - start_time=issue_times[0], - n_timesteps=n_timesteps, - ) - - precip_decomposed_nwp_second, v_nwp_second = load_NWP( - input_nc_path_decomp=os.path.join(decomp_file), - input_path_velocities=os.path.join(motion_file), - start_time=issue_times[1], - n_timesteps=n_timesteps, - ) + with xr.open_dataset(preprocessed_file) as nwp_file: + precip_decomposed_nwp_first, v_nwp_first = load_NWP( + input_nc_path_decomp=os.path.join(decomp_file), + input_path_velocities=os.path.join(preprocessed_file), + start_time=issue_times[0], + n_timesteps=n_timesteps, + ) - # Check if the output type and shapes are correct - assert isinstance(precip_decomposed_nwp_first, list) - assert isinstance(precip_decomposed_nwp_second, list) - assert isinstance(precip_decomposed_nwp_first[0], dict) - assert isinstance(precip_decomposed_nwp_second[0], dict) - - assert "domain" in precip_decomposed_nwp_first[0] - assert "normalized" in precip_decomposed_nwp_first[0] - assert "compact_output" in precip_decomposed_nwp_first[0] - assert "valid_times" in precip_decomposed_nwp_first[0] - assert "cascade_levels" in precip_decomposed_nwp_first[0] - assert "means" in precip_decomposed_nwp_first[0] - assert "stds" in precip_decomposed_nwp_first[0] - - assert precip_decomposed_nwp_first[0]["cascade_levels"].shape == ( - 8, - shape[0], - shape[1], - ) - assert precip_decomposed_nwp_first[0]["domain"] == "spatial" - assert precip_decomposed_nwp_first[0]["normalized"] == True - assert precip_decomposed_nwp_first[0]["compact_output"] == False - assert len(precip_decomposed_nwp_first) == n_timesteps + 1 - assert len(precip_decomposed_nwp_second) == n_timesteps + 1 - assert precip_decomposed_nwp_first[0]["means"].shape[0] == 8 - assert precip_decomposed_nwp_first[0]["stds"].shape[0] == 8 - - assert np.array(v_nwp_first).shape == (n_timesteps + 1, 2, shape[0], shape[1]) - assert np.array(v_nwp_second).shape == (n_timesteps + 1, 2, shape[0], shape[1]) - - # Check if the right times are loaded - assert ( - precip_decomposed_nwp_first[0]["valid_times"][0] == valid_times[0] - ), "Not the right valid times were loaded for the first forecast" - assert ( - precip_decomposed_nwp_second[0]["valid_times"][0] == valid_times[3] - ), "Not the right valid times were loaded for the second forecast" - - # Check, for a sample, if the stored motion fields are as expected - assert_array_almost_equal( - v_nwp_first[1], - oflow_method(precip_nwp_dataset[0:2, :, :]), - decimal=3, - err_msg="Stored motion field of first forecast not equal to expected motion field", - ) - assert_array_almost_equal( - v_nwp_second[1], - oflow_method(precip_nwp_dataset[3:5, :, :]), - decimal=3, - err_msg="Stored motion field of second forecast not equal to expected motion field", - ) + precip_decomposed_nwp_second, v_nwp_second = load_NWP( + input_nc_path_decomp=os.path.join(decomp_file), + input_path_velocities=os.path.join(preprocessed_file), + start_time=issue_times[1], + n_timesteps=n_timesteps, + ) - ### - # Stack the cascades - ### - precip_decomposed_first_stack, mu_first_stack, sigma_first_stack = stack_cascades( - R_d=precip_decomposed_nwp_first, donorm=False - ) + # Check, for a sample, if the stored motion fields are as expected + assert_array_almost_equal( + v_nwp_first[1], + oflow_method(precip_nwp_dataset[0:2, :, :]), + decimal=3, + err_msg="Stored motion field of first forecast not equal to expected motion field", + ) + assert_array_almost_equal( + v_nwp_second[1], + oflow_method(precip_nwp_dataset[3:5, :, :]), + decimal=3, + err_msg="Stored motion field of second forecast not equal to expected motion field", + ) - print(precip_decomposed_nwp_first) - print(precip_decomposed_first_stack) - print(mu_first_stack) + ### + # Stack the cascades + ### + precip_decomposed_first_stack, mu_first_stack, sigma_first_stack = ( + stack_cascades(R_d=precip_decomposed_nwp_first, donorm=False) + ) - ( - precip_decomposed_second_stack, - mu_second_stack, - sigma_second_stack, - ) = stack_cascades(R_d=precip_decomposed_nwp_second, donorm=False) - - # Check if the array shapes are still correct - assert precip_decomposed_first_stack.shape == ( - n_timesteps + 1, - 8, - shape[0], - shape[1], - ) - assert mu_first_stack.shape == (n_timesteps + 1, 8) - assert sigma_first_stack.shape == (n_timesteps + 1, 8) + print(precip_decomposed_nwp_first) + print(precip_decomposed_first_stack) + print(mu_first_stack) - ### - # Blend the cascades - ### - precip_decomposed_blended = blend_cascades( - cascades_norm=np.stack( - (precip_decomposed_first_stack[0], precip_decomposed_second_stack[0]) - ), - weights=weights, - ) + ( + precip_decomposed_second_stack, + mu_second_stack, + sigma_second_stack, + ) = stack_cascades(R_d=precip_decomposed_nwp_second, donorm=False) + + # Check if the array shapes are still correct + assert precip_decomposed_first_stack.shape == ( + n_timesteps + 1, + 8, + shape[0], + shape[1], + ) + assert mu_first_stack.shape == (n_timesteps + 1, 8) + assert sigma_first_stack.shape == (n_timesteps + 1, 8) + + ### + # Blend the cascades + ### + precip_decomposed_blended = blend_cascades( + cascades_norm=np.stack( + (precip_decomposed_first_stack[0], precip_decomposed_second_stack[0]) + ), + weights=weights, + ) - assert precip_decomposed_blended.shape == precip_decomposed_first_stack[0].shape + assert precip_decomposed_blended.shape == precip_decomposed_first_stack[0].shape - ### - # Blend the optical flow fields - ### - v_nwp_blended = blend_optical_flows( - flows=np.stack((v_nwp_first[1], v_nwp_second[1])), weights=weights[:, 1] - ) + ### + # Blend the optical flow fields + ### + v_nwp_blended = blend_optical_flows( + flows=np.stack((v_nwp_first[1], v_nwp_second[1])), weights=weights[:, 1] + ) - assert v_nwp_blended.shape == v_nwp_first[1].shape - assert_array_almost_equal( - v_nwp_blended, - ( - oflow_method(precip_nwp_dataset[0:2, :, :]) - + oflow_method(precip_nwp_dataset[3:5, :, :]) + assert v_nwp_blended.shape == v_nwp_first[1].shape + assert_array_almost_equal( + v_nwp_blended, + ( + oflow_method(precip_nwp_dataset[0:2, :, :]) + + oflow_method(precip_nwp_dataset[3:5, :, :]) + ) + / 2, + decimal=3, + err_msg="Blended motion field does not equal average of the two motion fields", ) - / 2, - decimal=3, - err_msg="Blended motion field does not equal average of the two motion fields", - ) - ### - # Recompose the fields (the non-blended fields are used for this here) - ### - precip_recomposed_first = recompose_cascade( - combined_cascade=precip_decomposed_first_stack[0], - combined_mean=mu_first_stack[0], - combined_sigma=sigma_first_stack[0], - ) - precip_recomposed_second = recompose_cascade( - combined_cascade=precip_decomposed_second_stack[0], - combined_mean=mu_second_stack[0], - combined_sigma=sigma_second_stack[0], - ) + ### + # Recompose the fields (the non-blended fields are used for this here) + ### + precip_recomposed_first = recompose_cascade( + combined_cascade=precip_decomposed_first_stack[0], + combined_mean=mu_first_stack[0], + combined_sigma=sigma_first_stack[0], + ) + precip_recomposed_second = recompose_cascade( + combined_cascade=precip_decomposed_second_stack[0], + combined_mean=mu_second_stack[0], + combined_sigma=sigma_second_stack[0], + ) - assert_array_almost_equal( - precip_recomposed_first, - precip_nwp_dataset[0, :, :], - decimal=3, - err_msg="Recomposed field of first forecast does not equal original field", - ) - assert_array_almost_equal( - precip_recomposed_second, - precip_nwp_dataset[3, :, :], - decimal=3, - err_msg="Recomposed field of second forecast does not equal original field", - ) + assert_array_almost_equal( + precip_recomposed_first, + precip_nwp_dataset[0, :, :], + decimal=3, + err_msg="Recomposed field of first forecast does not equal original field", + ) + assert_array_almost_equal( + precip_recomposed_second, + precip_nwp_dataset[3, :, :], + decimal=3, + err_msg="Recomposed field of second forecast does not equal original field", + ) - precip_arr = precip_nwp_dataset - # rainy fraction is 0.005847 - assert not check_norain(precip_arr, win_fun=None) - assert not check_norain( - precip_arr, precip_thr=nwp_metadata["threshold"], win_fun=None - ) - assert not check_norain( - precip_arr, precip_thr=nwp_metadata["threshold"], norain_thr=0.005, win_fun=None - ) - assert not check_norain(precip_arr, norain_thr=0.005, win_fun=None) - # so with norain_thr beyond this number it should report that there's no rain - assert check_norain(precip_arr, norain_thr=0.006, win_fun=None) - assert check_norain( - precip_arr, precip_thr=nwp_metadata["threshold"], norain_thr=0.006, win_fun=None - ) + precip_arr = precip_nwp_dataset + # rainy fraction is 0.005847 + assert not check_norain(precip_arr, win_fun=None) + assert not check_norain( + precip_arr, precip_thr=nwp_metadata["threshold"], win_fun=None + ) + assert not check_norain( + precip_arr, + precip_thr=nwp_metadata["threshold"], + norain_thr=0.005, + win_fun=None, + ) + assert not check_norain(precip_arr, norain_thr=0.005, win_fun=None) + # so with norain_thr beyond this number it should report that there's no rain + assert check_norain(precip_arr, norain_thr=0.006, win_fun=None) + assert check_norain( + precip_arr, + precip_thr=nwp_metadata["threshold"], + norain_thr=0.006, + win_fun=None, + ) - # also if we set the precipitation threshold sufficiently high, it should report there's no rain - # rainy fraction > 4mm/h is 0.004385 - assert not check_norain(precip_arr, precip_thr=4.0, norain_thr=0.004, win_fun=None) - assert check_norain(precip_arr, precip_thr=4.0, norain_thr=0.005, win_fun=None) + # also if we set the precipitation threshold sufficiently high, it should report there's no rain + # rainy fraction > 4mm/h is 0.004385 + assert not check_norain( + precip_arr, precip_thr=4.0, norain_thr=0.004, win_fun=None + ) + assert check_norain(precip_arr, precip_thr=4.0, norain_thr=0.005, win_fun=None) - # no rain above 100mm/h so it should give norain - assert check_norain(precip_arr, precip_thr=100, win_fun=None) + # no rain above 100mm/h so it should give norain + assert check_norain(precip_arr, precip_thr=100, win_fun=None) - # should always give norain if the threshold is set to 100% - assert check_norain(precip_arr, norain_thr=1.0, win_fun=None) + # should always give norain if the threshold is set to 100% + assert check_norain(precip_arr, norain_thr=1.0, win_fun=None) # Finally, also test the compute_smooth_dilated mask functionality diff --git a/pysteps/tests/test_io_bom_rf3.py b/pysteps/tests/test_io_bom_rf3.py index 8d6a3cec5..294751fe4 100644 --- a/pysteps/tests/test_io_bom_rf3.py +++ b/pysteps/tests/test_io_bom_rf3.py @@ -15,7 +15,6 @@ ) test_metadata_bom = [ - ("transform", None, None), ("zerovalue", 0.0, 0.1), ("projection", expected_proj1, None), ("unit", "mm", None), diff --git a/pysteps/tests/test_io_dwd_hdf5.py b/pysteps/tests/test_io_dwd_hdf5.py index 59de1692d..bf35916a0 100644 --- a/pysteps/tests/test_io_dwd_hdf5.py +++ b/pysteps/tests/test_io_dwd_hdf5.py @@ -56,7 +56,6 @@ def test_io_import_dwd_hdf5_ry_shape(): (precip_dataarray.attrs["accutime"], 5.0, 1e-10), (precip_dataset.time.attrs["stepsize"], 5.0, 1e-10), (precip_dataarray.attrs["units"], "mm/h", None), - (precip_dataarray.attrs["transform"], None, None), (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6), (precip_dataarray.attrs["threshold"], 0.12, 1e-6), ] diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py index 25b96acdc..858a9d8ed 100644 --- a/pysteps/tests/test_io_knmi_hdf5.py +++ b/pysteps/tests/test_io_knmi_hdf5.py @@ -48,7 +48,6 @@ def test_io_import_knmi_hdf5_shape(): "KNMI - Royal Netherlands Meteorological Institute", None, ), - (precip_dataarray.attrs["transform"], None, None), (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6), (precip_dataarray.attrs["threshold"], 0.01, 1e-6), (precip_dataarray.attrs["zr_a"], 200.0, None), diff --git a/pysteps/tests/test_io_mch_gif.py b/pysteps/tests/test_io_mch_gif.py index 13d71f9a5..e6096c3e5 100644 --- a/pysteps/tests/test_io_mch_gif.py +++ b/pysteps/tests/test_io_mch_gif.py @@ -38,7 +38,6 @@ def test_io_import_mch_gif_shape(): (precip_dataarray.attrs["accutime"], 5.0, 1e-10), (precip_dataset.time.attrs["stepsize"], 5.0, 1e-10), (precip_dataarray.attrs["units"], "mm", None), - (precip_dataarray.attrs["transform"], None, None), (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6), (precip_dataarray.attrs["threshold"], 0.0008258007600496956, 1e-19), (precip_dataarray.attrs["zr_a"], 316.0, None), diff --git a/pysteps/tests/test_io_mrms_grib.py b/pysteps/tests/test_io_mrms_grib.py index cc7e20a84..47e3c38c6 100644 --- a/pysteps/tests/test_io_mrms_grib.py +++ b/pysteps/tests/test_io_mrms_grib.py @@ -28,7 +28,6 @@ def test_io_import_mrms_grib(): "xpixelsize": 0.01, "ypixelsize": 0.01, "unit": "mm/h", - "transform": None, "zerovalue": 0, "projection": "+proj=longlat +ellps=IAU76", "yorigin": "upper", diff --git a/pysteps/tests/test_io_nowcast_importers.py b/pysteps/tests/test_io_nowcast_importers.py index e35a5c26b..a19388afe 100644 --- a/pysteps/tests/test_io_nowcast_importers.py +++ b/pysteps/tests/test_io_nowcast_importers.py @@ -14,10 +14,6 @@ upscale=2000, ) -# XR: The following lines are currently needed as we don't have a dedicated xarray writer, mayve we should fix the xarray reader to disallow non types in attrs and to store datetime units as encoding instead of attributes -precip_dataset.time.encoding = {"units": precip_dataset.time.attrs["units"]} -del precip_dataset.time.attrs["units"] - precip_var = precip_dataset.attrs["precip_var"] precip_dataarray = precip_dataset[precip_var] diff --git a/pysteps/tests/test_io_opera_hdf5.py b/pysteps/tests/test_io_opera_hdf5.py index dc62e637b..248fbe899 100644 --- a/pysteps/tests/test_io_opera_hdf5.py +++ b/pysteps/tests/test_io_opera_hdf5.py @@ -92,7 +92,6 @@ def test_io_import_opera_hdf5_nimbus_rain_accum_shape(): # ("yorigin", "upper", None), (precip_odyssey_dataarray.attrs["units"], "mm/h", None), (precip_odyssey.attrs["institution"], "Odyssey datacentre", None), - (precip_odyssey_dataarray.attrs["transform"], None, None), (precip_odyssey_dataarray.attrs["zerovalue"], 0.0, 1e-6), (precip_odyssey_dataarray.attrs["threshold"], 0.01, 1e-6), ] @@ -153,7 +152,6 @@ def test_io_import_opera_hdf5_cirrus_dataset_attrs(variable, expected, tolerance (precip_nimbus_rain_rate_dataarray.attrs["accutime"], 15.0, 1e-10), (precip_nimbus_rain_rate_dataarray.attrs["units"], "mm/h", None), (precip_nimbus_rain_rate.attrs["institution"], "Odyssey datacentre", None), - (precip_nimbus_rain_rate_dataarray.attrs["transform"], None, None), (precip_nimbus_rain_rate_dataarray.attrs["zerovalue"], 0.0, 1e-10), (precip_nimbus_rain_rate_dataarray.attrs["threshold"], 0.01, 1e-10), ] @@ -185,7 +183,6 @@ def test_io_import_opera_hdf5_nimbus_rain_rate_dataset_attrs( (precip_nimbus_rain_accum_dataarray.attrs["accutime"], 15.0, 1e-10), (precip_nimbus_rain_accum_dataarray.attrs["units"], "mm", None), (precip_nimbus_rain_accum.attrs["institution"], "Odyssey datacentre", None), - (precip_nimbus_rain_accum_dataarray.attrs["transform"], None, None), (precip_nimbus_rain_accum_dataarray.attrs["zerovalue"], 0.0, 1e-10), (precip_nimbus_rain_accum_dataarray.attrs["threshold"], 0.01, 1e-10), ] diff --git a/pysteps/tests/test_io_saf_crri.py b/pysteps/tests/test_io_saf_crri.py index 074e39f70..9ecedf135 100644 --- a/pysteps/tests/test_io_saf_crri.py +++ b/pysteps/tests/test_io_saf_crri.py @@ -86,7 +86,6 @@ def test_io_import_saf_crri_geodata(variable, expected, tolerance): ), (precip_dataarray.attrs["accutime"], None, None), (precip_dataarray.attrs["units"], "mm/h", None), - (precip_dataarray.attrs["transform"], None, None), (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6), ] diff --git a/pysteps/tests/test_nowcasts_anvil.py b/pysteps/tests/test_nowcasts_anvil.py index 706f588b1..f9259fc93 100644 --- a/pysteps/tests/test_nowcasts_anvil.py +++ b/pysteps/tests/test_nowcasts_anvil.py @@ -31,11 +31,8 @@ def test_default_anvil_norain(): "time": ( ["time"], np.arange(4.0) * 5.0, - { - "long_name": "forecast time", - "units": "seconds since 1970-01-01 00:00:00", - "stepsize": 5.0, - }, + {"long_name": "forecast time", "stepsize": 5.0}, + {"units": "seconds since 1970-01-01 00:00:00"}, ), "y": ( ["y"], diff --git a/pysteps/tests/test_nowcasts_sseps.py b/pysteps/tests/test_nowcasts_sseps.py index e8ddfedd7..ee7a6b885 100644 --- a/pysteps/tests/test_nowcasts_sseps.py +++ b/pysteps/tests/test_nowcasts_sseps.py @@ -45,11 +45,8 @@ def test_default_sseps_norain(): "time": ( ["time"], np.arange(3.0) * 5.0, - { - "long_name": "forecast time", - "units": "seconds since 1970-01-01 00:00:00", - "stepsize": 5.0, - }, + {"long_name": "forecast time", "stepsize": 5.0}, + {"units": "seconds since 1970-01-01 00:00:00"}, ), "y": ( ["y"], diff --git a/pysteps/tests/test_utils_conversion.py b/pysteps/tests/test_utils_conversion.py index bdf1fa42f..ecfc00ff2 100644 --- a/pysteps/tests/test_utils_conversion.py +++ b/pysteps/tests/test_utils_conversion.py @@ -17,7 +17,6 @@ np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -33,7 +32,6 @@ np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -51,7 +49,6 @@ np.array([1.0]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -67,7 +64,6 @@ np.array([12.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 12.0, "zerovalue": 12.0, @@ -101,7 +97,6 @@ np.array([1.25892541]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.25892541, "zerovalue": 0.0, @@ -135,7 +130,6 @@ np.array([15.10710494]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 15.10710494, "zerovalue": 0.0, @@ -169,7 +163,6 @@ np.array([0.04210719]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 0.04210719, "zerovalue": 0.0, @@ -203,7 +196,6 @@ np.array([2.71828183]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 2.71828183, "zerovalue": 0.0, @@ -237,7 +229,6 @@ np.array([32.61938194]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 32.61938194, "zerovalue": 0.0, @@ -271,7 +262,6 @@ np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -305,7 +295,6 @@ np.array([12.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 12.0, "zerovalue": 12.0, @@ -335,7 +324,6 @@ def test_to_rainrate(dataset, expected): np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -351,7 +339,6 @@ def test_to_rainrate(dataset, expected): np.array([0.08333333]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 0.08333333, "zerovalue": 0.08333333, @@ -369,7 +356,6 @@ def test_to_rainrate(dataset, expected): np.array([1.0]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -385,7 +371,6 @@ def test_to_rainrate(dataset, expected): np.array([1.0]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -419,7 +404,6 @@ def test_to_rainrate(dataset, expected): np.array([0.10491045]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 0.10491045, "zerovalue": 0.0, @@ -453,7 +437,6 @@ def test_to_rainrate(dataset, expected): np.array([1.25892541]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 1.25892541, "zerovalue": 0.0, @@ -487,7 +470,6 @@ def test_to_rainrate(dataset, expected): np.array([0.00350893]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 0.00350893, "zerovalue": 0.0, @@ -521,7 +503,6 @@ def test_to_rainrate(dataset, expected): np.array([0.22652349]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 0.22652349, "zerovalue": 0.0, @@ -555,7 +536,6 @@ def test_to_rainrate(dataset, expected): np.array([0.08333333]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 0.08333333, "zerovalue": 0.08333333, @@ -589,7 +569,6 @@ def test_to_rainrate(dataset, expected): np.array([1.0]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -619,7 +598,6 @@ def test_to_raindepth(dataset, expected): np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, @@ -653,7 +631,6 @@ def test_to_raindepth(dataset, expected): np.array([1.0]), { "units": "mm", - "transform": None, "accutime": 5, "threshold": 1.0, "zerovalue": 1.0, diff --git a/pysteps/tests/test_utils_dimension.py b/pysteps/tests/test_utils_dimension.py index 038b725a0..0be35eb12 100644 --- a/pysteps/tests/test_utils_dimension.py +++ b/pysteps/tests/test_utils_dimension.py @@ -19,7 +19,6 @@ "zerovalue": 0, "yorigin": "lower", "unit": "mm/h", - "transform": None, "accutime": 5, "threshold": 1.0, "projection": "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0", diff --git a/pysteps/tests/test_utils_reprojection.py b/pysteps/tests/test_utils_reprojection.py index 2fcb4f7aa..58a1231cf 100644 --- a/pysteps/tests/test_utils_reprojection.py +++ b/pysteps/tests/test_utils_reprojection.py @@ -57,12 +57,12 @@ def build_precip_dataset( name=precip_var_name, attrs={ "units": units, - "transform": transform_attr, "accutime": float(accutime_min), "threshold": float(threshold), "zerovalue": float(zerovalue), "zr_a": float(zr_a), "zr_b": float(zr_b), + **({"transform": transform_attr} if transform_attr is not None else {}), }, ) diff --git a/pysteps/tests/test_utils_transformation.py b/pysteps/tests/test_utils_transformation.py index 29e6e639c..0b73b0c72 100644 --- a/pysteps/tests/test_utils_transformation.py +++ b/pysteps/tests/test_utils_transformation.py @@ -16,7 +16,6 @@ np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": np.e, "zerovalue": 0, @@ -74,7 +73,6 @@ np.array([np.exp(1.0)]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": np.e, "zerovalue": 0, @@ -92,7 +90,6 @@ np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": np.e, "zerovalue": 0, @@ -150,7 +147,6 @@ np.array([0.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": np.e, "zerovalue": 0, @@ -185,7 +181,6 @@ def test_boxcox_transform(dataset, Lambda, threshold, zerovalue, inverse, expect np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1, "zerovalue": 1, @@ -241,7 +236,6 @@ def test_boxcox_transform(dataset, Lambda, threshold, zerovalue, inverse, expect np.array([1.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 1, "zerovalue": 0, @@ -273,7 +267,6 @@ def test_dB_transform(dataset, threshold, zerovalue, inverse, expected): np.array([1.0, 2.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 0, "zerovalue": 0, @@ -320,7 +313,6 @@ def test_NQ_transform(dataset, inverse, expected): np.array([1.0, 4.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 4, "zerovalue": 0, @@ -372,7 +364,6 @@ def test_NQ_transform(dataset, inverse, expected): np.array([1.0, 4.0]), { "units": "mm/h", - "transform": None, "accutime": 5, "threshold": 4, "zerovalue": 0, diff --git a/pysteps/utils/conversion.py b/pysteps/utils/conversion.py index 92e08a450..f4e77d173 100644 --- a/pysteps/utils/conversion.py +++ b/pysteps/utils/conversion.py @@ -92,7 +92,7 @@ def to_rainrate(dataset: xr.Dataset, zr_a=None, zr_b=None): precip_var = dataset.attrs["precip_var"] metadata = dataset[precip_var].attrs - if metadata["transform"] is not None: + if "transform" in metadata: if metadata["transform"] == "dB": dataset = transformation.dB_transform(dataset, inverse=True) @@ -182,7 +182,7 @@ def to_raindepth(dataset: xr.Dataset, zr_a=None, zr_b=None): precip_var = dataset.attrs["precip_var"] metadata = dataset[precip_var].attrs - if metadata["transform"] is not None: + if "transform" in metadata: if metadata["transform"] == "dB": dataset = transformation.dB_transform(dataset, inverse=True) @@ -202,7 +202,7 @@ def to_raindepth(dataset: xr.Dataset, zr_a=None, zr_b=None): metadata = dataset[precip_var].attrs precip_data = dataset[precip_var].values - if metadata["units"] == "mm" and metadata["transform"] is None: + if metadata["units"] == "mm" and "transform" not in metadata: pass elif metadata["units"] == "mm/h": @@ -272,7 +272,7 @@ def to_reflectivity(dataset: xr.Dataset, zr_a=None, zr_b=None): precip_var = dataset.attrs["precip_var"] metadata = dataset[precip_var].attrs - if metadata["transform"] is not None: + if "transform" in metadata: if metadata["transform"] == "dB": dataset = transformation.dB_transform(dataset, inverse=True) diff --git a/pysteps/utils/transformation.py b/pysteps/utils/transformation.py index 3e48fe0d8..49b4e4fe8 100644 --- a/pysteps/utils/transformation.py +++ b/pysteps/utils/transformation.py @@ -79,7 +79,7 @@ def boxcox_transform( precip_data = dataset[precip_var].values if not inverse: - if metadata["transform"] == "BoxCox": + if "transform" in metadata and metadata["transform"] == "BoxCox": return dataset if Lambda is None: @@ -131,7 +131,7 @@ def boxcox_transform( precip_data[precip_data < threshold] = zerovalue - metadata["transform"] = None + del metadata["transform"] metadata["zerovalue"] = zerovalue metadata["threshold"] = threshold @@ -174,7 +174,7 @@ def dB_transform( # to dB units if not inverse: - if metadata["transform"] == "dB": + if "transform" in metadata and metadata["transform"] == "dB": return dataset if threshold is None: @@ -209,7 +209,7 @@ def dB_transform( threshold = 10.0 ** (threshold / 10.0) precip_data[precip_data < threshold] = zerovalue - metadata["transform"] = None + del metadata["transform"] metadata["threshold"] = threshold metadata["zerovalue"] = zerovalue @@ -295,7 +295,7 @@ def NQ_transform(dataset: xr.Dataset, inverse: bool = False, **kwargs) -> xr.Dat else: f = metadata.pop("inqt") precip_data__ = f(precip_data_) - metadata["transform"] = None + del metadata["transform"] metadata["zerovalue"] = precip_data__.min() metadata["threshold"] = precip_data__[precip_data__ > precip_data__.min()].min() @@ -341,7 +341,7 @@ def sqrt_transform(dataset: xr.Dataset, inverse: bool = False, **kwargs) -> xr.D # inverse sqrt transform precip_data = precip_data**2 - metadata["transform"] = None + del metadata["transform"] metadata["zerovalue"] = metadata["zerovalue"] ** 2 metadata["threshold"] = metadata["threshold"] ** 2 diff --git a/pysteps/xarray_helpers.py b/pysteps/xarray_helpers.py index 2ad5b0557..403c77247 100644 --- a/pysteps/xarray_helpers.py +++ b/pysteps/xarray_helpers.py @@ -85,6 +85,7 @@ def convert_input_to_xarray_dataset( quality: np.ndarray | None, metadata: dict[str, str | float | None], startdate: datetime | None = None, + timestep: int | None = None, ) -> xr.Dataset: """ Read a precip, quality, metadata tuple as returned by the importers @@ -94,15 +95,18 @@ def convert_input_to_xarray_dataset( Parameters ---------- precip: array - 2D array containing imported precipitation data. + ND array containing imported precipitation data. quality: array, None - 2D array containing the quality values of the imported precipitation + ND array containing the quality values of the imported precipitation data, can be None. metadata: dict Metadata dictionary containing the attributes described in the documentation of :py:mod:`pysteps.io.importers`. startdate: datetime, None Datetime object containing the start date and time for the nowcast + timestep: int, None + The timestep in seconds between 2 consecutive fields, mandatory if + the precip has 3 or more dimensions Returns ------- @@ -122,6 +126,8 @@ def convert_input_to_xarray_dataset( if startdate is None: raise Exception("startdate missing") + if timestep is None: + raise Exception("timestep missing") elif precip.ndim == 3: timesteps, h, w = precip.shape @@ -129,6 +135,8 @@ def convert_input_to_xarray_dataset( if startdate is None: raise Exception("startdate missing") + if timestep is None: + raise Exception("timestep missing") elif precip.ndim == 2: h, w = precip.shape @@ -268,8 +276,12 @@ def convert_input_to_xarray_dataset( coords["time"] = ( ["time"], - list(range(1, timesteps + 1, 1)), - {"long_name": "forecast time", "units": "seconds since %s" % startdate_str}, + [ + startdate + timedelta(seconds=float(second)) + for second in np.arange(timesteps) * timestep + ], + {"long_name": "forecast time", "stepsize": timestep}, + {"units": "seconds since %s" % startdate_str}, ) if grid_mapping_var_name is not None: coords[grid_mapping_name] = ( From 8184755dc2f3724c25e1d918f38fd206aa6d0a35 Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Thu, 11 Sep 2025 14:54:30 +0200 Subject: [PATCH 2/8] fix test blending utils --- pysteps/blending/utils.py | 40 ------ pysteps/io/importers.py | 3 +- pysteps/tests/test_blending_utils.py | 181 +++++++++++++-------------- 3 files changed, 91 insertions(+), 133 deletions(-) diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py index 49e95619b..2e23d6299 100644 --- a/pysteps/blending/utils.py +++ b/pysteps/blending/utils.py @@ -47,46 +47,6 @@ CV2_IMPORTED = False -def stack_cascades(R_d, donorm=True): - """Stack the given cascades into a larger array. - - Parameters - ---------- - R_d : dict - Dictionary containing a list of cascades obtained by calling a method - implemented in pysteps.cascade.decomposition. - donorm : bool - If True, normalize the cascade levels before stacking. - - Returns - ------- - out : tuple - A three-element tuple containing a four-dimensional array of stacked - cascade levels and arrays of mean values and standard deviations for each - cascade level. - """ - R_c = [] - mu_c = [] - sigma_c = [] - - for cascade in R_d: - R_ = [] - R_i = cascade["cascade_levels"] - n_levels = R_i.shape[0] - mu_ = np.asarray(cascade["means"]) - sigma_ = np.asarray(cascade["stds"]) - if donorm: - for j in range(n_levels): - R__ = (R_i[j, :, :] - mu_[j]) / sigma_[j] - R_.append(R__) - else: - R_ = R_i - R_c.append(np.stack(R_)) - mu_c.append(mu_) - sigma_c.append(sigma_) - return np.stack(R_c), np.stack(mu_c), np.stack(sigma_c) - - def blend_cascades(cascades_norm, weights): """Calculate blended normalized cascades using STEPS weights following eq. 10 in :cite:`BPS2006`. diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index cd45df034..3706321e2 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -215,6 +215,7 @@ from pysteps.decorators import postprocess_import from pysteps.exceptions import DataModelError, MissingOptionalDependency from pysteps.utils import aggregate_fields +from pysteps.xarray_helpers import convert_input_to_xarray_dataset try: from osgeo import gdal, gdalconst, osr @@ -550,7 +551,7 @@ def _import_mrms_grib(filename, extent=None, window_size=4, **kwargs): cartesian_unit="degrees", ) - return precip, None, metadata + return convert_input_to_xarray_dataset(precip, None, metadata) @postprocess_import() diff --git a/pysteps/tests/test_blending_utils.py b/pysteps/tests/test_blending_utils.py index fccbe9067..c7a455340 100644 --- a/pysteps/tests/test_blending_utils.py +++ b/pysteps/tests/test_blending_utils.py @@ -13,10 +13,7 @@ blend_optical_flows, compute_smooth_dilated_mask, preprocess_and_store_nwp_data, - decompose_NWP, - load_NWP, recompose_cascade, - stack_cascades, ) from pysteps.utils.check_norain import check_norain from pysteps.xarray_helpers import convert_input_to_xarray_dataset @@ -100,27 +97,11 @@ # Set the testing arguments # Test function arguments -utils_arg_names = ( - "precip_nwp_dataset", - "nwp_model", - "issue_times", - "timestep", - "n_timesteps", - "shape", - "weights", -) +utils_arg_names = ("precip_nwp_dataset", "nwp_model", "issue_times", "weights") # Test function values utils_arg_values = [ - ( - precip_nwp_dataset, - "test", - [issue_time_first, issue_time_second], - 5.0, - 3, - precip_nwp_dataset[nwp_precip_var].values.shape[1:], - weights, - ) + (precip_nwp_dataset, "test", [issue_time_first, issue_time_second], weights) ] smoothing_arg_names = ( @@ -145,15 +126,7 @@ ### @pytest.mark.parametrize(utils_arg_names, utils_arg_values) # The test function to be used -def test_blending_utils( - precip_nwp_dataset, - nwp_model, - issue_times, - timestep, - n_timesteps, - shape, - weights, -): +def test_blending_utils(precip_nwp_dataset, nwp_model, issue_times, weights): """Tests if all blending utils functions behave correctly.""" # First, make the output path if it does not exist yet @@ -204,87 +177,110 @@ def test_blending_utils( ### # Now check if files load correctly for two different issue times ### - with xr.open_dataset(preprocessed_file) as nwp_file: - precip_decomposed_nwp_first, v_nwp_first = load_NWP( - input_nc_path_decomp=os.path.join(decomp_file), - input_path_velocities=os.path.join(preprocessed_file), - start_time=issue_times[0], - n_timesteps=n_timesteps, - ) - - precip_decomposed_nwp_second, v_nwp_second = load_NWP( - input_nc_path_decomp=os.path.join(decomp_file), - input_path_velocities=os.path.join(preprocessed_file), - start_time=issue_times[1], - n_timesteps=n_timesteps, - ) + with xr.open_dataset(preprocessed_file) as nwp_file_dataset: + nwp_file_dataset_first = nwp_file_dataset.sel(time=issue_times[0]) + nwp_file_dataset_second = nwp_file_dataset.sel(time=issue_times[1]) + precip_var = nwp_file_dataset.attrs["precip_var"] # Check, for a sample, if the stored motion fields are as expected assert_array_almost_equal( - v_nwp_first[1], - oflow_method(precip_nwp_dataset[0:2, :, :]), + nwp_file_dataset_first.velocity_x.values, + oflow_method(precip_nwp_dataset.isel(time=slice(0, 2))).velocity_x.values, decimal=3, err_msg="Stored motion field of first forecast not equal to expected motion field", ) assert_array_almost_equal( - v_nwp_second[1], - oflow_method(precip_nwp_dataset[3:5, :, :]), + nwp_file_dataset_first.velocity_y.values, + oflow_method(precip_nwp_dataset.isel(time=slice(0, 2))).velocity_y.values, decimal=3, - err_msg="Stored motion field of second forecast not equal to expected motion field", + err_msg="Stored motion field of first forecast not equal to expected motion field", ) - - ### - # Stack the cascades - ### - precip_decomposed_first_stack, mu_first_stack, sigma_first_stack = ( - stack_cascades(R_d=precip_decomposed_nwp_first, donorm=False) + assert_array_almost_equal( + nwp_file_dataset_second.velocity_x.values, + oflow_method(precip_nwp_dataset.isel(time=slice(3, 5))).velocity_x.values, + decimal=3, + err_msg="Stored motion field of second forecast not equal to expected motion field", ) - - print(precip_decomposed_nwp_first) - print(precip_decomposed_first_stack) - print(mu_first_stack) - - ( - precip_decomposed_second_stack, - mu_second_stack, - sigma_second_stack, - ) = stack_cascades(R_d=precip_decomposed_nwp_second, donorm=False) - - # Check if the array shapes are still correct - assert precip_decomposed_first_stack.shape == ( - n_timesteps + 1, - 8, - shape[0], - shape[1], + assert_array_almost_equal( + nwp_file_dataset_second.velocity_y.values, + oflow_method(precip_nwp_dataset.isel(time=slice(3, 5))).velocity_y.values, + decimal=3, + err_msg="Stored motion field of second forecast not equal to expected motion field", ) - assert mu_first_stack.shape == (n_timesteps + 1, 8) - assert sigma_first_stack.shape == (n_timesteps + 1, 8) ### # Blend the cascades ### precip_decomposed_blended = blend_cascades( cascades_norm=np.stack( - (precip_decomposed_first_stack[0], precip_decomposed_second_stack[0]) + ( + nwp_file_dataset_first[precip_var].values, + nwp_file_dataset_second[precip_var].values, + ) ), weights=weights, ) - assert precip_decomposed_blended.shape == precip_decomposed_first_stack[0].shape + assert ( + precip_decomposed_blended.shape + == nwp_file_dataset_first[precip_var].values.shape + ) ### # Blend the optical flow fields ### v_nwp_blended = blend_optical_flows( - flows=np.stack((v_nwp_first[1], v_nwp_second[1])), weights=weights[:, 1] + flows=np.stack( + ( + np.array( + [ + nwp_file_dataset_first.velocity_x.values, + nwp_file_dataset_first.velocity_y.values, + ] + ), + np.array( + [ + nwp_file_dataset_second.velocity_x.values, + nwp_file_dataset_second.velocity_y.values, + ] + ), + ) + ), + weights=weights[:, 1], ) - assert v_nwp_blended.shape == v_nwp_first[1].shape + assert ( + v_nwp_blended.shape + == np.array( + [ + nwp_file_dataset_first.velocity_x.values, + nwp_file_dataset_first.velocity_y.values, + ] + ).shape + ) + assert_array_almost_equal( + v_nwp_blended[0], + ( + oflow_method( + precip_nwp_dataset.isel(time=slice(0, 2)) + ).velocity_x.values + + oflow_method( + precip_nwp_dataset.isel(time=slice(3, 5)) + ).velocity_x.values + ) + / 2, + decimal=3, + err_msg="Blended motion field does not equal average of the two motion fields", + ) assert_array_almost_equal( - v_nwp_blended, + v_nwp_blended[1], ( - oflow_method(precip_nwp_dataset[0:2, :, :]) - + oflow_method(precip_nwp_dataset[3:5, :, :]) + oflow_method( + precip_nwp_dataset.isel(time=slice(0, 2)) + ).velocity_y.values + + oflow_method( + precip_nwp_dataset.isel(time=slice(3, 5)) + ).velocity_y.values ) / 2, decimal=3, @@ -295,30 +291,30 @@ def test_blending_utils( # Recompose the fields (the non-blended fields are used for this here) ### precip_recomposed_first = recompose_cascade( - combined_cascade=precip_decomposed_first_stack[0], - combined_mean=mu_first_stack[0], - combined_sigma=sigma_first_stack[0], + combined_cascade=nwp_file_dataset_first[precip_var].values, + combined_mean=nwp_file_dataset_first["means"].values, + combined_sigma=nwp_file_dataset_first["stds"].values, ) precip_recomposed_second = recompose_cascade( - combined_cascade=precip_decomposed_second_stack[0], - combined_mean=mu_second_stack[0], - combined_sigma=sigma_second_stack[0], + combined_cascade=nwp_file_dataset_second[precip_var].values, + combined_mean=nwp_file_dataset_second["means"].values, + combined_sigma=nwp_file_dataset_second["stds"].values, ) assert_array_almost_equal( precip_recomposed_first, - precip_nwp_dataset[0, :, :], + precip_nwp_dataset.isel(time=0)[nwp_precip_var].values, decimal=3, err_msg="Recomposed field of first forecast does not equal original field", ) assert_array_almost_equal( precip_recomposed_second, - precip_nwp_dataset[3, :, :], + precip_nwp_dataset.isel(time=3)[nwp_precip_var].values, decimal=3, err_msg="Recomposed field of second forecast does not equal original field", ) - precip_arr = precip_nwp_dataset + precip_arr = precip_nwp_dataset[nwp_precip_var].values # rainy fraction is 0.005847 assert not check_norain(precip_arr, win_fun=None) assert not check_norain( @@ -364,8 +360,9 @@ def test_blending_smoothing_utils( non_linear_growth_kernel_sizes, ): # First add some nans to indicate a mask - precip_nwp_dataset[:, 0:100, 0:100] = np.nan - nan_indices = np.isnan(precip_nwp_dataset[0]) + nwp_precip_var = precip_nwp_dataset.attrs["precip_var"] + precip_nwp_dataset[nwp_precip_var].data[:, 0:100, 0:100] = np.nan + nan_indices = np.isnan(precip_nwp_dataset[nwp_precip_var].values[0]) new_mask = compute_smooth_dilated_mask( nan_indices, max_padding_size_in_px=max_padding_size_in_px, From 2e356240d5610756812ed2d5c54fe450ac956a7e Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Mon, 15 Sep 2025 08:30:44 +0200 Subject: [PATCH 3/8] fix steps blending tests --- pysteps/blending/steps.py | 115 +++++++++--- pysteps/blending/utils.py | 256 +++++++++++++++------------ pysteps/cascade/decomposition.py | 2 +- pysteps/io/importers.py | 2 +- pysteps/tests/test_blending_steps.py | 136 +++++++------- pysteps/tests/test_blending_utils.py | 4 +- pysteps/tests/test_cascade.py | 12 +- pysteps/xarray_helpers.py | 3 +- 8 files changed, 307 insertions(+), 223 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 7a6a2543b..0ba7bbea3 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -42,6 +42,7 @@ calculate_weights_spn blend_means_sigmas """ +from datetime import datetime import math import time from copy import copy, deepcopy @@ -49,6 +50,7 @@ from multiprocessing.pool import ThreadPool import numpy as np +import xarray as xr from scipy.linalg import inv from scipy.ndimage import binary_dilation, generate_binary_structure, iterate_structure @@ -57,6 +59,7 @@ from pysteps.postprocessing import probmatching from pysteps.timeseries import autoregression, correlation from pysteps.utils.check_norain import check_norain +from pysteps.xarray_helpers import convert_output_to_xarray_dataset try: import dask @@ -412,20 +415,76 @@ class StepsBlendingState: class StepsBlendingNowcaster: def __init__( self, - precip, - precip_models, - velocity, - velocity_models, + radar_dataset: xr.Dataset, + model_dataset: xr.Dataset, time_steps, - issue_time, + issue_time: datetime, steps_blending_config: StepsBlendingConfig, ): """Initializes the StepsBlendingNowcaster with inputs and configurations.""" # Store inputs - self.__precip = precip - self.__precip_models = precip_models - self.__velocity = velocity - self.__velocity_models = velocity_models + radar_precip_var = radar_dataset.attrs["precip_var"] + model_precip_var = model_dataset.attrs["precip_var"] + if issue_time != radar_dataset.time.isel(time=-1).values.astype( + "datetime64[us]" + ).astype(datetime): + raise ValueError( + "Issue time should be equal to last timestep in radar dataset" + ) + time_stepsize_seconds = radar_dataset.time.attrs["stepsize"] + if isinstance(time_steps, list): + # XR: validates this works or just remove the subtimestep stuff + timesteps_seconds = ( + np.array(list(range(time_steps[-1] + 1))) * time_stepsize_seconds + ) + else: + timesteps_seconds = ( + np.array(list(range(time_steps + 1))) * time_stepsize_seconds + ) + model_times = radar_dataset.time.isel( + time=-1 + ).values + timesteps_seconds.astype("timedelta64[s]") + model_dataset = model_dataset.sel(time=model_times) + + self.__precip = radar_dataset[radar_precip_var].values + # XR: don't extract to dict but pass dataset + if model_dataset[model_precip_var].ndim == 5: + self.__precip_models = np.array( + [ + [ + { + "cascade_levels": model_dataset[model_precip_var] + .sel(time=time, ens_number=ens_number) + .values, + "means": model_dataset["means"] + .sel(time=time, ens_number=ens_number) + .values, + "stds": model_dataset["stds"] + .sel(time=time, ens_number=ens_number) + .values, + "domain": model_dataset[model_precip_var].attrs["domain"], + "normalized": model_dataset[model_precip_var].attrs[ + "normalized" + ], + "compact_output": model_dataset[model_precip_var].attrs[ + "compact_output" + ], + } + for time in model_dataset.time + ] + for ens_number in model_dataset.ens_number + ] + ) + else: + self.__precip_models = model_dataset[model_precip_var].values + self.__velocity = np.array( + [radar_dataset["velocity_x"].values, radar_dataset["velocity_y"].values] + ) + self.__velocity_models = np.array( + [model_dataset["velocity_x"].values, model_dataset["velocity_y"].values] + ).transpose(1, 2, 0, 3, 4) + self.__original_timesteps = time_steps + self.__input_radar_dataset = radar_dataset self.__timesteps = time_steps self.__issuetime = issue_time @@ -447,6 +506,7 @@ def compute_forecast(self): Parameters ---------- + # XR: fix docstring precip: array-like Array of shape (ar_order+1,m,n) containing the input precipitation fields ordered by timestamp from oldest to newest. The time steps between the @@ -545,7 +605,7 @@ def compute_forecast(self): # Determine if rain is present in both radar and NWP fields if self.__params.zero_precip_radar and self.__params.zero_precip_model_fields: - return self.__zero_precipitation_forecast() + result = self.__zero_precipitation_forecast() else: # Prepare the data for the zero precipitation radar case and initialize the noise correctly if self.__params.zero_precip_radar: @@ -572,16 +632,20 @@ def compute_forecast(self): for j in range(self.__config.n_ens_members) ] ) - if self.__config.measure_time: - return ( - self.__state.final_blended_forecast, - self.__init_time, - self.__mainloop_time, - ) - else: - return self.__state.final_blended_forecast + result = self.__state.final_blended_forecast else: return None + result_dataset = convert_output_to_xarray_dataset( + self.__input_radar_dataset, self.__original_timesteps, result + ) + if self.__config.measure_time: + return ( + result_dataset, + self.__init_time, + self.__mainloop_time, + ) + else: + return result_dataset def __blended_nowcast_main_loop(self): """ @@ -2817,10 +2881,8 @@ def __measure_time(self, label, start_time): def forecast( - precip, - precip_models, - velocity, - velocity_models, + radar_dataset, + model_dataset, timesteps, timestep, issuetime, @@ -2864,6 +2926,7 @@ def forecast( Parameters ---------- + # XR: fix docstring precip: array-like Array of shape (ar_order+1,m,n) containing the input precipitation fields ordered by timestamp from oldest to newest. The time steps between the @@ -3182,13 +3245,7 @@ def forecast( """ # Create an instance of the new class with all the provided arguments blended_nowcaster = StepsBlendingNowcaster( - precip, - precip_models, - velocity, - velocity_models, - timesteps, - issuetime, - blending_config, + radar_dataset, model_dataset, timesteps, issuetime, blending_config ) forecast_steps_nowcast = blended_nowcaster.compute_forecast() diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py index 2e23d6299..4cae9c28a 100644 --- a/pysteps/blending/utils.py +++ b/pysteps/blending/utils.py @@ -19,7 +19,7 @@ """ import datetime -from typing import Any +from typing import Any, Callable import warnings from pathlib import Path @@ -54,24 +54,24 @@ def blend_cascades(cascades_norm, weights): Parameters ---------- cascades_norm : array-like - Array of shape [number_components + 1, scale_level, ...] - with the cascade for each component (NWP, nowcasts, noise) and scale level, - obtained by calling a method implemented in pysteps.blending.utils.stack_cascades + Array of shape [number_components + 1, scale_level, ...] + with the cascade for each component (NWP, nowcasts, noise) and scale level, + obtained by calling a method implemented in pysteps.blending.utils.stack_cascades weights : array-like - An array of shape [number_components + 1, scale_level, ...] - containing the weights to be used in this routine - for each component plus noise, scale level, and optionally [y, x] - dimensions, obtained by calling a method implemented in - pysteps.blending.steps.calculate_weights + An array of shape [number_components + 1, scale_level, ...] + containing the weights to be used in this routine + for each component plus noise, scale level, and optionally [y, x] + dimensions, obtained by calling a method implemented in + pysteps.blending.steps.calculate_weights Returns ------- combined_cascade : array-like - An array of shape [scale_level, y, x] - containing per scale level (cascade) the weighted combination of - cascades from multiple components (NWP, nowcasts and noise) to be used - in STEPS blending. + An array of shape [scale_level, y, x] + containing per scale level (cascade) the weighted combination of + cascades from multiple components (NWP, nowcasts and noise) to be used + in STEPS blending. """ # check inputs if isinstance(cascades_norm, (list, tuple)): @@ -115,18 +115,18 @@ def recompose_cascade(combined_cascade, combined_mean, combined_sigma): Parameters ---------- combined_cascade : array-like - An array of shape [scale_level, y, x] - containing per scale level (cascade) the weighted combination of - cascades from multiple components (NWP, nowcasts and noise) to be used - in STEPS blending. + An array of shape [scale_level, y, x] + containing per scale level (cascade) the weighted combination of + cascades from multiple components (NWP, nowcasts and noise) to be used + in STEPS blending. combined_mean : array-like - An array of shape [scale_level, ...] - similar to combined_cascade, but containing the normalization parameter - mean. + An array of shape [scale_level, ...] + similar to combined_cascade, but containing the normalization parameter + mean. combined_sigma : array-like - An array of shape [scale_level, ...] - similar to combined_cascade, but containing the normalization parameter - standard deviation. + An array of shape [scale_level, ...] + similar to combined_cascade, but containing the normalization parameter + standard deviation. Returns ------- @@ -151,17 +151,17 @@ def blend_optical_flows(flows, weights): Parameters ---------- flows : array-like - A stack of multiple advenction fields having shape - (S, 2, m, n), where flows[N, :, :, :] contains the motion vectors - for source N. - Advection fields for each source can be obtanined by - calling any of the methods implemented in - pysteps.motion and then stack all together + A stack of multiple advenction fields having shape + (S, 2, m, n), where flows[N, :, :, :] contains the motion vectors + for source N. + Advection fields for each source can be obtanined by + calling any of the methods implemented in + pysteps.motion and then stack all together weights : array-like - An array of shape [number_sources] - containing the weights to be used to combine - the advection fields of each source. - weights are modified to make their sum equal to one. + An array of shape [number_sources] + containing the weights to be used to combine + the advection fields of each source. + weights are modified to make their sum equal to one. Returns ------- out: ndarray @@ -204,7 +204,7 @@ def blend_optical_flows(flows, weights): def decompose_NWP( precip_nwp_dataset: xr.Dataset, - num_cascade_levels=8, + num_cascade_levels=6, num_workers=1, decomp_method="fft", fft_method="numpy", @@ -219,57 +219,57 @@ def decompose_NWP( Parameters ---------- R_NWP: array-like - Array of dimension (n_timesteps, x, y) containing the precipitation forecast - from some NWP model. + Array of dimension (n_timesteps, x, y) containing the precipitation forecast + from some NWP model. NWP_model: str - The name of the NWP model + The name of the NWP model analysis_time: numpy.datetime64 - The analysis time of the NWP forecast. The analysis time is assumed to be a - numpy.datetime64 type as imported by the pysteps importer + The analysis time of the NWP forecast. The analysis time is assumed to be a + numpy.datetime64 type as imported by the pysteps importer output_path: str - The location where to save the file with the NWP cascade. Defaults to the - path_workdir specified in the rcparams file. + The location where to save the file with the NWP cascade. Defaults to the + path_workdir specified in the rcparams file. num_cascade_levels: int, optional - The number of frequency bands to use. Must be greater than 2. Defaults to 8. + The number of frequency bands to use. Must be greater than 2. Defaults to 8. num_workers: int, optional - The number of workers to use for parallel computation. Applicable if dask - is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it - is advisable to disable OpenMP by setting the environment variable - OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous - threads. + The number of workers to use for parallel computation. Applicable if dask + is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it + is advisable to disable OpenMP by setting the environment variable + OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous + threads. Other Parameters ---------------- decomp_method: str, optional - A string defining the decomposition method to use. Defaults to "fft". + A string defining the decomposition method to use. Defaults to "fft". fft_method: str or tuple, optional - A string or a (function,kwargs) tuple defining the FFT method to use - (see :py:func:`pysteps.utils.interface.get_method`). - Defaults to "numpy". This option is not used if input_domain and - output_domain are both set to "spectral". + A string or a (function,kwargs) tuple defining the FFT method to use + (see :py:func:`pysteps.utils.interface.get_method`). + Defaults to "numpy". This option is not used if input_domain and + output_domain are both set to "spectral". domain: {"spatial", "spectral"}, optional - If "spatial", the output cascade levels are transformed back to the - spatial domain by using the inverse FFT. If "spectral", the cascade is - kept in the spectral domain. Defaults to "spatial". + If "spatial", the output cascade levels are transformed back to the + spatial domain by using the inverse FFT. If "spectral", the cascade is + kept in the spectral domain. Defaults to "spatial". normalize: bool, optional - If True, normalize the cascade levels to zero mean and unit variance. - Requires that compute_stats is True. Implies that compute_stats is True. - Defaults to False. + If True, normalize the cascade levels to zero mean and unit variance. + Requires that compute_stats is True. Implies that compute_stats is True. + Defaults to False. compute_stats: bool, optional - If True, the output dictionary contains the keys "means" and "stds" - for the mean and standard deviation of each output cascade level. - Defaults to False. + If True, the output dictionary contains the keys "means" and "stds" + for the mean and standard deviation of each output cascade level. + Defaults to False. compact_output: bool, optional - Applicable if output_domain is "spectral". If set to True, only the - parts of the Fourier spectrum with non-negligible filter weights are - stored. Defaults to False. + Applicable if output_domain is "spectral". If set to True, only the + parts of the Fourier spectrum with non-negligible filter weights are + stored. Defaults to False. Returns ------- xarray.Dataset - The same dataset as was passed in but with the precip data replaced - with decomposed precip data and means and stds added + The same dataset as was passed in but with the precip data replaced + with decomposed precip data and means and stds added """ nwp_precip_var = precip_nwp_dataset.attrs["precip_var"] @@ -329,12 +329,48 @@ def decompose_NWP( ) precip_nwp_dataset["means"] = (["time", "cascade_level"], means) precip_nwp_dataset["stds"] = (["time", "cascade_level"], stds) + + precip_nwp_dataset[nwp_precip_var].attrs["domain"] = domain + precip_nwp_dataset[nwp_precip_var].attrs["normalized"] = int(normalize) + precip_nwp_dataset[nwp_precip_var].attrs["compact_output"] = int(compact_output) + + return precip_nwp_dataset + + +def _preprocess_nwp_data_single_member( + precip_nwp_dataset: xr.Dataset, + oflow_method: Callable[..., Any], + decompose_nwp: bool, + decompose_kwargs: dict[str, Any] = {}, +) -> xr.Dataset: + nwp_precip_var = precip_nwp_dataset.attrs["precip_var"] + precip_nwp = precip_nwp_dataset[nwp_precip_var].values + + # Get the velocity field per time step + v_nwp_x = np.zeros((precip_nwp.shape[0], precip_nwp.shape[1], precip_nwp.shape[2])) + v_nwp_y = np.zeros((precip_nwp.shape[0], precip_nwp.shape[1], precip_nwp.shape[2])) + # Loop through the timesteps. We need two images to construct a motion + # field, so we can start from timestep 1. + for t in range(1, precip_nwp.shape[0]): + v_nwp_dataset = oflow_method(precip_nwp_dataset.isel(time=slice(t - 1, t + 1))) + v_nwp_x[t] = v_nwp_dataset.velocity_x + v_nwp_y[t] = v_nwp_dataset.velocity_y + + # Make timestep 0 the same as timestep 1. + v_nwp_x[0] = v_nwp_x[1] + v_nwp_y[0] = v_nwp_y[1] + precip_nwp_dataset["velocity_x"] = (["time", "y", "x"], v_nwp_x) + precip_nwp_dataset["velocity_y"] = (["time", "y", "x"], v_nwp_y) + + if decompose_nwp: + precip_nwp_dataset = decompose_NWP(precip_nwp_dataset, **decompose_kwargs) + return precip_nwp_dataset -def preprocess_and_store_nwp_data( +def preprocess_nwp_data( precip_nwp_dataset: xr.Dataset, - oflow_method: str, + oflow_method: Callable[..., Any], nwp_model: str, output_path: str | None, decompose_nwp: bool, @@ -345,22 +381,22 @@ def preprocess_and_store_nwp_data( Parameters ---------- precip_nwp_dataset: xarray.Dataset - xarray Dataset containing the precipitation forecast - from some NWP model. + xarray Dataset containing the precipitation forecast + from some NWP model. oflow_method: {'constant', 'darts', 'lucaskanade', 'proesmans', 'vet'}, optional - An optical flow method from pysteps.motion.get_method. + An optical flow method from pysteps.motion.get_method. nwp_model: str - The name of the NWP model. + The name of the NWP model. output_path: str, optional - The location where to save the netcdf file with the NWP velocity fields. Defaults - to the path_workdir specified in the rcparams file. + The location where to save the netcdf file with the NWP velocity fields. Defaults + to the path_workdir specified in the rcparams file. decompose_nwp: bool - Defines wether or not the NWP needs to be decomposed before storing. This can - be beneficial for performance, because then the decomposition does not need - to happen during the blending anymore. It can however also be detrimental because - this increases the amount of storage and RAM required for the blending. + Defines wether or not the NWP needs to be decomposed before storing. This can + be beneficial for performance, because then the decomposition does not need + to happen during the blending anymore. It can however also be detrimental because + this increases the amount of storage and RAM required for the blending. decompose_kwargs: dict - Keyword arguments passed to the decompose_NWP method. + Keyword arguments passed to the decompose_NWP method. Returns ------- @@ -373,34 +409,32 @@ def preprocess_and_store_nwp_data( "but it is not installed" ) - # Set the output file - analysis_time = precip_nwp_dataset.time.values[0] - output_date = f"{analysis_time.astype('datetime64[us]').astype(datetime.datetime):%Y%m%d%H%M%S}" - outfn = Path(output_path) / f"preprocessed_{nwp_model}_{output_date}.nc" - nwp_precip_var = precip_nwp_dataset.attrs["precip_var"] - precip_nwp = precip_nwp_dataset[nwp_precip_var].values - - # Get the velocity field per time step - v_nwp_x = np.zeros((precip_nwp.shape[0], precip_nwp.shape[1], precip_nwp.shape[2])) - v_nwp_y = np.zeros((precip_nwp.shape[0], precip_nwp.shape[1], precip_nwp.shape[2])) - # Loop through the timesteps. We need two images to construct a motion - # field, so we can start from timestep 1. - for t in range(1, precip_nwp.shape[0]): - v_nwp_dataset = oflow_method(precip_nwp_dataset.isel(time=slice(t - 1, t + 1))) - v_nwp_x[t] = v_nwp_dataset.velocity_x - v_nwp_y[t] = v_nwp_dataset.velocity_y - - # Make timestep 0 the same as timestep 1. - v_nwp_x[0] = v_nwp_x[1] - v_nwp_y[0] = v_nwp_y[1] - precip_nwp_dataset["velocity_x"] = (["time", "y", "x"], v_nwp_x) - precip_nwp_dataset["velocity_y"] = (["time", "y", "x"], v_nwp_y) - - if decompose_nwp: - precip_nwp_dataset = decompose_NWP(precip_nwp_dataset, **decompose_kwargs) + if "ens_number" in precip_nwp_dataset.dims: + preprocessed_nwp_datasets = [] + for ens_number in precip_nwp_dataset["ens_number"]: + preprocessed_nwp_datasets.append( + _preprocess_nwp_data_single_member( + precip_nwp_dataset.sel(ens_number=ens_number), + oflow_method, + decompose_nwp, + decompose_kwargs, + ).expand_dims({"ens_number": [ens_number]}, axis=0) + ) + precip_nwp_dataset = xr.concat(preprocessed_nwp_datasets, "ens_number") + else: + precip_nwp_dataset = _preprocess_nwp_data_single_member( + precip_nwp_dataset, oflow_method, decompose_nwp, decompose_kwargs + ) # Save it as a numpy array - precip_nwp_dataset.to_netcdf(outfn) + if output_path: + analysis_time = precip_nwp_dataset.time.values[0] + output_date = f"{analysis_time.astype('datetime64[us]').astype(datetime.datetime):%Y%m%d%H%M%S}" + outfn = Path(output_path) / f"preprocessed_{nwp_model}_{output_date}.nc" + precip_nwp_dataset.to_netcdf(outfn) + return None + else: + return precip_nwp_dataset def check_norain(precip_arr, precip_thr=None, norain_thr=0.0): @@ -409,17 +443,17 @@ def check_norain(precip_arr, precip_thr=None, norain_thr=0.0): Parameters ---------- precip_arr: array-like - Array containing the input precipitation field + Array containing the input precipitation field precip_thr: float, optional - Specifies the threshold value for minimum observable precipitation intensity. If None, the - minimum value over the domain is taken. + Specifies the threshold value for minimum observable precipitation intensity. If None, the + minimum value over the domain is taken. norain_thr: float, optional - Specifies the threshold value for the fraction of rainy pixels in precip_arr below which we consider there to be - no rain. Standard set to 0.0 + Specifies the threshold value for the fraction of rainy pixels in precip_arr below which we consider there to be + no rain. Standard set to 0.0 Returns ------- norain: bool - Returns whether the fraction of rainy pixels is below the norain_thr threshold. + Returns whether the fraction of rainy pixels is below the norain_thr threshold. """ warnings.warn( diff --git a/pysteps/cascade/decomposition.py b/pysteps/cascade/decomposition.py index e3def7416..c96152c0e 100644 --- a/pysteps/cascade/decomposition.py +++ b/pysteps/cascade/decomposition.py @@ -286,7 +286,7 @@ def recompose_fft(decomp, **kwargs): ): result = np.sum(levels, axis=0) else: - if decomp["compact_output"]: + if decomp["domain"] == "spectral" and decomp["compact_output"]: weight_masks = decomp["weight_masks"] result = np.zeros(weight_masks.shape[1:], dtype=complex) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index 3706321e2..bab54e200 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -66,7 +66,7 @@ +------------------+----------------------------------------------------------+ # XR: Move this to appropriate place -# XR: Add means, stds cascade level +# XR: Add means, stds cascade level and decomposition attributes The data and metadata is then postprocessed into an xarray dataset. This dataset will always contain an x and y dimension, but can be extended with a time dimension and/or an ensemble member dimension over the course of the process. diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index 21064a45f..070cdda7a 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -1,12 +1,14 @@ # -*- coding: utf-8 -*- -import datetime +from datetime import datetime import numpy as np import pytest import pysteps from pysteps import blending, cascade +from pysteps.blending.utils import preprocess_nwp_data +from pysteps.xarray_helpers import convert_input_to_xarray_dataset # fmt:off steps_arg_values = [ @@ -159,13 +161,21 @@ def test_steps_blending( metadata = dict() metadata["unit"] = "mm" - metadata["transformation"] = "dB" + metadata["cartesian_unit"] = "km" metadata["accutime"] = 5.0 - metadata["transform"] = "dB" metadata["zerovalue"] = 0.0 metadata["threshold"] = 0.01 metadata["zr_a"] = 200.0 metadata["zr_b"] = 1.6 + metadata["x1"] = 0.0 + metadata["x2"] = 200.0 + metadata["y1"] = 0.0 + metadata["y2"] = 200.0 + metadata["yorigin"] = "lower" + metadata["institution"] = "test" + metadata["projection"] = ( + "+proj=lcc +lon_0=4.55 +lat_1=50.8 +lat_2=50.8 +a=6371229 +es=0 +lat_0=50.8 +x_0=365950 +y_0=-365950.000000001" + ) # Also set the outdir_path, clim_kwargs and mask_kwargs outdir_path_skill = "./tmp/" @@ -186,102 +196,78 @@ def test_steps_blending( radar_precip[radar_precip < metadata["threshold"]] = 0.0 nwp_precip[nwp_precip < metadata["threshold"]] = 0.0 + radar_dataset = convert_input_to_xarray_dataset( + radar_precip, + None, + metadata, + datetime.fromisoformat("2021-07-04T11:50:00.000000000"), + 300, + ) + model_dataset = convert_input_to_xarray_dataset( + nwp_precip, + None, + metadata, + datetime.fromisoformat("2021-07-04T12:00:00.000000000"), + 300, + ) # convert the data converter = pysteps.utils.get_method("mm/h") - radar_precip, _ = converter(radar_precip, metadata) - nwp_precip, metadata = converter(nwp_precip, metadata) + radar_dataset = converter(radar_dataset) + model_dataset = converter(model_dataset) # transform the data - transformer = pysteps.utils.get_method(metadata["transformation"]) - radar_precip, _ = transformer(radar_precip, metadata) - nwp_precip, metadata = transformer(nwp_precip, metadata) + transformer = pysteps.utils.get_method("dB") + radar_dataset = transformer(radar_dataset) + model_dataset = transformer(model_dataset) + + radar_precip_var = radar_dataset.attrs["precip_var"] + model_precip_var = model_dataset.attrs["precip_var"] # set NaN equal to zero - radar_precip[~np.isfinite(radar_precip)] = metadata["zerovalue"] - nwp_precip[~np.isfinite(nwp_precip)] = metadata["zerovalue"] + radar_dataset[radar_precip_var].data[ + ~np.isfinite(radar_dataset[radar_precip_var].values) + ] = radar_dataset[radar_precip_var].attrs["zerovalue"] + model_dataset[model_precip_var].data[ + ~np.isfinite(model_dataset[model_precip_var].values) + ] = model_dataset[model_precip_var].attrs["zerovalue"] assert ( - np.any(~np.isfinite(radar_precip)) == False + np.any(~np.isfinite(radar_dataset[radar_precip_var].values)) == False ), "There are still infinite values in the input radar data" assert ( - np.any(~np.isfinite(nwp_precip)) == False + np.any(~np.isfinite(model_dataset[radar_precip_var].values)) == False ), "There are still infinite values in the NWP data" ### # Decompose the R_NWP data ### - # Initial decomposition settings - decomp_method, _ = cascade.get_method("fft") - bandpass_filter_method = "gaussian" - precip_shape = radar_precip.shape[1:] - filter_method = cascade.get_method(bandpass_filter_method) - bp_filter = filter_method(precip_shape, n_cascade_levels) - - # If we only use one model: - if nwp_precip.ndim == 3: - nwp_precip = nwp_precip[None, :] - - if decomposed_nwp: - nwp_precip_decomp = [] - # Loop through the n_models - for i in range(nwp_precip.shape[0]): - R_d_models_ = [] - # Loop through the time steps - for j in range(nwp_precip.shape[1]): - R_ = decomp_method( - field=nwp_precip[i, j, :, :], - bp_filter=bp_filter, - normalize=True, - compute_stats=True, - compact_output=True, - ) - R_d_models_.append(R_) - nwp_precip_decomp.append(R_d_models_) - - nwp_precip_decomp = np.array(nwp_precip_decomp) - - assert nwp_precip_decomp.ndim == 2, "Wrong number of dimensions in R_d_models" + radar_precip = radar_dataset[radar_precip_var].values - else: - nwp_precip_decomp = nwp_precip.copy() - - assert nwp_precip_decomp.ndim == 4, "Wrong number of dimensions in R_d_models" + oflow_method = pysteps.motion.get_method("lucaskanade") + nwp_preproc_dataset = preprocess_nwp_data( + model_dataset, + oflow_method, + "test", + None, + decomposed_nwp, + {"num_cascade_levels": n_cascade_levels}, + ) ### # Determine the velocity fields ### - oflow_method = pysteps.motion.get_method("lucaskanade") - radar_velocity = oflow_method(radar_precip) - nwp_velocity = [] - # Loop through the models - for n_model in range(nwp_precip.shape[0]): - # Loop through the timesteps. We need two images to construct a motion - # field, so we can start from timestep 1. Timestep 0 will be the same - # as timestep 0. - _V_NWP_ = [] - for t in range(1, nwp_precip.shape[1]): - V_NWP_ = oflow_method(nwp_precip[n_model, t - 1 : t + 1, :]) - _V_NWP_.append(V_NWP_) - V_NWP_ = None - _V_NWP_ = np.insert(_V_NWP_, 0, _V_NWP_[0], axis=0) - nwp_velocity.append(_V_NWP_) - - nwp_velocity = np.stack(nwp_velocity) - - assert nwp_velocity.ndim == 5, "nwp_velocity must be a five-dimensional array" + radar_dataset_w_velocity = oflow_method(radar_dataset) ### # The nowcasting ### - precip_forecast = blending.steps.forecast( - precip=radar_precip, - precip_models=nwp_precip_decomp, - velocity=radar_velocity, - velocity_models=nwp_velocity, + precip_forecast_dataset = blending.steps.forecast( + radar_dataset=radar_dataset_w_velocity, + model_dataset=nwp_preproc_dataset, timesteps=timesteps, timestep=5.0, - issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), + issuetime=datetime.fromisoformat("2021-07-04T12:00:00.000000000"), n_ens_members=n_ens_members, n_cascade_levels=n_cascade_levels, blend_nwp_members=blend_nwp_members, @@ -315,6 +301,8 @@ def test_steps_blending( mask_kwargs=mask_kwargs, measure_time=False, ) + precip_var_forecast = precip_forecast_dataset.attrs["precip_var"] + precip_forecast = precip_forecast_dataset[precip_var_forecast].values assert precip_forecast.ndim == 4, "Wrong amount of dimensions in forecast output" assert ( @@ -325,7 +313,9 @@ def test_steps_blending( ), "Wrong amount of output time steps in forecast output" # Transform the data back into mm/h - precip_forecast, _ = converter(precip_forecast, metadata) + precip_forecast_dataset = converter(precip_forecast_dataset) + precip_var_forecast = precip_forecast_dataset.attrs["precip_var"] + precip_forecast = precip_forecast_dataset[precip_var_forecast].values assert ( precip_forecast.ndim == 4 diff --git a/pysteps/tests/test_blending_utils.py b/pysteps/tests/test_blending_utils.py index c7a455340..15312c36c 100644 --- a/pysteps/tests/test_blending_utils.py +++ b/pysteps/tests/test_blending_utils.py @@ -12,7 +12,7 @@ blend_cascades, blend_optical_flows, compute_smooth_dilated_mask, - preprocess_and_store_nwp_data, + preprocess_nwp_data, recompose_cascade, ) from pysteps.utils.check_norain import check_norain @@ -139,7 +139,7 @@ def test_blending_utils(precip_nwp_dataset, nwp_model, issue_times, weights): ### # Compute and store the motion ### - preprocess_and_store_nwp_data( + preprocess_nwp_data( precip_nwp_dataset=precip_nwp_dataset, oflow_method=oflow_method, nwp_model=nwp_model, diff --git a/pysteps/tests/test_cascade.py b/pysteps/tests/test_cascade.py index 4a428be99..d9f92b737 100644 --- a/pysteps/tests/test_cascade.py +++ b/pysteps/tests/test_cascade.py @@ -22,18 +22,20 @@ def test_decompose_recompose(): root_path = pysteps.rcparams.data_sources["bom"]["root_path"] rel_path = os.path.join("prcp-cscn", "2", "2018", "06", "16") filename = os.path.join(root_path, rel_path, "2_20180616_120000.prcp-cscn.nc") - precip, _, metadata = pysteps.io.import_bom_rf3(filename) + precip_dataset = pysteps.io.import_bom_rf3(filename) # Convert to rain rate from mm - precip, metadata = pysteps.utils.to_rainrate(precip, metadata) + precip_dataset = pysteps.utils.to_rainrate(precip_dataset) # Log-transform the data - precip, metadata = pysteps.utils.dB_transform( - precip, metadata, threshold=0.1, zerovalue=-15.0 + precip_dataset = pysteps.utils.dB_transform( + precip_dataset, threshold=0.1, zerovalue=-15.0 ) + precip_var = precip_dataset.attrs["precip_var"] + precip = precip_dataset[precip_var].values # Set Nans as the fill value - precip[~np.isfinite(precip)] = metadata["zerovalue"] + precip[~np.isfinite(precip)] = precip_dataset[precip_var].attrs["zerovalue"] # Set number of cascade levels num_cascade_levels = 9 diff --git a/pysteps/xarray_helpers.py b/pysteps/xarray_helpers.py index 403c77247..8e8757112 100644 --- a/pysteps/xarray_helpers.py +++ b/pysteps/xarray_helpers.py @@ -309,6 +309,7 @@ def convert_output_to_xarray_dataset( dataset["time"][-1].values.astype("datetime64[us]").astype(datetime) ) time_metadata = dataset["time"].attrs + time_encoding = dataset["time"].encoding timestep_seconds = dataset["time"].attrs["stepsize"] dataset = dataset.drop_vars([precip_var]).drop_dims(["time"]) if isinstance(timesteps, int): @@ -317,7 +318,7 @@ def convert_output_to_xarray_dataset( last_timestamp + timedelta(seconds=timestep_seconds * i) for i in timesteps ] dataset = dataset.assign_coords( - {"time": (["time"], next_timestamps, time_metadata)} + {"time": (["time"], next_timestamps, time_metadata, time_encoding)} ) if output.ndim == 4: From de057a0dc536fd88ae05c0c3c7bcb4e134f5aacf Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Mon, 15 Sep 2025 16:43:49 +0200 Subject: [PATCH 4/8] fix loose ends and make all relevant tests succeed --- pysteps/decorators.py | 73 ----- pysteps/extrapolation/eulerian_persistence.py | 26 +- pysteps/io/importers.py | 256 +++++++++++------- pysteps/io/interface.py | 3 - pysteps/io/nowcast_importers.py | 3 +- pysteps/nowcasts/interface.py | 6 +- .../tests/test_blending_linear_blending.py | 177 ++---------- pysteps/tests/test_downscaling_rainfarm.py | 2 +- pysteps/tests/test_exporters.py | 2 +- pysteps/tests/test_feature.py | 1 - pysteps/tests/test_interfaces.py | 15 +- pysteps/tests/test_io_bom_rf3.py | 2 - pysteps/tests/test_io_dwd_hdf5.py | 1 - pysteps/tests/test_io_fmi_geotiff.py | 2 - pysteps/tests/test_io_fmi_pgm.py | 1 - pysteps/tests/test_io_knmi_hdf5.py | 3 +- pysteps/tests/test_io_mch_gif.py | 3 +- pysteps/tests/test_io_mrms_grib.py | 8 +- pysteps/tests/test_io_nowcast_importers.py | 1 - pysteps/tests/test_io_readers.py | 1 - pysteps/tests/test_io_saf_crri.py | 2 - pysteps/tests/test_motion_lk.py | 1 - pysteps/tests/test_noise_fftgenerators.py | 1 - pysteps/tests/test_nowcasts_anvil.py | 1 - .../test_nowcasts_lagrangian_probability.py | 1 - pysteps/tests/test_nowcasts_linda.py | 2 - pysteps/tests/test_nowcasts_sprog.py | 1 - pysteps/tests/test_nowcasts_sseps.py | 1 - pysteps/tests/test_nowcasts_steps.py | 2 - pysteps/tests/test_nowcasts_utils.py | 2 - pysteps/tests/test_plt_animate.py | 1 - pysteps/tests/test_plugins_support.py | 14 +- pysteps/tests/test_utils_reprojection.py | 2 - pysteps/tests/test_verification_salscores.py | 12 +- 34 files changed, 220 insertions(+), 409 deletions(-) diff --git a/pysteps/decorators.py b/pysteps/decorators.py index 7c25f7866..8f6e42d80 100644 --- a/pysteps/decorators.py +++ b/pysteps/decorators.py @@ -9,7 +9,6 @@ .. autosummary:: :toctree: ../generated/ - postprocess_import check_input_frames prepare_interpolator memoize @@ -22,8 +21,6 @@ import numpy as np -from pysteps.xarray_helpers import convert_input_to_xarray_dataset - def _add_extra_kwrds_to_docstrings(target_func, extra_kwargs_doc_text): """ @@ -44,76 +41,6 @@ def _add_extra_kwrds_to_docstrings(target_func, extra_kwargs_doc_text): return target_func -def postprocess_import(fillna=np.nan, dtype="double"): - """ - Postprocess the imported precipitation data. - Operations: - - Allow type casting (dtype keyword) - - Set invalid or missing data to predefined value (fillna keyword) - This decorator replaces the text "{extra_kwargs}" in the function's - docstring with the documentation of the keywords used in the postprocessing. - The additional docstrings are added as "Other Parameters" in the importer function. - - Parameters - ---------- - dtype: str - Default data type for precipitation. Double precision by default. - fillna: float or np.nan - Default value used to represent the missing data ("No Coverage"). - By default, np.nan is used. - If the importer returns a MaskedArray, all the masked values are set to the - fillna value. If a numpy array is returned, all the invalid values (nan and inf) - are set to the fillna value. - - """ - - def _postprocess_import(importer): - @wraps(importer) - def _import_with_postprocessing(*args, **kwargs): - precip, quality, metadata = importer(*args, **kwargs) - - _dtype = kwargs.get("dtype", dtype) - - accepted_precisions = ["float32", "float64", "single", "double"] - if _dtype not in accepted_precisions: - raise ValueError( - "The selected precision does not correspond to a valid value." - "The accepted values are: " + str(accepted_precisions) - ) - - if isinstance(precip, np.ma.MaskedArray): - invalid_mask = np.ma.getmaskarray(precip) - precip.data[invalid_mask] = fillna - else: - # If plain numpy arrays are used, the importers should indicate - # the invalid values with np.nan. - _fillna = kwargs.get("fillna", fillna) - if _fillna is not np.nan: - mask = ~np.isfinite(precip) - precip[mask] = _fillna - - return convert_input_to_xarray_dataset( - precip.astype(_dtype), quality, metadata - ) - - extra_kwargs_doc = """ - Other Parameters - ---------------- - dtype: str - Data-type to which the array is cast. - Valid values: "float32", "float64", "single", and "double". - fillna: float or np.nan - Value used to represent the missing data ("No Coverage"). - By default, np.nan is used. - """ - - _add_extra_kwrds_to_docstrings(_import_with_postprocessing, extra_kwargs_doc) - - return _import_with_postprocessing - - return _postprocess_import - - def check_input_frames( minimum_input_frames=2, maximum_input_frames=np.inf, just_ndim=False ): diff --git a/pysteps/extrapolation/eulerian_persistence.py b/pysteps/extrapolation/eulerian_persistence.py index 7ada0e7e1..16e45ad75 100644 --- a/pysteps/extrapolation/eulerian_persistence.py +++ b/pysteps/extrapolation/eulerian_persistence.py @@ -1,47 +1,40 @@ import numpy as np -import xarray as xr -from pysteps.xarray_helpers import convert_output_to_xarray_dataset -def extrapolate(precip_dataset: xr.Dataset, timesteps, outval=np.nan, **kwargs): +def extrapolate(precip, velocity, timesteps, outval=np.nan, **kwargs): """ A dummy extrapolation method to apply Eulerian persistence to a two-dimensional precipitation field. The method returns the a sequence of the same initial field with no extrapolation applied (i.e. Eulerian persistence). - Parameters ---------- - precip_dataset : xarray.Dataset - xarray dataset containing the input precipitation field. All + precip : array-like + Array of shape (m,n) containing the input precipitation field. All values are required to be finite. + velocity : array-like + Not used by the method. timesteps : int or list of floats Number of time steps or a list of time steps. outval : float, optional Not used by the method. - Other Parameters ---------------- return_displacement : bool If True, return the total advection velocity (displacement) between the initial input field and the advected one integrated along the trajectory. Default : False - Returns ------- out : array or tuple If return_displacement=False, return a sequence of the same initial field of shape (num_timesteps,m,n). Otherwise, return a tuple containing the replicated fields and a (2,m,n) array of zeros. - References ---------- :cite:`GZ2002` - """ - del outval # Unused by _eulerian_persistence - precip_var = precip_dataset.attrs["precip_var"] - precip = precip_dataset[precip_var].values[-1] + del velocity, outval # Unused by _eulerian_persistence if isinstance(timesteps, int): num_timesteps = timesteps @@ -51,11 +44,8 @@ def extrapolate(precip_dataset: xr.Dataset, timesteps, outval=np.nan, **kwargs): return_displacement = kwargs.get("return_displacement", False) extrapolated_precip = np.repeat(precip[np.newaxis, :, :], num_timesteps, axis=0) - extrapolated_precip_dataset = convert_output_to_xarray_dataset( - precip_dataset, timesteps, extrapolated_precip - ) if not return_displacement: - return extrapolated_precip_dataset + return extrapolated_precip else: - return extrapolated_precip_dataset, np.zeros((2,) + extrapolated_precip.shape) + return extrapolated_precip, np.zeros((2,) + extrapolated_precip.shape) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index f7c7fb289..633bec3d4 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -208,11 +208,12 @@ import array import datetime from functools import partial +from typing import Any import numpy as np +import numpy.typing as npt from matplotlib.pyplot import imread -from pysteps.decorators import postprocess_import from pysteps.exceptions import DataModelError, MissingOptionalDependency from pysteps.utils import aggregate_fields from pysteps.xarray_helpers import convert_input_to_xarray_dataset @@ -267,6 +268,26 @@ PYGRIB_IMPORTED = False +def _postprocess_precip(precip, fillna=np.nan, dtype="double") -> npt.NDArray[Any]: + accepted_precisions = ["float32", "float64", "single", "double"] + if dtype not in accepted_precisions: + raise ValueError( + "The selected precision does not correspond to a valid value." + "The accepted values are: " + str(accepted_precisions) + ) + + if isinstance(precip, np.ma.MaskedArray): + invalid_mask = np.ma.getmaskarray(precip) + precip.data[invalid_mask] = fillna + else: + # If plain numpy arrays are used, the importers should indicate + # the invalid values with np.nan. + if fillna is not np.nan: + mask = ~np.isfinite(precip) + precip[mask] = fillna + return precip.astype(dtype) + + def _check_coords_range(selected_range, coordinate, full_range): """ Check that the coordinates range arguments follow the expected pattern in @@ -354,7 +375,9 @@ def _get_threshold_value(precip): return np.nan -def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): +def import_mrms_grib( + filename, extent=None, window_size=4, fillna=np.nan, dtype="float32" +): """ Importer for NSSL's Multi-Radar/Multi-Sensor System ([MRMS](https://www.nssl.noaa.gov/projects/mrms/)) rainrate product @@ -407,13 +430,10 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): If an integer value is given, the same block shape is used for all the image dimensions. Default: window_size=4. - - Other Parameters - ---------------- - dtype: str + dtype: str, optional Data-type to which the array is cast. Valid values: "float32", "float64", "single", and "double". - fillna: float or np.nan + fillna: float or np.nan, optional Value used to represent the missing data ("No Coverage"). By default, np.nan is used. @@ -428,33 +448,6 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): metadata: dict Associated metadata (pixel sizes, map projections, etc.). """ - dataset = _import_mrms_grib(filename, extent, window_size, **kwargs) - # Create a function with default arguments for aggregate_fields - block_reduce = partial(aggregate_fields, method="mean", trim=True) - - if window_size != (1, 1): - # Downscale data - precip_var = dataset.attrs["precip_var"] - # block_reduce does not handle nan values - if "fillna" in kwargs: - no_data_mask = dataset[precip_var].values == kwargs["fillna"] - else: - no_data_mask = np.isnan(dataset[precip_var].values) - dataset[precip_var].data[no_data_mask] = 0 - dataset["no_data_mask"] = (("y", "x"), no_data_mask) - dataset = block_reduce(dataset, window_size, dim=("y", "x")) - - # Consider that if a single invalid observation is located in the block, - # then mark that value as invalid. - no_data_mask = dataset.no_data_mask.values == 1.0 - dataset = dataset.drop_vars("no_data_mask") - - return dataset - - -@postprocess_import(dtype="float32") -def _import_mrms_grib(filename, extent=None, window_size=4, **kwargs): - del kwargs if not PYGRIB_IMPORTED: raise MissingOptionalDependency( @@ -550,12 +543,31 @@ def _import_mrms_grib(filename, extent=None, window_size=4, **kwargs): y2=y2 + ysize / 2, cartesian_unit="degrees", ) + precip = _postprocess_precip(precip, fillna, dtype) - return convert_input_to_xarray_dataset(precip, None, metadata) + precip_dataset = convert_input_to_xarray_dataset(precip, None, metadata) + + if window_size != (1, 1): + # Create a function with default arguments for aggregate_fields + block_reduce = partial(aggregate_fields, method="mean", trim=True) + # Downscale data + precip_var = precip_dataset.attrs["precip_var"] + # block_reduce does not handle nan values + no_data_mask = np.isnan(precip_dataset[precip_var].values) + precip_dataset[precip_var].data[no_data_mask] = 0 + precip_dataset["no_data_mask"] = (("y", "x"), no_data_mask) + precip_dataset = block_reduce(precip_dataset, window_size, dim=("y", "x")) + + # Consider that if a single invalid observation is located in the block, + # then mark that value as invalid. + no_data_mask = precip_dataset.no_data_mask.values == 1.0 + precip_dataset = precip_dataset.drop_vars("no_data_mask") + precip_dataset[precip_var].data[no_data_mask] = fillna + return precip_dataset -@postprocess_import() -def import_bom_rf3(filename, **kwargs): + +def import_bom_rf3(filename, gzipped=False, fillna=np.nan, dtype="double"): """ Import a NetCDF radar rainfall product from the BoM Rainfields3. @@ -563,8 +575,12 @@ def import_bom_rf3(filename, **kwargs): ---------- filename: str Name of the file to import. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -585,7 +601,9 @@ def import_bom_rf3(filename, **kwargs): metadata["zerovalue"] = np.nanmin(precip) metadata["threshold"] = _get_threshold_value(precip) - return precip, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _import_bom_rf3_data(filename): @@ -682,8 +700,7 @@ def _import_bom_rf3_geodata(ds_rainfall): return geodata -@postprocess_import() -def import_fmi_geotiff(filename, **kwargs): +def import_fmi_geotiff(filename, fillna=np.nan, dtype="double"): """ Import a reflectivity field (dBZ) from an FMI GeoTIFF file. @@ -691,8 +708,12 @@ def import_fmi_geotiff(filename, **kwargs): ---------- filename: str Name of the file to import. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -747,11 +768,12 @@ def import_fmi_geotiff(filename, **kwargs): metadata["zr_a"] = 223.0 metadata["zr_b"] = 1.53 - return precip, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) -@postprocess_import() -def import_fmi_pgm(filename, gzipped=False, **kwargs): +def import_fmi_pgm(filename, gzipped=False, fillna=np.nan, dtype="double"): """ Import a 8-bit PGM radar reflectivity composite from the FMI archive. @@ -761,8 +783,12 @@ def import_fmi_pgm(filename, gzipped=False, **kwargs): Name of the file to import. gzipped: bool If True, the input file is treated as a compressed gzip file. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -806,7 +832,9 @@ def import_fmi_pgm(filename, gzipped=False, **kwargs): metadata["zr_a"] = 223.0 metadata["zr_b"] = 1.53 - return precip, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _import_fmi_pgm_geodata(metadata): @@ -877,13 +905,8 @@ def _import_fmi_pgm_metadata(filename, gzipped=False): return metadata -@postprocess_import() def import_knmi_hdf5( - filename, - qty="ACRR", - accutime=5.0, - pixelsize=1000.0, - **kwargs, + filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, fillna=np.nan, dtype="double" ): """ Import a precipitation or reflectivity field (and optionally the quality @@ -905,8 +928,12 @@ def import_knmi_hdf5( The pixel size of a raster cell in meters. The default value for the KNMI datasets is a 1000 m grid cell size, but datasets with 2400 m pixel size are also available. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1030,11 +1057,12 @@ def import_knmi_hdf5( f.close() - return precip, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) -@postprocess_import() -def import_mch_gif(filename, product, unit, accutime, **kwargs): +def import_mch_gif(filename, product, unit, accutime, fillna=np.nan, dtype="double"): """ Import a 8-bit gif radar reflectivity composite from the MeteoSwiss archive. @@ -1063,8 +1091,12 @@ def import_mch_gif(filename, product, unit, accutime, **kwargs): the physical unit of the data accutime: float the accumulation time in minutes of the data - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1160,11 +1192,12 @@ def import_mch_gif(filename, product, unit, accutime, **kwargs): metadata["zr_a"] = 316.0 metadata["zr_b"] = 1.5 - return precip, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) -@postprocess_import() -def import_mch_hdf5(filename, qty="RATE", **kwargs): +def import_mch_hdf5(filename, qty="RATE", fillna=np.nan, dtype="double"): """ Import a precipitation field (and optionally the quality field) from a MeteoSwiss HDF5 file conforming to the ODIM specification. @@ -1178,8 +1211,12 @@ def import_mch_hdf5(filename, qty="RATE", **kwargs): are: 'RATE'=instantaneous rain rate (mm/h), 'ACRR'=hourly rainfall accumulation (mm) and 'DBZH'=max-reflectivity (dBZ). The default value is 'RATE'. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1295,7 +1332,9 @@ def import_mch_hdf5(filename, qty="RATE", **kwargs): f.close() - return precip, quality, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _read_mch_hdf5_what_group(whatgrp): @@ -1308,7 +1347,6 @@ def _read_mch_hdf5_what_group(whatgrp): return qty, gain, offset, nodata, undetect -@postprocess_import() def import_mch_metranet(filename, product, unit, accutime): """ Import a 8-bit bin radar reflectivity composite from the MeteoSwiss @@ -1338,8 +1376,12 @@ def import_mch_metranet(filename, product, unit, accutime): the physical unit of the data accutime: float the accumulation time in minutes of the data - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1370,7 +1412,9 @@ def import_mch_metranet(filename, product, unit, accutime): metadata["zr_a"] = 316.0 metadata["zr_b"] = 1.5 - return precip, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _import_mch_geodata(): @@ -1408,8 +1452,7 @@ def _import_mch_geodata(): return geodata -@postprocess_import() -def import_odim_hdf5(filename, qty="RATE", **kwargs): +def import_odim_hdf5(filename, qty="RATE", fillna=np.nan, dtype="double", **kwargs): """ Import a precipitation field (and optionally the quality field) from a HDF5 file conforming to the ODIM specification. @@ -1426,8 +1469,12 @@ def import_odim_hdf5(filename, qty="RATE", **kwargs): are: 'RATE'=instantaneous rain rate (mm/h), 'ACRR'=hourly rainfall accumulation (mm) and 'DBZH'=max-reflectivity (dBZ). The default value is 'RATE'. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1629,17 +1676,19 @@ def import_odim_hdf5(filename, qty="RATE", **kwargs): f.close() - return precip, quality, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) -def import_opera_hdf5(filename, qty="RATE", **kwargs): +def import_opera_hdf5(filename, qty="RATE", fillna=np.nan, dtype="double"): """ Wrapper to :py:func:`pysteps.io.importers.import_odim_hdf5` to maintain backward compatibility with previous pysteps versions. **Important:** Use :py:func:`~pysteps.io.importers.import_odim_hdf5` instead. """ - return import_odim_hdf5(filename, qty=qty, **kwargs) + return import_odim_hdf5(filename, qty, fillna, dtype) def _read_opera_hdf5_what_group(whatgrp): @@ -1652,8 +1701,9 @@ def _read_opera_hdf5_what_group(whatgrp): return qty, gain, offset, nodata, undetect -@postprocess_import() -def import_saf_crri(filename, extent=None, **kwargs): +def import_saf_crri( + filename, extent=None, gzipped=False, fillna=np.nan, dtype="double" +): """ Import a NetCDF radar rainfall product from the Convective Rainfall Rate Intensity (CRRI) product from the Satellite Application Facilities (SAF). @@ -1668,8 +1718,12 @@ def import_saf_crri(filename, extent=None, **kwargs): extent: scalars (left, right, bottom, top), optional The spatial extent specified in data coordinates. If None, the full extent is imported. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1725,7 +1779,9 @@ def import_saf_crri(filename, extent=None, **kwargs): metadata["zerovalue"] = np.nanmin(precip) metadata["threshold"] = _get_threshold_value(precip) - return precip, quality, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _import_saf_crri_data(filename, idx_x=None, idx_y=None): @@ -1786,8 +1842,7 @@ def _import_saf_crri_geodata(filename): return geodata -@postprocess_import() -def import_dwd_hdf5(filename, qty="RATE", **kwargs): +def import_dwd_hdf5(filename, qty="RATE", fillna=np.nan, dtype="double", **kwargs): """ Import a DWD precipitation product field (and optionally the quality field) from a HDF5 file conforming to the ODIM specification @@ -1801,8 +1856,12 @@ def import_dwd_hdf5(filename, qty="RATE", **kwargs): are: 'RATE'=instantaneous rain rate (mm/h), 'ACRR'=hourly rainfall accumulation (mm) and 'DBZH'=max-reflectivity (dBZ). The default value is 'RATE'. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -1984,7 +2043,9 @@ def import_dwd_hdf5(filename, qty="RATE", **kwargs): f.close() - return precip, quality, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _read_hdf5_cont(f, d): @@ -2062,8 +2123,7 @@ def _get_whatgrp(d, g): return -@postprocess_import() -def import_dwd_radolan(filename, product_name): +def import_dwd_radolan(filename, product_name, fillna=np.nan, dtype="double"): """ Import a RADOLAN precipitation product from a binary file. @@ -2076,8 +2136,12 @@ def import_dwd_radolan(filename, product_name): https://www.dwd.de/DE/leistungen/radolan/radolan_info/ radolan_radvor_op_komposit_format_pdf.pdf for a detailed description. - - {extra_kwargs_doc} + dtype: str, optional + Data-type to which the array is cast. + Valid values: "float32", "float64", "single", and "double". + fillna: float or np.nan, optional + Value used to represent the missing data ("No Coverage"). + By default, np.nan is used. Returns ------- @@ -2160,7 +2224,9 @@ def import_dwd_radolan(filename, product_name): geodata = _import_dwd_geodata(product_name, dims) metadata = geodata - return data, None, metadata + precip = _postprocess_precip(precip, fillna, dtype) + + return convert_input_to_xarray_dataset(precip, None, metadata) def _identify_info_bits(data): diff --git a/pysteps/io/interface.py b/pysteps/io/interface.py index 58b684ea3..828d791fa 100644 --- a/pysteps/io/interface.py +++ b/pysteps/io/interface.py @@ -14,7 +14,6 @@ """ from importlib.metadata import entry_points -from pysteps.decorators import postprocess_import from pysteps.io import importers, exporters, interface from pprint import pprint @@ -57,8 +56,6 @@ def discover_importers(): importer_function_name = _importer.__name__ importer_short_name = importer_function_name.replace("import_", "") - _postprocess_kws = getattr(_importer, "postprocess_kws", dict()) - _importer = postprocess_import(**_postprocess_kws)(_importer) if importer_short_name not in _importer_methods: _importer_methods[importer_short_name] = _importer else: diff --git a/pysteps/io/nowcast_importers.py b/pysteps/io/nowcast_importers.py index f353ec75b..cc14c7cd9 100755 --- a/pysteps/io/nowcast_importers.py +++ b/pysteps/io/nowcast_importers.py @@ -70,8 +70,7 @@ import numpy as np -from pysteps.decorators import postprocess_import -from pysteps.exceptions import MissingOptionalDependency, DataModelError +from pysteps.exceptions import MissingOptionalDependency import xarray as xr try: diff --git a/pysteps/nowcasts/interface.py b/pysteps/nowcasts/interface.py index 2af2048c3..7777f8310 100644 --- a/pysteps/nowcasts/interface.py +++ b/pysteps/nowcasts/interface.py @@ -30,7 +30,7 @@ get_method """ -from pysteps.extrapolation.interface import eulerian_persistence +from functools import partial from pysteps.nowcasts import ( anvil, extrapolation, @@ -41,9 +41,11 @@ ) from pysteps.nowcasts import lagrangian_probability +eulerian_persistence = partial(extrapolation.forecast, extrap_method="eulerian") + _nowcast_methods = dict() _nowcast_methods["anvil"] = anvil.forecast -_nowcast_methods["eulerian"] = eulerian_persistence.extrapolate +_nowcast_methods["eulerian"] = eulerian_persistence _nowcast_methods["extrapolation"] = extrapolation.forecast _nowcast_methods["lagrangian"] = extrapolation.forecast _nowcast_methods["lagrangian_probability"] = lagrangian_probability.forecast diff --git a/pysteps/tests/test_blending_linear_blending.py b/pysteps/tests/test_blending_linear_blending.py index 34d3f3875..2dc636cfd 100644 --- a/pysteps/tests/test_blending_linear_blending.py +++ b/pysteps/tests/test_blending_linear_blending.py @@ -10,155 +10,35 @@ # Test function arguments linear_arg_values = [ - (5, 30, 60, 20, 45, "eulerian", None, 1, False, True, False), - (5, 30, 60, 20, 45, "eulerian", None, 2, False, False, False), - (5, 30, 60, 20, 45, "eulerian", None, 0, False, False, False), - (4, 23, 33, 9, 28, "eulerian", None, 1, False, False, False), - (3, 18, 36, 13, 27, "eulerian", None, 1, False, False, False), - (7, 30, 68, 11, 49, "eulerian", None, 1, False, False, False), - (7, 30, 68, 11, 49, "eulerian", None, 1, False, False, True), - (10, 100, 160, 25, 130, "eulerian", None, 1, False, False, False), - (6, 60, 180, 22, 120, "eulerian", None, 1, False, False, False), - (5, 100, 200, 40, 150, "eulerian", None, 1, False, False, False), - ( - 5, - 30, - 60, - 20, - 45, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 4, - 23, - 33, - 9, - 28, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 3, - 18, - 36, - 13, - 27, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 7, - 30, - 68, - 11, - 49, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 10, - 100, - 160, - 25, - 130, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 6, - 60, - 180, - 22, - 120, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 5, - 100, - 200, - 40, - 150, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - False, - ), - ( - 5, - 100, - 200, - 40, - 150, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - False, - False, - True, - ), - (5, 30, 60, 20, 45, "eulerian", None, 1, True, True, False), - (5, 30, 60, 20, 45, "eulerian", None, 2, True, False, False), - (5, 30, 60, 20, 45, "eulerian", None, 0, True, False, False), - ( - 5, - 30, - 60, - 20, - 45, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - True, - False, - False, - ), - (4, 23, 33, 9, 28, "extrapolation", np.zeros((2, 200, 200)), 1, True, False, False), - ( - 3, - 18, - 36, - 13, - 27, - "extrapolation", - np.zeros((2, 200, 200)), - 1, - True, - False, - False, - ), + (5, 30, 60, 20, 45, "eulerian", 1, False, True, False), + (5, 30, 60, 20, 45, "eulerian", 2, False, False, False), + (5, 30, 60, 20, 45, "eulerian", 0, False, False, False), + (4, 23, 33, 9, 28, "eulerian", 1, False, False, False), + (3, 18, 36, 13, 27, "eulerian", 1, False, False, False), + (7, 30, 68, 11, 49, "eulerian", 1, False, False, False), + (7, 30, 68, 11, 49, "eulerian", 1, False, False, True), + (10, 100, 160, 25, 130, "eulerian", 1, False, False, False), + (6, 60, 180, 22, 120, "eulerian", 1, False, False, False), + (5, 100, 200, 40, 150, "eulerian", 1, False, False, False), + (5, 30, 60, 20, 45, "extrapolation", 1, False, False, False), + (4, 23, 33, 9, 28, "extrapolation", 1, False, False, False), + (3, 18, 36, 13, 27, "extrapolation", 1, False, False, False), + (7, 30, 68, 11, 49, "extrapolation", 1, False, False, False), + (10, 100, 160, 25, 130, "extrapolation", 1, False, False, False), + (6, 60, 180, 22, 120, "extrapolation", 1, False, False, False), + (5, 100, 200, 40, 150, "extrapolation", 1, False, False, False), + (5, 100, 200, 40, 150, "extrapolation", 1, False, False, True), + (5, 30, 60, 20, 45, "eulerian", 1, True, True, False), + (5, 30, 60, 20, 45, "eulerian", 2, True, False, False), + (5, 30, 60, 20, 45, "eulerian", 0, True, False, False), + (5, 30, 60, 20, 45, "extrapolation", 1, True, False, False), + (4, 23, 33, 9, 28, "extrapolation", 1, True, False, False), + (3, 18, 36, 13, 27, "extrapolation", 1, True, False, False), ] @pytest.mark.parametrize( - "timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V, n_models, salient_blending, squeeze_nwp_array, fill_nwp", + "timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, n_models, salient_blending, squeeze_nwp_array, fill_nwp", linear_arg_values, ) def test_linear_blending( @@ -168,7 +48,6 @@ def test_linear_blending( n_timesteps, controltime, nowcast_method, - V, n_models, salient_blending, squeeze_nwp_array, @@ -253,9 +132,9 @@ def test_linear_blending( radar_dataset = transformation.dB_transform( radar_dataset, threshold=0.1, zerovalue=-15.0 ) - if V is not None: - radar_dataset["velocity_x"] = (["y", "x"], V[0]) - radar_dataset["velocity_y"] = (["y", "x"], V[1]) + V = np.zeros((2, 200, 200)) + radar_dataset["velocity_x"] = (["y", "x"], V[0]) + radar_dataset["velocity_y"] = (["y", "x"], V[1]) if r_nwp is None: model_dataset = None diff --git a/pysteps/tests/test_downscaling_rainfarm.py b/pysteps/tests/test_downscaling_rainfarm.py index 3b60f48b3..aeb6deba5 100644 --- a/pysteps/tests/test_downscaling_rainfarm.py +++ b/pysteps/tests/test_downscaling_rainfarm.py @@ -10,7 +10,7 @@ @pytest.fixture(scope="module") def dataset(): precip_dataset = get_precipitation_fields( - num_prev_files=0, num_next_files=0, return_raw=False, metadata=True + num_prev_files=0, num_next_files=0, return_raw=False ) precip_dataset = square_domain(precip_dataset, "crop") return precip_dataset diff --git a/pysteps/tests/test_exporters.py b/pysteps/tests/test_exporters.py index 274390724..6208f3f7d 100644 --- a/pysteps/tests/test_exporters.py +++ b/pysteps/tests/test_exporters.py @@ -71,7 +71,7 @@ def test_io_export_netcdf_one_member_one_time_step( pytest.importorskip("pyproj") precip_dataset: xr.Dataset = get_precipitation_fields( - num_prev_files=2, return_raw=True, metadata=True, source="fmi" + num_prev_files=2, return_raw=True, source="fmi" ) precip_var = precip_dataset.attrs["precip_var"] diff --git a/pysteps/tests/test_feature.py b/pysteps/tests/test_feature.py index d1a257aff..084744b29 100644 --- a/pysteps/tests/test_feature.py +++ b/pysteps/tests/test_feature.py @@ -18,7 +18,6 @@ def test_feature(method, max_num_features): num_prev_files=0, num_next_files=0, return_raw=True, - metadata=False, upscale=None, source="mch", ) diff --git a/pysteps/tests/test_interfaces.py b/pysteps/tests/test_interfaces.py index bfe1d6683..1f1c2ef6d 100644 --- a/pysteps/tests/test_interfaces.py +++ b/pysteps/tests/test_interfaces.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- +from functools import partial import numpy import pytest @@ -61,14 +62,14 @@ def test_extrapolation_interface(): from pysteps import extrapolation from pysteps.extrapolation import semilagrangian - from pysteps.extrapolation.interface import eulerian_persistence as eulerian + from pysteps.extrapolation import eulerian_persistence from pysteps.extrapolation.interface import _do_nothing as do_nothing method_getter = extrapolation.interface.get_method valid_returned_objs = dict() valid_returned_objs["semilagrangian"] = semilagrangian.extrapolate - valid_returned_objs["eulerian"] = eulerian + valid_returned_objs["eulerian"] = eulerian_persistence.extrapolate valid_returned_objs[None] = do_nothing valid_returned_objs["None"] = do_nothing @@ -295,6 +296,7 @@ def test_nowcasts_interface(): valid_names_func_pair = [ ("anvil", anvil.forecast), ("extrapolation", extrapolation.forecast), + ("eulerian", pysteps.nowcasts.interface.eulerian_persistence), ("lagrangian", extrapolation.forecast), ("linda", linda.forecast), ("probability", lagrangian_probability.forecast), @@ -307,15 +309,6 @@ def test_nowcasts_interface(): invalid_names = ["extrap", "step", "s-prog", "pysteps"] _generic_interface_test(method_getter, valid_names_func_pair, invalid_names) - # Test eulerian persistence method - precip = numpy.random.rand(100, 100) - velocity = numpy.random.rand(100, 100) - num_timesteps = 10 - for name in ["eulerian", "EULERIAN"]: - forecast = method_getter(name)(precip, velocity, num_timesteps) - for i in range(num_timesteps): - assert numpy.all(forecast[i] == precip) - def test_utils_interface(): """Test utils module interface.""" diff --git a/pysteps/tests/test_io_bom_rf3.py b/pysteps/tests/test_io_bom_rf3.py index 66f075e0e..22535c4dd 100644 --- a/pysteps/tests/test_io_bom_rf3.py +++ b/pysteps/tests/test_io_bom_rf3.py @@ -11,7 +11,6 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="bom", log_transform=False, ) @@ -46,7 +45,6 @@ def test_io_import_bom_shape(): (precip_dataset.x.attrs["units"], "m", None), (precip_dataset.y.attrs["units"], "m", None), (precip_dataarray.attrs["accutime"], 6, 1e-4), - (precip_dataarray.attrs["transform"], None, None), (precip_dataarray.attrs["zerovalue"], 0.0, 1e-4), (precip_dataarray.attrs["units"], "mm", None), ] diff --git a/pysteps/tests/test_io_dwd_hdf5.py b/pysteps/tests/test_io_dwd_hdf5.py index 950da8d86..e69b2a6bb 100644 --- a/pysteps/tests/test_io_dwd_hdf5.py +++ b/pysteps/tests/test_io_dwd_hdf5.py @@ -9,7 +9,6 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="dwd", log_transform=False, ) diff --git a/pysteps/tests/test_io_fmi_geotiff.py b/pysteps/tests/test_io_fmi_geotiff.py index 2a07f03b0..db2d0a3c2 100644 --- a/pysteps/tests/test_io_fmi_geotiff.py +++ b/pysteps/tests/test_io_fmi_geotiff.py @@ -7,8 +7,6 @@ precip_dataset = get_precipitation_fields( num_prev_files=0, num_next_files=0, - return_raw=True, - metadata=True, source="fmi_geotiff", log_transform=False, ) diff --git a/pysteps/tests/test_io_fmi_pgm.py b/pysteps/tests/test_io_fmi_pgm.py index a704e0e50..48075bd21 100644 --- a/pysteps/tests/test_io_fmi_pgm.py +++ b/pysteps/tests/test_io_fmi_pgm.py @@ -7,7 +7,6 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="fmi", log_transform=False, ) diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py index e585055f9..ef65d8b02 100644 --- a/pysteps/tests/test_io_knmi_hdf5.py +++ b/pysteps/tests/test_io_knmi_hdf5.py @@ -8,10 +8,9 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="knmi", log_transform=False, - importer_kwargs=dict(qty="ACRR"), + qty="ACRR", ) precip_var = precip_dataset.attrs["precip_var"] diff --git a/pysteps/tests/test_io_mch_gif.py b/pysteps/tests/test_io_mch_gif.py index aa08bcb53..254eb526f 100644 --- a/pysteps/tests/test_io_mch_gif.py +++ b/pysteps/tests/test_io_mch_gif.py @@ -8,10 +8,9 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="mch", log_transform=False, - importer_kwargs=dict(qty="AQC"), + product="AQC", ) precip_var = precip_dataset.attrs["precip_var"] diff --git a/pysteps/tests/test_io_mrms_grib.py b/pysteps/tests/test_io_mrms_grib.py index 3b28fc5af..d70b8ffbd 100644 --- a/pysteps/tests/test_io_mrms_grib.py +++ b/pysteps/tests/test_io_mrms_grib.py @@ -8,13 +8,11 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="mrms", log_transform=False, window_size=1, ) -print(precip_dataset) precip_var = precip_dataset.attrs["precip_var"] precip_dataarray = precip_dataset[precip_var] @@ -36,7 +34,6 @@ def test_io_import_mrms_grib(): None, ), (precip_dataarray.attrs["units"], "mm/h", None), - (precip_dataarray.attrs["transform"], None, None), (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6), (precip_dataarray.attrs["threshold"], 0.1, 1e-10), (precip_dataset.x.isel(x=0).values, -129.995, 1e-10), @@ -63,7 +60,6 @@ def test_io_import_mrms_grib_dataset_extent(): num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="mrms", log_transform=False, extent=(230, 300, 20, 55), @@ -79,7 +75,6 @@ def test_io_import_mrms_grib_dataset_extent(): num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="mrms", log_transform=False, extent=(250, 260, 30, 35), @@ -91,6 +86,7 @@ def test_io_import_mrms_grib_dataset_extent(): smart_assert(precip_dataarray_even_smaller.shape, (1, 500, 1000), None) # XR: we had to change the selection of the original field since these is a flip happening in the way the data is read in. # XR: We had two ways to solve this: precip_dataarray[:,::-1, :][:, 2000:2500, 2000:3000][:,::-1, :] or switch the 2000:2500 to + # I think this is logical as both the extend selected data and the reference data are flipped. That is why we need the double flip assert_array_almost_equal( precip_dataarray.values[:, 1000:1500, 2000:3000], precip_dataarray_even_smaller.values, @@ -100,7 +96,6 @@ def test_io_import_mrms_grib_dataset_extent(): num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="mrms", log_transform=False, extent=(250, 260, 30, 35), @@ -116,7 +111,6 @@ def test_io_import_mrms_grib_dataset_extent(): num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="mrms", log_transform=False, extent=(250, 260, 30, 35), diff --git a/pysteps/tests/test_io_nowcast_importers.py b/pysteps/tests/test_io_nowcast_importers.py index a19388afe..e6f98f1ad 100644 --- a/pysteps/tests/test_io_nowcast_importers.py +++ b/pysteps/tests/test_io_nowcast_importers.py @@ -10,7 +10,6 @@ num_prev_files=1, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) diff --git a/pysteps/tests/test_io_readers.py b/pysteps/tests/test_io_readers.py index 692e2d9eb..9417e931e 100644 --- a/pysteps/tests/test_io_readers.py +++ b/pysteps/tests/test_io_readers.py @@ -11,7 +11,6 @@ def test_read_timeseries_mch(): num_prev_files=1, num_next_files=1, return_raw=True, - metadata=True, source="mch", log_transform=False, ) diff --git a/pysteps/tests/test_io_saf_crri.py b/pysteps/tests/test_io_saf_crri.py index 9ecedf135..b75273e4a 100644 --- a/pysteps/tests/test_io_saf_crri.py +++ b/pysteps/tests/test_io_saf_crri.py @@ -8,7 +8,6 @@ num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="saf", log_transform=False, ) @@ -37,7 +36,6 @@ def test_io_import_saf_crri_extent(extent, expected_extent, expected_shape, tole num_prev_files=0, num_next_files=0, return_raw=True, - metadata=True, source="saf", log_transform=False, extent=extent, diff --git a/pysteps/tests/test_motion_lk.py b/pysteps/tests/test_motion_lk.py index 3ce9e2059..7cffb6768 100644 --- a/pysteps/tests/test_motion_lk.py +++ b/pysteps/tests/test_motion_lk.py @@ -64,7 +64,6 @@ def test_lk( num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) precip_var = dataset.attrs["precip_var"] diff --git a/pysteps/tests/test_noise_fftgenerators.py b/pysteps/tests/test_noise_fftgenerators.py index cc27258e1..fafdedaa3 100644 --- a/pysteps/tests/test_noise_fftgenerators.py +++ b/pysteps/tests/test_noise_fftgenerators.py @@ -8,7 +8,6 @@ num_prev_files=0, num_next_files=0, return_raw=False, - metadata=False, upscale=2000, ) diff --git a/pysteps/tests/test_nowcasts_anvil.py b/pysteps/tests/test_nowcasts_anvil.py index f9259fc93..7a082ad65 100644 --- a/pysteps/tests/test_nowcasts_anvil.py +++ b/pysteps/tests/test_nowcasts_anvil.py @@ -92,7 +92,6 @@ def test_anvil_rainrate( num_prev_files=4, num_next_files=0, return_raw=False, - metadata=False, upscale=2000, ) diff --git a/pysteps/tests/test_nowcasts_lagrangian_probability.py b/pysteps/tests/test_nowcasts_lagrangian_probability.py index d75b29e87..a92204ac1 100644 --- a/pysteps/tests/test_nowcasts_lagrangian_probability.py +++ b/pysteps/tests/test_nowcasts_lagrangian_probability.py @@ -86,7 +86,6 @@ def test_real_case(): num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) precip_var = dataset_input.attrs["precip_var"] diff --git a/pysteps/tests/test_nowcasts_linda.py b/pysteps/tests/test_nowcasts_linda.py index 51d688644..8bc6f8913 100644 --- a/pysteps/tests/test_nowcasts_linda.py +++ b/pysteps/tests/test_nowcasts_linda.py @@ -76,7 +76,6 @@ def test_linda( dataset_input = get_precipitation_fields( num_prev_files=2, num_next_files=0, - metadata=True, clip=(354000, 866000, -96000, 416000), upscale=4000, log_transform=False, @@ -203,7 +202,6 @@ def test_linda_callback(tmp_path): num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) precip_input = precip_input.filled() diff --git a/pysteps/tests/test_nowcasts_sprog.py b/pysteps/tests/test_nowcasts_sprog.py index 383fc155e..63e9b795d 100644 --- a/pysteps/tests/test_nowcasts_sprog.py +++ b/pysteps/tests/test_nowcasts_sprog.py @@ -61,7 +61,6 @@ def test_sprog( num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) diff --git a/pysteps/tests/test_nowcasts_sseps.py b/pysteps/tests/test_nowcasts_sseps.py index ee7a6b885..bfe05a002 100644 --- a/pysteps/tests/test_nowcasts_sseps.py +++ b/pysteps/tests/test_nowcasts_sseps.py @@ -108,7 +108,6 @@ def test_sseps( num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) precip_var = dataset_input.attrs["precip_var"] diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py index a10760192..4e489dcfb 100644 --- a/pysteps/tests/test_nowcasts_steps.py +++ b/pysteps/tests/test_nowcasts_steps.py @@ -77,7 +77,6 @@ def test_steps_skill( num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) @@ -133,7 +132,6 @@ def test_steps_callback(tmp_path): num_prev_files=2, num_next_files=0, return_raw=False, - metadata=True, upscale=2000, ) precip_input = precip_input.filled() diff --git a/pysteps/tests/test_nowcasts_utils.py b/pysteps/tests/test_nowcasts_utils.py index 1dfeb27a9..a1a8d8133 100644 --- a/pysteps/tests/test_nowcasts_utils.py +++ b/pysteps/tests/test_nowcasts_utils.py @@ -29,8 +29,6 @@ def test_nowcast_main_loop( dataset = get_precipitation_fields( num_prev_files=2, num_next_files=0, - return_raw=False, - metadata=False, upscale=2000, ) diff --git a/pysteps/tests/test_plt_animate.py b/pysteps/tests/test_plt_animate.py index f6b0d25a1..422ac0057 100644 --- a/pysteps/tests/test_plt_animate.py +++ b/pysteps/tests/test_plt_animate.py @@ -14,7 +14,6 @@ num_prev_files=2, num_next_files=0, return_raw=True, - metadata=True, upscale=2000, ) diff --git a/pysteps/tests/test_plugins_support.py b/pysteps/tests/test_plugins_support.py index f4170492c..46bfd0e99 100644 --- a/pysteps/tests/test_plugins_support.py +++ b/pysteps/tests/test_plugins_support.py @@ -20,13 +20,6 @@ from pysteps import io, postprocessing -# BUG: -# XR: Cookie cutter makes two calls to importers, one -# which is not wrapped with postprocess_import resulting in the importer -# returning 3 values on the first call. On the second call the importer goes -# through the postprocess_import resulting in one dataset being returned. -# Should fix this issue first before fixing tests. -# Test has been therefore commented out. def _check_installed_importer_plugin(import_func_name): # reload the pysteps module to detect the installed plugin io.discover_importers() @@ -90,7 +83,12 @@ def _uninstall_plugin(project_name): ) -# XR: Commented out test for reason explained above +# BUG: +# XR: Cookie cutter tries to add an importer using the postprocess importer +# decorator. I removed this decorator and made it so that any importer +# just directly returns an xarray. This example plugin needs to be updated +# to reflect that before this test will work again. + # def test_importers_plugins(): # with _create_and_install_plugin("pysteps-importer-institution-fun", "importer"): # _check_installed_importer_plugin("importer_institution_fun") diff --git a/pysteps/tests/test_utils_reprojection.py b/pysteps/tests/test_utils_reprojection.py index 58a1231cf..c5f1fa032 100644 --- a/pysteps/tests/test_utils_reprojection.py +++ b/pysteps/tests/test_utils_reprojection.py @@ -1,6 +1,5 @@ # -*- coding: utf-8 -*- -import os import numpy as np import pytest import xarray as xr @@ -92,7 +91,6 @@ def build_precip_dataset( num_next_files=0, source="rmi", return_raw=True, - metadata=True, log_transform=False, ) diff --git a/pysteps/tests/test_verification_salscores.py b/pysteps/tests/test_verification_salscores.py index fbed23def..aa225c2bd 100644 --- a/pysteps/tests/test_verification_salscores.py +++ b/pysteps/tests/test_verification_salscores.py @@ -20,9 +20,7 @@ class TestSAL: def test_sal_zeros(self, converter, thr_factor): """Test the SAL verification method.""" - dataset_input = get_precipitation_fields( - num_prev_files=0, log_transform=False, metadata=True - ) + dataset_input = get_precipitation_fields(num_prev_files=0, log_transform=False) dataset_input = converter(dataset_input) precip_var = dataset_input.attrs["precip_var"] precip = dataset_input[precip_var].values[0] @@ -37,9 +35,7 @@ def test_sal_zeros(self, converter, thr_factor): def test_sal_same_image(self, converter, thr_factor): """Test the SAL verification method.""" - dataset_input = get_precipitation_fields( - num_prev_files=0, log_transform=False, metadata=True - ) + dataset_input = get_precipitation_fields(num_prev_files=0, log_transform=False) dataset_input = converter(dataset_input) precip_var = dataset_input.attrs["precip_var"] precip = dataset_input[precip_var].values[0] @@ -49,9 +45,7 @@ def test_sal_same_image(self, converter, thr_factor): assert np.allclose(result, [0, 0, 0]) def test_sal_translation(self, converter, thr_factor): - dataset_input = get_precipitation_fields( - num_prev_files=0, log_transform=False, metadata=True - ) + dataset_input = get_precipitation_fields(num_prev_files=0, log_transform=False) dataset_input = converter(dataset_input) precip_var = dataset_input.attrs["precip_var"] precip = dataset_input[precip_var].values[0] From e5dff56e059349f72c5b26815ee193596a93d35d Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Tue, 16 Sep 2025 09:06:39 +0200 Subject: [PATCH 5/8] fix donothing --- pysteps/utils/interface.py | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/pysteps/utils/interface.py b/pysteps/utils/interface.py index ded0c3094..eb91dacf3 100644 --- a/pysteps/utils/interface.py +++ b/pysteps/utils/interface.py @@ -11,17 +11,21 @@ get_method """ -from . import arrays -from . import cleansing -from . import conversion -from . import dimension -from . import fft -from . import images -from . import interpolate -from . import reprojection -from . import spectral -from . import tapering -from . import transformation +import xarray as xr + +from . import ( + arrays, + cleansing, + conversion, + dimension, + fft, + images, + interpolate, + reprojection, + spectral, + tapering, + transformation, +) def get_method(name, **kwargs): @@ -163,8 +167,8 @@ def get_method(name, **kwargs): name = name.lower() - def donothing(R, metadata=None, *args, **kwargs): - return R.copy(), {} if metadata is None else metadata.copy() + def donothing(dataset: xr.Dataset, *args, **kwargs): + return dataset methods_objects = dict() methods_objects["none"] = donothing From 111635f08615f0488c4f8ae8c5704800494b0724 Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Tue, 16 Sep 2025 16:20:21 +0200 Subject: [PATCH 6/8] Add forecast reference time and drop velocity in output dataset --- pysteps/xarray_helpers.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/pysteps/xarray_helpers.py b/pysteps/xarray_helpers.py index 821fd9298..31ee9a2b6 100644 --- a/pysteps/xarray_helpers.py +++ b/pysteps/xarray_helpers.py @@ -19,7 +19,6 @@ import pyproj import xarray as xr - # TODO(converters): Write methods for converting Proj.4 projection definitions # into CF grid mapping attributes. Currently this has been implemented for # the stereographic projection. @@ -330,17 +329,22 @@ def convert_output_to_xarray_dataset( precip_var = dataset.attrs["precip_var"] metadata = dataset[precip_var].attrs - last_timestamp = ( + forecast_reference_time = ( dataset["time"][-1].values.astype("datetime64[us]").astype(datetime) ) time_metadata = dataset["time"].attrs time_encoding = dataset["time"].encoding timestep_seconds = dataset["time"].attrs["stepsize"] dataset = dataset.drop_vars([precip_var]).drop_dims(["time"]) + if "velocity_x" in dataset: + dataset = dataset.drop_vars(["velocity_x"]) + if "velocity_y" in dataset: + dataset = dataset.drop_vars(["velocity_y"]) if isinstance(timesteps, int): timesteps = list(range(1, timesteps + 1)) next_timestamps = [ - last_timestamp + timedelta(seconds=timestep_seconds * i) for i in timesteps + forecast_reference_time + timedelta(seconds=timestep_seconds * i) + for i in timesteps ] dataset = dataset.assign_coords( {"time": (["time"], next_timestamps, time_metadata, time_encoding)} @@ -364,4 +368,15 @@ def convert_output_to_xarray_dataset( else: dataset[precip_var] = (["time", "y", "x"], output, metadata) + dataset = dataset.assign_coords( + { + "forecast_reference_time": ( + [], + forecast_reference_time, + {"long_name": "forecast reference time"}, + time_encoding, + ) + } + ) + return dataset From df21567da322bb694fab641b30aad9fe779bfc2c Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Tue, 16 Sep 2025 16:54:00 +0200 Subject: [PATCH 7/8] fix grid_mapping_name --- pysteps/xarray_helpers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysteps/xarray_helpers.py b/pysteps/xarray_helpers.py index 31ee9a2b6..a45ed27ba 100644 --- a/pysteps/xarray_helpers.py +++ b/pysteps/xarray_helpers.py @@ -214,7 +214,7 @@ def convert_input_to_xarray_dataset( "units": attrs["units"], "standard_name": attrs["standard_name"], "long_name": attrs["long_name"], - "grid_mapping": "projection", + "grid_mapping": grid_mapping_name, }, ) } @@ -238,7 +238,7 @@ def convert_input_to_xarray_dataset( { "units": "1", "standard_name": "quality_flag", - "grid_mapping": "projection", + "grid_mapping": grid_mapping_name, }, ) coords = { From d8d8dcd77defaa3d2b5961a2b8b97a62a25022c9 Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Tue, 16 Sep 2025 17:08:24 +0200 Subject: [PATCH 8/8] extrapolate from last value not first --- pysteps/nowcasts/extrapolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/nowcasts/extrapolation.py b/pysteps/nowcasts/extrapolation.py index a70b6985c..6e45092d2 100644 --- a/pysteps/nowcasts/extrapolation.py +++ b/pysteps/nowcasts/extrapolation.py @@ -71,7 +71,7 @@ def forecast( dataset = dataset.copy(deep=True) precip_var = dataset.attrs["precip_var"] - precip = dataset[precip_var].values[0] + precip = dataset[precip_var].values[-1] velocity = np.stack([dataset["velocity_x"], dataset["velocity_y"]]) _check_inputs(precip, velocity, timesteps)