Reprojection functionality

pysteps/io/importers.py

@@ -917,7 +917,7 @@ def _import_fmi_pgm_metadata(filename, gzipped=False):
 
 
 @postprocess_import()
-def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.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.
 
@@ -1013,7 +1013,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
     # The 'where' group of mch- and Opera-data, is called 'geographic' in the
     # KNMI data.
     geographic = f["geographic"]
-    proj4str = geographic["map_projection"].attrs["projection_proj4_params"].decode()
+    proj4str = "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0"
     pr = pyproj.Proj(proj4str)
     metadata["projection"] = proj4str
 
@@ -1044,7 +1044,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
     metadata["y2"] = y2
     metadata["xpixelsize"] = pixelsize
     metadata["ypixelsize"] = pixelsize
-    metadata["cartesian_unit"] = "km"
+    metadata["cartesian_unit"] = "m"
     metadata["yorigin"] = "upper"
     metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute"
     metadata["accutime"] = accutime

pysteps/io/nwp_importers.py

@@ -517,10 +517,21 @@ def _import_knmi_nwp_geodata_xr(
     xpixelsize = abs(ds_in.x[1] - ds_in.x[0])
     ypixelsize = abs(ds_in.y[1] - ds_in.y[0])
 
+    x_coords = (
+        np.arange(xmin, xmin + float(xpixelsize) * ds_in.x.shape[0], float(xpixelsize))
+        + float(xpixelsize) / 2
+    )
+    y_coords = (
+        np.arange(ymin, ymin + float(ypixelsize) * ds_in.y.shape[0], float(ypixelsize))
+        + float(ypixelsize) / 2
+    )
+
     cartesian_unit = ds_in.x.units
 
-    # Add metadata needed by pySTEPS as attrs in X and Y variables
+    # Assign new coordinates that are in the middle of the cell
+    ds_in = ds_in.assign_coords(x=x_coords, y=y_coords)
 
+    # Add metadata needed by pySTEPS as attrs in X and Y variables
     ds_in.x.attrs.update(
         {
             # TODO: Remove before final 2.0 version
@@ -541,6 +552,8 @@ def _import_knmi_nwp_geodata_xr(
 
     # Add metadata needed by pySTEPS as attrs in rainfall variable
     da_rainfall = ds_in[varname].isel({varname_time: 0})
+    # Set values below 0.0 to 0.0
+    da_rainfall = da_rainfall.where(da_rainfall >= 0.0, 0.0)
 
     ds_in[varname].attrs.update(
         {

pysteps/utils/__init__.py

@@ -11,3 +11,4 @@
 from .spectral import *
 from .tapering import *
 from .transformation import *
+from .reprojection import *

pysteps/utils/interface.py

@@ -18,6 +18,7 @@
 from . import fft
 from . import images
 from . import interpolate
+from . import reprojection
 from . import spectral
 from . import tapering
 from . import transformation
@@ -187,6 +188,9 @@ def donothing(R, metadata=None, *args, **kwargs):
     methods_objects["rbfinterp2d"] = interpolate.rbfinterp2d
     methods_objects["idwinterp2d"] = interpolate.idwinterp2d
 
+    # reprojection methods
+    methods_objects["reprojection"] = reprojection.reprojection
+
     # spectral methods
     methods_objects["rapsd"] = spectral.rapsd
     methods_objects["rm_rdisc"] = spectral.remove_rain_norain_discontinuity

pysteps/utils/reprojection.py

@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+
+from rasterio import Affine as A
+from rasterio.warp import reproject, Resampling
+import xarray as xr
+import numpy as np
+
+
+def reprojection(R_src, R_dst):
+    """Reprojects precipitation fields to the domain of another precipiation
+    field.
+
+    Parameters
+    ----------
+    R_src: xarray
+        Three-dimensional xarray with dimensions (t, x, y) containing a
+        time series of precipitation fields. These precipitaiton fields
+        will be reprojected.
+    R_dst: xarray
+        Xarray containing a precipitation field or a time series of precipitation
+        fields. The xarray R_src will be reprojected to the domain of R_dst.
+
+    Returns
+    -------
+    R_rprj: xarray
+        Three-dimensional xarray with dimensions (t, x, y) containing the
+        precipitation fields of R_src, but reprojected to the domain of
+        R_dst.
+    """
+
+    # Extract the grid info from R_src
+    src_crs = R_src.attrs["projection"]
+    x1_src = R_src.x.attrs["x1"]
+    y2_src = R_src.y.attrs["y2"]
+    xpixelsize_src = R_src.attrs["xpixelsize"]
+    ypixelsize_src = R_src.attrs["ypixelsize"]
+    src_transform = A.translation(float(x1_src), float(y2_src)) * A.scale(
+        float(xpixelsize_src), float(-ypixelsize_src)
+    )
+
+    # Extract the grid info from R_dst
+    dst_crs = R_dst.attrs["projection"]
+    x1_dst = R_dst.x.attrs["x1"]
+    y2_dst = R_dst.y.attrs["y2"]
+    xpixelsize_dst = R_dst.attrs["xpixelsize"]
+    ypixelsize_dst = R_dst.attrs["ypixelsize"]
+    dst_transform = A.translation(float(x1_dst), float(y2_dst)) * A.scale(
+        float(xpixelsize_dst), float(-ypixelsize_dst)
+    )
+
+    # Initialise the reprojected (x)array
+    R_rprj = np.zeros((R_src.shape[0], R_dst.shape[-2], R_dst.shape[-1]))
+
+    # For every timestep, reproject the precipitation field of R_src to
+    # the domain of R_dst
+    if R_src.attrs["yorigin"] != R_dst.attrs["yorigin"]:
+        R_src = R_src[:, ::-1, :]
+
+    for i in range(R_src.shape[0]):
+        reproject(
+            R_src.values[i, :, :],
+            R_rprj[i, :, :],
+            src_transform=src_transform,
+            src_crs=src_crs,
+            dst_transform=dst_transform,
+            dst_crs=dst_crs,
+            resampling=Resampling.nearest,
+            dst_nodata=np.nan,
+        )
+
+    # Assign the necessary attributes from R_src and R_dst to R_rprj
+    R_rprj = xr.DataArray(
+        data=R_rprj,
+        dims=("t", "y", "x"),
+        coords=dict(
+            t=("t", R_src.coords["t"].data),
+            x=("x", R_dst.coords["x"].data),
+            y=("y", R_dst.coords["y"].data),
+        ),
+    )
+    R_rprj.attrs.update(R_src.attrs)
+    R_rprj.x.attrs.update(R_dst.x.attrs)
+    R_rprj.y.attrs.update(R_dst.y.attrs)
+    for key in ["projection", "yorigin", "xpixelsize", "ypixelsize"]:
+        R_rprj.attrs[key] = R_dst.attrs[key]
+
+    return R_rprj