Assign GeoCoords#

This notebook demonstrates how to work with geocoordinates (longitude and latitude) instead of radar-centric x and y coordinates.

import cmweather  # noqa: F401
import matplotlib.pyplot as plt
import xarray as xr
from open_radar_data import DATASETS

import xradar as xd

Define Functions#

def get_geocoords(ds):
    """
    Converts Cartesian coordinates (x, y, z) in a radar dataset to geographic
    coordinates (longitude, latitude, altitude) using CRS transformation.

    Parameters
    ----------
    ds : xarray.Dataset
        Radar dataset with Cartesian coordinates.

    Returns
    -------
    xarray.Dataset
        Dataset with added 'lon', 'lat', and 'alt' coordinates and their attributes.
    """
    from pyproj import CRS, Transformer

    # Convert the dataset to georeferenced coordinates
    ds = ds.xradar.georeference()
    # Define source and target coordinate reference systems (CRS)
    src_crs = ds.xradar.get_crs()
    trg_crs = CRS.from_user_input(4326)  # EPSG:4326 (WGS 84)
    # Create a transformer for coordinate conversion
    transformer = Transformer.from_crs(src_crs, trg_crs)
    # Transform x, y, z coordinates to latitude, longitude, and altitude
    trg_y, trg_x, trg_z = transformer.transform(ds.x, ds.y, ds.z)
    # Assign new coordinates with appropriate attributes
    ds = ds.assign_coords(
        {
            "lon": (ds.x.dims, trg_x, xd.model.get_longitude_attrs()),
            "lat": (ds.y.dims, trg_y, xd.model.get_latitude_attrs()),
            "alt": (ds.z.dims, trg_z, xd.model.get_altitude_attrs()),
        }
    )
    return ds


def fix_sitecoords(ds):
    coords = ["longitude", "latitude", "altitude", "altitude_agl"]
    for coord in coords:
        if coord not in ds:
            continue
        # Skip coords that are already scalar
        if ds[coord].ndim == 0:
            continue
        # Compute median excluding NaN
        data = ds[coord].median(skipna=True).item()
        attrs = ds[coord].attrs
        ds = ds.assign_coords({coord: xr.DataArray(data=data, attrs=attrs)})
    return ds

Load Data#

file1 = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc")
file2 = DATASETS.fetch("cfrad.20211011_201557.188_to_20211011_201617.720_DOW8_PPI.nc")

Downloading file 'cfrad.20080604_002217_000_SPOL_v36_SUR.nc' from 'https://github.com/openradar/open-radar-data/raw/main/data/cfrad.20080604_002217_000_SPOL_v36_SUR.nc' to '/home/docs/.cache/open-radar-data'.

Downloading file 'cfrad.20211011_201557.188_to_20211011_201617.720_DOW8_PPI.nc' from 'https://github.com/openradar/open-radar-data/raw/main/data/cfrad.20211011_201557.188_to_20211011_201617.720_DOW8_PPI.nc' to '/home/docs/.cache/open-radar-data'.

Example #2#

dtree2 = xd.io.open_cfradial1_datatree(file2)

try:
    dtree2 = dtree2.xradar.georeference()
except Exception:
    print("Georeferencing failed!")

Georeferencing failed!

ds2 = dtree2["sweep_0"].to_dataset(inherit="all_coords")
print("Longitude", ds2["longitude"].shape)
print("Latitude", ds2["latitude"].shape)

Longitude (731,)
Latitude (731,)

Important Note:

This radar data is from a mobile research radar called Doppler on Wheels (DOW). In previous versions of xradar, site coordinates (latitude, longitude) were stored per-ray on each sweep, resulting in shape (731,). Now, station coordinates are stored as scalar values on the root node of the DataTree and inherited by sweeps via to_dataset(inherit="all_coords"). The fix_sitecoords helper gracefully handles both cases — it skips coordinates that are already scalar.

Assign GeoCoords#

Define Functions#

Load Data#

Example #1#

Assign lat, lon, and alt#

Example #2#

Fix Coords#