From eec8eed42e7387eb4244ef7cd4dcb379e3a50690 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff <31476760+RubenImhoff@users.noreply.github.com> Date: Mon, 23 Aug 2021 17:18:11 +0200 Subject: [PATCH] Update to steps_blending branch * First basic functions to implement STEPS blending * Add compute of blend means,sigmas and recompose * pysteps.io with xarray (#219) * Add xarray dependency * MCH importer returns an xarray Dataset * Remove plot lines * Remove import * Adapt readers to xarray format * Rewrite as more general decorator * Add missing metadata * Adapt io tests * Mrms bounding box (#222) * Fix bounding box coordinates * Add missing metadata * Import xarray DataArray Ignore quality field * Black * Do not hardcode metadata * Address review comments by ruben * Add a legacy option to the io functions A "legacy" options is added to revert back the importers and readers behavior to version 1. This is a temporary solution to allow the examples, and other functions, to run as usual (v1.*). Hopefully, this is will allow a easier transition into version 2 during the development process and will allow testing functions that were not updated to v2. * Fix compatibility problems with tests Many of the tests were updated to use the legacy data structures (v1). The tests that still contains issues, were tagged with a TODO message and they are skipped. This will allow incremental changes to be tested in the new v2. IMPORTANT: once the v2 branch is stable, we may remove the legacy compatibility from the code base and the tests. * Update dependencies * Ignore plugins test Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> * Add blend_optical_flow * changes to steps blending procedure - weights according to adjusted BPS2006 method * changes to blending procedures - adjust weights from original BPS2006 method * Determine spatial correlation of NWP model forecast * First attempt to make correlations and thus weights lead time dependent (in progress..) * Change back to original BPS2006 blending formulation and add regression of skill values to climatological values for weights determination * Reformat code with Black * Skill score script imports climatological correlation-values from file now * Small changes to skill score script * Add skill score tests and an interface * Add skill score tests and an interface * Small change to docstring * Bom import xarray (#228) * Add import_bom_rf3 using xarray * Add tests to xarray version * Fix mrms importer tests * Pass **kwargs to internal functions * Add nwp_importers to read bom nwp sample data * Add bom nwp data to source file * Add tests for bom_nwp reader * Fix pystepsrc Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> * Functions to store and compute climatological weights (#231) * Implement the functions get_default_weights, save_weights, calc_clim_weights. These functions are used to evolve the weights in the scale- and skill-dependent blending with NWP in the STEPS blending algorithm. The current weights, based on the correlations per cascade level, are regressed towards these climatological weights in the course of the forecast. These functions save the current and compute the climatological weights (a running mean of the weights of the past n days, where typically n=30). First daily averages are stored and these are then averaged over the running window of n days. * Add tests for pysteps climatological weight io and calculations. * Add path_workdir to outputs section in pystepsrc file and use it as a default path to store/retrieve blending weights. * Minor changes to docstrings, changes to skill scores and testing scripts * Completed documentation for blending clim module, cleanup. Co-authored-by: RubenImhoff Co-authored-by: Carlos Velasco Co-authored-by: ned Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> Co-authored-by: Ruben Imhoff Co-authored-by: Carlos Velasco Co-authored-by: Lesley De Cruz --- .github/workflows/test_pysteps.yml | 8 +- .pre-commit-config.yaml | 2 +- environment.yml | 1 + environment_dev.yml | 1 + examples/LK_buffer_mask.py | 4 +- examples/advection_correction.py | 2 +- examples/anvil_nowcast.py | 6 +- examples/data_transformations.py | 2 +- examples/optical_flow_methods_convergence.py | 2 +- examples/plot_cascade_decomposition.py | 2 +- examples/plot_ensemble_verification.py | 4 +- examples/plot_extrapolation_nowcast.py | 4 +- examples/plot_noise_generators.py | 4 +- examples/plot_optical_flow.py | 2 +- examples/plot_steps_nowcast.py | 2 +- examples/probability_forecast.py | 4 +- examples/rainfarm_downscale.py | 2 +- .../thunderstorm_detection_and_tracking.py | 2 +- pysteps/blending/__init__.py | 4 + pysteps/blending/clim.py | 188 ++++++++++++ pysteps/blending/interface.py | 56 ++++ pysteps/blending/skill_scores.py | 250 +++++++++++++++ pysteps/blending/steps.py | 248 +++++++++++++++ pysteps/blending/utils.py | 138 +++++++++ pysteps/decorators.py | 107 ++++++- pysteps/io/__init__.py | 1 + pysteps/io/importers.py | 165 +++++++++- pysteps/io/nwp_importers.py | 277 +++++++++++++++++ pysteps/io/readers.py | 66 ++-- pysteps/pystepsrc | 14 +- pysteps/tests/helpers.py | 118 +++---- pysteps/tests/test_blending_clim.py | 148 +++++++++ pysteps/tests/test_blending_skill_scores.py | 289 ++++++++++++++++++ pysteps/tests/test_cascade.py | 3 + pysteps/tests/test_datasets.py | 3 + pysteps/tests/test_detcatscores.py | 3 + pysteps/tests/test_exporters.py | 5 + pysteps/tests/test_io_archive.py | 35 +++ pysteps/tests/test_io_bom_nwp.py | 82 +++++ pysteps/tests/test_io_bom_rf3.py | 114 +++++-- pysteps/tests/test_io_fmi_pgm.py | 95 +++--- pysteps/tests/test_io_knmi_hdf5.py | 47 ++- pysteps/tests/test_io_mch_gif.py | 54 ++-- pysteps/tests/test_io_mrms_grib.py | 65 ++-- pysteps/tests/test_io_opera_hdf5.py | 45 +-- pysteps/tests/test_io_readers.py | 38 +++ pysteps/tests/test_io_saf_crri.py | 58 ++-- pysteps/tests/test_nowcasts_steps.py | 5 + pysteps/tests/test_plt_animate.py | 3 + pysteps/tests/test_plugins_support.py | 3 + pysteps/tests/test_probscores.py | 3 +- pysteps/tests/test_spatialscores.py | 3 +- .../tests/test_timeseries_autoregression.py | 4 +- requirements.txt | 2 + requirements_dev.txt | 2 + setup.py | 2 + 56 files changed, 2443 insertions(+), 354 deletions(-) create mode 100644 pysteps/blending/__init__.py create mode 100644 pysteps/blending/clim.py create mode 100644 pysteps/blending/interface.py create mode 100644 pysteps/blending/skill_scores.py create mode 100644 pysteps/blending/steps.py create mode 100644 pysteps/blending/utils.py create mode 100644 pysteps/io/nwp_importers.py create mode 100644 pysteps/tests/test_blending_clim.py create mode 100644 pysteps/tests/test_blending_skill_scores.py create mode 100644 pysteps/tests/test_io_archive.py create mode 100644 pysteps/tests/test_io_bom_nwp.py create mode 100644 pysteps/tests/test_io_readers.py diff --git a/.github/workflows/test_pysteps.yml b/.github/workflows/test_pysteps.yml index 41d8c8426..ecb0db7b9 100644 --- a/.github/workflows/test_pysteps.yml +++ b/.github/workflows/test_pysteps.yml @@ -8,9 +8,13 @@ env: on: # Triggers the workflow on push or pull request events but only for the master branch push: - branches: [ master ] + branches: + - master + - pysteps-v2 pull_request: - branches: [ master ] + branches: + - master + - pysteps-v2 jobs: unit_tests: diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c80571951..f35ad636b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,6 +1,6 @@ repos: - repo: https://github.com/psf/black - rev: 21.6b0 + rev: 21.7b0 hooks: - id: black language_version: python3 \ No newline at end of file diff --git a/environment.yml b/environment.yml index b71cd8de3..20609195b 100644 --- a/environment.yml +++ b/environment.yml @@ -13,3 +13,4 @@ dependencies: - pillow - pyproj - scipy + - xarray diff --git a/environment_dev.yml b/environment_dev.yml index b29e347fe..34b894cd4 100644 --- a/environment_dev.yml +++ b/environment_dev.yml @@ -30,3 +30,4 @@ dependencies: - cartopy>=0.18 - scikit-image - pandas + - xarray diff --git a/examples/LK_buffer_mask.py b/examples/LK_buffer_mask.py index 35efd5297..821f54f19 100644 --- a/examples/LK_buffer_mask.py +++ b/examples/LK_buffer_mask.py @@ -55,7 +55,7 @@ # Read the radar composites importer = io.get_method(importer_name, "importer") -R, quality, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, quality, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) del quality # Not used @@ -210,7 +210,7 @@ ) # Read and convert the radar composites -R_o, _, metadata_o = io.read_timeseries(fns, importer, **importer_kwargs) +R_o, _, metadata_o = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) R_o, metadata_o = conversion.to_rainrate(R_o, metadata_o) # Compute Spearman correlation diff --git a/examples/advection_correction.py b/examples/advection_correction.py index 5bab1d9bd..706807dbb 100644 --- a/examples/advection_correction.py +++ b/examples/advection_correction.py @@ -60,7 +60,7 @@ # Read the radar composites importer = io.get_method(importer_name, "importer") -R, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, __, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to mm/h R, metadata = conversion.to_rainrate(R, metadata) diff --git a/examples/anvil_nowcast.py b/examples/anvil_nowcast.py index 57802ad4e..fe0f6ad0a 100644 --- a/examples/anvil_nowcast.py +++ b/examples/anvil_nowcast.py @@ -48,7 +48,7 @@ # Read the input time series importer = io.get_method(importer_name, "importer") rainrate_field, quality, metadata = io.read_timeseries( - filenames, importer, **importer_kwargs + filenames, importer, legacy=True, **importer_kwargs ) # Convert to rain rate (mm/h) @@ -118,8 +118,8 @@ date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=3 ) -refobs_field, quality, metadata = io.read_timeseries( - filenames, importer, **importer_kwargs +refobs_field, _, metadata = io.read_timeseries( + filenames, importer, legacy=True, **importer_kwargs ) refobs_field, metadata = utils.to_rainrate(refobs_field[-1], metadata) diff --git a/examples/data_transformations.py b/examples/data_transformations.py index 24962ec6f..8b601ac14 100644 --- a/examples/data_transformations.py +++ b/examples/data_transformations.py @@ -58,7 +58,7 @@ # Read the radar composites importer = io.get_method(importer_name, "importer") -Z, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +Z, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Keep only positive rainfall values Z = Z[Z > metadata["zerovalue"]].flatten() diff --git a/examples/optical_flow_methods_convergence.py b/examples/optical_flow_methods_convergence.py index 8dd7522fc..f6289450a 100644 --- a/examples/optical_flow_methods_convergence.py +++ b/examples/optical_flow_methods_convergence.py @@ -62,7 +62,7 @@ # Read the reference radar composite importer = io.get_method(importer_name, "importer") reference_field, quality, metadata = io.read_timeseries( - fns, importer, **importer_kwargs + fns, importer, legacy=True, **importer_kwargs ) del quality # Not used diff --git a/examples/plot_cascade_decomposition.py b/examples/plot_cascade_decomposition.py index 06e884156..49a3798d8 100644 --- a/examples/plot_cascade_decomposition.py +++ b/examples/plot_cascade_decomposition.py @@ -30,7 +30,7 @@ filename = os.path.join( root_path, "20160928", "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz" ) -R, _, metadata = io.import_fmi_pgm(filename, gzipped=True) +R, _, metadata = io.import_fmi_pgm(filename, gzipped=True, legacy=True) # Convert to rain rate R, metadata = conversion.to_rainrate(R, metadata) diff --git a/examples/plot_ensemble_verification.py b/examples/plot_ensemble_verification.py index 54d2500c8..799c43899 100644 --- a/examples/plot_ensemble_verification.py +++ b/examples/plot_ensemble_verification.py @@ -56,7 +56,7 @@ # Read the data from the archive importer = io.get_method(importer_name, "importer") -R, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to rain rate R, metadata = conversion.to_rainrate(R, metadata) @@ -142,7 +142,7 @@ ) # Read the observations -R_o, _, metadata_o = io.read_timeseries(fns, importer, **importer_kwargs) +R_o, _, metadata_o = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to mm/h R_o, metadata_o = conversion.to_rainrate(R_o, metadata_o) diff --git a/examples/plot_extrapolation_nowcast.py b/examples/plot_extrapolation_nowcast.py index 203a5909d..2e76eeeef 100644 --- a/examples/plot_extrapolation_nowcast.py +++ b/examples/plot_extrapolation_nowcast.py @@ -48,7 +48,7 @@ # Read the radar composites importer = io.get_method(importer_name, "importer") -Z, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +Z, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to rain rate R, metadata = conversion.to_rainrate(Z, metadata) @@ -113,7 +113,7 @@ num_next_files=n_leadtimes, ) # Read the radar composites -R_o, _, metadata_o = io.read_timeseries(fns, importer, **importer_kwargs) +R_o, _, metadata_o = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) R_o, metadata_o = conversion.to_rainrate(R_o, metadata_o, 223.0, 1.53) # Compute fractions skill score (FSS) for all lead times, a set of scales and 1 mm/h diff --git a/examples/plot_noise_generators.py b/examples/plot_noise_generators.py index dd392ca4f..a809eefa1 100644 --- a/examples/plot_noise_generators.py +++ b/examples/plot_noise_generators.py @@ -34,7 +34,9 @@ # Import the example radar composite root_path = rcparams.data_sources["mch"]["root_path"] filename = os.path.join(root_path, "20160711", "AQC161932100V_00005.801.gif") -R, _, metadata = io.import_mch_gif(filename, product="AQC", unit="mm", accutime=5.0) +R, _, metadata = io.import_mch_gif( + filename, product="AQC", unit="mm", accutime=5.0, legacy=True +) # Convert to mm/h R, metadata = conversion.to_rainrate(R, metadata) diff --git a/examples/plot_optical_flow.py b/examples/plot_optical_flow.py index e1bc6035b..89ed3b37e 100644 --- a/examples/plot_optical_flow.py +++ b/examples/plot_optical_flow.py @@ -47,7 +47,7 @@ # Read the radar composites importer = io.get_method(importer_name, "importer") -R, quality, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, quality, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) del quality # Not used diff --git a/examples/plot_steps_nowcast.py b/examples/plot_steps_nowcast.py index 43b01a6b4..38139b848 100644 --- a/examples/plot_steps_nowcast.py +++ b/examples/plot_steps_nowcast.py @@ -51,7 +51,7 @@ # Read the data from the archive importer = io.get_method(importer_name, "importer") -R, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to rain rate R, metadata = conversion.to_rainrate(R, metadata) diff --git a/examples/probability_forecast.py b/examples/probability_forecast.py index ee996c9ea..059115bc1 100644 --- a/examples/probability_forecast.py +++ b/examples/probability_forecast.py @@ -74,7 +74,7 @@ date = datetime.strptime("201607112100", "%Y%m%d%H%M") fns = io.find_by_date(date, root, fmt, pattern, ext, timestep, num_prev_files=2) importer = io.get_method(importer_name, "importer") -precip, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +precip, __, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) precip, metadata = utils.to_rainrate(precip, metadata) # precip[np.isnan(precip)] = 0 @@ -124,7 +124,7 @@ fns = io.find_by_date( date, root, fmt, pattern, ext, timestep, num_next_files=nleadtimes ) -obs, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +obs, __, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) obs, metadata = utils.to_rainrate(obs, metadata) obs[np.isnan(obs)] = 0 diff --git a/examples/rainfarm_downscale.py b/examples/rainfarm_downscale.py index 1f6a635f1..be4033989 100644 --- a/examples/rainfarm_downscale.py +++ b/examples/rainfarm_downscale.py @@ -33,7 +33,7 @@ root_path = rcparams.data_sources["mch"]["root_path"] filename = os.path.join(root_path, "20160711", "AQC161932100V_00005.801.gif") precip, _, metadata = io.import_mch_gif( - filename, product="AQC", unit="mm", accutime=5.0 + filename, product="AQC", unit="mm", accutime=5.0, legacy=True ) # Convert to mm/h diff --git a/examples/thunderstorm_detection_and_tracking.py b/examples/thunderstorm_detection_and_tracking.py index 1c82ded82..de3cc39e9 100644 --- a/examples/thunderstorm_detection_and_tracking.py +++ b/examples/thunderstorm_detection_and_tracking.py @@ -63,7 +63,7 @@ date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_next_files=20 ) importer = io.get_method(importer_name, "importer") -R, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs) +R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs) # Convert to reflectivity (it is possible to give the a- and b- parameters of the # Marshall-Palmer relationship here: zr_a = and zr_b =). diff --git a/pysteps/blending/__init__.py b/pysteps/blending/__init__.py new file mode 100644 index 000000000..420a11bd0 --- /dev/null +++ b/pysteps/blending/__init__.py @@ -0,0 +1,4 @@ +# -*- coding: utf-8 -*- +"""Methods for blending NWP model(s) with nowcasts.""" + +from pysteps.blending.interface import get_method diff --git a/pysteps/blending/clim.py b/pysteps/blending/clim.py new file mode 100644 index 000000000..fbb03e1ff --- /dev/null +++ b/pysteps/blending/clim.py @@ -0,0 +1,188 @@ +""" +pysteps.blending.clim +===================== + +Module with methods to read, write and compute past and climatological model weights. + +.. autosummary:: + :toctree: ../generated/ + + get_default_weights + save_weights + calc_clim_weights +""" + +import numpy as np +from os.path import exists, join +import pickle +from pysteps import rcparams + + +def get_default_weights(n_cascade_levels=8, n_models=1): + """ + Get the default weights as given in BPS2004. + Take subset of n_cascade_levels or add entries with small values (1e-4) if + n_cascade_levels differs from 8. + + Parameters + ---------- + n_cascade_levels: int, optional + Number of cascade levels. Defaults to 8. + n_models: int, optional + Number of NWP models. Defaults to 1. + + Returns + ------- + default_weights: array-like + Array of shape [model, scale_level] containing the climatological weights. + + """ + + default_weights = np.array( + [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040] + ) + n_weights = default_weights.shape[0] + if n_cascade_levels < n_weights: + default_weights = default_weights[0:n_cascade_levels] + elif n_cascade_levels > n_weights: + default_weights = np.append( + default_weights, np.repeat(1e-4, n_cascade_levels - n_weights) + ) + return np.resize(default_weights, (n_models, n_cascade_levels)) + + +def save_weights( + current_weights, + validtime, + outdir_path=rcparams.outputs["path_workdir"], + window_length=30, +): + """ + Add the current NWP weights to update today's daily average weight. If the + day is over, update the list of daily average weights covering a rolling + window. + + Parameters + ---------- + current_weights: array-like + Array of shape [model, scale_level, ...] + containing the current weights of the different NWP models per cascade + level. + validtime: datetime + Datetime object containing the date and time at for which the current + weights are valid. + outdir_path: string, optional + Path to folder where the historical weights are stored. Defaults to + path_workdir from rcparams. + window_length: int, optional + Length of window (in days) of daily weights that should be retained. + Defaults to 30. + + Returns + ------- + None + + """ + + n_cascade_levels = current_weights.shape[1] + + # Load weights_today, a dictionary containing {mean_weights, n, last_validtime} + weights_today_file = join(outdir_path, "NWP_weights_today.pkl") + if exists(weights_today_file): + with open(weights_today_file, "rb") as f: + weights_today = pickle.load(f) + else: + weights_today = { + "mean_weights": np.copy(current_weights), + "n": 0, + "last_validtime": validtime, + } + + # Load the past weights which is an array with dimensions day x model x scale_level + past_weights_file = join(outdir_path, "NWP_weights_window.npy") + past_weights = None + if exists(past_weights_file): + past_weights = np.load(past_weights_file) + # First check if we have started a new day wrt the last written weights, in which + # case we should update the daily weights file and reset daily statistics. + if weights_today["last_validtime"].date() < validtime.date(): + # Append weights to the list of the past X daily averages. + if past_weights is not None: + past_weights = np.append( + past_weights, [weights_today["mean_weights"]], axis=0 + ) + else: + past_weights = np.array([weights_today["mean_weights"]]) + print(past_weights.shape) + # Remove oldest if the number of entries exceeds the window length. + if past_weights.shape[0] > window_length: + past_weights = np.delete(past_weights, 0, axis=0) + # FIXME also write out last_validtime.date() in this file? + # In that case it will need pickling or netcdf. + # Write out the past weights within the rolling window. + np.save(past_weights_file, past_weights) + # Reset statistics for today. + weights_today["n"] = 0 + weights_today["mean_weights"] = 0 + + # Reset today's weights if needed and/or compute online average from the + # current weights using numerically stable algorithm + weights_today["n"] += 1 + weights_today["mean_weights"] += ( + current_weights - weights_today["mean_weights"] + ) / weights_today["n"] + weights_today["last_validtime"] = validtime + with open(weights_today_file, "wb") as f: + pickle.dump(weights_today, f) + + return None + + +def calc_clim_weights( + n_cascade_levels=8, + outdir_path=rcparams.outputs["path_workdir"], + n_models=1, + window_length=30, +): + """ + Return the climatological weights based on the daily average weights in the + rolling window. This is done using a geometric mean. + + Parameters + ---------- + n_cascade_levels: int, optional + Number of cascade levels. + outdir_path: string, optional + Path to folder where the historical weights are stored. Defaults to + path_workdir from rcparams. + n_models: int, optional + Number of NWP models. Defaults to 1. + window_length: int, optional + Length of window (in days) over which to compute the climatological + weights. Defaults to 30. + + Returns + ------- + climatological_mean_weights: array-like + Array of shape [model, scale_level, ...] containing the climatological + (geometric) mean weights. + + """ + past_weights_file = join(outdir_path, "NWP_weights_window.npy") + # past_weights has dimensions date x model x scale_level x .... + if exists(past_weights_file): + past_weights = np.load(past_weights_file) + else: + past_weights = None + # check if there's enough data to compute the climatological skill + if not past_weights or past_weights.shape[0] < window_length: + return get_default_weights(n_cascade_levels, n_models) + # reduce window if necessary + else: + past_weights = past_weights[-window_length:] + + # Calculate climatological weights from the past_weights using the + # geometric mean. + geomean_weights = np.exp(np.log(past_weights).mean(axis=0)) + + return geomean_weights diff --git a/pysteps/blending/interface.py b/pysteps/blending/interface.py new file mode 100644 index 000000000..0486188c7 --- /dev/null +++ b/pysteps/blending/interface.py @@ -0,0 +1,56 @@ +# -*- coding: utf-8 -*- +""" +pysteps.blending.interface +========================== + +Interface for the blending module. It returns a callable function for computing +blended nowcasts with NWP models. + +.. autosummary:: + :toctree: ../generated/ + + get_method +""" + +from pysteps.blending import steps + +_blending_methods = dict() +_blending_methods["steps"] = steps.forecast + + +def get_method(name): + """Return a callable function for computing nowcasts. + + Description: + Return a callable function for computing deterministic or ensemble + precipitation nowcasts. + + Implemented methods: + + +-----------------+-------------------------------------------------------+ + | Name | Description | + +=================+=======================================================+ + +-----------------+-------------------------------------------------------+ + | steps | the STEPS stochastic nowcasting blending method | + | | described in :cite:`Seed2003`, :cite:`BPS2006` and | + | | :cite:`SPN2013`. The blending weights approach | + | | currently follows :cite:`BPS2006`. | + +-----------------+-------------------------------------------------------+ + """ + if isinstance(name, str): + name = name.lower() + else: + raise TypeError( + "Only strings supported for the method's names.\n" + + "Available names:" + + str(list(_blending_methods.keys())) + ) from None + + try: + return _blending_methods[name] + except KeyError: + raise ValueError( + "Unknown blending method {}\n".format(name) + + "The available methods are:" + + str(list(_blending_methods.keys())) + ) from None diff --git a/pysteps/blending/skill_scores.py b/pysteps/blending/skill_scores.py new file mode 100644 index 000000000..d53ac3548 --- /dev/null +++ b/pysteps/blending/skill_scores.py @@ -0,0 +1,250 @@ +# -*- coding: utf-8 -*- +""" +pysteps.blending.skill_scores +============================== + +Methods for computing skill scores, needed for the blending weights, of two- +dimensional model fields with the latest observation field. + +.. autosummary:: + :toctree: ../generated/ + + spatial_correlation + lt_dependent_cor_nwp + lt_dependent_cor_extrapolation + clim_regr_values +""" + +import numpy as np +from pysteps.blending import clim + + +def spatial_correlation(obs, mod): + """Determine the spatial correlation between the cascade of the latest + available observed (radar) rainfall field and a time-synchronous cascade + derived from a model (generally NWP) field. Both fields are assumed to use + the same grid. + + + Parameters + ---------- + obs : array-like + Array of shape [cascade_level, y, x] with per cascade_level the + normalized cascade of the observed (radar) rainfall field. + mod : array-like + Array of shape [cascade_level, y, x] with per cascade_level the + normalized cascade of the model field. + + Returns + ------- + rho : array-like + Array of shape [n_cascade_levels] containing per cascade_level the + correlation between the normalized cascade of the observed (radar) + rainfall field and the normalized cascade of the model field. + + References + ---------- + :cite:`BPS2006` + :cite:`SPN2013` + + """ + rho = [] + # Fill rho per cascade level, so loop through the cascade levels + for cascade_level in range(0, obs.shape[0]): + # Flatten both arrays + obs_1d = obs[cascade_level, :, :].flatten() + mod_1d = mod[cascade_level, :, :].flatten() + # Calculate the correlation between the two + cov = np.sum( + (mod_1d - np.mean(mod_1d)) * (obs_1d - np.mean(obs_1d)) + ) # Without 1/n, as this cancels out (same for stdx and -y) + std_obs = np.sqrt(np.sum((obs_1d - np.mean(obs_1d)) ** 2.0)) + std_mod = np.sqrt(np.sum((mod_1d - np.mean(mod_1d)) ** 2.0)) + rho.append(cov / (std_mod * std_obs)) + + return rho + + +def lt_dependent_cor_nwp(lt, correlations, skill_kwargs=dict()): + """Determine the correlation of a model field for lead time lt and + cascade k, by assuming that the correlation determined at t=0 regresses + towards the climatological values. + + + Parameters + ---------- + lt : int + The lead time of the forecast in minutes. + correlations : array-like + Array of shape [n_cascade_levels] containing per cascade_level the + correlation between the normalized cascade of the observed (radar) + rainfall field and the normalized cascade of the model field. + skill_kwargs : dict, optional + Dictionary containing e.g. the outdir_path, nmodels and window_length + parameters. + + Returns + ------- + rho : array-like + Array of shape [n_cascade_levels] containing, for lead time lt, per + cascade_level the correlation between the normalized cascade of the + observed (radar) rainfall field and the normalized cascade of the + model field. + + References + ---------- + :cite:`BPS2004` + :cite:`BPS2006` + """ + # Obtain the climatological values towards which the correlations will + # regress + clim_cor_values, regr_pars = clim_regr_values(n_cascade_levels=len(correlations)) + # Determine the speed of the regression (eq. 24 in BPS2004) + qm = np.exp(-lt / regr_pars[0, :]) * (2 - np.exp(-lt / regr_pars[1, :])) + # Determine the correlation for lead time lt + rho = qm * correlations + (1 - qm) * clim_cor_values + + return rho + + +def lt_dependent_cor_extrapolation(PHI, correlations=None, correlations_prev=None): + """Determine the correlation of the extrapolation (nowcast) component for + lead time lt and cascade k, by assuming that the correlation determined at + t=0 regresses towards the climatological values. + + + Parameters + ---------- + PHI : array-like + Array of shape [n_cascade_levels, ar_order + 1] containing per + cascade level the autoregression parameters. + correlations : array-like, optional + Array of shape [n_cascade_levels] containing per cascade_level the + latest available correlation from the extrapolation component that can + be found from the AR-2 model. + correlations_prev : array-like, optional + Similar to correlations, but from the timestep before that. + + Returns + ------- + rho : array-like + Array of shape [n_cascade_levels] containing, for lead time lt, per + cascade_level the correlation of the extrapolation component. + + References + ---------- + :cite:`BPS2004` + :cite:`BPS2006` + + """ + # Check if correlations_prev exists, if not, we set it to 1.0 + if correlations_prev is None: + correlations_prev = np.repeat(1.0, PHI.shape[0]) + # Same for correlations at first time step, we set it to + # phi1 / (1 - phi2), see BPS2004 + if correlations is None: + correlations = PHI[:, 0] / (1.0 - PHI[:, 1]) + + # Calculate the correlation for lead time lt + rho = PHI[:, 0] * correlations + PHI[:, 1] * correlations_prev + + # Finally, set the current correlations array as the previous one for the + # next time step + rho_prev = correlations + + return rho, rho_prev + + +# TODO: Make sure the method can also handle multiple model inputs.. +def clim_regr_values(n_cascade_levels, skill_kwargs=dict()): + """Obtains the climatological correlation values and regression parameters + from a file called NWP_weights_window.bin in the outdir_path. If this file + is not present yet, the values from :cite:`BPS2004` are used. + + + Parameters + ---------- + n_cascade_levels : int + The number of cascade levels to use. + skill_kwargs : dict, optional + Dictionary containing e.g. the outdir_path, nmodels and window_length + parameters. + + Returns + ------- + clim_cor_values : array-like + Array of shape [n_cascade_levels] containing the + climatological values of the lag 1 and lag 2 auto-correlation + coefficients, obtained by calling a method implemented in + pysteps.blending.skill_scores.get_clim_skill_scores. + regr_pars : array-like + Array of shape [2, n_cascade_levels] containing the regression + parameters. These are fixed values that should be hard-coded in this + function. Unless changed by the user, the standard values from + `BPS2004` are used. + + Notes + ----- + The literature climatological values assume 8 cascade levels. In case + less or more cascade levels are used, the script will handle this by taking + the first n values or extending the array with a small value. This is not + ideal, but will be fixed once the clim_regr_file is made. Hence, this + requires a warm-up period of the model. + In addition, the regression parameter values (eq. 24 in BPS2004) are hard- + coded and can only be optimized by the user after (re)fitting of the + equation. + + """ + # First, obtain climatological skill values + try: + clim_cor_values = clim.calc_clim_weights( + n_cascade_levels=n_cascade_levels, **skill_kwargs + ) + except FileNotFoundError: + # The climatological skill values file does not exist yet, so we'll + # use the default values from BPS2004. + clim_cor_values = np.array( + [[0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040]] + ) + + clim_cor_values = clim_cor_values[0] + + # Check if clim_cor_values has only one model, otherwise it has + # returned the skill values for multiple models + if clim_cor_values.ndim != 1: + raise IndexError( + "The climatological cor. values of multiple models were returned, but only one model should be specified. Make sure to pass the argument nmodels in the function" + ) + + # Also check whether the number of cascade_levels is correct + if clim_cor_values.shape[0] > n_cascade_levels: + clim_cor_values = clim_cor_values[0:n_cascade_levels] + elif clim_cor_values.shape[0] < n_cascade_levels: + # Get the number of cascade levels that is missing + n_extra_lev = n_cascade_levels - clim_cor_values.shape[0] + # Append the array with correlation values of 10e-4 + clim_cor_values = np.append(clim_cor_values, np.repeat(1e-4, n_extra_lev)) + + # Get the regression parameters (L in eq. 24 in BPS2004) - Hard coded, + # change to own values when present. + regr_pars = np.array( + [ + [130.0, 165.0, 120.0, 55.0, 50.0, 15.0, 15.0, 10.0], + [155.0, 220.0, 200.0, 75.0, 10e4, 10e4, 10e4, 10e4], + ] + ) + + # Check whether the number of cascade_levels is correct + if regr_pars.shape[1] > n_cascade_levels: + regr_pars = regr_pars[:, 0:n_cascade_levels] + elif regr_pars.shape[1] < n_cascade_levels: + # Get the number of cascade levels that is missing + n_extra_lev = n_cascade_levels - regr_pars.shape[1] + # Append the array with correlation values of 10e-4 + regr_pars = np.append( + regr_pars, + [np.repeat(10.0, n_extra_lev), np.repeat(10e4, n_extra_lev)], + axis=1, + ) + + return clim_cor_values, regr_pars diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py new file mode 100644 index 000000000..2c053dc14 --- /dev/null +++ b/pysteps/blending/steps.py @@ -0,0 +1,248 @@ +# -*- coding: utf-8 -*- +""" +pysteps.blending.steps +====================== + +Implementation of the STEPS stochastic blending method as described in +:cite:`BPS2006`. The method assumes the presence of one NWP model or ensemble +member to be blended with one nowcast. More models, such as in :cite:`SPN2013` +is possible with this code, but we recommend the use of just two models. + +.. autosummary:: + :toctree: ../generated/ + + forecast + calculate_ratios + calculate_weights + blend_cascades + blend_means_sigmas + recompose_cascade + +References +---------- +:cite:`BPS2004` +:cite:`BPS2006` +:cite:`SPN2013` +""" + +import numpy as np +from pysteps.blending import utils + + +def forecast(): + # TODO: .. + + return ... + + +def calculate_ratios(correlations): + """Calculate explained variance ratios from correlation. + + Parameters + ---------- + Array of shape [component, scale_level, ...] + containing correlation (skills) for each component (NWP and nowcast), + scale level, and optionally along [y, x] dimensions. + + Returns + ------- + out : numpy array + An array containing the ratios of explain variance for each + component, scale level, ... + """ + # correlations: [component, scale, ...] + square_corrs = np.square(correlations) + # Calculate the ratio of the explained variance to the unexplained + # variance of the nowcast and NWP model components + out = square_corrs / (1 - square_corrs) + # out: [component, scale, ...] + return out + + +def calculate_weights(correlations): + """Calculate blending weights for STEPS blending from correlation. + + Parameters + ---------- + correlations : array-like + Array of shape [component, scale_level, ...] + containing correlation (skills) for each component (NWP and nowcast), + scale level, and optionally along [y, x] dimensions. + + Returns + ------- + weights : array-like + Array of shape [component+1, scale_level, ...] + containing the weights to be used in STEPS blending for + each original component plus an addtional noise component, scale level, + and optionally along [y, x] dimensions. + + Notes + ----- + The weights in the BPS method can sum op to more than 1.0. Hence, the + blended cascade has the be (re-)normalized (mu = 0, sigma = 1.0) first + before the blended cascade can be recomposed. + """ + # correlations: [component, scale, ...] + # Check if the correlations are positive, otherwise rho = 10e-5 + correlations = np.where(correlations < 10e-5, 10e-5, correlations) + # Calculate weights for each source + ratios = calculate_ratios(correlations) + # ratios: [component, scale, ...] + total_ratios = np.sum(ratios, axis=0) + # total_ratios: [scale, ...] - the denominator of eq. 11 & 12 in BPS2006 + weights = correlations * np.sqrt(ratios / total_ratios) + # weights: [component, scale, ...] + # Calculate the weight of the noise component. + # Original BPS2006 method in the following two lines (eq. 13) + total_square_weights = np.sum(np.square(weights), axis=0) + noise_weight = np.sqrt(1.0 - total_square_weights) + # TODO: determine the weights method and/or add different functions + + # Finally, add the noise_weights to the weights variable. + weights = np.concatenate((weights, noise_weight[None, ...]), axis=0) + return weights + + +# TODO: Make sure that where the radar rainfall data has no data, the NWP +# data is used. +def blend_cascades(cascades_norm, weights): + """Calculate blended normalized cascades using STEPS weights following eq. + 10 in :cite:`BPS2006`. + + Parameters + ---------- + cascades_norm : array-like + Array of shape [number_components + 1, scale_level, ...] + with normalized cascades components 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 + + 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. + """ + # cascade_norm component, scales, y, x + # weights component, scales, .... + all_c_wn = weights * cascades_norm + combined_cascade = np.sum(all_c_wn, axis=0) + # combined_cascade [scale, ...] + return combined_cascade + + +def blend_means_sigmas(means, sigmas, weights): + """Calculate the blended means and sigmas, the normalization parameters + needed to recompose the cascade. This procedure uses the weights of the + blending of the normalized cascades and follows eq. 32 and 33 in BPS2004. + + + Parameters + ---------- + means : array-like + Array of shape [number_components, scale_level, ...] + with the mean for each component (NWP, nowcasts, noise). + sigmas : array-like + Array of shape [number_components, scale_level, ...] + with the standard deviation for each component. + 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 + + Returns + ------- + combined_means : array-like + An array of shape [scale_level, ...] + containing per scale level (cascade) the weighted combination of + means from multiple components (NWP, nowcasts and noise). + combined_sigmas : array-like + An array of shape [scale_level, ...] + similar to combined_means, but containing the standard deviations. + + """ + # Check if the dimensions are the same + diff_dims = weights.ndim - means.ndim + if diff_dims: + for i in range(diff_dims): + means = np.expand_dims(means, axis=means.ndim) + diff_dims = weights.ndim - sigmas.ndim + if diff_dims: + for i in range(diff_dims): + sigmas = np.expand_dims(sigmas, axis=sigmas.ndim) + # Weight should have one component more (the noise component) than the + # means and sigmas. Check this + if ( + weights.shape[0] - means.shape[0] != 1 + or weights.shape[0] - sigmas.shape[0] != 1 + ): + raise ValueError( + "The weights array does not have one (noise) component more than mu and sigma" + ) + else: + # Throw away the last component, which is the noise component + weights = weights[:-1] + + # Combine (blend) the means and sigmas + combined_means = np.zeros(weights.shape[1]) + combined_sigmas = np.zeros(weights.shape[1]) + total_weight = np.sum((weights), axis=0) + for i in range(weights.shape[0]): + combined_means += (weights[i] / total_weight) * means[i] + combined_sigmas += (weights[i] / total_weight) * sigmas[i] + # TODO: substract covarainces to weigthed sigmas - still necessary? + + return combined_means, combined_sigmas + + +def recompose_cascade(combined_cascade, combined_mean, combined_sigma): + """Recompose the cascades into a transformed rain rate field. + + + 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. + combined_mean : array-like + 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. + + Returns + ------- + out: array-like + A two-dimensional array containing the recomposed cascade. + + Notes + ----- + The combined_cascade is made with weights that do not have to sum up to + 1.0. Therefore, the combined_cascade is first normalized again using + pysteps.blending.steps.normalize_cascade. + """ + # First, normalize the combined_cascade again + combined_cascade = utils.normalize_cascade(combined_cascade) + # Now, renormalize it with the blended sigma and mean values + renorm = (combined_cascade * combined_sigma) + combined_mean + # print(renorm.shape) + out = np.sum(renorm, axis=0) + # print(out.shape) + return out diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py new file mode 100644 index 000000000..041e9e1b6 --- /dev/null +++ b/pysteps/blending/utils.py @@ -0,0 +1,138 @@ +# -*- coding: utf-8 -*- +""" +pysteps.blending.utils +====================== + +Module with common utilities used by blending methods. + +.. autosummary:: + :toctree: ../generated/ + + stack_cascades + normalize_cascade + blend_optical_flows +""" + +import numpy as np + + +def stack_cascades(R_d, donorm=True): + """Stack the given cascades into a larger array. + + Parameters + ---------- + R_d : list + 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.ones(n_levels) * 0.0 + sigma_ = np.ones(n_levels) * 1.0 + if donorm: + mu_ = np.asarray(cascade["means"]) + sigma_ = np.asarray(cascade["stds"]) + for j in range(n_levels): + R__ = (R_i[j, :, :] - mu_[j]) / sigma_[j] + R_.append(R__) + 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 normalize_cascade(cascade): + """Normalizes a cascade (again). + + Parameters + ---------- + cascade : array-like + An array of shape [scale_level, y, x] + containing per scale level a cascade that has to be normalized (again). + + Returns + ------- + out : array-like + The normalized cascade with the same shape as cascade. + + """ + # Determine the mean and standard dev. of the combined cascade + mu = np.mean(cascade, axis=(1, 2)) + sigma = np.std(cascade, axis=(1, 2)) + # Normalize the cascade + out = [(cascade[i] - mu[i]) / sigma[i] for i in range(cascade.shape[0])] + out = np.stack(out) + + return out + + +def blend_optical_flows(flows, weights): + """Combine advection fields using given weights. Following :cite:`BPS2006` + the second level of the cascade is used for the 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 + 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. + Returns + ------- + out: ndarray_ + Return the blended advection field having shape + (2, m, n), where out[0, :, :] contains the x-components of + the blended motion vectors and out[1, :, :] contains the y-components. + The velocities are in units of pixels / timestep. + """ + + # check inputs + if isinstance(flows, (list, tuple)): + flows = np.stack(flows) + + if isinstance(weights, (list, tuple)): + weights = np.asarray(weights) + + # check weights dimensions match number of sources + num_sources = flows.shape[0] + num_weights = weights.shape[0] + + if num_weights != num_sources: + raise ValueError( + "dimension mismatch between flows and weights.\n" + "weights dimension must match the number of flows.\n" + f"number of flows={num_sources}, number of weights={num_weights}" + ) + # normalize weigths + weights = weights / np.sum(weights) + + # flows dimension sources, 2, m, n + # weights dimension sources + # move source axis to last to allow broadcasting + all_c_wn = weights * np.moveaxis(flows, 0, -1) + # sum uses last axis + combined_flows = np.sum(all_c_wn, axis=-1) + # combined_flows [2, m, n] + return combined_flows diff --git a/pysteps/decorators.py b/pysteps/decorators.py index 7afefe846..b9cb401fb 100644 --- a/pysteps/decorators.py +++ b/pysteps/decorators.py @@ -16,6 +16,7 @@ """ import inspect import uuid +import xarray as xr from collections import defaultdict from functools import wraps @@ -67,7 +68,11 @@ def _postprocess_import(importer): @wraps(importer) def _import_with_postprocessing(*args, **kwargs): - precip, *other_args = importer(*args, **kwargs) + data_array = importer(*args, **kwargs) + + if not isinstance(data_array, xr.DataArray): + array, _, metadata = data_array + data_array = _to_xarray(array, metadata) _dtype = kwargs.get("dtype", dtype) @@ -78,18 +83,15 @@ def _import_with_postprocessing(*args, **kwargs): "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 + _fillna = kwargs.get("fillna", fillna) + if _fillna is not np.nan: + data_array = data_array.fillna(_fillna) + + data_array = data_array.astype(_dtype) - return (precip.astype(_dtype),) + tuple(other_args) + if kwargs.get("legacy", False): + return _xarray2legacy(data_array) + return data_array extra_kwargs_doc = """ Other Parameters @@ -283,3 +285,84 @@ def _func_with_cache(*args, **kwargs): return _func_with_cache return _memoize + + +def _to_xarray(array, metadata): + """Convert to xarray DataArray.""" + x1, x2 = metadata["x1"], metadata["x2"] + y1, y2 = metadata["y1"], metadata["y2"] + xsize, ysize = metadata["xpixelsize"], metadata["ypixelsize"] + # x_coords = np.arange(x1, x2, xsize) + xsize / 2 + # y_coords = np.arange(y1, y2, ysize) + ysize / 2 + x_coords = np.arange(x1, x1 + xsize * array.shape[1], xsize) + xsize / 2 + y_coords = np.arange(y1, y1 + ysize * array.shape[0], ysize) + ysize / 2 + + data_array = xr.DataArray( + data=array, + dims=("y", "x"), + coords=dict( + x=("x", x_coords), + y=("y", y_coords), + ), + ) + + data_array.attrs.update( + { + "unit": metadata["unit"], + "accutime": metadata["accutime"], + "transform": metadata["transform"], + "zerovalue": metadata["zerovalue"], + "threshold": metadata["threshold"], + "zr_a": metadata.get("zr_a", None), + "zr_b": metadata.get("zr_b", None), + "institution": metadata["institution"], + "projection": metadata["projection"], + # + # TODO: Remove before final 2.0 version + "yorigin": metadata["yorigin"], + "xpixelsize": metadata["xpixelsize"], + "ypixelsize": metadata["ypixelsize"], + } + ) + + data_array.x.attrs.update( + { + "standard_name": "projection_x_coordinate", + "units": metadata["cartesian_unit"], + # + # TODO: Remove before final 2.0 version + "x1": x1, + "x2": x2, + } + ) + + data_array.y.attrs.update( + { + "standard_name": "projection_y_coordinate", + "units": metadata["cartesian_unit"], + # + # TODO: Remove before final 2.0 version + "y1": y1, + "y2": y2, + } + ) + + return data_array + + +# TODO: Remove before final 2.0 version +def _xarray2legacy(data_array): + """ + Convert the new DataArrays to the legacy format used in pysteps v1.* + """ + _array = data_array.values + + metadata = data_array.x.attrs.copy() + metadata.update(**data_array.y.attrs) + metadata.update(**data_array.attrs) + + if "t" in data_array.coords: + print(data_array["t"]) + metadata["timestamps"] = data_array["t"] + + return _array, None, metadata diff --git a/pysteps/io/__init__.py b/pysteps/io/__init__.py index 6f5fd1677..b37a2f5d4 100644 --- a/pysteps/io/__init__.py +++ b/pysteps/io/__init__.py @@ -8,3 +8,4 @@ from .importers import * from .nowcast_importers import * from .readers import * +from .nwp_importers import * diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index b55319a14..40d865c77 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -90,6 +90,7 @@ from functools import partial import numpy as np +import xarray as xr from matplotlib.pyplot import imread from pysteps.decorators import postprocess_import @@ -456,6 +457,164 @@ def import_bom_rf3(filename, **kwargs): return precip, None, metadata +@postprocess_import() +def import_bom_rf3_xr(filename, **kwargs): + """Import a NetCDF radar rainfall product from the BoM Rainfields3 + using xarray. + + Parameters + ---------- + filename: str + Name of the file to import. + + {extra_kwargs_doc} + + Returns + ------- + out_da : xr.DataArray + A xarray DataArray containing the rainfall field in mm/h imported + from the Bureau RF3 netcdf, the quality field and the metadata. The + quality field is currently set to None. + """ + + if not NETCDF4_IMPORTED: + raise MissingOptionalDependency( + "netCDF4 package is required to import BoM Rainfields3 products " + "but it is not installed" + ) + + ds = _import_bom_rf3_data_xr(filename, **kwargs,) + ds_meta = _import_bom_rf3_geodata_xr(ds, **kwargs,) + + # rename valid_time to t if exists + if 'valid_time' in ds_meta: + ds_meta = ds_meta.rename({'valid_time': 't'}) + + return ds_meta.precipitation + + +def _import_bom_rf3_data_xr(filename, **kwargs,): + + varname_time = kwargs.get('varname_time', 'valid_time') + # Tested in python3.6 and chunks did not work properly + # commenting next line until find the reason + # chunks = kwargs.get('chunks', {varname_time: 1}) + + ds_rainfall = xr.open_mfdataset( + filename, + combine='nested', + concat_dim=varname_time, + # chunks=chunks, + lock=False, + parallel=True, + ) + + return ds_rainfall + + +def _import_bom_rf3_geodata_xr(ds_in, + **kwargs, + ): + + # select a default varname if none is passed + varname = kwargs.get('varname', 'precipitation') + + # extract useful information + # projection + projdef = None + if "proj" in ds_in: + projection = ds_in.proj + if projection.grid_mapping_name == "albers_conical_equal_area": + projdef = "+proj=aea " + lon_0 = projection.longitude_of_central_meridian + projdef += f" +lon_0={lon_0:.3f}" + lat_0 = projection.latitude_of_projection_origin + projdef += f" +lat_0={lat_0:.3f}" + standard_parallel = projection.standard_parallel + projdef += f" +lat_1={standard_parallel[0]:.3f}" + projdef += f" +lat_2={standard_parallel[1]:.3f}" + + # get the accumulation period + valid_time = None + if "valid_time" in ds_in: + valid_time = ds_in.valid_time + + start_time = None + if "start_time" in ds_in: + start_time = ds_in.start_time + + time_step = None + if start_time is not None: + if valid_time is not None: + time_step = (valid_time - start_time).isel(valid_time=0) + time_step = time_step.values.astype('timedelta64[m]') + + # get the units of precipitation + units = None + if "units" in ds_in[varname].attrs: + units = ds_in[varname].units + if units in ("kg m-2", "mm"): + units = "mm" + ds_in[varname].attrs.update({'units': units}) + # get spatial boundaries and pixelsize + # move to meters if coordiantes in kilometers + if "units" in ds_in.x.attrs: + if ds_in.x.units == "km": + ds_in['x'] = ds_in.x*1000. + ds_in.x.attrs.update({'units': 'm'}) + ds_in['y'] = ds_in.y*1000. + ds_in.y.attrs.update({'units': 'm'}) + + xmin = ds_in.x.min().values + xmax = ds_in.x.max().values + ymin = ds_in.y.min().values + ymax = ds_in.y.max().values + xpixelsize = abs(ds_in.x[1] - ds_in.x[0]) + ypixelsize = abs(ds_in.y[1] - ds_in.y[0]) + + cartesian_unit = ds_in.x.units + + # Add metadata needed by pySTEPS as attrs in X and Y variables + + ds_in.x.attrs.update({ + # TODO: Remove before final 2.0 version + "x1": xmin, + "x2": xmax, + "cartesian_unit": cartesian_unit, + } + ) + + ds_in.y.attrs.update({ + # TODO: Remove before final 2.0 version + "y1": ymin, + "y2": ymax, + "cartesian_unit": cartesian_unit, + } + ) + + # Add metadata needed by pySTEPS as attrs in rainfall variable + da_rainfall = ds_in[varname].isel(valid_time=0) + + ds_in[varname].attrs.update( + {"transform": None, + "unit": units, # copy 'units' in 'unit' for legacy reasons + "projection": projdef, + "accutime": time_step, + "zr_a": None, + "zr_b": None, + "zerovalue": np.nanmin(da_rainfall), + "institution": "Commonwealth of Australia, Bureau of Meteorology", + "threshold": _get_threshold_value(da_rainfall.values), + # TODO(_import_bom_rf3_geodata_xr): Remove before final 2.0 version + "yorigin": "upper", + "xpixelsize": xpixelsize.values, + "ypixelsize": ypixelsize.values, + } + ) + + return ds_in + + def _import_bom_rf3_data(filename): ds_rainfall = netCDF4.Dataset(filename) if "precipitation" in ds_rainfall.variables.keys(): @@ -522,7 +681,8 @@ def _import_bom_rf3_geodata(filename): calendar = "standard" if "calendar" in times.ncattrs(): calendar = times.calendar - valid_time = netCDF4.num2date(times[:], units=times.units, calendar=calendar) + valid_time = netCDF4.num2date(times[:], units=times.units, + calendar=calendar) start_time = None if "start_time" in ds_rainfall.variables.keys(): @@ -530,7 +690,8 @@ def _import_bom_rf3_geodata(filename): calendar = "standard" if "calendar" in times.ncattrs(): calendar = times.calendar - start_time = netCDF4.num2date(times[:], units=times.units, calendar=calendar) + start_time = netCDF4.num2date(times[:], units=times.units, + calendar=calendar) time_step = None diff --git a/pysteps/io/nwp_importers.py b/pysteps/io/nwp_importers.py new file mode 100644 index 000000000..680bea081 --- /dev/null +++ b/pysteps/io/nwp_importers.py @@ -0,0 +1,277 @@ +""" +pysteps.io.nwp_importers +==================== + +Methods for importing files containing two-dimensional NWP mosaics. + +The methods in this module implement the following interface:: + + nwp_import_xxx(filename, optional arguments) + +where **xxx** is the name (or abbreviation) of the file format and filename +is the name of the input file. + +The output of each method is a xarray DataArray containing +forecast rainfall fields in mm/timestep and the metadata as attributes. + +The metadata should contain the following recommended key-value pairs: + +.. tabularcolumns:: |p{2cm}|L| + ++------------------+----------------------------------------------------------+ +| Key | Value | ++==================+==========================================================+ +| projection | PROJ.4-compatible projection definition | ++------------------+----------------------------------------------------------+ +| x1 | x-coordinate of the lower-left corner of the data raster | ++------------------+----------------------------------------------------------+ +| y1 | y-coordinate of the lower-left corner of the data raster | ++------------------+----------------------------------------------------------+ +| x2 | x-coordinate of the upper-right corner of the data raster| ++------------------+----------------------------------------------------------+ +| y2 | y-coordinate of the upper-right corner of the data raster| ++------------------+----------------------------------------------------------+ +| xpixelsize | grid resolution in x-direction | ++------------------+----------------------------------------------------------+ +| ypixelsize | grid resolution in y-direction | ++------------------+----------------------------------------------------------+ +| cartesian_unit | the physical unit of the cartesian x- and y-coordinates: | +| | e.g. 'm' or 'km' | ++------------------+----------------------------------------------------------+ +| yorigin | a string specifying the location of the first element in | +| | the data raster w.r.t. y-axis: | +| | 'upper' = upper border | +| | 'lower' = lower border | ++------------------+----------------------------------------------------------+ +| institution | name of the institution who provides the data | ++------------------+----------------------------------------------------------+ +| unit | the physical unit of the data: 'mm/h', 'mm' or 'dBZ' | ++------------------+----------------------------------------------------------+ +| transform | the transformation of the data: None, 'dB', 'Box-Cox' or | +| | others | ++------------------+----------------------------------------------------------+ +| accutime | the accumulation time in minutes of the data, float | ++------------------+----------------------------------------------------------+ +| threshold | the rain/no rain threshold with the same unit, | +| | transformation and accutime of the data. | ++------------------+----------------------------------------------------------+ +| zerovalue | the value assigned to the no rain pixels with the same | +| | unit, transformation and accutime of the data. | ++------------------+----------------------------------------------------------+ +| zr_a | the Z-R constant a in Z = a*R**b | ++------------------+----------------------------------------------------------+ +| zr_b | the Z-R exponent b in Z = a*R**b | ++------------------+----------------------------------------------------------+ + +Available Importers +------------------- + +.. autosummary:: + :toctree: ../generated/ + + import_bom_nwp +""" + +import numpy as np +import xarray as xr + +from pysteps.decorators import postprocess_import +from pysteps.exceptions import MissingOptionalDependency + +try: + import netCDF4 + + NETCDF4_IMPORTED = True +except ImportError: + NETCDF4_IMPORTED = False + + +@postprocess_import() +def import_bom_nwp_xr(filename, **kwargs): + """Import a NetCDF with NWP rainfall forecasts regrided to a BoM Rainfields3 + using xarray. + + Parameters + ---------- + filename: str + Name of the file to import. + + {extra_kwargs_doc} + + Returns + ------- + out_da : xr.DataArray + A xarray DataArray containing forecast rainfall fields in mm/timestep + imported from a netcdf, and the metadata. + """ + + if not NETCDF4_IMPORTED: + raise MissingOptionalDependency( + "netCDF4 package is required to import BoM NWP regrided rainfall " + "products but it is not installed" + ) + + ds = _import_bom_nwp_data_xr(filename, **kwargs) + ds_meta = _import_bom_nwp_geodata_xr(ds, **kwargs) + + # rename varname_time (def: time) to t + varname_time = kwargs.get('varname_time', 'time') + ds_meta = ds_meta.rename({varname_time: 't'}) + varname_time = 't' + + # if data variable is named accum_prcp + # it is assumed that NWP rainfall data is accumulated + # so it needs to be disagregated by time step + varname = kwargs.get('varname', 'accum_prcp') + if varname == 'accum_prcp': + print("Rainfall values are accumulated. Disagreagating by time step") + accum_prcp = ds_meta[varname] + precipitation = (accum_prcp - accum_prcp.shift({varname_time: 1})) + precipitation = precipitation.dropna(varname_time, 'all') + # update/copy attributes + precipitation.name = 'precipitation' + # copy attributes + precipitation.attrs.update({**accum_prcp.attrs}) + # update attributes + precipitation.attrs.update({'standard_name': "precipitation_amount"}) + else: + precipitation = ds_meta[varname] + + return precipitation + + +def _import_bom_nwp_data_xr(filename, **kwargs): + + varname_time = kwargs.get('varname_time', 'time') + chunks = kwargs.get('chunks', {varname_time: 1}) + + ds_rainfall = xr.open_mfdataset( + filename, + combine='nested', + concat_dim=varname_time, + chunks=chunks, + lock=False, + parallel=True, + ) + + return ds_rainfall + + +def _import_bom_nwp_geodata_xr(ds_in, + **kwargs, + ): + + varname = kwargs.get('varname', 'accum_prcp') + varname_time = kwargs.get('varname_time', 'time') + + # extract useful information + # projection + projdef = None + if "proj" in ds_in: + projection = ds_in.proj + if projection.grid_mapping_name == "albers_conical_equal_area": + projdef = "+proj=aea " + lon_0 = projection.longitude_of_central_meridian + projdef += f" +lon_0={lon_0:.3f}" + lat_0 = projection.latitude_of_projection_origin + projdef += f" +lat_0={lat_0:.3f}" + standard_parallel = projection.standard_parallel + projdef += f" +lat_1={standard_parallel[0]:.3f}" + projdef += f" +lat_2={standard_parallel[1]:.3f}" + + # get the accumulation period + time = ds_in[varname_time] + # shift the array to calculate the time step + delta_time = time - time.shift({varname_time: 1}) + # assuming first valid delta_time is representative of all time steps + time_step = delta_time[1] + time_step = time_step.values.astype('timedelta64[m]') + + # get the units of precipitation + units = None + if "units" in ds_in[varname].attrs: + units = ds_in[varname].units + if units in ("kg m-2", "mm"): + units = "mm" + ds_in[varname].attrs.update({'units': units}) + # get spatial boundaries and pixelsize + # move to meters if coordiantes in kilometers + if "units" in ds_in.x.attrs: + if ds_in.x.units == "km": + ds_in['x'] = ds_in.x*1000. + ds_in.x.attrs.update({'units': 'm'}) + ds_in['y'] = ds_in.y*1000. + ds_in.y.attrs.update({'units': 'm'}) + + xmin = ds_in.x.min().values + xmax = ds_in.x.max().values + ymin = ds_in.y.min().values + ymax = ds_in.y.max().values + xpixelsize = abs(ds_in.x[1] - ds_in.x[0]) + ypixelsize = abs(ds_in.y[1] - ds_in.y[0]) + + cartesian_unit = ds_in.x.units + + # Add metadata needed by pySTEPS as attrs in X and Y variables + + ds_in.x.attrs.update({ + # TODO: Remove before final 2.0 version + "x1": xmin, + "x2": xmax, + "cartesian_unit": cartesian_unit, + } + ) + + ds_in.y.attrs.update({ + # TODO: Remove before final 2.0 version + "y1": ymin, + "y2": ymax, + "cartesian_unit": cartesian_unit, + } + ) + + # Add metadata needed by pySTEPS as attrs in rainfall variable + da_rainfall = ds_in[varname].isel({varname_time: 0}) + + ds_in[varname].attrs.update( + {"transform": None, + "unit": units, # copy 'units' in 'unit' for legacy reasons + "projection": projdef, + "accutime": time_step, + "zr_a": None, + "zr_b": None, + "zerovalue": np.nanmin(da_rainfall), + "institution": "Commonwealth of Australia, Bureau of Meteorology", + "threshold": _get_threshold_value(da_rainfall.values), + # TODO(_import_bom_rf3_geodata_xr): Remove before final 2.0 version + "yorigin": "upper", + "xpixelsize": xpixelsize.values, + "ypixelsize": ypixelsize.values, + } + ) + + return ds_in + + +def _get_threshold_value(precip): + """ + Get the the rain/no rain threshold with the same unit, transformation and + accutime of the data. + If all the values are NaNs, the returned value is `np.nan`. + Otherwise, np.min(precip[precip > precip.min()]) is returned. + + Returns + ------- + threshold: float + """ + valid_mask = np.isfinite(precip) + if valid_mask.any(): + _precip = precip[valid_mask] + min_precip = _precip.min() + above_min_mask = _precip > min_precip + if above_min_mask.any(): + return np.min(_precip[above_min_mask]) + else: + return min_precip + else: + return np.nan diff --git a/pysteps/io/readers.py b/pysteps/io/readers.py index a5ef9a661..92599be11 100644 --- a/pysteps/io/readers.py +++ b/pysteps/io/readers.py @@ -12,16 +12,19 @@ """ import numpy as np +import xarray as xr +from pysteps.decorators import _xarray2legacy -def read_timeseries(inputfns, importer, **kwargs): + +def read_timeseries(input_fns, importer, **kwargs): """Read a time series of input files using the methods implemented in the - :py:mod:`pysteps.io.importers` module and stack them into a 3d array of - shape (num_timesteps, height, width). + :py:mod:`pysteps.io.importers` module and concatenate them into a 3d array + of shape (t, y, x). Parameters ---------- - inputfns: tuple + input_fns: tuple Input files returned by a function implemented in the :py:mod:`pysteps.io.archive` module. importer: function @@ -31,50 +34,45 @@ def read_timeseries(inputfns, importer, **kwargs): Returns ------- - out: tuple - A three-element tuple containing the read data and quality rasters and + out: xr.DataArray + A xarray DataArray containing the data rasters and associated metadata. If an input file name is None, the corresponding - precipitation and quality fields are filled with nan values. If all - input file names are None or if the length of the file name list is - zero, a three-element tuple containing None values is returned. - + fields are filled with nan values. + If all input file names are None or if the length of the file name list + is zero, None is returned. """ - + legacy = kwargs.get("legacy", False) + kwargs["legacy"] = False # check for missing data precip_ref = None - if all(ifn is None for ifn in inputfns): - return None, None, None + if all(ifn is None for ifn in input_fns): + return None else: - if len(inputfns[0]) == 0: - return None, None, None - for ifn in inputfns[0]: + if len(input_fns[0]) == 0: + return None + for ifn in input_fns[0]: if ifn is not None: - precip_ref, quality_ref, metadata = importer(ifn, **kwargs) + precip_ref = importer(ifn, **kwargs) break if precip_ref is None: - return None, None, None + return None precip = [] - quality = [] timestamps = [] - for i, ifn in enumerate(inputfns[0]): + for i, ifn in enumerate(input_fns[0]): if ifn is not None: - precip_, quality_, _ = importer(ifn, **kwargs) + precip_ = importer(ifn, **kwargs) precip.append(precip_) - quality.append(quality_) - timestamps.append(inputfns[1][i]) + timestamps.append(input_fns[1][i]) else: - precip.append(precip_ref * np.nan) - if quality_ref is not None: - quality.append(quality_ref * np.nan) - else: - quality.append(None) - timestamps.append(inputfns[1][i]) + precip.append(xr.full_like(precip_ref, np.nan)) + timestamps.append(input_fns[1][i]) + + precip = xr.concat(precip, "t") + precip = precip.assign_coords({"t": timestamps}) - # Replace this with stack? - precip = np.concatenate([precip_[None, :, :] for precip_ in precip]) - # TODO: Q should be organized as R, but this is not trivial as Q_ can be also None or a scalar - metadata["timestamps"] = np.array(timestamps) + if legacy: + return _xarray2legacy(precip) - return precip, quality, metadata + return precip diff --git a/pysteps/pystepsrc b/pysteps/pystepsrc index 7e2b2c9d7..4d7f2f878 100644 --- a/pysteps/pystepsrc +++ b/pysteps/pystepsrc @@ -4,7 +4,8 @@ "silent_import": false, "outputs": { // path_outputs : path where to save results (figures, forecasts, etc) - "path_outputs": "./" + "path_outputs": "./", + "path_workdir": "./tmp/" }, "plot": { // "motion_plot" : "streamplot" or "quiver" @@ -89,6 +90,17 @@ "importer_kwargs": { "gzipped": true } + }, + "bom_nwp": { + "root_path": "./nwp/bom", + "path_fmt": "%Y/%m/%d", + "fn_pattern": "%Y%m%d_%H00_regrid_short", + "fn_ext": "nc", + "importer": "bom_nwp", + "timestep": 10, + "importer_kwargs": { + "gzipped": true + } } } } diff --git a/pysteps/tests/helpers.py b/pysteps/tests/helpers.py index 76a07bf38..86c59b663 100644 --- a/pysteps/tests/helpers.py +++ b/pysteps/tests/helpers.py @@ -11,6 +11,7 @@ import pysteps as stp from pysteps import io, rcparams +from pysteps.decorators import _to_xarray, _xarray2legacy from pysteps.utils import aggregate_fields_space _reference_dates = dict() @@ -30,80 +31,82 @@ def get_precipitation_fields( metadata=False, upscale=None, source="mch", + legacy=True, log_transform=True, clip=None, **importer_kwargs, ): """ - Get a precipitation field from the archive to be used as reference. + Get a precipitation field from the archive to be used as reference. - Source: bom - Reference time: 2018/06/16 10000 UTC + Source: bom + Reference time: 2018/06/16 10000 UTC - Source: fmi - Reference time: 2016/09/28 1600 UTC + Source: fmi + Reference time: 2016/09/28 1600 UTC - Source: knmi - Reference time: 2010/08/26 0000 UTC + Source: knmi + Reference time: 2010/08/26 0000 UTC - Source: mch - Reference time: 2015/05/15 1630 UTC + Source: mch + Reference time: 2015/05/15 1630 UTC - Source: opera - Reference time: 2018/08/24 1800 UTC + Source: opera + Reference time: 2018/08/24 1800 UTC - Source: saf - Reference time: 2018/06/01 0700 UTC + Source: saf + Reference time: 2018/06/01 0700 UTC - Source: mrms - Reference time: 2019/06/10 0000 UTC + Source: mrms + Reference time: 2019/06/10 0000 UTC - Parameters - ---------- + Parameters + ---------- - num_prev_files: int, optional - Number of previous times (files) to return with respect to the - reference time. + num_prev_files: int, optional + Number of previous times (files) to return with respect to the + reference time. - num_next_files: int, optional - Number of future times (files) to return with respect to the - reference time. + num_next_files: int, optional + Number of future times (files) to return with respect to the + reference time. - return_raw: bool, optional - Do not preprocess the precipitation fields. False by default. - The pre-processing steps are: 1) Convert to mm/h, - 2) Mask invalid values, 3) Log-transform the data [dBR]. + return_raw: bool, optional + Do not preprocess the precipitation fields. False by default. + The pre-processing steps are: 1) Convert to mm/h, + 2) Mask invalid values, 3) Log-transform the data [dBR]. - metadata: bool, optional - If True, also return file metadata. + <<<<<<< HEAD + ======= + metadata: bool, optional + If True, also return file metadata. - clip: scalars (left, right, bottom, top), optional - The extent of the bounding box in data coordinates to be used to clip - the data. + clip: scalars (left, right, bottom, top), optional + The extent of the bounding box in data coordinates to be used to clip + the data. - upscale: float or None, optional - Upscale fields in space during the pre-processing steps. - If it is None, the precipitation field is not modified. - If it is a float, represents the length of the space window that is - used to upscale the fields. + >>>>>>> master + upscale: float or None, optional + Upscale fields in space during the pre-processing steps. + If it is None, the precipitation field is not modified. + If it is a float, represents the length of the space window that is + used to upscale the fields. - source: {"bom", "fmi" , "knmi", "mch", "opera", "saf", "mrms"}, optional - Name of the data source to be used. + source: {"bom", "fmi" , "knmi", "mch", "opera", "saf", "mrms"}, optional + Name of the data source to be used. - log_transform: bool - Whether to transform the output to dB. + log_transform: bool + Whether to transform the output to dB. - Other Parameters - ---------------- + Other Parameters + ---------------- - importer_kwargs : dict - Additional keyword arguments passed to the importer. + importer_kwargs : dict + Additional keyword arguments passed to the importer. - Returns - ------- - reference_field : array - - metadata : dict + Returns + ------- + data_array : xr.DataArray """ if source == "bom": @@ -161,9 +164,9 @@ def get_precipitation_fields( # Read the radar composites importer = io.get_method(importer_name, "importer") - reference_field, __, ref_metadata = io.read_timeseries( - fns, importer, **_importer_kwargs - ) + reference_field = io.read_timeseries(fns, importer, **_importer_kwargs) + + reference_field, _, ref_metadata = _xarray2legacy(reference_field) if not return_raw: @@ -199,10 +202,13 @@ def get_precipitation_fields( np.ma.set_fill_value(reference_field, -15.0) reference_field.data[reference_field.mask] = -15.0 - if metadata: - return reference_field, ref_metadata + if legacy or metadata: + if metadata: + return reference_field, ref_metadata + else: + return reference_field - return reference_field + return _to_xarray(reference_field, ref_metadata) def smart_assert(actual_value, expected, tolerance=None): diff --git a/pysteps/tests/test_blending_clim.py b/pysteps/tests/test_blending_clim.py new file mode 100644 index 000000000..e790984c7 --- /dev/null +++ b/pysteps/tests/test_blending_clim.py @@ -0,0 +1,148 @@ +# -*- coding: utf-8 -*- +import numpy as np +import pytest + +from pysteps.blending.clim import save_weights, calc_clim_weights +import random +from datetime import datetime, timedelta +from os.path import join, exists +import pickle +from numpy.testing import assert_array_equal + +random.seed(12356) +n_cascade_levels = 7 +model_names = ["alaro13", "arome13"] +default_start_weights = [0.8, 0.5] +""" Helper functions """ + + +def generate_random_weights(n_cascade_levels, n_models=1): + """ + Generate random weights which decay exponentially with scale. + """ + start_weights = np.array([random.uniform(0.5, 0.99) for i in range(n_models)]) + powers = np.arange(1, n_cascade_levels + 1) + return pow(start_weights[:, np.newaxis], powers) + + +def generate_fixed_weights(n_cascade_levels, n_models=1): + """ + Generate weights starting at default_start_weights which decay exponentially with scale. + """ + start_weights = np.resize(default_start_weights, n_models) + powers = np.arange(1, n_cascade_levels + 1) + return pow(start_weights[:, np.newaxis], powers) + + +""" Test arguments """ + +clim_arg_names = ("startdatestr", "enddatestr", "n_models", "expected_weights_today") + +test_enddates = ["20210701235500", "20210702000000", "20200930235500"] + +clim_arg_values = [ + ( + "20210701230000", + "20210701235500", + 1, + { + "mean_weights": generate_fixed_weights(n_cascade_levels), + "n": 12, + "last_validtime": datetime.strptime(test_enddates[0], "%Y%m%d%H%M%S"), + }, + ), + ( + "20210701235500", + "20210702000000", + 1, + { + "mean_weights": generate_fixed_weights(n_cascade_levels), + "n": 1, + "last_validtime": datetime.strptime(test_enddates[1], "%Y%m%d%H%M%S"), + }, + ), + ( + "20200801000000", + "20200930235500", + 1, + { + "mean_weights": generate_fixed_weights(n_cascade_levels), + "n": 288, + "last_validtime": datetime.strptime(test_enddates[2], "%Y%m%d%H%M%S"), + }, + ), + ( + "20210701230000", + "20210701235500", + 2, + { + "mean_weights": generate_fixed_weights(n_cascade_levels, 2), + "n": 12, + "last_validtime": datetime.strptime(test_enddates[0], "%Y%m%d%H%M%S"), + }, + ), + ( + "20210701230000", + "20210702000000", + 2, + { + "mean_weights": generate_fixed_weights(n_cascade_levels, 2), + "n": 1, + "last_validtime": datetime.strptime(test_enddates[1], "%Y%m%d%H%M%S"), + }, + ), + ( + "20200801000000", + "20200930235500", + 2, + { + "mean_weights": generate_fixed_weights(n_cascade_levels, 2), + "n": 288, + "last_validtime": datetime.strptime(test_enddates[2], "%Y%m%d%H%M%S"), + }, + ), +] + + +@pytest.mark.parametrize(clim_arg_names, clim_arg_values) +def test_save_weights( + startdatestr, enddatestr, n_models, expected_weights_today, tmpdir +): + """Test if the weights are saved correctly and the daily average is computed""" + + # get validtime + currentdate = datetime.strptime(startdatestr, "%Y%m%d%H%M%S") + enddate = datetime.strptime(enddatestr, "%Y%m%d%H%M%S") + timestep = timedelta(minutes=5) + + outdir_path = tmpdir + + while currentdate <= enddate: + current_weights = generate_fixed_weights(n_cascade_levels, n_models) + print("Saving weights: ", current_weights, currentdate, outdir_path) + save_weights(current_weights, currentdate, outdir_path, window_length=2) + currentdate += timestep + + weights_today_file = join(outdir_path, "NWP_weights_today.pkl") + assert exists(weights_today_file) + with open(weights_today_file, "rb") as f: + weights_today = pickle.load(f) + + # Check type + assert type(weights_today) == type({}) + assert "mean_weights" in weights_today + assert "n" in weights_today + assert "last_validtime" in weights_today + assert_array_equal( + weights_today["mean_weights"], expected_weights_today["mean_weights"] + ) + assert weights_today["n"] == expected_weights_today["n"] + assert weights_today["last_validtime"] == expected_weights_today["last_validtime"] + + +if __name__ == "__main__": + save_weights( + generate_fixed_weights(n_cascade_levels, 1), + datetime.strptime("20200801000000", "%Y%m%d%H%M%S"), + "/tmp/", + ) diff --git a/pysteps/tests/test_blending_skill_scores.py b/pysteps/tests/test_blending_skill_scores.py new file mode 100644 index 000000000..0bf4f1a98 --- /dev/null +++ b/pysteps/tests/test_blending_skill_scores.py @@ -0,0 +1,289 @@ +# -*- coding: utf-8 -*- + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal + +from pysteps.blending.skill_scores import ( + spatial_correlation, + lt_dependent_cor_nwp, + lt_dependent_cor_extrapolation, + clim_regr_values, +) + +# TODO: Fix tests for xarray fields + +# Set the climatological correlations values +clim_cor_values_8lev = np.array( + [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040] +) +clim_cor_values_6lev = np.array([0.848, 0.537, 0.237, 0.065, 0.020, 0.0044]) +clim_cor_values_9lev = np.array( + [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040, 1e-4] +) + +# Set the regression values +regr_pars_8lev = np.array( + [ + [130.0, 165.0, 120.0, 55.0, 50.0, 15.0, 15.0, 10.0], + [155.0, 220.0, 200.0, 75.0, 10e4, 10e4, 10e4, 10e4], + ] +) +regr_pars_6lev = np.array( + [ + [130.0, 165.0, 120.0, 55.0, 50.0, 15.0], + [155.0, 220.0, 200.0, 75.0, 10e4, 10e4], + ] +) +regr_pars_9lev = np.array( + [ + [130.0, 165.0, 120.0, 55.0, 50.0, 15.0, 15.0, 10.0, 10.0], + [155.0, 220.0, 200.0, 75.0, 10e4, 10e4, 10e4, 10e4, 10e4], + ] +) + +# Set the dummy observation and model values +dummy_2d_array = np.array([[1, 2], [3, 4]]) +obs_8lev = np.repeat(dummy_2d_array[None, :, :], 8, axis=0) +obs_6lev = np.repeat(dummy_2d_array[None, :, :], 6, axis=0) +obs_9lev = np.repeat(dummy_2d_array[None, :, :], 9, axis=0) +mod_8lev = np.repeat(dummy_2d_array[None, :, :], 8, axis=0) +mod_6lev = np.repeat(dummy_2d_array[None, :, :], 6, axis=0) +mod_9lev = np.repeat(dummy_2d_array[None, :, :], 9, axis=0) + +# Gives some dummy values to PHI +dummy_phi = np.array([0.472650, 0.523825, 0.103454]) +PHI_8lev = np.repeat(dummy_phi[None, :], 8, axis=0) +PHI_6lev = np.repeat(dummy_phi[None, :], 6, axis=0) +PHI_9lev = np.repeat(dummy_phi[None, :], 9, axis=0) + +# Test function arguments +skill_scores_arg_names = ( + "obs", + "mod", + "lt", + "PHI", + "cor_prev", + "clim_cor_values", + "regr_pars", + "n_cascade_levels", + "expected_cor_t0", + "expected_cor_nwp_lt", + "expected_cor_nowcast_lt", +) + +# Test function values +skill_scores_arg_values = [ + ( + obs_8lev, + mod_8lev, + 60, + PHI_8lev, + None, + clim_cor_values_8lev, + regr_pars_8lev, + 8, + np.repeat(1.0, 8), + np.array( + [ + 0.97455941, + 0.9356775, + 0.81972779, + 0.55202975, + 0.31534738, + 0.02264599, + 0.02343133, + 0.00647032, + ] + ), + np.array( + [ + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + ] + ), + ), + ( + obs_6lev, + mod_6lev, + 60, + PHI_6lev, + None, + clim_cor_values_6lev, + regr_pars_6lev, + 6, + np.repeat(1.0, 6), + np.array( + [0.97455941, 0.9356775, 0.81972779, 0.55202975, 0.31534738, 0.02264599] + ), + np.array([0.996475, 0.996475, 0.996475, 0.996475, 0.996475, 0.996475]), + ), + ( + obs_9lev, + mod_9lev, + 60, + PHI_9lev, + None, + clim_cor_values_9lev, + regr_pars_9lev, + 9, + np.repeat(1.0, 9), + np.array( + [ + 0.97455941, + 0.9356775, + 0.81972779, + 0.55202975, + 0.31534738, + 0.02264599, + 0.02343133, + 0.00647032, + 0.00347776, + ] + ), + np.array( + [ + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + ] + ), + ), + ( + obs_8lev, + mod_8lev, + 0, + PHI_8lev, + None, + clim_cor_values_8lev, + regr_pars_8lev, + 8, + np.repeat(1.0, 8), + np.repeat(1.0, 8), + np.array( + [ + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + 0.996475, + ] + ), + ), +] + + +# The test +@pytest.mark.parametrize(skill_scores_arg_names, skill_scores_arg_values) + +# The test function to be used +def test_blending_skill_scores( + obs, + mod, + lt, + PHI, + cor_prev, + clim_cor_values, + regr_pars, + n_cascade_levels, + expected_cor_t0, + expected_cor_nwp_lt, + expected_cor_nowcast_lt, +): + """Tests if the skill_score functions behave correctly. A dummy gridded + model and observation field should be given for n_cascade_levels, which + leads to a given spatial correlation per cascade level. Then, the function + tests if the correlation regresses towards the climatological values given + lead time lt for the NWP fields or given the PHI-values for the + extrapolation field. + + """ + # Calculate the spatial correlation of the given model field + correlations_t0 = np.array(spatial_correlation(obs, mod)) + + # Check if the field has the same number of cascade levels as the model + # field and as the given n_cascade_levels + assert ( + correlations_t0.shape[0] == mod.shape[0] + ), "Number of cascade levels should be the same as in the model field" + assert ( + correlations_t0.shape[0] == n_cascade_levels + ), "Number of cascade levels should be the same as n_cascade_levels" + + # Check if the returned values are as expected + assert_array_almost_equal( + correlations_t0, + expected_cor_t0, + decimal=3, + err_msg="Returned spatial correlation is not the same as the expected value", + ) + + # Test if the NWP correlation regresses towards the correct value given + # a lead time in minutes + # First, check if the climatological values are returned correctly + correlations_clim, regr_clim = clim_regr_values(n_cascade_levels=n_cascade_levels) + assert ( + correlations_clim.shape[0] == n_cascade_levels + ), "Number of cascade levels should be the same as n_cascade_levels" + assert_array_almost_equal( + correlations_clim, + clim_cor_values, + decimal=3, + err_msg="Not the correct climatological correlations were returned", + ) + assert_array_almost_equal( + regr_clim, + regr_pars, + decimal=3, + err_msg="Not the correct regression parameters were returned", + ) + + # Then, check the regression of the correlation values + correlations_nwp_lt = lt_dependent_cor_nwp(lt=lt, correlations=correlations_t0) + assert ( + correlations_nwp_lt.shape[0] == mod.shape[0] + ), "Number of cascade levels should be the same as in the model field" + assert ( + correlations_nwp_lt.shape[0] == n_cascade_levels + ), "Number of cascade levels should be the same as n_cascade_levels" + assert_array_almost_equal( + correlations_nwp_lt, + expected_cor_nwp_lt, + decimal=3, + err_msg="Correlations of NWP not equal to the expected correlations", + ) + + # Finally, make sure nowcast correlation regresses towards the correct + # value given some PHI-values. + correlations_nowcast_lt, __ = lt_dependent_cor_extrapolation( + PHI, correlations_t0, cor_prev + ) + + print(correlations_nowcast_lt) + assert ( + correlations_nowcast_lt.shape[0] == mod.shape[0] + ), "Number of cascade levels should be the same as in the model field" + assert ( + correlations_nowcast_lt.shape[0] == n_cascade_levels + ), "Number of cascade levels should be the same as n_cascade_levels" + assert_array_almost_equal( + correlations_nowcast_lt, + expected_cor_nowcast_lt, + decimal=3, + err_msg="Correlations of nowcast not equal to the expected correlations", + ) diff --git a/pysteps/tests/test_cascade.py b/pysteps/tests/test_cascade.py index 4a428be99..77951922c 100644 --- a/pysteps/tests/test_cascade.py +++ b/pysteps/tests/test_cascade.py @@ -13,6 +13,9 @@ from pysteps.cascade.decomposition import decomposition_fft, recompose_fft from pysteps.tests.helpers import smart_assert +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip("Needs migration to xarray") + def test_decompose_recompose(): """Tests cascade decomposition.""" diff --git a/pysteps/tests/test_datasets.py b/pysteps/tests/test_datasets.py index 3c6f9581b..e524099e7 100644 --- a/pysteps/tests/test_datasets.py +++ b/pysteps/tests/test_datasets.py @@ -13,6 +13,9 @@ ) from pysteps.exceptions import DirectoryNotEmpty +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip("Needs migration to xarray") + _datasets_opt_deps = dict( fmi=["pyproj"], mch=["PIL"], diff --git a/pysteps/tests/test_detcatscores.py b/pysteps/tests/test_detcatscores.py index 2126df6bb..e1dc3de4b 100644 --- a/pysteps/tests/test_detcatscores.py +++ b/pysteps/tests/test_detcatscores.py @@ -6,6 +6,9 @@ from pysteps.verification import det_cat_fct +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip("Needs migration to xarray") + # CREATE A LARGE DATASET TO MATCH # EXAMPLES IN # http://www.cawcr.gov.au/projects/verification/ diff --git a/pysteps/tests/test_exporters.py b/pysteps/tests/test_exporters.py index 95776df23..9515828b7 100644 --- a/pysteps/tests/test_exporters.py +++ b/pysteps/tests/test_exporters.py @@ -15,6 +15,11 @@ from pysteps.io.exporters import initialize_forecast_exporter_netcdf from pysteps.tests.helpers import get_precipitation_fields, get_invalid_mask +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip( + "Needs migration to xarray. Problem with cartesian_unit attribute." +) + def test_get_geotiff_filename(): """Test the geotif name generator.""" diff --git a/pysteps/tests/test_io_archive.py b/pysteps/tests/test_io_archive.py new file mode 100644 index 000000000..74a0567ec --- /dev/null +++ b/pysteps/tests/test_io_archive.py @@ -0,0 +1,35 @@ +from datetime import datetime + +import pytest + +import pysteps + + +def test_find_by_date_mch(): + + pytest.importorskip("PIL") + + date = datetime.strptime("201505151630", "%Y%m%d%H%M") + data_source = pysteps.rcparams.data_sources["mch"] + root_path = data_source["root_path"] + path_fmt = data_source["path_fmt"] + fn_pattern = data_source["fn_pattern"] + fn_ext = data_source["fn_ext"] + timestep = data_source["timestep"] + + fns = pysteps.io.archive.find_by_date( + date, + root_path, + path_fmt, + fn_pattern, + fn_ext, + timestep=timestep, + num_prev_files=1, + num_next_files=1, + ) + + assert len(fns) == 2 + assert len(fns[0]) == 3 + assert len(fns[1]) == 3 + assert isinstance(fns[0][0], str) + assert isinstance(fns[1][0], datetime) diff --git a/pysteps/tests/test_io_bom_nwp.py b/pysteps/tests/test_io_bom_nwp.py new file mode 100644 index 000000000..e2d5e5d1d --- /dev/null +++ b/pysteps/tests/test_io_bom_nwp.py @@ -0,0 +1,82 @@ +# -*- coding: utf-8 -*- + +import os +import numpy as np + +import xarray as xr +import pytest + +import pysteps +from pysteps.tests.helpers import smart_assert + +pytest.importorskip("netCDF4") + +root_path = pysteps.rcparams.data_sources["bom_nwp"]["root_path"] +rel_path = os.path.join("2020", "10", "31") +filename = os.path.join(root_path, rel_path, "20201031_0000_regrid_short.nc") +data_array_xr = pysteps.io.import_bom_nwp_xr(filename) + +expected_proj = ( + "+proj=aea +lon_0=153.240 +lat_0=-27.718 +lat_1=-26.200 +lat_2=-29.300" +) + + +def test_io_import_bom_nwp_xarray(): + """Test the importer Bom NWP.""" + assert isinstance(data_array_xr, xr.DataArray) + + +def test_io_import_bom_nwp_xarray_shape(): + """Test the importer Bom NWP shape.""" + assert isinstance(data_array_xr, xr.DataArray) + assert data_array_xr.shape == (144, 512, 512) + + +test_attrs_xr = [ + ("projection", expected_proj, None), + ("institution", "Commonwealth of Australia, Bureau of Meteorology", None), + ("transform", None, None), + ("zerovalue", 0.0, 0.1), + ("unit", "mm", None), + ("accutime", np.timedelta64(10, 'm'), None), + ("zr_a", None, None), + ("zr_b", None, None), + ("xpixelsize", 500.0, 0.1), + ("ypixelsize", 500.0, 0.1), + ("units", "mm", None), + ("yorigin", "upper", None), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs_xr) +def test_io_import_bom_nwp_xarray_attrs(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_x = [ + ("units", "m", None), + ("x1", -127750.0, 0.1), + ("x2", 127750.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_x) +def test_io_import_bom_nwp_xarray_attrs_coordx(variable, expected, tolerance): + """Test the importer Bom NWP.""" + smart_assert(data_array_xr.x.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_y = [ + ("units", "m", None), + ("y1", -127750.0, 0.1), + ("y2", 127750.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_y) +def test_io_import_bom_nwp_xarray_attrs_coordy(variable, expected, tolerance): + """Test the importer Bom NWP.""" + smart_assert(data_array_xr.y.attrs[variable], expected, tolerance) diff --git a/pysteps/tests/test_io_bom_rf3.py b/pysteps/tests/test_io_bom_rf3.py index 6c2f17062..b276c3a61 100644 --- a/pysteps/tests/test_io_bom_rf3.py +++ b/pysteps/tests/test_io_bom_rf3.py @@ -1,7 +1,9 @@ # -*- coding: utf-8 -*- import os +import numpy as np +import xarray as xr import pytest import pysteps @@ -9,55 +11,41 @@ pytest.importorskip("netCDF4") +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_100000.prcp-cscn.nc") +data_array = pysteps.io.import_bom_rf3(filename) + def test_io_import_bom_rf3_shape(): """Test the importer Bom RF3.""" - 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_100000.prcp-cscn.nc") - precip, _, _ = pysteps.io.import_bom_rf3(filename) - assert precip.shape == (512, 512) + assert isinstance(data_array, xr.DataArray) + assert data_array.shape == (512, 512) -expected_proj1 = ( +expected_proj = ( "+proj=aea +lon_0=144.752 +lat_0=-37.852 " "+lat_1=-18.000 +lat_2=-36.000" ) -# test_metadata_bom: list of (variable,expected,tolerance) tuples -test_metadata_bom = [ +test_attrs = [ + ("projection", expected_proj, None), + ("institution", "Commonwealth of Australia, Bureau of Meteorology", None), ("transform", None, None), ("zerovalue", 0.0, 0.1), - ("projection", expected_proj1, None), ("unit", "mm", None), ("accutime", 6, 0.1), - ("x1", -128000.0, 0.1), - ("x2", 127500.0, 0.1), - ("y1", -127500.0, 0.1), - ("y2", 128000.0, 0.1), - ("xpixelsize", 500.0, 0.1), - ("ypixelsize", 500.0, 0.1), - ("cartesian_unit", "m", None), - ("yorigin", "upper", None), - ("institution", "Commonwealth of Australia, Bureau of Meteorology", None), ] -@pytest.mark.parametrize("variable, expected, tolerance", test_metadata_bom) -def test_io_import_bom_rf3_metadata(variable, expected, tolerance): +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs) +def test_io_import_bom_rf3_dataset_attrs(variable, expected, tolerance): """Test the importer Bom RF3.""" - 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_100000.prcp-cscn.nc") - _, _, metadata = pysteps.io.import_bom_rf3(filename) - smart_assert(metadata[variable], expected, tolerance) + smart_assert(data_array.attrs[variable], expected, tolerance) -expected_proj2 = ( - "+proj=aea +lon_0=144.752 +lat_0=-37.852 " "+lat_1=-18.000 +lat_2=-36.000" -) # test_geodata: list of (variable,expected,tolerance) tuples test_geodata_bom = [ - ("projection", expected_proj2, None), + ("projection", expected_proj, None), ("unit", "mm", None), ("accutime", 6, 0.1), ("x1", -128000.0, 0.1), @@ -77,6 +65,72 @@ def test_io_import_bom_rf3_geodata(variable, expected, tolerance): """Test the importer Bom RF3.""" 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_100000.prcp-cscn.nc") + filename = os.path.join(root_path, rel_path, + "2_20180616_100000.prcp-cscn.nc") geodata = pysteps.io.importers._import_bom_rf3_geodata(filename) smart_assert(geodata[variable], expected, tolerance) + + +# TEST XARRAY IMPLEMENTATION +data_array_xr = pysteps.io.import_bom_rf3_xr(filename) + + +def test_io_import_bom_rf3_xarray(): + """Test the importer Bom RF3.""" + assert isinstance(data_array_xr, xr.DataArray) + + +def test_io_import_bom_rf3_xarray_shape(): + """Test the importer Bom RF3.""" + assert isinstance(data_array_xr, xr.DataArray) + assert data_array_xr.shape == (1, 512, 512) + + +test_attrs_xr = [ + ("projection", expected_proj, None), + ("institution", "Commonwealth of Australia, Bureau of Meteorology", None), + ("transform", None, None), + ("zerovalue", 0.0, 0.1), + ("unit", "mm", None), + ("accutime", np.timedelta64(6,'m'), None), + ("zr_a", None, None), + ("zr_b", None, None), + ("xpixelsize", 500.0, 0.1), + ("ypixelsize", 500.0, 0.1), + ("units", "mm", None), + ("yorigin", "upper", None), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs_xr) +def test_io_import_bom_rf3_xarray_attrs(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_x = [ + ("units", "m", None), + ("x1", -128000.0, 0.1), + ("x2", 127500.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_x) +def test_io_import_bom_rf3_xarray_attrs_coordx(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.x.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_y = [ + ("units", "m", None), + ("y1", -127500.0, 0.1), + ("y2", 128000.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_y) +def test_io_import_bom_rf3_xarray_attrs_coordy(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.y.attrs[variable], expected, tolerance) diff --git a/pysteps/tests/test_io_fmi_pgm.py b/pysteps/tests/test_io_fmi_pgm.py index 58dfd5137..a07ebb9ee 100644 --- a/pysteps/tests/test_io_fmi_pgm.py +++ b/pysteps/tests/test_io_fmi_pgm.py @@ -2,6 +2,7 @@ import os +import xarray as xr import pytest import pysteps @@ -10,62 +11,66 @@ pytest.importorskip("pyproj") -def test_io_import_fmi_pgm_shape(): - """Test the importer FMI PGM.""" - root_path = pysteps.rcparams.data_sources["fmi"]["root_path"] - filename = os.path.join( - root_path, - "20160928", - "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz", - ) - R, _, _ = pysteps.io.import_fmi_pgm(filename, gzipped=True) - assert R.shape == (1226, 760) - - -# test_metadata: list of (variable,expected) tuples -test_metadata = [ - ("composite_area", ["FIN"]), - ("projection_name", ["SUOMI1"]), - ("radar", ["LUO", "1", "26.9008", "67.1386"]), - ("obstime", ["201609281600"]), - ("producttype", ["CAPPI"]), - ("productname", ["LOWEST"]), - ("param", ["CorrectedReflectivity"]), - ("metersperpixel_x", ["999.674053"]), - ("metersperpixel_y", ["999.62859"]), - ("projection", ["radar", "{"]), - ("type", ["stereographic"]), - ("centrallongitude", ["25"]), - ("centrallatitude", ["90"]), - ("truelatitude", ["60"]), - ("bottomleft", ["18.600000", "57.930000"]), - ("topright", ["34.903000", "69.005000"]), - ("missingval", 255), -] +root_path = pysteps.rcparams.data_sources["fmi"]["root_path"] +filename = os.path.join( + root_path, + "20160928", + "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz", +) +data_array = pysteps.io.import_fmi_pgm(filename, gzipped=True) -@pytest.mark.parametrize("variable,expected", test_metadata) -def test_io_import_fmi_pmg_metadata(variable, expected): - """Test the importer FMI PMG.""" - root_path = pysteps.rcparams.data_sources["fmi"]["root_path"] - filename = os.path.join( - root_path, - "20160928", - "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz", - ) - metadata = pysteps.io.importers._import_fmi_pgm_metadata(filename, gzipped=True) - assert metadata[variable] == expected +def test_io_import_fmi_pgm_shape(): + """Test the importer FMI PGM.""" + assert isinstance(data_array, xr.DataArray) + assert data_array.shape == (1226, 760) -expected_proj1 = ( +expected_proj = ( "+proj=stere +lon_0=25E +lat_0=90N " "+lat_ts=60 +a=6371288 +x_0=380886.310 " "+y_0=3395677.920 +no_defs" ) +test_attrs = [ + ("projection", expected_proj, None), + ("institution", "Finnish Meteorological Institute", None), + # ("composite_area", ["FIN"]), + # ("projection_name", ["SUOMI1"]), + # ("radar", ["LUO", "1", "26.9008", "67.1386"]), + # ("obstime", ["201609281600"]), + # ("producttype", ["CAPPI"]), + # ("productname", ["LOWEST"]), + # ("param", ["CorrectedReflectivity"]), + # ("metersperpixel_x", ["999.674053"]), + # ("metersperpixel_y", ["999.62859"]), + # ("projection", ["radar", "{"]), + # ("type", ["stereographic"]), + # ("centrallongitude", ["25"]), + # ("centrallatitude", ["90"]), + # ("truelatitude", ["60"]), + # ("bottomleft", ["18.600000", "57.930000"]), + # ("topright", ["34.903000", "69.005000"]), + # ("missingval", 255), + ("accutime", 5.0, 0.1), + ("unit", "dBZ", None), + ("transform", "dB", None), + ("zerovalue", -32.0, 0.1), + ("threshold", -31.5, 0.1), + ("zr_a", 223.0, 0.1), + ("zr_b", 1.53, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs) +def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance): + """Test the importer FMI PMG.""" + smart_assert(data_array.attrs[variable], expected, tolerance) + + # test_geodata: list of (variable,expected,tolerance) tuples test_geodata = [ - ("projection", expected_proj1, None), + ("projection", expected_proj, None), ("x1", 0.0049823258887045085, 1e-20), ("x2", 759752.2852757066, 1e-10), ("y1", 0.009731985162943602, 1e-18), diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py index 44e777a11..657e522bf 100644 --- a/pysteps/tests/test_io_knmi_hdf5.py +++ b/pysteps/tests/test_io_knmi_hdf5.py @@ -2,6 +2,7 @@ import os +import xarray as xr import pytest import pysteps @@ -10,54 +11,40 @@ pytest.importorskip("h5py") -def test_io_import_knmi_hdf5_shape(): - """Test the importer KNMI HDF5.""" - - root_path = pysteps.rcparams.data_sources["knmi"]["root_path"] +root_path = pysteps.rcparams.data_sources["knmi"]["root_path"] +filename = os.path.join(root_path, "2010/08", "RAD_NL25_RAP_5min_201008260000.h5") +data_array = pysteps.io.import_knmi_hdf5(filename) - filename = os.path.join(root_path, "2010/08", "RAD_NL25_RAP_5min_201008260000.h5") - precip, __, __ = pysteps.io.import_knmi_hdf5(filename) - - assert precip.shape == (765, 700) +def test_io_import_knmi_hdf5_shape(): + """Test the importer KNMI HDF5.""" + assert isinstance(data_array, xr.DataArray) + assert data_array.shape == (765, 700) # test_metadata: list of (variable,expected, tolerance) tuples -expected_proj1 = ( +expected_proj = ( "+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" ) -test_metadata = [ - ("projection", expected_proj1, None), - ("x1", 0.0, 1e-10), - ("y1", -4415.038179210632, 1e-10), - ("x2", 699.9842646331593, 1e-10), - ("y2", -3649.950360247753, 1e-10), - ("xpixelsize", 1.0, 1e-10), - ("xpixelsize", 1.0, 1e-10), - ("cartesian_unit", "km", None), +# list of (variable,expected,tolerance) tuples +test_attrs = [ + ("projection", expected_proj, None), + ("institution", "KNMI - Royal Netherlands Meteorological Institute", None), ("accutime", 5.0, 1e-10), - ("yorigin", "upper", None), ("unit", "mm", None), - ("institution", "KNMI - Royal Netherlands Meteorological Institute", None), ("transform", None, None), ("zerovalue", 0.0, 1e-10), ("threshold", 0.01, 1e-10), - ("zr_a", 200.0, None), - ("zr_b", 1.6, None), + ("zr_a", 200.0, 0.1), + ("zr_b", 1.6, 0.1), ] -@pytest.mark.parametrize("variable,expected,tolerance", test_metadata) +@pytest.mark.parametrize("variable,expected,tolerance", test_attrs) def test_io_import_knmi_hdf5_metadata(variable, expected, tolerance): """Test the importer KNMI HDF5.""" - root_path = pysteps.rcparams.data_sources["knmi"]["root_path"] - - filename = os.path.join(root_path, "2010/08", "RAD_NL25_RAP_5min_201008260000.h5") - - __, __, metadata = pysteps.io.import_knmi_hdf5(filename) - - smart_assert(metadata[variable], expected, tolerance) + smart_assert(data_array.attrs[variable], expected, tolerance) diff --git a/pysteps/tests/test_io_mch_gif.py b/pysteps/tests/test_io_mch_gif.py index f569dd4a2..d0eeae2b1 100644 --- a/pysteps/tests/test_io_mch_gif.py +++ b/pysteps/tests/test_io_mch_gif.py @@ -3,22 +3,25 @@ import os import pytest +import xarray as xr import pysteps from pysteps.tests.helpers import smart_assert pytest.importorskip("PIL") +root_path = pysteps.rcparams.data_sources["mch"]["root_path"] +filename = os.path.join(root_path, "20170131", "AQC170310945F_00005.801.gif") +data_array = pysteps.io.import_mch_gif(filename, "AQC", "mm", 5.0) + def test_io_import_mch_gif_shape(): """Test the importer MCH GIF.""" - root_path = pysteps.rcparams.data_sources["mch"]["root_path"] - filename = os.path.join(root_path, "20170131", "AQC170310945F_00005.801.gif") - precip, _, metadata = pysteps.io.import_mch_gif(filename, "AQC", "mm", 5.0) - assert precip.shape == (640, 710) + assert isinstance(data_array, xr.DataArray) + assert data_array.shape == (640, 710) -expected_proj1 = ( +expected_proj = ( "+proj=somerc +lon_0=7.43958333333333 " "+lat_0=46.9524055555556 +k_0=1 " "+x_0=600000 +y_0=200000 +ellps=bessel " @@ -26,48 +29,29 @@ def test_io_import_mch_gif_shape(): "+units=m +no_defs" ) -# test_metadata: list of (variable,expected,tolerance) tuples -test_metadata = [ - ("projection", expected_proj1, None), - ("x1", 255000.0, 0.1), - ("y1", -160000.0, 0.1), - ("x2", 965000.0, 0.1), - ("y2", 480000.0, 0.1), - ("xpixelsize", 1000.0, 0.1), - ("ypixelsize", 1000.0, 0.1), - ("cartesian_unit", "m", None), - ("yorigin", "upper", None), +# list of (variable,expected,tolerance) tuples +test_attrs = [ + ("projection", expected_proj, None), + ("institution", "MeteoSwiss", None), ("accutime", 5.0, 0.1), ("unit", "mm", None), ("transform", None, None), - ("zerovalue", 0.0, None), + ("zerovalue", 0.0, 0.1), ("threshold", 0.0009628129986471908, 1e-19), - ("institution", "MeteoSwiss", None), - ("product", "AQC", None), - ("zr_a", 316.0, None), - ("zr_b", 1.5, None), + ("zr_a", 316.0, 0.1), + ("zr_b", 1.5, 0.1), ] -@pytest.mark.parametrize("variable, expected, tolerance", test_metadata) -def test_io_import_mch_gif_metadata(variable, expected, tolerance): +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs) +def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance): """Test the importer MCH GIF.""" - root_path = pysteps.rcparams.data_sources["mch"]["root_path"] - filename = os.path.join(root_path, "20170131", "AQC170310945F_00005.801.gif") - _, _, metadata = pysteps.io.import_mch_gif(filename, "AQC", "mm", 5.0) - smart_assert(metadata[variable], expected, tolerance) - + smart_assert(data_array.attrs[variable], expected, tolerance) -expected_proj2 = ( - "+proj=somerc +lon_0=7.43958333333333 +lat_0=46.9524055555556 " - "+k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel " - "+towgs84=674.374,15.056,405.346,0,0,0,0 " - "+units=m +no_defs" -) # test_geodata: list of (variable,expected,tolerance) tuples test_geodata = [ - ("projection", expected_proj2, None), + ("projection", expected_proj, None), ("x1", 255000.0, 0.1), ("y1", -160000.0, 0.1), ("x2", 965000.0, 0.1), diff --git a/pysteps/tests/test_io_mrms_grib.py b/pysteps/tests/test_io_mrms_grib.py index 65dd0caeb..afebb7b37 100644 --- a/pysteps/tests/test_io_mrms_grib.py +++ b/pysteps/tests/test_io_mrms_grib.py @@ -4,6 +4,7 @@ import numpy as np import pytest +import xarray as xr from numpy.testing import assert_array_almost_equal import pysteps @@ -15,33 +16,32 @@ def test_io_import_mrms_grib(): """Test the importer for NSSL data.""" root_path = pysteps.rcparams.data_sources["mrms"]["root_path"] - filename = os.path.join( root_path, "2019/06/10/", "PrecipRate_00.00_20190610-000000.grib2" ) + data_Array = pysteps.io.import_mrms_grib(filename, fillna=0, window_size=1) - precip_full, _, metadata = pysteps.io.import_mrms_grib( - filename, fillna=0, window_size=1 - ) - assert precip_full.shape == (3500, 7000) - assert precip_full.dtype == "single" + assert isinstance(data_Array, xr.DataArray) + assert data_Array.shape == (3500, 7000) + assert data_Array.dtype == "single" expected_metadata = { - "xpixelsize": 0.01, - "ypixelsize": 0.01, + "institution": "NOAA National Severe Storms Laboratory", + "projection": "+proj=longlat +ellps=IAU76", "unit": "mm/h", + "accutime": 2.0, "transform": None, "zerovalue": 0, - "projection": "+proj=longlat +ellps=IAU76", - "yorigin": "upper", "threshold": 0.1, "x1": -129.99999999999997, "x2": -60.00000199999991, "y1": 20.000001, "y2": 55.00000000000001, - "cartesian_unit": "degrees", + "units": "degrees", } - + metadata = data_Array.attrs + metadata.update(**data_Array.x.attrs) + metadata.update(**data_Array.y.attrs) for key, value in expected_metadata.items(): if isinstance(value, float): assert_array_almost_equal(metadata[key], expected_metadata[key]) @@ -50,6 +50,7 @@ def test_io_import_mrms_grib(): x = np.arange(metadata["x1"], metadata["x2"], metadata["xpixelsize"]) y = np.arange(metadata["y1"], metadata["y2"], metadata["ypixelsize"]) + precip_full = data_Array.values assert y.size == precip_full.shape[0] assert x.size == precip_full.shape[1] @@ -59,36 +60,38 @@ def test_io_import_mrms_grib(): # Test that if the bounding box is larger than the domain, all the points are returned. precip_full2 = pysteps.io.import_mrms_grib( filename, fillna=0, extent=(220, 300, 20, 55), window_size=1 - )[0] + ) assert precip_full2.shape == (3500, 7000) - assert_array_almost_equal(precip_full, precip_full2) + assert_array_almost_equal(data_Array, precip_full2) del precip_full2 # Test that a portion of the domain is returned correctly precip_clipped = pysteps.io.import_mrms_grib( filename, fillna=0, extent=(250, 260, 30, 35), window_size=1 - )[0] + ) assert precip_clipped.shape == (500, 1000) - assert_array_almost_equal(precip_clipped, precip_full[2000:2500, 2000:3000]) + assert_array_almost_equal( + precip_clipped, data_Array.isel(y=slice(2000, 2500), x=slice(2000, 3000)) + ) del precip_clipped - precip_single = pysteps.io.import_mrms_grib(filename, dtype="double", fillna=0)[0] + precip_single = pysteps.io.import_mrms_grib(filename, dtype="double", fillna=0) assert precip_single.dtype == "double" del precip_single - precip_single = pysteps.io.import_mrms_grib(filename, dtype="single", fillna=0)[0] + precip_single = pysteps.io.import_mrms_grib(filename, dtype="single", fillna=0) assert precip_single.dtype == "single" del precip_single precip_donwscaled = pysteps.io.import_mrms_grib( filename, dtype="single", fillna=0, window_size=2 - )[0] + ) assert precip_donwscaled.shape == (3500 / 2, 7000 / 2) - precip_donwscaled, _, metadata = pysteps.io.import_mrms_grib( + precip_donwscaled = pysteps.io.import_mrms_grib( filename, dtype="single", fillna=0, window_size=3 ) print(metadata) @@ -108,6 +111,22 @@ def test_io_import_mrms_grib(): # precipitation pattern (i.e. the file being read). expected_metadata.pop("threshold", None) - for key in expected_metadata.keys(): - assert metadata[key] == expected_metadata[key] - assert precip_donwscaled.shape == (3500 // 3, 7000 // 3) + # expected_metadata.update( + # xpixelsize=0.03, + # ypixelsize=0.03, + # x1=-129.98500000028577, + # x2=-60.02500199942843, + # y1=20.03500099914261, + # y2=54.985000000285794, + # ) + + # # Remove the threshold keyword from the test when the window_size>1 is used. + # # The threshold is computed automatically from the minimum valid precipitation values. + # # The minimum value is affected by the the block_averaging (window_size=3 keyword) + # # of the precipitation fields. Hence, the "threshold" value will depend on the + # # precipitation pattern (i.e. the file being read). + # expected_metadata.pop("threshold", None) + + # for key in expected_metadata.keys(): + # assert metadata[key] == expected_metadata[key] + # assert precip_donwscaled.shape == (3500 // 3, 7000 // 3) diff --git a/pysteps/tests/test_io_opera_hdf5.py b/pysteps/tests/test_io_opera_hdf5.py index 6411bcd33..71032bbd3 100644 --- a/pysteps/tests/test_io_opera_hdf5.py +++ b/pysteps/tests/test_io_opera_hdf5.py @@ -2,6 +2,7 @@ import os +import xarray as xr import pytest import pysteps @@ -10,53 +11,39 @@ pytest.importorskip("h5py") +root_path = pysteps.rcparams.data_sources["opera"]["root_path"] +filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf") +data_array = pysteps.io.import_opera_hdf5(filename) + + def test_io_import_opera_hdf5_shape(): """Test the importer OPERA HDF5.""" - - root_path = pysteps.rcparams.data_sources["opera"]["root_path"] - filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf") - R, _, _ = pysteps.io.import_opera_hdf5(filename) - assert R.shape == (2200, 1900) + assert isinstance(data_array, xr.DataArray) + assert data_array.shape == (2200, 1900) # test_metadata: list of (variable,expected, tolerance) tuples -expected_proj1 = ( +expected_proj = ( "+proj=laea +lat_0=55.0 +lon_0=10.0 " "+x_0=1950000.0 " "+y_0=-2100000.0 " "+units=m +ellps=WGS84" ) -test_metadata = [ - ("projection", expected_proj1, None), - ("ll_lon", -10.434576838640398, 1e-10), - ("ll_lat", 31.746215319325056, 1e-10), - ("ur_lon", 57.81196475014995, 1e-10), - ("ur_lat", 67.62103710275053, 1e-10), - ("x1", -0.0004161088727414608, 1e-6), - ("y1", -4400000.001057557, 1e-10), - ("x2", 3800000.0004256153, 1e-10), - ("y2", -0.0004262728616595268, 1e-6), - ("xpixelsize", 2000.0, 1e-10), - ("xpixelsize", 2000.0, 1e-10), - ("cartesian_unit", "m", None), +# list of (variable,expected,tolerance) tuples +test_attrs = [ + ("projection", expected_proj, None), + ("institution", "Odyssey datacentre", None), ("accutime", 15.0, 1e-10), - ("yorigin", "upper", None), ("unit", "mm/h", None), - ("institution", "Odyssey datacentre", None), ("transform", None, None), ("zerovalue", 0.0, 1e-10), ("threshold", 0.01, 1e-10), ] -@pytest.mark.parametrize("variable,expected,tolerance", test_metadata) -def test_io_import_opera_hdf5_metadata(variable, expected, tolerance): +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs) +def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance): """Test the importer OPERA HDF5.""" - root_path = pysteps.rcparams.data_sources["opera"]["root_path"] - filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf") - _, _, metadata = pysteps.io.import_opera_hdf5(filename) - print(metadata) - # assert metadata[variable] == expected - smart_assert(metadata[variable], expected, tolerance) + smart_assert(data_array.attrs[variable], expected, tolerance) diff --git a/pysteps/tests/test_io_readers.py b/pysteps/tests/test_io_readers.py new file mode 100644 index 000000000..4feaf87c9 --- /dev/null +++ b/pysteps/tests/test_io_readers.py @@ -0,0 +1,38 @@ +from datetime import datetime + +import pytest +import xarray as xr + +import pysteps + + +def test_read_timeseries_mch(): + + pytest.importorskip("PIL") + + date = datetime.strptime("201505151630", "%Y%m%d%H%M") + data_source = pysteps.rcparams.data_sources["mch"] + root_path = data_source["root_path"] + path_fmt = data_source["path_fmt"] + fn_pattern = data_source["fn_pattern"] + fn_ext = data_source["fn_ext"] + importer_name = data_source["importer"] + importer_kwargs = data_source["importer_kwargs"] + timestep = data_source["timestep"] + + fns = pysteps.io.archive.find_by_date( + date, + root_path, + path_fmt, + fn_pattern, + fn_ext, + timestep=timestep, + num_prev_files=1, + num_next_files=1, + ) + + importer = pysteps.io.get_method(importer_name, "importer") + data_array = pysteps.io.read_timeseries(fns, importer, **importer_kwargs) + + assert isinstance(data_array, xr.DataArray) + assert data_array.t.size == 3 diff --git a/pysteps/tests/test_io_saf_crri.py b/pysteps/tests/test_io_saf_crri.py index 66faa29d2..6e4e5dce9 100644 --- a/pysteps/tests/test_io_saf_crri.py +++ b/pysteps/tests/test_io_saf_crri.py @@ -2,6 +2,8 @@ import os +import numpy as np +import xarray as xr import pytest import pysteps @@ -36,30 +38,31 @@ def test_io_import_saf_crri_geodata(variable, expected, tolerance): root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc" ) geodata = pysteps.io.importers._import_saf_crri_geodata(filename) - print(variable, geodata[variable], expected, tolerance) smart_assert(geodata[variable], expected, tolerance) -test_metadata_crri = [ +root_path = pysteps.rcparams.data_sources["saf"]["root_path"] +rel_path = "20180601/CRR" +filename = os.path.join( + root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc" +) +data_array = pysteps.io.import_saf_crri(filename) + +# list of (variable,expected,tolerance) tuples +test_attrs = [ + ("projection", expected_proj, None), + ("institution", "Agencia Estatal de Meteorología (AEMET)", None), ("transform", None, None), ("zerovalue", 0.0, 0.1), ("unit", "mm/h", None), ("accutime", None, None), - ("institution", "Agencia Estatal de Meteorología (AEMET)", None), ] -@pytest.mark.parametrize("variable, expected, tolerance", test_metadata_crri) -def test_io_import_saf_crri_metadata(variable, expected, tolerance): +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs) +def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance): """Test the importer SAF CRRI.""" - root_path = pysteps.rcparams.data_sources["saf"]["root_path"] - rel_path = "20180601/CRR" - filename = os.path.join( - root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc" - ) - _, _, metadata = pysteps.io.import_saf_crri(filename) - print(variable, metadata[variable], expected, tolerance) - smart_assert(metadata[variable], expected, tolerance) + smart_assert(data_array.attrs[variable], expected, tolerance) test_extent_crri = [ @@ -83,22 +86,17 @@ def test_io_import_saf_crri_extent(extent, expected_extent, expected_shape, tole filename = os.path.join( root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc" ) - precip, _, metadata = pysteps.io.import_saf_crri(filename, extent=extent) - extent_out = (metadata["x1"], metadata["x2"], metadata["y1"], metadata["y2"]) - - print(extent, extent_out, expected_extent) + data_array = pysteps.io.import_saf_crri(filename, extent=extent) + + xgridsize = np.diff(data_array.x)[0] + ygridsize = np.diff(data_array.y)[0] + extent_out = ( + data_array.x.min() - xgridsize / 2, + data_array.x.max() + xgridsize / 2, + data_array.y.min() - ygridsize / 2, + data_array.y.max() + ygridsize / 2, + ) + print([int(out) for out in extent_out]) smart_assert(extent_out, expected_extent, tolerance) - print(precip.shape, expected_shape) - smart_assert(precip.shape, expected_shape, tolerance) - - -if __name__ == "__main__": - - for i, args in enumerate(test_geodata_crri): - test_io_import_saf_crri_geodata(args[0], args[1], args[2]) - - for i, args in enumerate(test_metadata_crri): - test_io_import_saf_crri_metadata(args[0], args[1], args[2]) - for i, args in enumerate(test_extent_crri): - test_io_import_saf_crri_extent(args[0], args[1], args[2], args[3]) + smart_assert(data_array.shape, expected_shape, tolerance) diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py index ec01998db..f432f1d49 100644 --- a/pysteps/tests/test_nowcasts_steps.py +++ b/pysteps/tests/test_nowcasts_steps.py @@ -9,6 +9,11 @@ from pysteps import io, motion, nowcasts, verification from pysteps.tests.helpers import get_precipitation_fields +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip( + "Needs migration to xarray. Problem with cartesian_unit attribute." +) + steps_arg_names = ( "n_ens_members", "n_cascade_levels", diff --git a/pysteps/tests/test_plt_animate.py b/pysteps/tests/test_plt_animate.py index 80cc45f8f..1021080c5 100644 --- a/pysteps/tests/test_plt_animate.py +++ b/pysteps/tests/test_plt_animate.py @@ -8,6 +8,9 @@ from pysteps.tests.helpers import get_precipitation_fields from pysteps.visualization.animations import animate +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip("Needs migration to xarray") + def test_animate(tmp_path): diff --git a/pysteps/tests/test_plugins_support.py b/pysteps/tests/test_plugins_support.py index 91f9f8de5..d64446178 100644 --- a/pysteps/tests/test_plugins_support.py +++ b/pysteps/tests/test_plugins_support.py @@ -10,6 +10,9 @@ import sys import tempfile +# TODO: Fix tests for xarray fields +pytestmark = pytest.mark.skip("Needs migration to xarray") + __ = pytest.importorskip("cookiecutter") from cookiecutter.main import cookiecutter diff --git a/pysteps/tests/test_probscores.py b/pysteps/tests/test_probscores.py index c7f9990b8..b1af325a7 100644 --- a/pysteps/tests/test_probscores.py +++ b/pysteps/tests/test_probscores.py @@ -8,7 +8,8 @@ from pysteps.tests.helpers import get_precipitation_fields from pysteps.verification import probscores -precip = get_precipitation_fields(num_next_files=10, return_raw=True) +# TODO: Fix tests for xarray fields +precip = get_precipitation_fields(num_next_files=10, return_raw=True, legacy=True) # CRPS test_data = [(precip[:10], precip[-1], 0.01470871)] diff --git a/pysteps/tests/test_spatialscores.py b/pysteps/tests/test_spatialscores.py index ac6a2b63c..a2d6714c4 100644 --- a/pysteps/tests/test_spatialscores.py +++ b/pysteps/tests/test_spatialscores.py @@ -6,7 +6,8 @@ from pysteps.tests.helpers import get_precipitation_fields from pysteps.verification import spatialscores -R = get_precipitation_fields(num_prev_files=1, return_raw=True) +# TODO: Fix tests for xarray fields +R = get_precipitation_fields(num_prev_files=1, return_raw=True, legacy=True) test_data = [ (R[0], R[1], "FSS", [1], [10], None, 0.85161531), (R[0], R[1], "BMSE", [1], None, "Haar", 0.99989651), diff --git a/pysteps/tests/test_timeseries_autoregression.py b/pysteps/tests/test_timeseries_autoregression.py index f1cc76816..ff3326db4 100644 --- a/pysteps/tests/test_timeseries_autoregression.py +++ b/pysteps/tests/test_timeseries_autoregression.py @@ -211,7 +211,7 @@ def _create_data_multivariate(): R = [] for fn in filenames: filename = os.path.join(root_path, "20160928", fn) - R_, _, _ = pysteps.io.import_fmi_pgm(filename, gzipped=True) + R_, _, _ = pysteps.io.import_fmi_pgm(filename, gzipped=True, legacy=True) R_[~np.isfinite(R_)] = 0.0 R.append(np.stack([R_, np.roll(R_, 5, axis=0)])) @@ -235,7 +235,7 @@ def _create_data_univariate(): R = [] for fn in filenames: filename = os.path.join(root_path, "20160928", fn) - R_, _, _ = pysteps.io.import_fmi_pgm(filename, gzipped=True) + R_, _, _ = pysteps.io.import_fmi_pgm(filename, gzipped=True, legacy=True) R_[~np.isfinite(R_)] = 0.0 R.append(R_) diff --git a/requirements.txt b/requirements.txt index 1804df1d9..b8aea1788 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,3 +7,5 @@ matplotlib jsmin jsonschema netCDF4 +xarray +bottleneck diff --git a/requirements_dev.txt b/requirements_dev.txt index bf85f7d86..89ea53ca7 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -9,6 +9,8 @@ matplotlib jsmin jsonschema netCDF4 +xarray +bottleneck # Optional dependencies dask diff --git a/setup.py b/setup.py index 61d0d19ae..5542775c6 100644 --- a/setup.py +++ b/setup.py @@ -66,6 +66,8 @@ "scipy", "matplotlib", "jsonschema", + "xarray", + "bottleneck", ] setup(