fix in odim h5 reader and writer

pyart/aux_io/odim_h5.py

@@ -452,7 +452,6 @@ def read_odim_grid_h5(filename, field_names=None, additional_metadata=None,
         for odim_field, h_field_key, dset in zip(
                 odim_fields, h_field_keys, dsets):
             field_name = filemetadata.get_field_name(_to_str(odim_field))
-
             if field_name is None:
                 # warn(f'field {odim_field} not available in {filename}')
                 continue
@@ -679,9 +678,9 @@ def read_odim_h5(filename, field_names=None, additional_metadata=None,
         altitude = filemetadata('altitude')
 
         h_where = hfile['where'].attrs
-        latitude['data'] = np.array([h_where['lat']], dtype='float64')
-        longitude['data'] = np.array([h_where['lon']], dtype='float64')
-        altitude['data'] = np.array([h_where['height']], dtype='float64')
+        latitude['data'] = np.array([h_where['lat']], dtype='float64').flatten()
+        longitude['data'] = np.array([h_where['lon']], dtype='float64').flatten()
+        altitude['data'] = np.array([h_where['height']], dtype='float64').flatten()
 
         # metadata
         metadata = filemetadata('metadata')
@@ -741,7 +740,7 @@ def read_odim_h5(filename, field_names=None, additional_metadata=None,
             sweep_el = [hfile[d]['where'].attrs['az_angle'] for d in datasets]
         else:
             sweep_el = [hfile[d]['where'].attrs['elangle'] for d in datasets]
-        fixed_angle['data'] = np.array(sweep_el, dtype='float32')
+        fixed_angle['data'] = np.array(sweep_el, dtype='float32').flatten()
 
         # elevation
         elevation = filemetadata('elevation')
@@ -882,6 +881,8 @@ def read_odim_h5(filename, field_names=None, additional_metadata=None,
         # fields
         fields = {}
         h_field_keys = [k for k in hfile['dataset1'] if k.startswith('data')]
+        # reorder field names to match correct order 1 to N
+        h_field_keys.sort(key=lambda x: int(x[4:]))
         odim_fields = [hfile['dataset1'][d]['what'].attrs['quantity'] for d in
                        h_field_keys]
         for odim_field, h_field_key in zip(odim_fields, h_field_keys):
@@ -1012,9 +1013,9 @@ def read_odim_vp_h5(filename, field_names=None, additional_metadata=None,
         altitude = filemetadata('altitude')
 
         h_where = hfile['where'].attrs
-        latitude['data'] = np.array([h_where['lat']], dtype='float64')
-        longitude['data'] = np.array([h_where['lon']], dtype='float64')
-        altitude['data'] = np.array([h_where['height']], dtype='float64')
+        latitude['data'] = np.array([h_where['lat']], dtype='float64').flatten()
+        longitude['data'] = np.array([h_where['lon']], dtype='float64').flatten()
+        altitude['data'] = np.array([h_where['height']], dtype='float64').flatten()
 
         # metadata
         metadata = filemetadata('metadata')

pyart/aux_io/odim_h5_writer.py

@@ -413,6 +413,16 @@ def write_odim_h5(filename, radar, field_names=None, physical=True,
     # Initialize hdf5 file
     hdf_id = _create_odim_h5_file(filename)
 
+    # Check existence of ray_angle_res
+    if not radar.ray_angle_res:
+        radar.ray_angle_res = {'data': []}
+        for i in range(radar.nsweeps):
+            start = radar.sweep_start_ray_index['data'][i]
+            end = radar.sweep_end_ray_index['data'][i]
+            azdiff = (radar.azimuth['data'][start+1:end] -
+                        radar.azimuth['data'][start:end-1]).mean()
+            radar.ray_angle_res['data'].append(np.around(azdiff, 1))
+
     # Determine odim object, number of datasets (PPIs, RHIs) and
     # check if sweep mode does not change between different sweeps
     if radar.scan_type == 'rhi':
@@ -427,7 +437,9 @@ def write_odim_h5(filename, radar, field_names=None, physical=True,
                 if radar.sweep_mode is not None:
                     if len(set(radar.sweep_mode['data'])) <= 1:
                         nray_in_sweep = radar.rays_per_sweep['data'][0]
-                        ray_angle_res = radar.ray_angle_res['data'][0]
+                        if radar.ray_angle_res:
+                            ray_angle_res = radar.ray_angle_res['data'][0]
+                        print(ray_angle_res)
                         if nray_in_sweep * ray_angle_res < 360:
                             odim_object = 'AZIM'
                         else:
@@ -457,6 +469,7 @@ def write_odim_h5(filename, radar, field_names=None, physical=True,
 
     else:
         field_names = list(radar.fields.keys())
+
     n_datatypes = len(field_names)
 
     # Create level 1 group structure
@@ -593,8 +606,12 @@ def write_odim_h5(filename, radar, field_names=None, physical=True,
         where2_dict = {}
 
         where2_dict['nbins'] = np.int64(np.repeat(radar.ngates, n_datasets))
-        where2_dict['rstart'] = np.repeat(
-            np.double((radar.range['data'][0]) / (1000.)), n_datasets)  # [km]
+        if 'meters_to_center_of_first_gate' in radar.range:
+            rstart = radar.range['meters_to_center_of_first_gate']
+        else:
+            rstart = 1.5 * radar.range['data'][0] - 0.5 * radar.range['data'][1]
+
+        where2_dict['rstart'] = np.repeat(np.double(rstart / 1000.), n_datasets)  # [km]
         where2_dict['nrays'] = np.int64(radar.rays_per_sweep['data'])
 
         if len(set(radar.range['data'][1:] - radar.range['data'][0:-1])) <= 1:

pyart/testing/sample_files.py

@@ -36,6 +36,7 @@
 CFRADIAL_PPI_FILE = os.path.join(DATA_PATH, 'example_cfradial_ppi.nc')
 CFRADIAL_RHI_FILE = os.path.join(DATA_PATH, 'example_cfradial_rhi.nc')
 CFRADIAL_CR_RASTER_FILE = os.path.join(DATA_PATH, "example_cfradial_cr_raster.nc")
+ODIM_H5_FILE = os.path.join(DATA_PATH, "example_radar.polar.fikor.h5")
 CHL_RHI_FILE = os.path.join(DATA_PATH, 'example_chl_rhi.chl')
 SIGMET_PPI_FILE = os.path.join(DATA_PATH, 'example_sigmet_ppi.sigmet')
 SIGMET_RHI_FILE = os.path.join(DATA_PATH, 'example_sigmet_rhi.sigmet')

tests/io/test_odim_h5.py

@@ -0,0 +1,93 @@
+""" Unit Tests for Py-ART's io/cfradial.py module. """
+
+import h5py
+import numpy 
+from numpy.testing import assert_almost_equal
+
+import pyart
+
+#################################################
+# read_cfradial tests (verify radar attributes) #
+#################################################
+
+# read in the sample file and create a a Radar object
+radar = pyart.aux_io.read_odim_h5(pyart.testing.ODIM_H5_FILE)
+
+
+# time attribute
+def test_time():
+    assert "comment" in radar.time.keys()
+    assert "long_name" in radar.time.keys()
+    assert "standard_name" in radar.time.keys()
+    assert "units" in radar.time.keys()
+    assert "calendar" in radar.time.keys()
+    assert "data" in radar.time.keys()
+    assert radar.time["units"] == "seconds since 2023-08-07T16:10:08Z"
+    assert radar.time["data"].shape == (720,)
+    assert_almost_equal(radar.time["data"][38], 2, 0)
+
+
+# range attribute
+def test_range():
+    assert "long_name" in radar.range
+    assert "standard_name" in radar.range
+    assert "meters_to_center_of_first_gate" in radar.range
+    assert "meters_between_gates" in radar.range
+    assert "units" in radar.range
+    assert "data" in radar.range
+    assert "spacing_is_constant" in radar.range
+    assert radar.range["data"].shape == (500,)
+    assert_almost_equal(radar.range["data"][0], 250, 0)
+
+
+# metadata attribute
+def test_metadata():
+    assert "Conventions" in radar.metadata
+    assert "source" in radar.metadata
+    assert "version" in radar.metadata
+    assert "odim_conventions" in radar.metadata
+
+# scan_type attribute
+def test_scan_type():
+    assert radar.scan_type == "ppi"
+
+########################################################################
+# write_odim_h5 tests (verify data in written netCDF matches original #
+########################################################################
+
+def assert_hdf5_files_equal(file1, file2):
+    # Open the HDF5 files
+
+    def check_group(group1, group2):
+        # Check attributes of the group
+        assert group1.attrs.keys() == group2.attrs.keys()
+        for attr_name in group1.attrs:
+            if type(group1.attrs[attr_name]) in [list, numpy.ndarray]:
+                    assert_almost_equal(group1.attrs[attr_name], 
+                        group1.attrs[attr_name], 3)
+            else:
+                assert (group1.attrs[attr_name] == 
+                        group2.attrs[attr_name])
+
+        # Check variables in the group
+        assert set(group1.keys()) == set(group2.keys())
+        for var_name in group1.keys():
+            # Check attributes of the variable
+            if var_name == 'data':
+                assert_almost_equal(group1[var_name][:], group2[var_name][:])
+            else:
+                check_group(group1[var_name], group2[var_name])
+
+    # Recursively check groups
+    check_group(file1, file2)
+
+def test_write_ppi():
+    # CF/Radial example file -> Radar object -> netCDF file
+    with pyart.testing.InTemporaryDirectory():
+        tmpfile = "tmp_ppi.h5"
+        radar = pyart.aux_io.read_odim_h5(pyart.testing.ODIM_H5_FILE)
+        pyart.aux_io.write_odim_h5(tmpfile, radar)
+        ref = h5py.File(pyart.testing.ODIM_H5_FILE)
+        dset = h5py.File(tmpfile)
+        assert_hdf5_files_equal(dset, ref)
+        dset.close()