CfRadial1 - Export#

Imports#

[1]:
import cmweather  # noqa
import xarray as xr
from open_radar_data import DATASETS

import xradar as xd

Download#

Fetching CfRadial1 radar data file from open-radar-data repository.

[2]:
filename = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc")
[3]:
radar = xd.io.open_cfradial1_datatree(filename, first_dim="auto")
display(radar)
/home/docs/checkouts/readthedocs.org/user_builds/xradar/conda/latest/lib/python3.12/site-packages/xradar/io/backends/cfradial1.py:349: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  ds = open_dataset(filename_or_obj, engine=engine, **kwargs)
/home/docs/checkouts/readthedocs.org/user_builds/xradar/conda/latest/lib/python3.12/site-packages/xradar/io/backends/cfradial1.py:349: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  ds = open_dataset(filename_or_obj, engine=engine, **kwargs)
<xarray.DatasetView> Size: 477B
Dimensions:              (sweep: 9)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int32 4B ...
    platform_type        |S32 32B ...
    primary_axis         |S32 32B ...
    status_str           |S1 1B ...
    instrument_type      |S32 32B ...
    time_coverage_start  |S32 32B ...
    time_coverage_end    |S32 32B ...
    latitude             float64 8B ...
    longitude            float64 8B ...
    altitude             float64 8B ...
    sweep_group_name     (sweep) <U7 252B 'sweep_0' 'sweep_1' ... 'sweep_8'
    sweep_fixed_angle    (sweep) float32 36B ...
Attributes: (12/13)
    Conventions:         CF/Radial instrument_parameters radar_parameters rad...
    version:             1.2
    title:               TIMREX
    institution:
    references:
    source:
    ...                  ...
    comment:
    instrument_name:     SPOLRVP8
    site_name:
    scan_name:
    scan_id:             0
    platform_is_mobile:  false

Plot Azimuth vs. Range#

[4]:
radar.sweep_0.DBZ.plot(cmap="ChaseSpectral", vmin=-10, vmax=70)
[4]:
<matplotlib.collections.QuadMesh at 0x70f9fc7beff0>
../_images/notebooks_CfRadial1_Export_7_1.png

Plot Time vs. Range#

[5]:
radar.sweep_0.DBZ.swap_dims({"azimuth": "time"}).sortby("time").plot(
    cmap="ChaseSpectral", vmin=-10, vmax=70
)
[5]:
<matplotlib.collections.QuadMesh at 0x70f9fc6e7bc0>
../_images/notebooks_CfRadial1_Export_9_1.png

Georeference#

[6]:
radar = radar.xradar.georeference()
display(radar)
<xarray.DatasetView> Size: 477B
Dimensions:              (sweep: 9)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int32 4B ...
    platform_type        |S32 32B ...
    primary_axis         |S32 32B ...
    status_str           |S1 1B ...
    instrument_type      |S32 32B ...
    time_coverage_start  |S32 32B ...
    time_coverage_end    |S32 32B ...
    latitude             float64 8B ...
    longitude            float64 8B ...
    altitude             float64 8B ...
    sweep_group_name     (sweep) <U7 252B 'sweep_0' 'sweep_1' ... 'sweep_8'
    sweep_fixed_angle    (sweep) float32 36B ...
Attributes: (12/13)
    Conventions:         CF/Radial instrument_parameters radar_parameters rad...
    version:             1.2
    title:               TIMREX
    institution:
    references:
    source:
    ...                  ...
    comment:
    instrument_name:     SPOLRVP8
    site_name:
    scan_name:
    scan_id:             0
    platform_is_mobile:  false

Plot PPI#

[7]:
radar["sweep_0"]["DBZ"].plot(x="x", y="y", cmap="ChaseSpectral", vmin=-10, vmax=70)
[7]:
<matplotlib.collections.QuadMesh at 0x70f9fc5efce0>
../_images/notebooks_CfRadial1_Export_13_1.png

Filter#

Apply basic reflectivity filter. This is just a demonstration.

[8]:
def ref_filter(dtree, sweep="sweep_0", field="DBZ"):
    ds = dtree[sweep].ds
    ds = ds.where((ds[field] >= -10) & (ds[field] <= 70))
    red_patch = ds.where(
        (
            (ds[field] >= ds[field].max().values - 0.5)
            & (ds[field] <= ds[field].max().values + 0.5)
        ),
        drop=True,
    )
    rmin, rmax = int(red_patch.range.min().values - 150), int(
        red_patch.range.max().values + 150
    )
    out_of_range_mask = (ds.range < rmin) | (ds.range > rmax)
    ds[field] = ds[field].where(out_of_range_mask)
    # Interpolate missing values using the slinear method along the 'range' dimension
    ds[field] = ds[field].interpolate_na(dim="range", method="slinear")
    dtree[sweep][f"corr_{field}"] = ds[field].copy()
    return dtree[sweep]
[9]:
swp0 = ref_filter(radar, sweep="sweep_0", field="DBZ")
[10]:
swp0.corr_DBZ.plot(x="x", y="y", cmap="ChaseSpectral", vmin=-10, vmax=70)
[10]:
<matplotlib.collections.QuadMesh at 0x70f9fc49dc70>
../_images/notebooks_CfRadial1_Export_17_1.png

Filter full volume#

[11]:
# Initialize an empty DataTree
result_tree = xr.DataTree()

for sweep in radar.sweep_group_name.values:
    corrected_data = ref_filter(radar, sweep, field="DBZ")

    # Convert the xarray Dataset to a DataTree and add it to the result_tree
    data_tree = xr.DataTree.from_dict(corrected_data.to_dict())

    # Copy the contents of data_tree into result_tree
    for key, value in data_tree.items():
        result_tree[key] = value
[12]:
radar.sweep_6.corr_DBZ.plot(x="x", y="y", cmap="ChaseSpectral", vmin=-10, vmax=70)
[12]:
<matplotlib.collections.QuadMesh at 0x70f9fc2f2ff0>
../_images/notebooks_CfRadial1_Export_20_1.png

Export#

Export to CfRadial1

[13]:
xd.io.to_cfradial1(dtree=radar, filename="cfradial1_qced.nc", calibs=True)
[14]:
?xd.io.to_cfradial1

Note#

If filename is None in the xd.io.to_cfradial1 function, it will automatically generate a filename using the instrument name and the first available timestamp from the data.