Access S-Pol PRECIP Radar Data through the NCAR GDEX and Convert CfRadial1 to CfRadial2¶
Required Packages¶
The following packages are needed for this example. Please make sure they are installed before moving forward.
xarrayxradarh5netcdf>=1.8.0netCDF4zarrnumpymatplotlibcmweather
Why CfRadial1 to CfRadial2¶
Radar data are often distributed in the CfRadial1 format: one flat NetCDF file per volume, where all sweeps are stacked along the ray (time) dimension.
The newer CfRadial2 data model (based on the WMO FM 301 standard) stores each sweep in its own group, together with dedicated groups for radar parameters, calibration and georeferencing. That hierarchical structure maps naturally onto an xarray.DataTree, and onto cloud-optimized formats such as Zarr.
In this notebook we will:
Open a CfRadial1 NetCDF file from GDEX.
Convert it to a CfRadial2
DataTreewith xradar’s.xradar.to_cf2()method (see the xradar Transform notebook).Write the CfRadial2
DataTreeto a local Zarr store and to a NetCDF4 file.Do a round trip, CfRadial2 → CfRadial1, and check that the radar data are unchanged compared with the original NetCDF file (see the xradar CfRadial1 to CfRadial2 Data Model Transformation notebook).
1. Open CfRadial1¶
The example file is available on GDEX both through the GDEX POSIX path (on NCAR HPC systems such as Derecho/Casper) and through an HTTPS data URL.
In this example we use the GDEX POSIX path. If you are not on an NCAR system, set use_posix = False to read the same file over HTTPS.
# Please specify your preferred data access method: the Data URL or the GDEX POSIX path.
# filename = "cfrad.20220525_021216.343_to_20220525_021248.129_SPOL_PrecipRhi1_RHI.nc"
filename = "cfrad.20220525_015443.665_to_20220525_015709.985_SPOL_PrecipRhi2_RHI.nc"
data_url = f"https://data.gdex.ucar.edu/special_projects/ard_radar/tutorial/{filename}"
data_posix_path = f"/gdex/data/special_projects/ard_radar/tutorial/{filename}"
# True on NCAR HPC (POSIX access), False to stream the file over HTTPS
use_posix = TrueWe open the CfRadial1 file as a plain xarray.Dataset. Note that import xradar is required to make the .xradar accessor available.
import numpy as np
import xarray as xr
import xradar as xd
if use_posix:
source = data_posix_path
else:
source = data_url
ds_cf1 = xr.open_dataset(source, engine="h5netcdf")
ds_cf1The CfRadial1 file is flat: all sweeps share the time dimension, and because n_gates_vary = "true" the radar moments (e.g. DBZ) are 1-D arrays along n_points.
print("n_gates_vary :", ds_cf1.attrs.get("n_gates_vary"))
print("sweeps :", ds_cf1.sizes["sweep"])
print("rays (time) :", ds_cf1.sizes["time"])
print("n_points :", ds_cf1.sizes["n_points"])
print("DBZ dims :", ds_cf1["DBZ"].dims)n_gates_vary : true
sweeps : 11
rays (time) : 2866
n_points : 5729134
DBZ dims : ('n_points',)
2. Convert CfRadial1 to CfRadial2¶
The .xradar.to_cf2() method (an alias of .xradar.to_cfradial2() / .xradar.to_cfradial2_datatree()) converts the flat CfRadial1 Dataset
into a CfRadial2 DataTree, following the xradar Transform notebook.
The ragged n_points fields are unpacked onto a regular (azimuth, range) grid per sweep, with NaN padding where a ray has fewer gates.
The method returns a new object; it does not modify ds_cf1 in place.
dtree = ds_cf1.xradar.to_cf2()
dtree# Check the CfRadial2 hierarchy: root metadata, metadata subgroups and one group per sweep
for node in dtree.subtree:
print(f"{node.path:25s} dims={dict(node.dataset.sizes)}")/ dims={'sweep': 11, 'frequency': 1}
/radar_parameters dims={'frequency': 1}
/radar_calibration dims={'frequency': 1}
/georeferencing_correction dims={'frequency': 1}
/sweep_0 dims={'elevation': 252, 'range': 1999, 'frequency': 1}
/sweep_1 dims={'elevation': 259, 'range': 1999, 'frequency': 1}
/sweep_2 dims={'elevation': 259, 'range': 1999, 'frequency': 1}
/sweep_3 dims={'elevation': 255, 'range': 1999, 'frequency': 1}
/sweep_4 dims={'elevation': 256, 'range': 1999, 'frequency': 1}
/sweep_5 dims={'elevation': 260, 'range': 1999, 'frequency': 1}
/sweep_6 dims={'elevation': 266, 'range': 1999, 'frequency': 1}
/sweep_7 dims={'elevation': 265, 'range': 1999, 'frequency': 1}
/sweep_8 dims={'elevation': 266, 'range': 1999, 'frequency': 1}
/sweep_9 dims={'elevation': 265, 'range': 1999, 'frequency': 1}
/sweep_10 dims={'elevation': 263, 'range': 1999, 'frequency': 1}
A quick look at the data¶
Before writing anything out, we plot the reflectivity of the first RHI sweep. georeference() adds Cartesian x, y, z coordinates
(in meters, relative to the radar), and for an RHI we plot height against the horizontal distance from the radar.
import cmweather # registers radar colormaps such as "ChaseSpectral"
import matplotlib.pyplot as plt
ds = dtree["sweep_0"].to_dataset()
# bring in the radar location information as coordinates to the Dataset object for the georeferencing
ds['longitude']=dtree["longitude"].values
ds['latitude']=dtree["latitude"].values
ds['altitude']=dtree['altitude'].values
ds = ds.set_coords(["longitude", "latitude", "altitude"])
sweep = ds.xradar.georeference()
horizontal_distance = np.sqrt(sweep["x"] ** 2 + sweep["y"] ** 2) / 1000.0 # in kilometers
sweep = sweep.assign_coords(
horizontal_distance=horizontal_distance.assign_attrs(long_name="Horizontal distance", units="km"),
height=(sweep["z"] / 1000.0).assign_attrs(long_name="Height", units="km"),
)
fig, ax = plt.subplots(figsize=(10, 4))
sweep["DBZ"].plot(
ax=ax,
x="horizontal_distance",
y="height",
cmap="ChaseSpectral",
vmin=-10,
vmax=60,
)
ax.set_ylim(0, 15)
ax.set_title(f"S-Pol RHI reflectivity, fixed azimuth {float(sweep.sweep_fixed_angle):.1f}°")
Write the CfRadial2 DataTree to a local Zarr store¶
Some variables decoded from the CfRadial1 character arrays (e.g. time_coverage_start, instrument_name) are fixed-length string arrays.
These are not reliably supported by Zarr, so we convert them to variable-length (object) strings before writing.
Only the string variables are copied; all other arrays are shared by reference.
from pathlib import Path
output_dir = Path("cf1_to_cf2_tutorial")
output_dir.mkdir(exist_ok=True)
zarr_path = output_dir / filename.replace(".nc", "_cfradial2.zarr")
nc4_path = output_dir / filename.replace(".nc", "_cfradial2.nc")
def fix_string_encoding(ds):
"""Convert fixed-length string variables to object dtype for Zarr compatibility."""
string_vars = {}
for var in ds.variables:
if ds[var].dtype.kind in ("U", "S"):
string_vars[var] = ds[var].astype(str).astype(object)
if not string_vars:
return ds
return ds.assign(string_vars)
def exact_time_encoding(tree):
"""Encode each sweep's ``time`` as int64 nanoseconds so it round-trips exactly.
The default float64 "seconds since ..." encoding can shift a few timestamps
by ~1 ns after a write/read cycle.
"""
encoding = {}
for node in tree.subtree:
if "time" not in node.variables:
continue
# Skip groups whose `time` is not a datetime (some metadata groups store it as a string)
if node["time"].dtype.kind != "M":
continue
reference = str(node["time"].values.min().astype("datetime64[s]"))
encoding[node.path] = {
"time": {
"units": f"nanoseconds since {reference}",
"calendar": "gregorian",
"dtype": "int64",
}
}
return encoding%%time
# Apply the string encoding fix to the entire dataset tree and save it as a Zarr store
dtree_for_zarr = dtree.map_over_datasets(fix_string_encoding)
# Save the fixed dataset tree as a Zarr store
# Store time as exact int64 nanoseconds (see exact_time_encoding above)
dtree_for_zarr.to_zarr(
zarr_path,
mode="w",
encoding=exact_time_encoding(dtree_for_zarr),
)/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
CPU times: user 7.91 s, sys: 2.16 s, total: 10.1 s
Wall time: 10.2 s
/glade/u/home/chiaweih/test_xradar/.venv/lib/python3.11/site-packages/zarr/api/asynchronous.py:247: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
<xarray.backends.zarr.ZarrStore at 0x7f3d141b8ee0>Reopen the Zarr store lazily as a DataTree. The groups are the same CfRadial2 groups we wrote.
dtree_zarr = xr.open_datatree(zarr_path, engine="zarr")
dtree_zarrWrite the same CfRadial2 DataTree to NetCDF4¶
The CfRadial2 model is not tied to Zarr. The same DataTree can be written to a grouped NetCDF4 file with xradar’s
xd.io.to_cfradial2() exporter (or directly with dtree.to_netcdf()).
%%time
# The h5netcdf engine is required to write this DataTree to NetCDF4.
# engine="netcdf4" fails with "RuntimeError: NetCDF: HDF error".
# Use the same exact time encoding as the Zarr store so the two outputs match
dtree_zarr.to_netcdf(
nc4_path,
engine="h5netcdf",
encoding=exact_time_encoding(dtree_zarr),
)
# xd.io.to_cfradial2(dtree, nc4_path, engine="h5netcdf")CPU times: user 2.53 s, sys: 398 ms, total: 2.93 s
Wall time: 2.5 s
# read the NetCDF4 file into a DataTree
dtree_nc4 = xr.open_datatree(nc4_path, engine="h5netcdf")# Both stores should hold the same groups and the same radar data
print("Same group structure:", dtree_zarr.isomorphic(dtree_nc4))
for key in list(dtree_zarr.keys()):
ds_zarr = dtree_zarr[key].to_dataset()
ds_nc4 = dtree_nc4[key].to_dataset()
same_data = ds_zarr.equals(ds_nc4)
print(f"{key:25s} data equal: {same_data}")Same group structure: True
altitude_agl data equal: True
instrument_type data equal: True
platform_type data equal: True
primary_axis data equal: True
status_str data equal: True
sweep_fixed_angle data equal: True
sweep_group_name data equal: True
time_coverage_end data equal: True
time_coverage_start data equal: True
volume_number data equal: True
georeferencing_correction data equal: True
radar_calibration data equal: True
radar_parameters data equal: True
sweep_0 data equal: True
sweep_1 data equal: True
sweep_10 data equal: True
sweep_2 data equal: True
sweep_3 data equal: True
sweep_4 data equal: True
sweep_5 data equal: True
sweep_6 data equal: True
sweep_7 data equal: True
sweep_8 data equal: True
sweep_9 data equal: True
# Size on disk of the original CfRadial1 file and the two CfRadial2 outputs
def size_mb(path):
path = Path(path)
if path.is_dir():
total = 0
for file in path.rglob("*"):
if file.is_file():
total += file.stat().st_size
return total / 1e6
return path.stat().st_size / 1e6
if use_posix:
print(f"CfRadial1 NetCDF : {size_mb(data_posix_path):7.2f} MB")
print(f"CfRadial2 Zarr : {size_mb(zarr_path):7.2f} MB")
print(f"CfRadial2 NetCDF4: {size_mb(nc4_path):7.2f} MB")CfRadial1 NetCDF : 65.17 MB
CfRadial2 Zarr : 62.26 MB
CfRadial2 NetCDF4: 597.48 MB
Round trip: CfRadial2 → CfRadial1, compared with the original file¶
To check that nothing is lost, we convert the CfRadial2 tree read back from the Zarr store into a CfRadial1 Dataset
with .xradar.to_cfradial1_dataset() (alias .xradar.to_cfradial1()). This is the inverse of the transformation described in the xradar
CfRadial1 to CfRadial2 Data Model Transformation notebook.
The full round trip is therefore: CfRadial1 NetCDF → CfRadial2 DataTree → Zarr → CfRadial2 DataTree → CfRadial1 Dataset.
ds_roundtrip = dtree_zarr.xradar.to_cfradial1_dataset()
ds_roundtripxradar always exports the padded CfRadial1 layout: fields are (time, range) arrays with NaN beyond each ray’s last gate, and
n_gates_vary is set to "false". The original file uses the ragged n_points layout, so the two datasets cannot be compared with a
single identical() call. We compare them in three steps instead:
which variables exist on each side,
the per-ray metadata (
time,azimuth,elevation, ...),the radar moments, after unpacking the original ragged arrays onto the padded grid.
original_vars = set(ds_cf1.variables)
roundtrip_vars = set(ds_roundtrip.variables)
print("Only in original :", sorted(original_vars - roundtrip_vars))
print("Only in round trip:", sorted(roundtrip_vars - original_vars))Only in original : ['grid_mapping', 'r_calib_base_dbz_1km_hc', 'r_calib_base_dbz_1km_hx', 'r_calib_base_dbz_1km_vc', 'r_calib_base_dbz_1km_vx', 'radar_rx_bandwidth', 'ray_angle_res', 'ray_n_gates', 'ray_start_index', 'rays_are_indexed', 'status_xml', 'target_scan_rate']
Only in round trip: ['r_calib_base_1km_hc', 'r_calib_base_1km_hx', 'r_calib_base_1km_vc', 'r_calib_base_1km_vx', 'radar_receiver_bandwidth', 'status_str']
def compare_values(a, b):
"""Return a short status string comparing two arrays (NaNs in the same place count as equal)."""
a = np.asarray(a)
b = np.asarray(b)
if a.shape != b.shape:
return f"shape differs {a.shape} vs {b.shape}"
if a.dtype.kind in "fc" or b.dtype.kind in "fc":
if np.array_equal(a, b, equal_nan=True):
return "identical"
if np.allclose(a, b, equal_nan=True):
max_diff = np.nanmax(np.abs(a.astype(float) - b.astype(float)))
return f"close (max abs diff {max_diff:.3g})"
return "DIFFERS"
if np.array_equal(a, b):
return "identical"
return "DIFFERS"# Per-ray metadata: every variable on the `time` dimension only, present on both sides
ray_vars = []
for var in sorted(original_vars & roundtrip_vars):
if ds_cf1[var].dims == ("time",):
ray_vars.append(var)
for var in ray_vars:
status = compare_values(ds_cf1[var].values, ds_roundtrip[var].values)
print(f"{var:28s} {status}")antenna_transition identical
azimuth identical
elevation identical
measured_transmit_power_h identical
measured_transmit_power_v identical
n_samples identical
nyquist_velocity identical
prt identical
prt_ratio identical
pulse_width identical
r_calib_index identical
ray_gate_spacing identical
ray_start_range identical
scan_rate identical
time identical
unambiguous_range identical
For the radar moments we rebuild a boolean mask of the valid gates from ray_n_gates. Selecting the padded round-trip array with that mask
returns the gates ray by ray, in exactly the order they are stored in the original n_points array, so the two can be compared directly.
We also check that the padding outside the mask is all NaN.
ray_n_gates = ds_cf1["ray_n_gates"].values
ray_start_index = ds_cf1["ray_start_index"].values
# The ragged layout must be contiguous for the mask comparison to be valid
assert np.array_equal(ray_start_index, np.cumsum(ray_n_gates) - ray_n_gates)
n_range = ds_roundtrip.sizes["range"]
valid_gates = np.arange(n_range)[np.newaxis, :] < ray_n_gates[:, np.newaxis]
moment_vars = []
for var in ds_cf1.data_vars:
if ds_cf1[var].dims == ("n_points",):
moment_vars.append(var)
all_identical = True
for var in moment_vars:
if var not in ds_roundtrip:
print(f"{var:14s} MISSING from round trip")
all_identical = False
continue
padded = ds_roundtrip[var].transpose("time", "range").values
status = compare_values(ds_cf1[var].values, padded[valid_gates])
padding_is_nan = bool(np.isnan(padded[~valid_gates]).all())
print(f"{var:14s} {status:12s} padding all NaN: {padding_is_nan}")
if status != "identical" or not padding_is_nan:
all_identical = False
print()
print("All radar moments identical after the round trip:", all_identical)RATE_ZH identical padding all NaN: True
RATE_HYBRID identical padding all NaN: True
PID identical padding all NaN: True
TEMP_FOR_PID identical padding all NaN: True
DBZ identical padding all NaN: True
DBZ_F identical padding all NaN: True
VEL identical padding all NaN: True
VEL_F identical padding all NaN: True
WIDTH identical padding all NaN: True
WIDTH_F identical padding all NaN: True
NCP identical padding all NaN: True
NCP_F identical padding all NaN: True
ZDR_F identical padding all NaN: True
PHIDP_F identical padding all NaN: True
KDP_F identical padding all NaN: True
RHOHV_F identical padding all NaN: True
RHOHV_NNC_F identical padding all NaN: True
SNRHC identical padding all NaN: True
SNRHC_F identical padding all NaN: True
SNRVC identical padding all NaN: True
SNRVC_F identical padding all NaN: True
DBMHC identical padding all NaN: True
DBMHC_F identical padding all NaN: True
DBMVC identical padding all NaN: True
DBMVC_F identical padding all NaN: True
CMD_FLAG identical padding all NaN: True
All radar moments identical after the round trip: True
# Global attributes that changed during the round trip
for key in sorted(set(ds_cf1.attrs) | set(ds_roundtrip.attrs)):
original_value = ds_cf1.attrs.get(key, "<missing>")
roundtrip_value = ds_roundtrip.attrs.get(key, "<missing>")
if str(original_value) != str(roundtrip_value):
print(f"{key}:\n original : {original_value!r}\n round trip: {roundtrip_value!r}")Conventions:
original : 'CF-1.7'
round trip: 'Cf/Radial'
Sub_conventions:
original : 'CF-Radial instrument_parameters radar_parameters radar_calibration'
round trip: '<missing>'
created:
original : '2024/02/20 20:39:37.094'
round trip: '<missing>'
driver:
original : 'RadxConvert(NCAR)'
round trip: '<missing>'
end_datetime:
original : '2022-05-25T01:57:09Z'
round trip: '<missing>'
end_time:
original : '2022-05-25 01:57:09.985'
round trip: '<missing>'
history:
original : ''
round trip: 'xradar v0.12.1.dev36+g00aba7834 CfRadial1 export'
n_gates_vary:
original : 'true'
round trip: '<missing>'
original_format:
original : 'CFRADIAL'
round trip: '<missing>'
start_datetime:
original : '2022-05-25T01:54:43Z'
round trip: '<missing>'
start_time:
original : '2022-05-25 01:54:43.665'
round trip: '<missing>'
time_coverage_end:
original : '2022-05-25T01:57:09Z'
round trip: '<missing>'
time_coverage_start:
original : '2022-05-25T01:54:43Z'
round trip: '<missing>'
version:
original : 'CF-Radial-1.4'
round trip: '1.2'