diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index 9cfaa17f..cb2cb365 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -82,10 +82,14 @@ import_odim_hdf5 import_opera_hdf5 import_saf_crri + import_dwd_hdf5 + import_dwd_radolan """ import gzip import os +import array +import datetime from functools import partial import numpy as np @@ -191,7 +195,10 @@ def _get_grib_projection(grib_msg): _grib_shapes_of_earth[4] = {"ellps": "GRS80"} _grib_shapes_of_earth[5] = {"ellps": "WGS84"} _grib_shapes_of_earth[6] = {"R": 6371229} - _grib_shapes_of_earth[8] = {"datum": "WGS84", "R": 6371200} + _grib_shapes_of_earth[8] = { + "datum": "WGS84", + "R": 6371200, + } _grib_shapes_of_earth[9] = {"datum": "OSGB36"} # pygrib defines the ellipsoids using "a" and "b" only. @@ -354,7 +361,10 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): lons = block_reduce(lons, window_size[1]) # Update the limits - ul_lat, lr_lat = lats[0], lats[-1] # Lat from North to south! + ul_lat, lr_lat = ( + lats[0], + lats[-1], + ) # Lat from North to south! ul_lon, lr_lon = lons[0], lons[-1] precip[no_data_mask] = 0 # block_reduce does not handle nan values @@ -363,7 +373,9 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): # Consider that if a single invalid observation is located in the block, # then mark that value as invalid. no_data_mask = block_reduce( - no_data_mask.astype("int"), window_size, axis=(0, 1) + no_data_mask.astype("int"), + window_size, + axis=(0, 1), ).astype(bool) lons, lats = np.meshgrid(lons, lats) @@ -372,11 +384,15 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): if extent is not None: # clip domain ul_lon, lr_lon = _check_coords_range( - (extent[0], extent[1]), "longitude", (ul_lon, lr_lon) + (extent[0], extent[1]), + "longitude", + (ul_lon, lr_lon), ) lr_lat, ul_lat = _check_coords_range( - (extent[2], extent[3]), "latitude", (ul_lat, lr_lat) + (extent[2], extent[3]), + "latitude", + (ul_lat, lr_lat), ) mask_lat = (lats >= lr_lat) & (lats <= ul_lat) @@ -743,7 +759,13 @@ def _import_fmi_pgm_metadata(filename, gzipped=False): @postprocess_import() -def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwargs): +def import_knmi_hdf5( + filename, + qty="ACRR", + accutime=5.0, + pixelsize=1000.0, + **kwargs, +): """ Import a precipitation or reflectivity field (and optionally the quality field) from a HDF5 file conforming to the KNMI Data Centre specification. @@ -810,7 +832,9 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwa # the no data value. The precision of the data is two decimals (0.01 mm). if qty == "ACRR": precip = np.where( - precip_intermediate == 65535, np.nan, precip_intermediate / 100.0 + precip_intermediate == 65535, + np.nan, + precip_intermediate / 100.0, ) # In case reflectivities are imported, the no data value is 255. Values are @@ -818,7 +842,9 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwa # as: dBZ = 0.5 * pixel_value - 32.0 (this used to be 31.5). if qty == "DBZH": precip = np.where( - precip_intermediate == 255, np.nan, precip_intermediate * 0.5 - 32.0 + precip_intermediate == 255, + np.nan, + precip_intermediate * 0.5 - 32.0, ) if precip is None: @@ -948,14 +974,21 @@ def import_mch_gif(filename, product, unit, accutime, **kwargs): # load lookup table if product.lower() == "azc": lut_filename = os.path.join( - os.path.dirname(__file__), "mch_lut_8bit_Metranet_AZC_V104.txt" + os.path.dirname(__file__), + "mch_lut_8bit_Metranet_AZC_V104.txt", ) else: lut_filename = os.path.join( - os.path.dirname(__file__), "mch_lut_8bit_Metranet_v103.txt" + os.path.dirname(__file__), + "mch_lut_8bit_Metranet_v103.txt", ) lut = np.genfromtxt(lut_filename, skip_header=1) - lut = dict(zip(zip(lut[:, 1], lut[:, 2], lut[:, 3]), lut[:, -1])) + lut = dict( + zip( + zip(lut[:, 1], lut[:, 2], lut[:, 3]), + lut[:, -1], + ) + ) # apply lookup table conversion precip = np.zeros(len(img_rgb.getdata())) @@ -971,7 +1004,12 @@ def import_mch_gif(filename, product, unit, accutime, **kwargs): precip[precip < 0] = 0 precip[precip > 9999] = np.nan - elif product.lower() in ["aqc", "cpc", "acquire ", "combiprecip"]: + elif product.lower() in [ + "aqc", + "cpc", + "acquire ", + "combiprecip", + ]: # convert digital numbers to physical values img = np.array(img).astype(int) @@ -1305,9 +1343,13 @@ def import_odim_hdf5(filename, qty="RATE", **kwargs): if "what" in list(dsg[1].keys()): if "quantity" in dsg[1]["what"].attrs.keys(): try: - qty_, gain, offset, nodata, undetect = ( - _read_opera_hdf5_what_group(dsg[1]["what"]) - ) + ( + qty_, + gain, + offset, + nodata, + undetect, + ) = _read_opera_hdf5_what_group(dsg[1]["what"]) what_grp_found = True except KeyError: pass @@ -1531,11 +1573,19 @@ def import_saf_crri(filename, extent=None, **kwargs): if extent: xcoord = ( - np.arange(metadata["x1"], metadata["x2"], metadata["xpixelsize"]) + np.arange( + metadata["x1"], + metadata["x2"], + metadata["xpixelsize"], + ) + metadata["xpixelsize"] / 2 ) ycoord = ( - np.arange(metadata["y1"], metadata["y2"], metadata["ypixelsize"]) + np.arange( + metadata["y1"], + metadata["y2"], + metadata["ypixelsize"], + ) + metadata["ypixelsize"] / 2 ) ycoord = ycoord[::-1] # yorigin = "upper" @@ -1617,3 +1667,502 @@ def _import_saf_crri_geodata(filename): ds_rainfall.close() return geodata + + +@postprocess_import() +def import_dwd_hdf5(filename, qty="RATE", **kwargs): + """ + Import a DWD precipitation product field (and optionally the quality + field) from a HDF5 file conforming to the ODIM specification + + Parameters + ---------- + filename: str + Name of the file to import. + qty: {'RATE', 'ACRR', 'DBZH'} + The quantity to read from the file. The currently supported identitiers + are: 'RATE'=instantaneous rain rate (mm/h), 'ACRR'=hourly rainfall + accumulation (mm) and 'DBZH'=max-reflectivity (dBZ). The default value + is 'RATE'. + + {extra_kwargs_doc} + + Returns + ------- + tuple + A tuple containing: + - data : np.ndarray + The requested precipitation product imported from a HDF5 file + - quality : None + - metadata : dict + Dictionary containing geospatial metadata such as: + - 'projection': PROJ.4 string defining the stereographic projection. + - 'll_lon', 'll_lat': Coordinates of the lower-left corner. + - 'ur_lon', 'ur_lat': Coordinates of the upper-right corner. + - 'x1', 'y1': Carthesian coordinates of the lower-left corner. + - 'x2', 'y2': Carthesian coordinates of the upper-right corner. + - 'xpixelsize', 'ypixelsize': Pixel size in meters. + - 'cartesian_unit': Unit of the coordinate system (meters). + - 'yorigin': Origin of the y-axis ('lower'). + - 'institution': Originating institution ('DWD' or 'DWD Radolan') + - 'accutime': Accumulation period of the requested precipitation product + - 'unit': Unit. + - 'transform': Logarithmic transformation. + - 'zerovalue': No echo value. + - 'threshold': Precipitation threshold. + """ + if not H5PY_IMPORTED: + raise MissingOptionalDependency( + "h5py package is required to import " + "radar reflectivity composites using ODIM HDF5 specification " + "but it is not installed" + ) + + if not PYPROJ_IMPORTED: + raise MissingOptionalDependency( + "pyproj package is required to import " + "DWD's radar reflectivity composite " + "but it is not installed" + ) + + if qty not in ["ACRR", "DBZH", "RATE"]: + raise ValueError( + "unknown quantity %s: the available options are 'ACRR', 'DBZH' and 'RATE'" + ) + + # Open file + f = h5py.File(filename, "r") + precip = None + quality = None + + # Read data recursively + file_content = {} + _read_hdf5_cont(f, file_content) + f.close() + + # Read attributes + data_prop = {} + _get_whatgrp(file_content, data_prop) + + # Get data as well as no data and no echo masks + arr = file_content["dataset1"]["data1"]["data"] + mask_n = arr == data_prop["nodata"] + mask_u = arr == data_prop["undetect"] + mask = np.logical_and(~mask_u, ~mask_n) + + # If the requested quantity in the file + # Transform precipitation data by gain and offset + if data_prop["quantity"] == qty: + precip = np.empty(arr.shape) + precip[mask] = arr[mask] * data_prop["gain"] + data_prop["offset"] + if qty != "DBZH": + precip[mask_u] = data_prop["offset"] + else: + # Set the no echo value manually to -32.5 + # if the file contains horizontal reflectivity + precip[mask_u] = -32.5 + precip[mask_n] = np.nan + # Get possible information about data quality + elif data_prop["quantity"] == "QIND": + quality = np.empty(arr.shape, dtype=float) + quality[mask] = arr[mask] + quality[~mask] = np.nan + + if precip is None: + raise IOError("requested quantity %s not found" % qty) + + # Get the projection and grid information from the HDF5 file + pr = pyproj.Proj(file_content["where"]["projdef"]) + ll_x, ll_y = pr( + file_content["where"]["LL_lon"], + file_content["where"]["LL_lat"], + ) + ur_x, ur_y = pr( + file_content["where"]["UR_lon"], + file_content["where"]["UR_lat"], + ) + + # Determine domain corners in geographic and carthesian coordinates + if len([k for k in file_content["where"].keys() if "_lat" in k]) == 4: + lr_x, lr_y = pr( + file_content["where"]["LR_lon"], + file_content["where"]["LR_lat"], + ) + ul_x, ul_y = pr( + file_content["where"]["UL_lon"], + file_content["where"]["UL_lat"], + ) + x1 = min(ll_x, ul_x) + y1 = min(ll_y, lr_y) + x2 = max(lr_x, ur_x) + y2 = max(ul_y, ur_y) + else: + x1 = ll_x + y1 = ll_y + x2 = ur_x + y2 = ur_y + + # Get the grid cell size + if ( + "where" in file_content["dataset1"].keys() + and "xscale" in file_content["dataset1"]["where"].keys() + ): + xpixelsize = file_content["dataset1"]["where"]["xscale"] + ypixelsize = file_content["dataset1"]["where"]["yscale"] + elif "xscale" in file_content["where"].keys(): + xpixelsize = file_content["where"]["xscale"] + ypixelsize = file_content["where"]["yscale"] + else: + xpixelsize = None + ypixelsize = None + + # Get the unit and transform + if qty == "ACRR": + unit = "mm" + transform = None + elif qty == "DBZH": + unit = "dBZ" + transform = "dB" + else: + unit = "mm/h" + transform = None + + # Extract the time step + startdate = datetime.datetime.strptime( + file_content["dataset1"]["what"]["startdate"] + + file_content["dataset1"]["what"]["starttime"], + "%Y%m%d%H%M%S", + ) + enddate = datetime.datetime.strptime( + file_content["dataset1"]["what"]["enddate"] + + file_content["dataset1"]["what"]["endtime"], + "%Y%m%d%H%M%S", + ) + accutime = (enddate - startdate).total_seconds() / 60.0 + + # Finally, fill out the metadata + metadata = { + "projection": file_content["where"]["projdef"], + "ll_lon": file_content["where"]["LL_lon"], + "ll_lat": file_content["where"]["LL_lat"], + "ur_lon": file_content["where"]["UR_lon"], + "ur_lat": file_content["where"]["UR_lat"], + "x1": x1, + "y1": y1, + "x2": x2, + "y2": y2, + "xpixelsize": xpixelsize, + "ypixelsize": ypixelsize, + "cartesian_unit": "m", + "yorigin": "upper", + "institution": file_content["what"]["source"], + "accutime": accutime, + "unit": unit, + "transform": transform, + "zerovalue": np.nanmin(precip), + "threshold": _get_threshold_value(precip), + } + + metadata.update(kwargs) + + f.close() + + return precip, quality, metadata + + +def _read_hdf5_cont(f, d): + """ + Recursively read nested dictionaries from a HDF5 file. + + + Parameters: + ----------- + f : h5py.Group or h5py.File + The current group or file object from which to read data. + d : dict + The dictionary to populate with the contents of the HDF5 group. + + Returns: + -------- + None. + """ + # Set simple types of hdf content + group_type = h5py._hl.group.Group + + for key, value in f.items(): + if isinstance(value, group_type): + d[key] = {} + if len(list(value.items())) > 0: + # Recurse into non-empty group + _read_hdf5_cont(value, d[key]) + else: + # Handle empty group with attributes + d[key] = {attr: value.attrs[attr] for attr in value.attrs} + d[key] = { + k: (v.decode() if isinstance(v, np.bytes_) else v) + for k, v in d[key].items() + } + + else: + + # Save h5py.Dataset by group name + d[key] = np.array(value) + + return + + +def _get_whatgrp(d, g): + """ + Recursively get attributes of the what group containing + the scaling properties. + + Parameters: + ----------- + d : dict + Dictionary including content of an ODIM compliant + HDF5 file. + g : dict + Dictionary containing attributes of what group. + + + Returns: + -------- + None. + """ + if "what" in d.keys(): + # Searching for the corresponding what group + # that contains the scaling properties + if "gain" in d["what"].keys(): + g.update(d["what"]) + else: + k = [k for k in d.keys() if "data" in k][0] + _get_whatgrp(d[k], g) + else: + raise DataModelError( + "Non ODIM compliant file: " + "no what group found from {} " + "or its subgroups".format(d.keys()[0]) + ) + return + + +@postprocess_import() +def import_dwd_radolan(filename, product_name): + """ + Import a RADOLAN precipitation product from a binary file. + + Parameters + ---------- + filename: str + Name of the file to import. + product_name: {'WX','RX','EX','RY','RW','AY','RS','YW','WN'} + The specific product to read from the file. Please see + https://www.dwd.de/DE/leistungen/radolan/radolan_info/ + radolan_radvor_op_komposit_format_pdf.pdf + for a detailed description. + + {extra_kwargs_doc} + + Returns + ------- + tuple + A tuple containing: + - data : np.ndarray + The desired precipitation product in mm/h imported from a RADOLAN file + - quality : None + - metadata : dict + Dictionary containing geospatial metadata such as: + - 'projection': PROJ.4 string defining the stereographic projection. + - 'xpixelsize', 'ypixelsize': Pixel size in meters. + - 'cartesian_unit': Unit of the coordinate system (meters). + - 'yorigin': Origin of the y-axis ('upper'). + - 'x1', 'y1': Coordinates of the lower-left corner. + - 'x2', 'y2': Coordinates of the upper-right corner. + """ + # Determine file size and header size + size_file = os.path.getsize(filename) + size_data = np.round(size_file, -3) + size_header = size_file - size_data + + # Open file and read header + f = open(filename, "rb") + header = f.read(size_header).decode("utf-8") + + # Get product name from header + product = header[:2] + + # Check if its the desired product + assert product == product_name, "Product not in File!" + + # Distinguish between products saved with 8 or 16bit + product_cat1 = np.array(["WX", "RX"]) + product_cat2 = np.array(["RY", "RW", "YW"]) + + # Determine byte size and data type + nbyte = 1 if product in product_cat1 else 2 + signed = "B" if product in product_cat1 else "H" + + # Extract the scaling factor and grid dimensions + fac = int(header.split("E-")[1].split("INT")[0]) + dimsplit = header.split("x") + dims = np.array((dimsplit[0][-4:], dimsplit[1][:4]), dtype=int)[::-1] + + # Read binary data + data = array.array(signed) + data.fromfile(f, size_data // nbyte) + f.close() + + # Reshape and transpose data to match grid layout + data = np.array(np.reshape(data, dims, order="F"), dtype=float).T + + # Define no-echo values based on product type + if product == "SF": + no_echo_value = 0.0 + elif product in product_cat2: + no_echo_value = -0.01 + else: + no_echo_value = -32.5 + + # Apply scaling and handle missing data + if product in product_cat1: + data[data >= 249] = np.nan + data = data / 2.0 + no_echo_value + elif product in product_cat2: + data, no_data_mask = _identify_info_bits(data) + if product == "AY": + data = (10 ** (-fac)) * data / 2.0 + no_echo_value + else: + data = (10 ** (-fac)) * data + no_echo_value + else: + data, no_data_mask = _identify_info_bits(data) + data = (10 ** (-fac)) * data / 2.0 + no_echo_value + + # Mask out no-data values + data[no_data_mask] = np.nan + + # Load geospatial metadata + geodata = _import_dwd_geodata(product_name, dims) + metadata = geodata + + return data, None, metadata + + +def _identify_info_bits(data): + """ + Identifies and processes information bits embedded in RADOLAN data values. + + This function decodes metadata flags embedded in the RADOLAN data array, such as: + - Clutter (bit 16) + - Negative values (bit 15) + - No data (bit 14) + - Secondary data (bit 13) + + Parameters + ---------- + data : np.ndarray + The raw RADOLAN data array containing encoded information bits. + + Returns + ------- + tuple + A tuple containing: + - data : np.ndarray + The cleaned and decoded data array. + - no_data_mask : np.ndarray + A boolean mask indicating positions of no-data values. + """ + # Identify and remove clutter (bit 16) + clutter_mask = data - 2**15 >= 0.0 + data[clutter_mask] = 0 + + # Identify and convert negative values (bit 15) + mask = data - 2**14 >= 0.0 + data[mask] -= 2**14 + + # Identify no-data values (bit 14) + no_data_mask = data - 2**13 == 2500.0 + if np.sum(no_data_mask) == 0.0: + no_data_mask = data - 2**13 == 0.0 + data[no_data_mask] = 0.0 + + # Identify and remove secondary data flag (bit 13) + data[data - 2**12 > 0.0] -= 2**12 + + # Apply negative sign to previously marked negative values + data[mask] *= -1 + + return data, no_data_mask + + +def _import_dwd_geodata(product_name, dims): + """ + Generate geospatial metadata for RADOLAN precipitation products. + + Since RADOLAN binary files contain only limited projection metadata, + this function provides hard-coded geospatial definitions and calculates + the bounding box of the data grid based on the product type and dimensions. + + Parameters + ---------- + product : str + The RADOLAN product code (e.g., 'RX', 'WX', 'WN', etc.). + dims : tuple of int + The dimensions of the data grid (rows, columns). + + Returns + ------- + geodata : dict + A dictionary containing: + - 'projection': PROJ.4 string defining the stereographic projection. + - 'xpixelsize', 'ypixelsize': Pixel size in meters. + - 'cartesian_unit': Unit of the coordinate system (meters). + - 'yorigin': Origin of the y-axis ('upper'). + - 'x1', 'y1': Coordinates of the lower-left corner. + - 'x2', 'y2': Coordinates of the upper-right corner. + """ + geodata = {} + + # Define stereographic projection used by RADOLAN + projdef = ( + "+a=6378137.0 +b=6356752.0 +proj=stere +lat_ts=60.0 " + "+lat_0=90.0 +lon_0=10.0 +x_0=0 +y_0=0" + ) + geodata["projection"] = projdef + # Spatial resolution of 1km + geodata["xpixelsize"] = 1000.0 + geodata["ypixelsize"] = 1000.0 + geodata["cartesian_unit"] = "m" + geodata["yorigin"] = "upper" + + # Define product categories + product_cat1 = ["RX", "RY", "RW"] + product_cat2 = ["WN"] + product_cat3 = ["WX", "YW"] + + # Assign reference coordinates based on product type + if product_name in product_cat1: + lon, lat = 3.604382995, 46.95361536 + elif product_name in product_cat2: + lon, lat = 3.566994635, 45.69642538 + elif product_name in product_cat3: + lon, lat = 9.0, 51.0 + + # Project reference coordinates to Cartesian system + pr = pyproj.Proj(projdef) + x1, y1 = pr(lon, lat) + + # Adjust origin for center-based products + if product_name in product_cat3: + x1 -= dims[0] * 1000 // 2 + y1 -= dims[1] * 1000 // 2 - 80000 + + # Calculate bounding box + x2 = x1 + dims[0] * 1000 + y2 = y1 + dims[1] * 1000 + + geodata["x1"] = x1 + geodata["y1"] = y1 + geodata["x2"] = x2 + geodata["y2"] = y2 + + return geodata diff --git a/pysteps/io/interface.py b/pysteps/io/interface.py index b64ddd8e..8a8cb3f9 100644 --- a/pysteps/io/interface.py +++ b/pysteps/io/interface.py @@ -20,15 +20,17 @@ _importer_methods = dict( bom_rf3=importers.import_bom_rf3, + dwd_hdf5=importers.import_dwd_hdf5, + dwd_radolan=importers.import_dwd_radolan, fmi_geotiff=importers.import_fmi_geotiff, fmi_pgm=importers.import_fmi_pgm, + knmi_hdf5=importers.import_knmi_hdf5, mch_gif=importers.import_mch_gif, mch_hdf5=importers.import_mch_hdf5, mch_metranet=importers.import_mch_metranet, mrms_grib=importers.import_mrms_grib, odim_hdf5=importers.import_odim_hdf5, opera_hdf5=importers.import_opera_hdf5, - knmi_hdf5=importers.import_knmi_hdf5, saf_crri=importers.import_saf_crri, ) diff --git a/pysteps/pystepsrc b/pysteps/pystepsrc index 9b06ce75..91088e00 100644 --- a/pysteps/pystepsrc +++ b/pysteps/pystepsrc @@ -24,6 +24,17 @@ "gzipped": true } }, + "dwd": { + "root_path": "./radar/dwd", + "path_fmt": "%Y/%m/%d", + "fn_pattern": "%Y%m%d_%H%M_RY", + "fn_ext": "h5", + "importer": "dwd_hdf5", + "timestep": 5, + "importer_kwargs": { + "qty": "RATE" + } + }, "fmi": { "root_path": "./radar/fmi/pgm", "path_fmt": "%Y%m%d", @@ -121,16 +132,18 @@ "gzipped": true } }, - "rmi_nwp": { - "root_path": "./nwp/rmi", + "dwd_nwp": { + "root_path": "./nwp/dwd", "path_fmt": "%Y/%m/%d", - "fn_pattern": "ao13_%Y%m%d%H_native_5min", - "fn_ext": "nc", - "importer": "rmi_nwp", + "fn_ext": ".grib2", + "fn_pattern": "%Y%m%d_%H%M_PR_GSP_060_120", + "importer": "dwd_nwp", "timestep": 5, "importer_kwargs": { - "gzipped": true + "varname": "PR_GSP", + "grid_file_path": "./aux/grid_files/dwd/icon/R19B07/icon_grid_0047_R19B07_L.nc" } + }, "knmi_nwp": { "root_path": "./nwp/knmi", @@ -142,6 +155,17 @@ "importer_kwargs": { "gzipped": true } + }, + "rmi_nwp": { + "root_path": "./nwp/rmi", + "path_fmt": "%Y/%m/%d", + "fn_pattern": "ao13_%Y%m%d%H_native_5min", + "fn_ext": "nc", + "importer": "rmi_nwp", + "timestep": 5, + "importer_kwargs": { + "gzipped": true + } } } } diff --git a/pysteps/tests/test_io_dwd_hdf5.py b/pysteps/tests/test_io_dwd_hdf5.py new file mode 100644 index 00000000..22626888 --- /dev/null +++ b/pysteps/tests/test_io_dwd_hdf5.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8 -*- + +import os + +import pytest + +import pysteps +from pysteps.tests.helpers import smart_assert + +pytest.importorskip("h5py") + +# Test for RADOLAN RY product + +root_path = pysteps.rcparams.data_sources["dwd"]["root_path"] + +filename = os.path.join(root_path, "RY", "2025", "06", "04", "20250604_1700_RY.h5") +precip_ry, _, metadata_ry = pysteps.io.import_dwd_hdf5(filename, qty="RATE") + + +def test_io_import_dwd_hdf5_ry_shape(): + """Test the importer DWD HDF5.""" + assert precip_ry.shape == (1200, 1100) + + +# Test_metadata +# Expected projection definition +expected_proj = ( + "+proj=stere +lat_0=90 +lat_ts=60 " + "'+lon_0=10 +a=6378137 +b=6356752.3142451802" + "+no_defs +x_0=543196.83521776402 " + "+y_0=3622588.8619310018 +units=m" +) + +# List of (variable,expected,tolerance) tuples +test_ry_attrs = [ + ("projection", expected_proj, None), + ("ll_lon", 3.566994635, 1e-10), + ("ll_lat", 45.69642538, 1e-10), + ("ur_lon", 18.73161645, 1e-10), + ("ur_lat", 55.84543856, 1e-10), + ("x1", -500.0, 1e-6), + ("y1", -1199500.0, 1e-10), + ("x2", 1099500.0, 1e-10), + ("y2", 500.0, 1e-6), + ("xpixelsize", 1000.0, 1e-10), + ("xpixelsize", 1000.0, 1e-10), + ("cartesian_unit", "m", None), + ("yorigin", "upper", None), + ("institution", "ORG:78,CTY:616,CMT:Deutscher Wetterdienst radolan@dwd.de", None), + ("accutime", 5.0, 1e-10), + ("unit", "mm/h", None), + ("transform", None, None), + ("zerovalue", 0.0, 1e-10), + ("threshold", 0.12, 1e-10), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", test_ry_attrs) +def test_io_import_opera_hdf5_odyssey_dataset_attrs(variable, expected, tolerance): + """Test the importer OPERA HDF5.""" + smart_assert(metadata_ry[variable], expected, tolerance) diff --git a/pysteps/tests/test_utils_reprojection.py b/pysteps/tests/test_utils_reprojection.py index e6eb4fe5..9e0c8520 100644 --- a/pysteps/tests/test_utils_reprojection.py +++ b/pysteps/tests/test_utils_reprojection.py @@ -5,8 +5,11 @@ import pytest import pysteps from pysteps.utils import reprojection as rpj +import pysteps_nwp_importers pytest.importorskip("rasterio") +pytest.importorskip("pyproj") +# pytest.importorskip("pysteps_nwp_importers") root_path_radar = pysteps.rcparams.data_sources["rmi"]["root_path"] @@ -109,4 +112,66 @@ def test_utils_reproject_grids( assert ( metadata_reproj["projection"] == metadata_dst["projection"] - ), "projection is different than destionation projection" + ), "projection is different than destination projection" + + +root_path_radar = pysteps.rcparams.data_sources["dwd"]["root_path"] +root_path_nwp = pysteps.rcparams.data_sources["dwd_nwp"]["root_path"] + +filename_radar = os.path.join( + root_path_radar, "RY", "2025", "06", "04", "20250604_1700_RY.h5" +) +_, _, metadata_dst = pysteps.io.import_dwd_hdf5(filename_radar, qty="RATE") + +filename_nwp = os.path.join( + root_path_nwp, "2025", "06", "04", "20250604_1600_PR_GSP_060_120.grib2" +) +kwargs = { + "varname": "PR_GSP", + "grid_file_path": "./aux/grid_files/dwd/icon/R19B07/icon_grid_0047_R19B07_L.nc", +} +array_src, _, metadata_src = pysteps_nwp_importers.importer_dwd_nwp.import_dwd_nwp( + filename_nwp, **kwargs +) + +# Since output of NWP importer is based on xarray and the restructure function is +# rewritten for np.ndarray, there's is some improvisation +metadata_src["clon"] = array_src["longitude"].values +metadata_src["clat"] = array_src["latitude"].values +array_src = array_src.values + +restructure_arg_names = ("array_src", "metadata_src", "metadata_dst") +restructure_arg_values = [(array_src, metadata_src, metadata_dst)] + + +@pytest.mark.parametrize(restructure_arg_names, restructure_arg_values) +def test_utils_unstructured2regular(array_src, metadata_src, metadata_dst): + # Run unstructured2regular + array_rprj, metadata_rprj = rpj.unstructured2regular( + array_src, metadata_src, metadata_dst + ) + + # The tests + assert ( + array_rprj.shape[0] == array_src.shape[0] + ), "Time dimension has not the same length as source" + assert ( + array_rprj.shape[1] == array_src.shape[0] + ), "Ensemble member dimension has not the same length as source" + + assert ( + metadata_rprj["x1"] == metadata_dst["x1"] + ), "x-value lower left corner is not equal to radar composite" + assert ( + metadata_rprj["x2"] == metadata_dst["x2"] + ), "x-value upper right corner is not equal to radar composite" + assert ( + metadata_rprj["y1"] == metadata_dst["y1"] + ), "y-value lower left corner is not equal to radar composite" + assert ( + metadata_rprj["y2"] == metadata_dst["y2"] + ), "y-value upper right corner is not equal to radar composite" + + assert ( + metadata_rprj["projection"] == metadata_dst["projection"] + ), "projection is different than destination projection" diff --git a/pysteps/utils/reprojection.py b/pysteps/utils/reprojection.py index 6144a52c..1d070550 100644 --- a/pysteps/utils/reprojection.py +++ b/pysteps/utils/reprojection.py @@ -12,6 +12,7 @@ reproject_grids """ from pysteps.exceptions import MissingOptionalDependency +from scipy.interpolate import griddata import numpy as np @@ -23,6 +24,13 @@ except ImportError: RASTERIO_IMPORTED = False +try: + import pyproj + + PYPROJ_IMPORTED = True +except ImportError: + PYPROJ_IMPORTED = False + def reproject_grids(src_array, dst_array, metadata_src, metadata_dst): """ @@ -118,3 +126,115 @@ def reproject_grids(src_array, dst_array, metadata_src, metadata_dst): metadata[key] = metadata_dst[key] return r_rprj, metadata + + +def unstructured2regular(src_array, metadata_src, metadata_dst): + """ + Reproject unstructured data onto a regular grid on the assumption that + both src data and dst grid have the same projection. + + Parameters + ---------- + src_array: np.ndarray + Three-dimensional array of shape (t, n_ens, n_gridcells) containing a + time series of precipitation ensemble forecasts. These precipitation + fields will be reprojected. + metadata_src: dict + Metadata dictionary containing the projection, clon, clat, and ngridcells + and attributes of the src_array as described in the documentation of + :py:mod:`pysteps.io.importers`. + metadata_dst: dict + Metadata dictionary containing the projection, x- and ypixelsize, x1 and + y2 attributes of the dst_array. + + Returns + ------- + tuple + A tuple containing: + - r_rprj: np.ndarray + Four dimensional array of shape (t, n_ens, x, y) containing the + precipitation fields of src_array, but reprojected to the grid + of dst_array. + - metadata: dict + Dictionary containing geospatial metadat such as: + - 'projection' : PROJ.4 string defining the stereographic projection. + - 'xpixelsize', 'ypixelsize': Pixel size in meters. + - 'x1', 'y1': Carthesian coordinates of the lower-left corner. + - 'x2', 'y2': Carthesian coordinates of the upper-right corner. + - 'cartesian_unit': Unit of the coordinate system (meters). + """ + + if not PYPROJ_IMPORTED: + raise MissingOptionalDependency( + "pyproj package is required to reproject DWD's NWP data" + "but it is not installed" + ) + + if not "clon" in metadata_src.keys(): + raise KeyError("Center longitude (clon) is missing in metadata_src") + if not "clat" in metadata_src.keys(): + raise KeyError("Center latitude (clat) is missing in metadata_src") + + # Get number of grid cells + Nc = metadata_src["clon"].shape[0] + ic_in = np.arange(Nc) + + # Get cartesian coordinates of destination grid + x_dst = np.arange( + np.float32(metadata_dst["x1"]), + np.float32(metadata_dst["x2"]), + metadata_dst["xpixelsize"], + ) + + y_dst = np.arange( + np.float32(metadata_dst["y1"]), + np.float32(metadata_dst["y2"]), + metadata_dst["ypixelsize"], + ) + + # Create destination grid + if metadata_dst["yorigin"] == "upper": + y_dst = y_dst[::-1] + xx_dst, yy_dst = np.meshgrid(x_dst, y_dst) + s_out = yy_dst.shape + + # Extract the grid info of src_array assuming the same projection of src and dst + pr = pyproj.Proj(metadata_dst["projection"]) + x_src, y_src = pr(metadata_src["clon"], metadata_src["clat"]) + + # Create array of x-y pairs for interpolation + P_in = np.stack((x_src, y_src)).T + P_out = np.array((xx_dst.flatten(), yy_dst.flatten())).T + + # Nearest neighbor interpolation of x-y pairs + ic_out = ( + griddata(P_in, ic_in.flatten(), P_out, method="nearest") + .reshape(s_out) + .astype(int) + ) + + # Apply interpolation on all time steps and ensemble members + r_rprj = np.array( + [ + [src_array[i, j][ic_out] for j in range(src_array.shape[1])] + for i in range(src_array.shape[0]) + ] + ) + + # Update the src metadata + metadata = metadata_src.copy() + + for key in [ + "projection", + "yorigin", + "xpixelsize", + "ypixelsize", + "x1", + "x2", + "y1", + "y2", + "cartesian_unit", + ]: + metadata[key] = metadata_dst[key] + + return r_rprj, metadata