diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index ad687c53..d51cf9a5 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,6 @@ def import_knmi_hdf5( filename, qty="ACRR", accutime=5.0, - pixelsize=1000.0, **kwargs, ): """ @@ -786,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} @@ -870,39 +863,61 @@ 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" - pr = pyproj.Proj(proj4str) - metadata["projection"] = proj4str + proj4str = geographic["map_projection"].attrs["projection_proj4_params"].decode() - # 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] + # 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" + metadata["projection"] = proj4str - 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]) + 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"] = x1 - metadata["y1"] = y1 - metadata["x2"] = x2 - metadata["y2"] = y2 - metadata["xpixelsize"] = pixelsize - metadata["ypixelsize"] = pixelsize - metadata["cartesian_unit"] = "m" - metadata["yorigin"] = "upper" + 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]) + ) + metadata["ypixelsize"] = abs(ypixelsize) + dim_pixel = geographic.attrs["geo_dim_pixel"].decode().split(",")[0] + if fix_metadata: + metadata["cartesian_unit"] = "m" + elif 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 ypixelsize < 0 else "lower" metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute" metadata["accutime"] = accutime metadata["unit"] = unit @@ -1925,7 +1940,6 @@ def _read_hdf5_cont(f, d): } else: - # Save h5py.Dataset by group name d[key] = np.array(value) diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py index 3e30cb57..a8263b0f 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.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), ("accutime", 5.0, 1e-10), ("yorigin", "upper", None), ("unit", "mm", None),