Reprojection functionality (#236)

* Added Lesley's reprojection module to this branch

* Added compatibility for three-dimensional xarrays

* Add commentary to reprojection util

* Changes to make reprojection of KNMI data possible

* Changes after Daniele's review

* Add dependencies

* Changes to importers, see issue #215

* Add tests

* Fix some issues

* documentation

* Fixes for tests

* Set requirements again

* Some fixes

* Changes to nwp_importers after Carlos' response

* Remove wrong example script

* Remove rasterio dependencies from lists

* First try to prevent testing error

* Changes Daniele and fix knmi nwp importer

* Add rasterio to tox.ini

* Aesthetics

* rasterio import test

* Add rasterio to the test dependencies

* Reset try-except functionality for rasterio import

* Fix for failing test on windows python 3.6

* add importerskip rasterio

Co-authored-by: Wout Dewettinck <wout.dewettinck@ugent.be>

.github/workflows/test_pysteps.yml

@@ -2,7 +2,7 @@ name: Test Pysteps
 
 env:
   MINIMAL_DEPENDENCIES: cython numpy jsmin jsonschema matplotlib netCDF4 opencv pillow pyproj scipy dask
-  OPTIONAL_DEPENDENCIES: pyfftw cartopy h5py PyWavelets pandas scikit-image
+  OPTIONAL_DEPENDENCIES: pyfftw cartopy h5py PyWavelets pandas scikit-image rasterio
   TEST_DEPENDENCIES: pytest pytest-cov
 
 on:

doc/requirements.txt

@@ -1,4 +1,4 @@
-# Additional requeriments related to the documentation build only
+# Additional requirements related to the documentation build only
 sphinx
 numpydoc
 sphinxcontrib.bibtex

doc/source/pysteps_reference/utils.rst

@@ -16,3 +16,4 @@ Implementation of miscellaneous utility functions.
 .. automodule:: pysteps.utils.spectral
 .. automodule:: pysteps.utils.tapering
 .. automodule:: pysteps.utils.transformation
+.. automodule:: pysteps.utils.reprojection

environment_dev.yml

@@ -31,3 +31,4 @@ dependencies:
   - scikit-image
   - pandas
   - xarray
+  - rasterio

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.
 
@@ -934,8 +934,8 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
         is used as default, but hourly, daily and monthly accumulations
         are also available.
     pixelsize: float
-        The pixel size of a raster cell in kilometers. The default value for the
-        KNMI datasets is a 1 km grid cell size, but datasets with 2.4 km pixel
+        The pixel size of a raster cell in meters. The default value for the
+        KNMI datasets is a 1000 m grid cell size, but datasets with 2400 m pixel
         size are also available.
 
     {extra_kwargs_doc}
@@ -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

@@ -371,8 +371,20 @@ def _import_rmi_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
 
+    # 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(
@@ -517,16 +529,28 @@ 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])
 
-    cartesian_unit = ds_in.x.units
+    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
+    )
 
-    # Add metadata needed by pySTEPS as attrs in X and Y variables
+    cartesian_unit_x = ds_in.x.units
+    cartesian_unit_y = ds_in.y.units
+
+    # 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
             "x1": xmin,
             "x2": xmax,
-            "cartesian_unit": cartesian_unit,
+            "cartesian_unit": cartesian_unit_x,
         }
     )
 
@@ -535,12 +559,14 @@ def _import_knmi_nwp_geodata_xr(
             # TODO: Remove before final 2.0 version
             "y1": ymin,
             "y2": ymax,
-            "cartesian_unit": cartesian_unit,
+            "cartesian_unit": cartesian_unit_y,
         }
     )
 
     # 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/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"

pysteps/tests/test_io_bom_nwp.py

