"""
pyart.aux_io.knmi_h5
====================
Routines for reading ODIM_H5 files.
.. autosummary::
:toctree: generated/
read_knmi_grid_h5
_get_knmi_data
"""
import datetime
import numpy as np
try:
import h5py
_H5PY_AVAILABLE = True
except ImportError:
_H5PY_AVAILABLE = False
try:
import pyproj
_PYPROJ_AVAILABLE = True
except ImportError:
_PYPROJ_AVAILABLE = False
from ..aux_io.odim_h5 import _to_str, proj4_to_dict
from ..config import FileMetadata, get_fillvalue
from ..core.grid import Grid
from ..exceptions import MissingOptionalDependency
from ..io.common import _test_arguments, make_time_unit_str
from ..util import ma_broadcast_to
KNMI_H5_FIELD_NAMES = {
"RAINFALL_RATE_[MM/H]": "radar_estimated_rain_rate",
"ACCUMULATION_[MM]": "rainfall_accumulation",
"QUALITY_[-]": "signal_quality_index",
"ADJUSTMENT_FACTOR_[DB]": "adjustment_factor",
}
[docs]def read_knmi_grid_h5(
filename,
field_names=None,
additional_metadata=None,
file_field_names=False,
exclude_fields=None,
include_fields=None,
time_ref="end",
**kwargs
):
"""
Read a KNMI_H5 grid file.
Parameters
----------
filename : str
Name of the field_name file to read.
field_names : dict, optional
Dictionary mapping field_name field names to radar field names. If a
data type found in the file does not appear in this dictionary or has
a value of None it will not be placed in the radar.fields dictionary.
A value of None, the default, will use the mapping defined in the
Py-ART configuration file.
additional_metadata : dict of dicts, optional
Dictionary of dictionaries to retrieve metadata from during this read.
This metadata is not used during any successive file reads unless
explicitly included. A value of None, the default, will not
introduct any addition metadata and the file specific or default
metadata as specified by the Py-ART configuration file will be used.
file_field_names : bool, optional
True to use the MDV data type names for the field names. If this
case the field_names parameter is ignored. The field dictionary will
likely only have a 'data' key, unless the fields are defined in
`additional_metadata`.
exclude_fields : list or None, optional
List of fields to exclude from the radar object. This is applied
after the `file_field_names` and `field_names` parameters. Set
to None to include all fields specified by include_fields.
include_fields : list or None, optional
List of fields to include from the radar object. This is applied
after the `file_field_names` and `field_names` parameters. Set
to None to include all fields not specified by exclude_fields.
time_ref : str
Time reference in the /what/time attribute. Can be either 'start',
'mid' or 'end'. If 'start' the attribute is expected to be the
starttime of the scan, if 'mid', the middle time, if 'end' the
endtime.
Returns
-------
grid : Grid
Grid object containing data from ODIM_H5 file.
"""
# check that h5py is available
if not _H5PY_AVAILABLE:
raise MissingOptionalDependency(
"h5py is required to use read_odim_h5 but is not installed"
)
# test for non empty kwargs
_test_arguments(kwargs)
# create metadata retrieval object
if field_names is None:
field_names = KNMI_H5_FIELD_NAMES
filemetadata = FileMetadata(
"odim_h5",
field_names,
additional_metadata,
file_field_names,
exclude_fields,
include_fields,
)
with h5py.File(filename, "r") as hfile:
# determine the number of datasets by the number of groups which
# begin with image
datasets = [k for k in hfile if k.startswith("image")]
datasets.sort(key=lambda x: int(x[5:]))
# latitude, longitude and altitude
x = filemetadata("x")
y = filemetadata("y")
z = filemetadata("z")
geo_group = hfile["geographic"].attrs
map_proj = hfile["geographic"]["map_projection"].attrs
# projection definition
projection = map_proj["projection_proj4_params"]
# 0, 49.3621, 0, 55.9736, 10.8565, 55.389, 9.0093, 48.8953
geo_product_corners = geo_group["geo_product_corners"]
if _PYPROJ_AVAILABLE:
wgs84 = pyproj.CRS.from_epsg(4326)
try: # pyproj doens't like bytearrays
projection = projection.decode("utf-8")
except Exception:
pass
# projection = f'{projection} +units=km'
projection = (
"+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0"
" +a=6378137 +b=6356752 +x_0=0 +y_0=0"
)
coordTrans = pyproj.Transformer.from_crs(wgs84, projection)
xstart, ystart = coordTrans.transform(
geo_product_corners[1], geo_product_corners[0]
)
xend, yend = coordTrans.transform(
geo_product_corners[5], geo_product_corners[4]
)
xvec = np.linspace(xstart, xend, geo_group["geo_number_columns"][0])
yvec = np.linspace(ystart, yend, geo_group["geo_number_rows"][0])
else:
xvec = np.linspace(
geo_product_corners[5],
geo_product_corners[1],
geo_group["geo_number_columns"][0],
)
yvec = np.linspace(
geo_product_corners[0],
geo_product_corners[4],
geo_group["geo_number_rows"][0],
)
x["data"] = xvec
y["data"] = yvec
z["data"] = np.array([0], dtype="float64")
# metadata
metadata = filemetadata("metadata")
# metadata['source'] = _to_str(hfile['what'].attrs['source'])
metadata["original_container"] = "knmi_h5"
# metadata['odim_conventions'] = _to_str(hfile.attrs['Conventions'])
# h_what = hfile['what'].attrs
# metadata['version'] = _to_str(h_what['version'])
# metadata['source'] = _to_str(h_what['source'])
# Get the MeteoSwiss-specific data
# try:
# h_how2 = hfile['how']['MeteoSwiss'].attrs
# except KeyError:
# # if no how group exists mock it with an empty dictionary
# h_how2 = {}
# if 'radar' in h_how2:
# metadata['radar'] = h_how2['radar']
# try:
# ds1_how = hfile[datasets[0]]['how'].attrs
# except KeyError:
# # if no how group exists mock it with an empty dictionary
# ds1_how = {}
# if 'system' in ds1_how:
# metadata['system'] = ds1_how['system']
# if 'software' in ds1_how:
# metadata['software'] = ds1_how['software']
# if 'sw_version' in ds1_how:
# metadata['sw_version'] = ds1_how['sw_version']
dset_list = []
field_list = []
knmi_list = []
for knmi_field in field_names.keys():
for dset in datasets:
dset_field = _to_str(hfile[dset].attrs["image_geo_parameter"])
if knmi_field == dset_field:
dset_list.append(dset)
field_list.append(field_names[knmi_field])
knmi_list.append(knmi_field)
break
if (
knmi_field == "RAINFALL_RATE_[MM/H]"
and dset_field == "ACCUMULATION_[MM]"
):
dset_list.append(dset)
field_list.append(field_names[knmi_field])
knmi_list.append(knmi_field)
break
fields = {}
for knmi_field, field_name, dset in zip(knmi_list, field_list, dset_list):
fdata = _get_knmi_data(
hfile[dset], knmi_field, hfile[dset]["calibration"].attrs
)
field_dic = filemetadata(field_name)
# grid object expects a 3D field
ny = geo_group["geo_number_rows"][0]
nx = geo_group["geo_number_columns"][0]
field_dic["data"] = ma_broadcast_to(fdata[::-1, :], (1, ny, nx))
field_dic["_FillValue"] = get_fillvalue()
fields[field_name] = field_dic
if not fields:
# warn(f'No fields could be retrieved from file')
return None
_time = filemetadata("time")
origin_latitude = filemetadata("origin_latitude")
origin_longitude = filemetadata("origin_longitude")
origin_altitude = filemetadata("origin_altitude")
# format 12-OCT-2023;08:00:04.000
start_time = hfile["overview"].attrs["product_datetime_start"]
start_time = datetime.datetime.strptime(
_to_str(start_time), "%d-%b-%Y;%H:%M:%S.000"
)
end_time = hfile["overview"].attrs["product_datetime_end"]
end_time = datetime.datetime.strptime(
_to_str(end_time), "%d-%b-%Y;%H:%M:%S.000"
)
_time["data"] = [0]
if time_ref == "mid":
mid_delta = (end_time - start_time) / 2
mid_ts = start_time + mid_delta
_time["units"] = make_time_unit_str(mid_ts)
elif time_ref == "end":
_time["units"] = make_time_unit_str(end_time)
else:
_time["units"] = make_time_unit_str(start_time)
projection = proj4_to_dict(projection)
if "lat_0" in projection:
origin_latitude["data"] = np.array([projection["lat_0"]])
else:
origin_latitude["data"] = np.array([0.0])
if "lon_0" in projection:
origin_longitude["data"] = np.array([projection["lat_0"]])
else:
origin_longitude["data"] = np.array([0.0])
origin_altitude["data"] = np.array([0.0])
# radar variables
radar_latitude = None
radar_longitude = None
radar_altitude = None
radar_name = None
radar_time = None
return Grid(
_time,
fields,
metadata,
origin_latitude,
origin_longitude,
origin_altitude,
x,
y,
z,
projection=projection,
radar_latitude=radar_latitude,
radar_longitude=radar_longitude,
radar_altitude=radar_altitude,
radar_name=radar_name,
radar_time=radar_time,
)
def _get_knmi_data(group, knmi_field, calibration):
"""Get KNMI data from an HDF5 group."""
nodata = calibration["calibration_missing_data"]
undetect = calibration["calibration_out_of_image"]
formula = _to_str(calibration["calibration_formulas"]).split("=")[1]
# b'GEO = 0.500000 * PV + -32.000000'
gain = float(formula.split("*")[0])
offset = float(formula.split("+")[1])
# mask raw data
raw_data = group["image_data"][:]
data = np.ma.masked_values(raw_data, nodata)
data = np.ma.masked_values(data, undetect)
data = data * gain + offset
if knmi_field == "RAINFALL_RATE_[MM/H]":
data *= 12.0
return data