From ff50359c0cbffffdfc282a0723b6a69b6d74f9db Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Wed, 22 Jul 2026 12:41:19 +0200 Subject: [PATCH 1/6] fix: extend knmi importer with support for other projections --- pysteps/io/importers.py | 41 ++++++++++++++++++++++++++++------------- 1 file changed, 28 insertions(+), 13 deletions(-) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index ad687c53..ca9ec986 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -86,19 +86,17 @@ import_dwd_radolan """ -import gzip -import os import array import datetime +import gzip +import os from functools import partial import numpy as np - from matplotlib.pyplot import imread from pysteps.decorators import postprocess_import -from pysteps.exceptions import DataModelError -from pysteps.exceptions import MissingOptionalDependency +from pysteps.exceptions import DataModelError, MissingOptionalDependency from pysteps.utils import aggregate_fields try: @@ -321,7 +319,7 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs): try: grib_file = pygrib.open(filename) except OSError: - raise OSError(f"Error opening NCEP's MRMS file. " f"File Not Found: {filename}") + raise OSError(f"Error opening NCEP's MRMS file. File Not Found: {filename}") if isinstance(window_size, int): window_size = (window_size, window_size) @@ -767,7 +765,7 @@ def import_knmi_hdf5( filename, qty="ACRR", accutime=5.0, - pixelsize=1000.0, + pixelsize=None, **kwargs, ): """ @@ -870,7 +868,7 @@ def import_knmi_hdf5( # The 'where' group of mch- and Opera-data, is called 'geographic' in the # KNMI data. geographic = f["geographic"] - proj4str = "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0" + proj4str = geographic["map_projection"].attrs["projection_proj4_params"].decode() pr = pyproj.Proj(proj4str) metadata["projection"] = proj4str @@ -899,10 +897,28 @@ def import_knmi_hdf5( metadata["y1"] = y1 metadata["x2"] = x2 metadata["y2"] = y2 - metadata["xpixelsize"] = pixelsize - metadata["ypixelsize"] = pixelsize - metadata["cartesian_unit"] = "m" - metadata["yorigin"] = "upper" + metadata["xpixelsize"] = ( + float(geographic.attrs["geo_pixel_size_x"][0]) + if pixelsize is None + else pixelsize + ) + metadata["ypixelsize"] = ( + float(geographic.attrs["geo_pixel_size_y"][0]) + if pixelsize is None + else pixelsize + ) + dim_pixel = geographic.attrs["geo_dim_pixel"].decode().split(",")[0] + if dim_pixel == "KM": + metadata["cartesian_unit"] = "km" + elif dim_pixel == "M": + metadata["cartesian_unit"] = "m" + elif dim_pixel == "DEG": + metadata["cartesian_unit"] = "degrees" + else: + metadata["cartesian_unit"] = "km" + metadata["yorigin"] = ( + "upper" if geographic.attrs["geo_pixel_def"].decode()[1] == "U" else "lower" + ) metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute" metadata["accutime"] = accutime metadata["unit"] = unit @@ -1925,7 +1941,6 @@ def _read_hdf5_cont(f, d): } else: - # Save h5py.Dataset by group name d[key] = np.array(value) From 4680eda81346ce18075a029624f7859666c0d05a Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Wed, 22 Jul 2026 13:41:28 +0200 Subject: [PATCH 2/6] fix: absolute ypixelsize and determine upper based on that --- pysteps/io/importers.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index ca9ec986..a163a520 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -902,11 +902,12 @@ def import_knmi_hdf5( if pixelsize is None else pixelsize ) - metadata["ypixelsize"] = ( + ypixelsize = ( float(geographic.attrs["geo_pixel_size_y"][0]) if pixelsize is None else pixelsize ) + metadata["ypixelsize"] = abs(ypixelsize) dim_pixel = geographic.attrs["geo_dim_pixel"].decode().split(",")[0] if dim_pixel == "KM": metadata["cartesian_unit"] = "km" @@ -916,9 +917,7 @@ def import_knmi_hdf5( metadata["cartesian_unit"] = "degrees" else: metadata["cartesian_unit"] = "km" - metadata["yorigin"] = ( - "upper" if geographic.attrs["geo_pixel_def"].decode()[1] == "U" else "lower" - ) + metadata["yorigin"] = "upper" if ypixelsize < 0 else "lower" metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute" metadata["accutime"] = accutime metadata["unit"] = unit From f8cfdc9b2dfdeeeeaa166f893062e0d2376f539f Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Wed, 22 Jul 2026 14:17:47 +0200 Subject: [PATCH 3/6] fix: support files with broken proj string --- pysteps/io/importers.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index a163a520..b8beb8fc 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -765,7 +765,6 @@ def import_knmi_hdf5( filename, qty="ACRR", accutime=5.0, - pixelsize=None, **kwargs, ): """ @@ -784,10 +783,6 @@ def import_knmi_hdf5( The accumulation time of the dataset in minutes. A 5 min accumulation is used as default, but hourly, daily and monthly accumulations are also available. - pixelsize: float - 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} @@ -869,6 +864,15 @@ def import_knmi_hdf5( # KNMI data. geographic = f["geographic"] proj4str = geographic["map_projection"].attrs["projection_proj4_params"].decode() + + # There are a bunch of knmi hdf5 files out there with incorrect projection string, fix those + fix_metadata = False + if ( + proj4str + == "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" + ): + fix_metadata = True + 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 @@ -898,18 +902,16 @@ def import_knmi_hdf5( metadata["x2"] = x2 metadata["y2"] = y2 metadata["xpixelsize"] = ( - float(geographic.attrs["geo_pixel_size_x"][0]) - if pixelsize is None - else pixelsize + 1000.0 if fix_metadata else float(geographic.attrs["geo_pixel_size_x"][0]) ) ypixelsize = ( - float(geographic.attrs["geo_pixel_size_y"][0]) - if pixelsize is None - else pixelsize + -1000.0 if fix_metadata else float(geographic.attrs["geo_pixel_size_y"][0]) ) metadata["ypixelsize"] = abs(ypixelsize) dim_pixel = geographic.attrs["geo_dim_pixel"].decode().split(",")[0] - if dim_pixel == "KM": + if fix_metadata: + metadata["cartesian_unit"] = "m" + elif dim_pixel == "KM": metadata["cartesian_unit"] = "km" elif dim_pixel == "M": metadata["cartesian_unit"] = "m" From 1c13c8a2afe6682ce8716c94d116f08cd88fbfd1 Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Wed, 22 Jul 2026 15:17:40 +0200 Subject: [PATCH 4/6] fix: test --- pysteps/tests/test_io_knmi_hdf5.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py index 3e30cb57..0fc674bc 100644 --- a/pysteps/tests/test_io_knmi_hdf5.py +++ b/pysteps/tests/test_io_knmi_hdf5.py @@ -22,20 +22,18 @@ 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=6378137 +b=6356752 +x_0=0 +y_0=0" -) +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" # list of (variable,expected,tolerance) tuples test_attrs = [ ("projection", expected_proj, None), ("x1", 0.0, 1e-10), - ("y1", -4415038.179210632, 1e-10), - ("x2", 699984.2646331593, 1e-10), - ("y2", -3649950.360247753, 1e-10), - ("xpixelsize", 1000.0, 1e-10), - ("xpixelsize", 1000.0, 1e-10), - ("cartesian_unit", "m", None), + ("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), ("accutime", 5.0, 1e-10), ("yorigin", "upper", None), ("unit", "mm", None), From 179c88116c09df6f422e5ece82fcdc678b70bfeb Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Mon, 27 Jul 2026 17:38:50 +0200 Subject: [PATCH 5/6] fix x1 x2 y1 y2 --- pysteps/io/importers.py | 42 ++++++++++++++++++----------------------- 1 file changed, 18 insertions(+), 24 deletions(-) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index b8beb8fc..e7d3bfbb 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -873,34 +873,28 @@ def import_knmi_hdf5( ): fix_metadata = True 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 - # Get coordinates - latlon_corners = geographic.attrs["geo_product_corners"] - ll_lat = latlon_corners[1] - ll_lon = latlon_corners[0] - ur_lat = latlon_corners[5] - ur_lon = latlon_corners[4] - lr_lat = latlon_corners[7] - lr_lon = latlon_corners[6] - ul_lat = latlon_corners[3] - ul_lon = latlon_corners[2] - - ll_x, ll_y = pr(ll_lon, ll_lat) - ur_x, ur_y = pr(ur_lon, ur_lat) - 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) + x1 = float(geographic.attrs["geo_column_offset"][0]) * float( + geographic.attrs["geo_pixel_size_x"][0] + ) + y1 = float(geographic.attrs["geo_row_offset"][0]) * float( + geographic.attrs["geo_pixel_size_y"][0] + ) + x2 = ( + float(geographic.attrs["geo_column_offset"][0]) + + float(geographic.attrs["geo_number_columns"][0]) + ) * float(geographic.attrs["geo_pixel_size_x"][0]) + y2 = ( + float(geographic.attrs["geo_row_offset"][0]) + + float(geographic.attrs["geo_number_rows"][0]) + ) * float(geographic.attrs["geo_pixel_size_y"][0]) # Fill in the metadata - metadata["x1"] = x1 - metadata["y1"] = y1 - metadata["x2"] = x2 - metadata["y2"] = y2 + metadata["x1"] = 0.0 if fix_metadata else x1 + metadata["y1"] = -4415000.0 if fix_metadata else y1 + metadata["x2"] = 700000.0 if fix_metadata else x2 + metadata["y2"] = -3650000.0 if fix_metadata else y2 metadata["xpixelsize"] = ( 1000.0 if fix_metadata else float(geographic.attrs["geo_pixel_size_x"][0]) ) From 7b14d56cc60f1f193ff832c9e2649f4b102530e9 Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Tue, 28 Jul 2026 11:02:34 +0200 Subject: [PATCH 6/6] fix: tests --- pysteps/io/importers.py | 10 +++++++--- pysteps/tests/test_io_knmi_hdf5.py | 6 +++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index e7d3bfbb..d51cf9a5 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -889,6 +889,13 @@ def import_knmi_hdf5( float(geographic.attrs["geo_row_offset"][0]) + float(geographic.attrs["geo_number_rows"][0]) ) * float(geographic.attrs["geo_pixel_size_y"][0]) + ypixelsize = ( + -1000.0 if fix_metadata else float(geographic.attrs["geo_pixel_size_y"][0]) + ) + if ypixelsize < 0: + y_temp = y1 + y1 = y2 + y2 = y_temp # Fill in the metadata metadata["x1"] = 0.0 if fix_metadata else x1 @@ -898,9 +905,6 @@ def import_knmi_hdf5( metadata["xpixelsize"] = ( 1000.0 if fix_metadata else float(geographic.attrs["geo_pixel_size_x"][0]) ) - ypixelsize = ( - -1000.0 if fix_metadata else float(geographic.attrs["geo_pixel_size_y"][0]) - ) metadata["ypixelsize"] = abs(ypixelsize) dim_pixel = geographic.attrs["geo_dim_pixel"].decode().split(",")[0] if fix_metadata: diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py index 0fc674bc..a8263b0f 100644 --- a/pysteps/tests/test_io_knmi_hdf5.py +++ b/pysteps/tests/test_io_knmi_hdf5.py @@ -28,9 +28,9 @@ def test_io_import_knmi_hdf5_shape(): test_attrs = [ ("projection", expected_proj, None), ("x1", 0.0, 1e-10), - ("y1", -4415.038179210632, 1e-10), - ("x2", 699.9842646331593, 1e-10), - ("y2", -3649.950360247753, 1e-10), + ("y1", -4415.0, 1e-10), + ("x2", 700.0, 1e-10), + ("y2", -3650.0, 1e-10), ("xpixelsize", 1.0, 1e-10), ("xpixelsize", 1.0, 1e-10), ("cartesian_unit", "km", None),