FIX: implement more robust angle_res retrieval in extract_angle_param…

…eters (#118)

docs/history.md

@@ -5,7 +5,8 @@
 * FIX: Correct DB_DBTE8/DB_DBZE8 and DB_DBTE16/DB_DBZE16 decoding for iris-backend ({pull}`110`) by [@kmuehlbauer](https://github.com/kmuehlbauer)
 * FIX: Cast boolean string to int in rainbow dictionary ({pull}`113`) by [@egouden](https://github.com/egouden)
 * MNT: switch to mamba-org/setup-micromamba, split CI tests ({issue}`115`), ({pull}`116`) by [@kmuehlbauer](https://github.com/kmuehlbauer)
-* FIX: time interpolation
+* FIX: time interpolation ({pull}`117`) by [@kmuehlbauer](https://github.com/kmuehlbauer)
+* FIX: robust ``angle_res`` retrieval in ``extract_angle_parameters`` ({issue}`112`), ({pull}`118`) by [@kmuehlbauer](https://github.com/kmuehlbauer)
 
 ## 0.2.0 (2023-03-24)
 

tests/test_util.py

@@ -81,6 +81,104 @@ def test_extract_angle_parameters():
     assert angle_dict == expected_dict
 
 
+@pytest.mark.parametrize("low", [-0.25, -0.1, 0])
+@pytest.mark.parametrize("high", [0.25, 0.1, 0])
+def test_extract_angle_parameters2(low, high):
+    seed = 12345
+    rng = np.random.default_rng(seed)
+    noise = rng.uniform(low=low, high=high, size=360)
+    azimuth = np.arange(0.5, 360.0, 1.0) + noise
+    ds = model.create_sweep_dataset(a1gate=0, direction=1, azimuth=azimuth)
+    ds = ds.swap_dims({"time": "azimuth"})
+    ds = ds.sortby("azimuth")
+    ds = ds.assign_coords(sweep_mode="azimuth_surveillance")
+
+    # remove block of rays
+    ds_in = xr.concat(
+        [
+            ds.isel(azimuth=slice(None, 100)),
+            ds.isel(azimuth=slice(110, None)),
+        ],
+        "azimuth",
+        data_vars="minimal",
+    )
+    angle_dict = util.extract_angle_parameters(ds_in)
+    expected_dict = {
+        "first_angle": "azimuth",
+        "second_angle": "elevation",
+        "min_angle": np.array(0.36366801123358483),
+        "max_angle": np.array(359.6511639563352),
+        "min_time": np.datetime64("2022-08-27T10:00:00.000000000"),
+        "max_time": np.datetime64("2022-08-27T10:01:29.750000000"),
+        "angles_are_unique": True,
+        "times_are_unique": True,
+        "a1gate_idx": np.array(0),
+        "a1gate_val": np.array(0.36366801123358483),
+        "uniform_angle_spacing": False,
+        "ascending": np.array(True),
+        "direction": 1,
+        "angle_res": np.array(1.0),
+        "start_angle": 0,
+        "stop_angle": 360,
+        "expected_angle_span": 360,
+        "missing_rays": True,
+        "excess_rays": False,
+        "expected_number_rays": 360,
+    }
+    expected = {
+        (0.25, -0.25): {
+            "min_angle": np.array(0.36366801123358483),
+            "max_angle": np.array(359.6511639563352),
+            "a1gate_val": np.array(0.36366801123358483),
+        },
+        (0.25, -0.1): {
+            "min_angle": np.array(0.4795676078635094),
+            "max_angle": np.array(359.6808147694347),
+            "a1gate_val": np.array(0.4795676078635094),
+        },
+        (0.25, 0.0): {
+            "min_angle": np.array(0.5568340056167924),
+            "max_angle": np.array(359.70058197816763),
+            "a1gate_val": np.array(0.5568340056167924),
+        },
+        (0.1, -0.25): {
+            "min_angle": np.array(0.32956760786350936),
+            "max_angle": np.array(359.53081476943464),
+            "a1gate_val": np.array(0.32956760786350936),
+        },
+        (0.1, -0.1): {
+            "min_angle": np.array(0.4454672044934339),
+            "max_angle": np.array(359.5604655825341),
+            "a1gate_val": np.array(0.4454672044934339),
+        },
+        (0.1, 0.0): {
+            "min_angle": np.array(0.5227336022467169),
+            "max_angle": np.array(359.58023279126706),
+            "a1gate_val": np.array(0.5227336022467169),
+        },
+        (0.0, -0.25): {
+            "min_angle": np.array(0.3068340056167924),
+            "max_angle": np.array(359.45058197816763),
+            "a1gate_val": np.array(0.3068340056167924),
+        },
+        (0.0, -0.1): {
+            "min_angle": np.array(0.42273360224671697),
+            "max_angle": np.array(359.48023279126704),
+            "a1gate_val": np.array(0.42273360224671697),
+        },
+        (0.0, 0.0): {
+            "min_angle": np.array(0.5),
+            "max_angle": np.array(359.5),
+            "a1gate_val": np.array(0.5),
+        },
+    }
+    for k in angle_dict:
+        if k in ["min_angle", "max_angle", "a1gate_val"]:
+            assert angle_dict[k] == expected[(high, low)][k]
+            continue
+        assert angle_dict[k] == expected_dict[k]
+
+
 def test_ipol_time():
     filename = DATASETS.fetch("DWD-Vol-2_99999_20180601054047_00.h5")
     ds = xr.open_dataset(filename, group="sweep_7", engine="gamic", first_dim="auto")

xradar/util.py

@@ -273,12 +273,20 @@ def extract_angle_parameters(ds):
     # and other erroneous ray alignments which result in different first_angle values
 
     # 6. this finds the median of those diff-values which are over some quantile
-    # infact this removes angle differences which are too small
-    median_diff = diff.where(diff >= diff.quantile(0.30)).median(skipna=True)
+    # in fact this removes angle differences which are too small
+    # calculate std/median
+    # see GH112, if we get further issues, we might add a warning here in the future
+    std = diff.std(skipna=True)
+    median = diff.median(skipna=True)
+    # remove values above and below std (centered around median), calculate median
+    median_diff = diff.where(
+        (diff >= (median - std)) & (diff <= (median + std))
+    ).median(skipna=True)
+
     # unique differences
     diffset = set(diff.values)
     # if there are more than 1 unique angle differences that means
-    # non uniform angle spacing
+    # non-uniform angle spacing
     uniform_angle_spacing = len(diffset) == 1
     angle_dict.update(uniform_angle_spacing=uniform_angle_spacing)