refactoring : added proposed import_DPC_hdf5

pysteps/io/importers.py

@@ -1452,6 +1452,183 @@ def import_opera_hdf5(filename, qty="RATE", **kwargs):
 
     return precip, quality, metadata
 
+@postprocess_import()
+def import_DPC_hdf5(filename, qty="RATE", **kwargs):
+    """Import a precipitation field (and optionally the quality field) from a DPC
+    OPERA HDF5 file conforming to the ODIM specification 2.1.
+
+    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
+    -------
+
+    out : tuple
+        A three-element tuple containing the OPERA product for the requested
+        quantity and the associated quality field and metadata. The quality
+        field is read from the file if it contains a dataset whose quantity
+        identifier is 'QIND'. The quality dataset belongs to the quality subgroup
+        under the group of the main variable, having quantity attribute of the
+        dataset in {'RATE','ACRR','DBZH'} .
+    """
+    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 qty not in ["ACRR", "DBZH", "RATE"]:
+        raise ValueError(
+            "unknown quantity %s: the available options are 'ACRR', 'DBZH' and 'RATE'"
+        )
+
+    f = h5py.File(filename, "r")
+
+    precip = None
+    quality = None
+
+    for dsg in f.items():
+        if dsg[0].startswith("dataset"):
+            what_grp_found = False
+            # check if the "what" group is in the "dataset" group
+            if "what" in list(dsg[1].keys()):
+                if "quantity" in dsg[1]["what"].attrs.keys():
+                    qty_, gain, offset, nodata, undetect = _read_opera_hdf5_what_group(
+                        dsg[1]["what"]
+                    )
+                    what_grp_found = True
+
+            for dg in dsg[1].items():
+                if "what" in list(dg[1].keys()):
+                    qty_, gain, offset, nodata, undetect = _read_opera_hdf5_what_group( dg[1]['what'] ) 
+
+                    if( qty_.decode()==qty ):
+                    # In DPC files precip data belong to the dataset whose quantity attribute is RATE.
+                    # We look for eventual quality data inside the same group.
+                        dg_items = list(dg[1].items())
+
+                        for dg2 in dg_items:
+                            if('data' in dg2[0] ):
+                                arr = dg[1]["data"][...]
+                                mask_n = arr == nodata
+                                mask_u = arr == undetect
+                                mask = np.logical_and(~mask_u, ~mask_n)
+                                precip = np.empty(arr.shape)
+                                precip[mask] = arr[mask] * gain + offset
+                                precip[mask_u] = 0.0
+                                precip[mask_n] = np.nan
+                            elif not what_grp_found:
+                                raise DataModelError(
+                                    "Non ODIM compilant file: "
+                                    "no what group found from {} "
+                                    "or its subgroups".format(dg[0])
+                                )
+                            elif('quality' in dg2[0]):
+                                if('what' in list(dg2[1].keys())):                                
+                                     qty_, gain, offset, nodata, undetect = _read_opera_hdf5_what_group( dg2[1]['what'] )
+                                    if(qty_.decode()=='QIND'):
+                                        arr = dg2[1]["data"][...]
+                                        mask_n = arr == nodata
+                                        mask_u = arr == undetect
+                                        mask = np.logical_and(~mask_u, ~mask_n)
+                                        quality = np.empty(arr.shape)#, dtype=float)
+                                        quality[mask] = arr[mask] * gain + offset
+                                        quality[~mask] = np.nan
+
+    if precip is None:
+        raise IOError("requested quantity %s not found" % qty)
+
+    where = f["where"]
+    proj4str = where.attrs["projdef"].decode()
+    pr = pyproj.Proj(proj4str)
+
+    ll_lat = where.attrs["LL_lat"]
+    ll_lon = where.attrs["LL_lon"]
+    ur_lat = where.attrs["UR_lat"]
+    ur_lon = where.attrs["UR_lon"]
+    if (
+        "LR_lat" in where.attrs.keys()
+        and "LR_lon" in where.attrs.keys()
+        and "UL_lat" in where.attrs.keys()
+        and "UL_lon" in where.attrs.keys()
+    ):
+        lr_lat = float(where.attrs["LR_lat"])
+        lr_lon = float(where.attrs["LR_lon"])
+        ul_lat = float(where.attrs["UL_lat"])
+        ul_lon = float(where.attrs["UL_lon"])
+        full_cornerpts = True
+    else:
+        full_cornerpts = False
+
+    ll_x, ll_y = pr(ll_lon, ll_lat)
+    ur_x, ur_y = pr(ur_lon, ur_lat)
+
+    if full_cornerpts:
+        lr_x, lr_y = pr(lr_lon, lr_lat)
+        ul_x, ul_y = pr(ul_lon, 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
+
+    if "xscale" in where.attrs.keys() and "yscale" in where.attrs.keys():
+        xpixelsize = where.attrs["xscale"]
+        ypixelsize = where.attrs["yscale"]
+    else:
+        xpixelsize = None
+        ypixelsize = None
+
+    if qty == "ACRR":
+        unit = "mm"
+        transform = None
+    elif qty == "DBZH":
+        unit = "dBZ"
+        transform = "dB"
+    else:
+        unit = "mm/h"
+        transform = None
+
+    metadata = {
+        "projection": proj4str,
+        "ll_lon": ll_lon,
+        "ll_lat": ll_lat,
+        "ur_lon": ur_lon,
+        "ur_lat": ur_lat,
+        "x1": x1,
+        "y1": y1,
+        "x2": x2,
+        "y2": y2,
+        "xpixelsize": xpixelsize,
+        "ypixelsize": ypixelsize,
+        "yorigin": "upper",
+        "institution": "Odyssey datacentre",
+        "accutime": 15.0,
+        "unit": unit,
+        "transform": transform,
+        "zerovalue": np.nanmin(precip),
+        "threshold": _get_threshold_value(precip),
+    }
+
+    f.close()
+
+    return precip, quality, metadata
 
 def _read_opera_hdf5_what_group(whatgrp):
     qty = whatgrp.attrs["quantity"]