best way to estimate tolerace parameter in angle reindexing?

[aladinor]

• xradar version: 0.1.0
• Python version: 3.9
• Operating System: Ubuntu 22.04

Hi everyone,

I am working with angle reindexing, and after doing some tests, I noticed that I got empty rays.

import fsspec
import xradar as xd
import matplotlib.pyplot as plt
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 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

def main():
    date_query = datetime(2023, 4, 7)
    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}*"))
    task_files = [fsspec.open_local(f'simplecache::s3://{i}', s3={'anon': True}, filecache={'cache_storage': '.'})
                 for i in radar_files[0:2]]
    radar_raw = xd.io.open_iris_datatree(task_files[1])
    radar_new = fix_angle(radar_raw['sweep_0'].ds)

    fig, (ax, ax1) = plt.subplots(1, 2, figsize=(9, 4), sharey=True)
    sc = radar_raw['sweep_0'].DBZH.sel(azimuth=slice(30, 60)).plot(vmin=-10, vmax=50, cmap="Spectral_r", ax=ax,
                                                                   add_colorbar=False)
    radar_new.sel(azimuth=slice(30, 60)).DBZH.plot(vmin=-10, vmax=50, cmap="Spectral_r", ax=ax1, add_colorbar=False)

    m2km = lambda x, _: f"{x / 1000:g}"
    fig.colorbar(sc, ax=[ax, ax1], label='Reflectivity [dBZ]')
    # set new ticks
    ax.xaxis.set_major_formatter(m2km)
    ax.set_xlabel("Range [Km]")
    ax.set_ylabel("Azimuth from true north [Deg]")

    ax.set_title("Raw data")
    ax1.xaxis.set_major_formatter(m2km)
    ax1.set_xlabel("Range [Km]")
    ax1.set_title("Reindexed data")
    _ = ax1.set_ylabel("")
    plt.show()

The output is this

Diggin around, I realized that in xd.util.reindex_angle method, we can specify the tolerance that by default is angle_res / 2.0. For this specific example, the default tolerance is 0.25. If I change the tolerance to 0.26 it will handle the empty rays as follows.

    ds = xd.util.reindex_angle(
        ds, start_ang, stop_ang, angle_res, direction, method="nearest", 
        tolerance=0.26
    )

Please let me know your comments.

Alfonso

[kmuehlbauer]

@aladinor Empty rays might happen if the original data is further away than the tolerance.

Should we make the tolerance more prominent in the reindex_angle-docstring or examples?

[aladinor]

@kmuehlbauer The way that I think we can handle this is to make tolerance more prominent, as you mentioned. In utils.py #L189 tolerance is defined as resolution / 2.0. I would make it resolution / 1.75 either by default value or by pointing it in the docstring or documentation.

This is how I solved it without any source code modification

def fix_angle(ds, tol=None):
    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"]
    tol = angle_res / 1.75   ### <--------------------- here
    # 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", 
        tolerance=tol
    )
    return ds

[kmuehlbauer]

@aladinor I'm hesitant to change it to values > angle_res/2. From the assumption that the angle_res is more or less comparable to actual beam width it means that there is at least the beam edge on the new angle.

Instead I'm suggesting to check if we get NaN rays and issue a warning. Then the user can take measures to resolve.

Another solution would be to check actual beam width and angle_res and try to calculate a decent value from that to make sure the above assumptions are fulfilled.

Maybe we should set up some test case with artificial values and try to work along until we are happy.

[aladinor]

Sorry @kmuehlbauer for my late reply. I think we don't need to change any default value since we can pass the tolerance parameter to xd.util.reindex_angle method. Thanks for your time!