@@ -53,7 +53,7 @@ def test_io_import_bom_nwp_xarray_attrs(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_x = [
-    ("units", "m", None),
+    ("cartesian_unit", "m", None),
     ("x1", -127750.0, 0.1),
     ("x2", 127750.0, 0.1),
 ]
@@ -66,7 +66,7 @@ def test_io_import_bom_nwp_xarray_attrs_coordx(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_y = [
-    ("units", "m", None),
+    ("cartesian_unit", "m", None),
     ("y1", -127750.0, 0.1),
     ("y2", 127750.0, 0.1),
 ]

pysteps/tests/test_io_knmi_hdf5.py

@@ -25,9 +25,7 @@ def test_io_import_knmi_hdf5_shape():
 # test_metadata: list of (variable,expected, tolerance) tuples
 
 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"
+    "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0"
 )
 
 # list of (variable,expected,tolerance) tuples

pysteps/tests/test_io_knmi_nwp.py

@@ -53,7 +53,7 @@ def test_io_import_knmi_nwp_xarray_attrs(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_x = [
-    ("units", "degrees_east", None),
+    ("cartesian_unit", "degrees_east", None),
     ("x1", 0.0, 0.0001),
     ("x2", 11.063, 0.0001),
 ]
@@ -66,7 +66,7 @@ def test_io_import_knmi_nwp_xarray_attrs_coordx(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_y = [
-    ("units", "degrees_north", None),
+    ("cartesian_unit", "degrees_north", None),
     ("y1", 49.0, 0.0001),
     ("y2", 55.877, 0.0001),
 ]

pysteps/tests/test_io_rmi_nwp.py

@@ -56,7 +56,7 @@ def test_io_import_rmi_nwp_xarray_attrs(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_x = [
-    ("units", "m", None),
+    ("cartesian_unit", "m", None),
     ("x1", 0, 0.1),
     ("x2", 731900.0, 0.1),
 ]
@@ -69,7 +69,7 @@ def test_io_import_rmi_nwp_xarray_attrs_coordx(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_y = [
-    ("units", "m", None),
+    ("cartesian_unit", "m", None),
     ("y1", -731900.0, 0.1),
     ("y2", 0.0, 0.1),
 ]

pysteps/tests/test_utils_reprojection.py

@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+
+import os
+import xarray as xr
+import pytest
+import pysteps
+from pysteps.utils import reprojection
+
+pytest.importorskip("rasterio")
+
+root_path_radar = pysteps.rcparams.data_sources["knmi"]["root_path"]
+root_path_nwp = pysteps.rcparams.data_sources["knmi_nwp"]["root_path"]
+
+rel_path_nwp = os.path.join("2018", "09", "05")
+rel_path_radar = os.path.join(
+    "2010", "08"
+)  # Different day, but that does not matter for the tester
+
+filename_nwp = os.path.join(
+    root_path_nwp, rel_path_nwp, "20180905_0600_Pforecast_Harmonie.nc"
+)
+filename_radar = os.path.join(
+    root_path_radar, rel_path_radar, "RAD_NL25_RAP_5min_201008260230.h5"
+)
+
+# Open the radar and NWP data
+radar_array_xr = pysteps.io.importers.import_knmi_hdf5(filename_radar)
+nwp_array_xr = pysteps.io.import_knmi_nwp_xr(filename_nwp)
+
+steps_arg_names = (
+    "radar_array_xr",
+    "nwp_array_xr",
+    "t",
+)
+
+steps_arg_values = [
+    (radar_array_xr, nwp_array_xr, 1),
+    (radar_array_xr, nwp_array_xr, 4),
+    (radar_array_xr, nwp_array_xr, 8),
+]
+
+
+@pytest.mark.parametrize(steps_arg_names, steps_arg_values)
+def test_utils_reprojection(
+    radar_array_xr,
+    nwp_array_xr,
+    t,
+):
+
+    # Reproject
+    nwp_array_xr_reproj = reprojection(nwp_array_xr, radar_array_xr)
+
+    # The tests
+    assert (
+        nwp_array_xr_reproj["t"].shape[0] == nwp_array_xr["t"].shape[0]
+    ), "Time dimension has not the same length as source"
+    assert (
+        nwp_array_xr_reproj["y"].shape[0] == radar_array_xr["y"].shape[0]
+    ), "y dimension has not the same length as radar composite"
+    assert (
+        nwp_array_xr_reproj["x"].shape[0] == radar_array_xr["x"].shape[0]
+    ), "x dimension has not the same length as radar composite"
+
+    assert (
+        nwp_array_xr_reproj.x.attrs["x1"] == radar_array_xr.x.attrs["x1"]
+    ), "x-value lower left corner is not equal to radar composite"
+    assert (
+        nwp_array_xr_reproj.x.attrs["x2"] == radar_array_xr.x.attrs["x2"]
+    ), "x-value upper right corner is not equal to radar composite"
+    assert (
+        nwp_array_xr_reproj.y.attrs["y1"] == radar_array_xr.y.attrs["y1"]
+    ), "y-value lower left corner is not equal to radar composite"
+    assert (
+        nwp_array_xr_reproj.y.attrs["y2"] == radar_array_xr.y.attrs["y2"]
+    ), "y-value upper right corner is not equal to radar composite"
+
+    assert (
+        nwp_array_xr_reproj.x.min().values == radar_array_xr.x.min().values
+    ), "First x-coordinate does not equal first x-coordinate radar composite"
+    assert (
+        nwp_array_xr_reproj.y.min().values == radar_array_xr.y.min().values
+    ), "First y-coordinate does not equal first y-coordinate radar composite"
+    assert (
+        nwp_array_xr_reproj.x.max().values == radar_array_xr.x.max().values
+    ), "Last x-coordinate does not equal last x-coordinate radar composite"
+    assert (
+        nwp_array_xr_reproj.y.max().values == radar_array_xr.y.max().values
+    ), "Last y-coordinate does not equal last y-coordinate radar composite"
+
+    assert (
+        nwp_array_xr_reproj.attrs["projection"] == radar_array_xr.attrs["projection"]
+    ), "projection is different than destionation projection"

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