Repository: openradar/xradar

Issue:
possible bug with angle reindexing (?)
* xradar version: 0.1.0
* Python version: 3.9
* Operating System: Ubuntu 22.04

Hi everyone, 

I am trying to apply the angle reindexing function to my dataset as is explained [here](https://docs.openradarscience.org/projects/xradar/en/stable/notebooks/angle_reindexing.html). the radar data I am analyzing is 360 degrees in azimuth as follows

``` python 
import fsspec
import xradar as xd
import xarray as xr
from pandas import to_datetime
from datetime import datetime

def create_query(date, radar_site):
    """
    Creates a string for quering the IDEAM radar files stored in AWS bucket
    :param date: date to be queried. e.g datetime(2021, 10, 3, 12). Datetime python object
    :param radar_site: radar site e.g. Guaviare
    :return: string with a IDEAM radar bucket format
    """
    if (date.hour != 0) and (date.hour != 0):
        return f'l2_data/{date:%Y}/{date:%m}/{date:%d}/{radar_site}/{radar_site[:3].upper()}{date:%y%m%d%H}'
    elif (date.hour != 0) and (date.hour == 0):
        return f'l2_data/{date:%Y}/{date:%m}/{date:%d}/{radar_site}/{radar_site[:3].upper()}{date:%y%m%d}'
    else:
        return f'l2_data/{date:%Y}/{date:%m}/{date:%d}/{radar_site}/{radar_site[:3].upper()}{date:%y%m%d}'
    
    
def data_accessor(file):
    """
    Open AWS S3 file(s), which can be resolved locally by file caching
    """
    return fsspec.open_local(f'simplecache::s3://{file}', s3={'anon': True}, filecache={'cache_storage': '/tmp/radar/'})

date_query = datetime(2023, 4, 7, 3)
radar_name = "Barrancabermeja" 
query = create_query(date=date_query, radar_site=radar_name)
str_bucket = 's3://s3-radaresideam/'
fs = fsspec.filesystem("s3", anon=True)
radar_files = sorted(fs.glob(f"{str_bucket}{query}*"))
radar_dt = xd.io.open_iris_datatree(data_accessor(radar_files[1])).xradar.georeference()
print(radar_dt.sweep_0)
DataTree('sweep_0', parent="root")
    Dimensions:            (azimuth: 360, range: 747)
    Coordinates:
      * azimuth            (azimuth) float64 0.05493 1.057 2.038 ... 358.1 359.0
        elevation          (azimuth) float64 1.472 1.472 1.472 ... 1.472 1.472 1.472
        time               (azimuth) datetime64
```

The datatree object is `360` angles in azimuth. Then, I try to reindex the azimuth angles using the `fix_angle` function

```python 
def fix_angle(ds):
    angle_dict = xd.util.extract_angle_parameters(ds)
    display(angle_dict)
    start_ang = angle_dict["start_angle"]
    stop_ang = angle_dict["stop_angle"]
    angle_res = angle_dict["angle_res"]
    direction = angle_dict["direction"]

    # first find exact duplicates and remove
    ds = xd.util.remove_duplicate_rays(ds)

    # second reindex according to retrieved parameters
    ds = xd.util.reindex_angle(
        ds, start_ang, stop_ang, angle_res, direction, method="nearest"
    )
    return ds
```
 The result looks like this, 

```python

ds2 = fix_angle(radar_dt.sweep_0)

{'first_angle': 'azimuth',
 'second_angle': 'elevation',
 'min_angle': array(0.05493164),
 'max_angle': array(359.04968262),
 'min_time': numpy.datetime64('2023-04-07T03:01:08.569000000'),
 'max_time': numpy.datetime64('2023-04-07T03:01:28.569000000'),
 'angles_are_unique': True,
 'times_are_unique': False,
 'a1gate_idx': array(247),
 'a1gate_val': array(247.03308105),
 'uniform_angle_spacing': False,
 'ascending': array(True),
 'direction': 1,
 'angle_res': array(1.01),
 'start_angle': 0,
 'stop_angle': 360,
 'expected_angle_span': 360,
 'missing_rays': array(False),
 'excess_rays': array(True),
 'expected_number_rays': 356}
```
Now our new xradar object has 356 angles in the azimuth dimension

```python
DataTree('sweep_0', parent=None)
    Dimensions:            (azimuth: 356, range: 747)
    Coordinates:
      * azimuth            (azimuth) float64 0.505 1.515 2.525 ... 357.0 358.0 359.1
      * range              (range) float32 1e+03 1.3e+03 ... 2.245e+05 2.248e+05
        elevation          (azimuth) float64 1.472 1.472 1.472 ... 1.472 1.472 1.472
        time               (azimuth) datetime64[ns] 2023-04-07T03:01:14.569000 .....
        longitude          float64 -73.76
        latitude           float64 6.933
        altitude           float64 105.0
```

Digging into [util.py](https://github.com/openradar/xradar/blob/main/xradar/util.py) seems like in the line [#L291](https://github.com/openradar/xradar/blob/6fe2ea094e0b68dc9817d3276ec2673bc17bc6f4/xradar/util.py#L291)  the angle resolution is affected by the number of decimals of the  `np.round` method. 

```python 
'angle_res': array(1.01)
```

If we change the number of decimals to 1 the output will look like this:

```python 
{'first_angle': 'azimuth',
 'second_angle': 'elevation',
 'min_angle': array(0.05493164),
 'max_angle': array(359.04968262),
 'min_time': numpy.datetime64('2023-04-07T03:01:08.569000000'),
 'max_time': numpy.datetime64('2023-04-07T03:01:28.569000000'),
 'angles_are_unique': True,
 'times_are_unique': False,
 'a1gate_idx': array(247),
 'a1gate_val': array(247.03308105),
 'uniform_angle_spacing': False,
 'ascending': array(True),
 'direction': 1,
 'angle_res': array(1.),
 'start_angle': 0,
 'stop_angle': 360,
 'expected_angle_span': 360,
 'missing_rays': array(False),
 'excess_rays': array(False),
 'expected_number_rays': 360}

DataTree('sweep_0', parent=None)
    Dimensions:            (azimuth: 360, range: 747)
    Coordinates:
      * azimuth            (azimuth) float64 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
      * range              (range) float32 1e+03 1.3e+03 ... 2.245e+05 2.248e+05
        elevation          (azimuth) float64
``` 

Please let me know your thoughts. 

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/openradar/xradar/pull/118?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=openradar) Report
> Merging [#118](https://app.codecov.io/gh/openradar/xradar/pull/118?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=openradar) (150ee85) into [main](https://app.codecov.io/gh/openradar/xradar/commit/e464ac7d8de09bd70c2c7b719175c4c93648feb7?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=openradar) (e464ac7) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

> :exclamation: Current head 150ee85 differs from pull request most recent head 4660b50. Consider uploading reports for the commit 4660b50 to get more accurate results

```diff
@@           Coverage Diff           @@
##             main     #118   +/-   ##
=======================================
  Coverage   88.17%   88.18%           
=======================================
  Files          19       19           
  Lines        3333     3335    +2     
=======================================
+ Hits         2939     2941    +2     
  Misses        394      394           
```

| Flag | Coverage Δ | |
|---|---|---|
| unittests | `88.18% <100.00%> (+<0.01%)` | :arrow_up: |

Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=openradar#carryforward-flags-in-the-pull-request-comment) to find out more.

| [Impacted Files](https://app.codecov.io/gh/openradar/xradar/pull/118?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=openradar) | Coverage Δ | |
|---|---|---|
| [xradar/util.py](https://app.codecov.io/gh/openradar/xradar/pull/118?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=openradar#diff-eHJhZGFyL3V0aWwucHk=) | `97.31% <100.00%> (+0.03%)` | :arrow_up: |



Patch:
diff --git a/docs/history.md b/docs/history.md
index 1abd0ad..b0456f8 100644
--- a/docs/history.md
+++ b/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)
 
diff --git a/xradar/util.py b/xradar/util.py
index 6ae1aec..c4441a1 100644
--- a/xradar/util.py
+++ b/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)
 


Test patch:
diff --git a/tests/test_util.py b/tests/test_util.py
index e10918e..f740dfe 100644
--- a/tests/test_util.py
+++ b/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")