Reproject Radar Coordinates#

This notebook demonstrates the target_crs option now exposed on .xradar.georeference() for xarray.Dataset and xarray.DataTree.

It lets you reproject radar x, y coordinates from the radar-native Azimuthal Equidistant (AEQD) projection into any pyproj-compatible CRS in a single accessor call:

radar = radar.xradar.georeference(target_crs=4326)

Imports#

import cartopy.crs as ccrs
import cmweather  # noqa
import matplotlib.pyplot as plt
import pyproj
from open_radar_data import DATASETS

import xradar as xd

1. Default georeferencing (AEQD, meters)#

Without target_crs, the accessor adds x, y, z in the radar-native AEQD projection (distances in meters from the radar site).

radar_aeqd = radar.xradar.georeference()
sweep = radar_aeqd["sweep_0"].to_dataset(inherit="all_coords")
crs = sweep.xradar.get_crs()
print("projection:", crs.coordinate_operation.method_name if crs.coordinate_operation else crs.name)
print("is_projected / is_geographic:", crs.is_projected, "/", crs.is_geographic)
print("x units:", sweep.x.attrs.get("units"))
print("x range:", float(sweep.x.min()), float(sweep.x.max()))

projection: Azimuthal Equidistant
is_projected / is_geographic: True / False
x units: meters
x range: -149353.7199056801 149354.1123565658

fig, ax = plt.subplots(figsize=(7, 6))
radar_aeqd["sweep_0"]["DBZ"].plot(x="x", y="y", cmap="ChaseSpectral", ax=ax)
ax.set_title("AEQD (default) — x, y in meters from radar")
ax.set_aspect("equal")

../_images/4abe39f19fd032a71a7b737b1d995645ff86722fae0b15853dbe9259b038343c.png

2. Reproject to geographic lon/lat (EPSG:4326)#

Pass target_crs=4326 to get x as longitude and y as latitude. Attributes are updated automatically (standard_name, units).

radar_geo = radar.xradar.georeference(target_crs=4326)
sweep = radar_geo["sweep_0"].to_dataset(inherit="all_coords")
crs = sweep.xradar.get_crs()
print("CRS name:", crs.name, "| EPSG:", crs.to_epsg())
print("x attrs:", dict(sweep.x.attrs))
print("y attrs:", dict(sweep.y.attrs))

CRS name: WGS 84 | EPSG: 4326
x attrs: {'standard_name': 'longitude', 'long_name': 'Geographic Longitude', 'units': 'degrees_east', 'axis': 'X'}
y attrs: {'standard_name': 'latitude', 'long_name': 'Geographic Latitude', 'units': 'degrees_north', 'axis': 'Y'}

fig, ax = plt.subplots(figsize=(8, 6))
radar_geo["sweep_0"]["DBZ"].plot(x="x", y="y", cmap="ChaseSpectral", ax=ax)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("EPSG:4326 — x=lon, y=lat")

Text(0.5, 1.0, 'EPSG:4326 — x=lon, y=lat')

../_images/84721baf60008747afdf3446a45fadb684710823f1d6a79dd27edddf379cf0e3.png

3. Reproject to a projected CRS (UTM, Web Mercator)#

Any pyproj-compatible CRS input works: EPSG code, WKT string, or a pyproj.CRS instance.

# Web Mercator (meters)
radar_wm = radar.xradar.georeference(target_crs="EPSG:3857")
sweep_wm = radar_wm["sweep_0"].to_dataset(inherit="all_coords")
crs = sweep_wm.xradar.get_crs()
print("CRS:", crs.name, "| EPSG:", crs.to_epsg())

CRS: WGS 84 / Pseudo-Mercator | EPSG: 3857

utm_crs = pyproj.CRS("EPSG:32633")  # UTM zone 33N (tweak for your radar site)
radar_utm = radar.xradar.georeference(target_crs=utm_crs)
sweep_utm = radar_utm["sweep_0"].to_dataset(inherit="all_coords")
crs = sweep_utm.xradar.get_crs()
print("CRS:", crs.name, "| EPSG:", crs.to_epsg())

CRS: WGS 84 / UTM zone 33N | EPSG: 32633

4. Plot on a cartopy map#

Once data is in a known CRS, hand it off to cartopy for map plotting.

radar_geo = radar.xradar.georeference(target_crs=4326)

fig = plt.figure(figsize=(9, 8))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
radar_geo["sweep_0"]["DBZ"].plot(
    x="x",
    y="y",
    cmap="ChaseSpectral",
    transform=ccrs.PlateCarree(),
    ax=ax,
    cbar_kwargs=dict(pad=0.05, shrink=0.7),
)
ax.coastlines()
ax.gridlines(draw_labels=True)
ax.set_title("Georeferenced PPI on PlateCarree map")

Text(0.5, 1.0, 'Georeferenced PPI on PlateCarree map')

/home/docs/checkouts/readthedocs.org/user_builds/xradar/conda/latest/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

../_images/b27e738f6f33dd901afbe1d05a66492db72ad857ec21f843fd43ee584fdfebd4.png

5. Side-by-side comparison of projections#

Compare the same sweep in AEQD (meters), geographic (lon/lat), and Lambert Conformal Conic — a projection commonly used for mid-latitude weather maps.

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# AEQD (default)
radar_aeqd = radar.xradar.georeference()
radar_aeqd["sweep_0"]["DBZ"].plot(
    x="x", y="y", cmap="ChaseSpectral", ax=axes[0], add_colorbar=False,
)
axes[0].set_title("AEQD (meters)")
axes[0].set_aspect("equal")

# EPSG:4326 (lon/lat)
radar_geo = radar.xradar.georeference(target_crs=4326)
radar_geo["sweep_0"]["DBZ"].plot(
    x="x", y="y", cmap="ChaseSpectral", ax=axes[1], add_colorbar=False,
)
axes[1].set_title("EPSG:4326 (lon/lat)")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")

# Lambert Conformal Conic
lcc_crs = pyproj.CRS.from_proj4(
    "+proj=lcc +lat_1=30 +lat_2=60 +lat_0=39.5 +lon_0=-105 +datum=WGS84"
)
radar_lcc = radar.xradar.georeference(target_crs=lcc_crs)
radar_lcc["sweep_0"]["DBZ"].plot(
    x="x", y="y", cmap="ChaseSpectral", ax=axes[2], add_colorbar=False,
)
axes[2].set_title("Lambert Conformal Conic")
axes[2].set_aspect("equal")

fig.tight_layout()

../_images/0e5c78545d1362318ebd0190cb23096d46d08af9c2a007c531b6e01bca77c59c.png

6. Radar data on different map projections#

The same georeferenced radar data rendered on four different cartopy map projections — coastlines and gridlines make the projection distortion visible.

import cartopy.feature as cfeature

radar_geo = radar.xradar.georeference(target_crs=4326)
lon0 = float(radar_geo.ds.coords["longitude"].values)
lat0 = float(radar_geo.ds.coords["latitude"].values)

projections = [
    ("PlateCarree", ccrs.PlateCarree()),
    ("Lambert Conformal", ccrs.LambertConformal(
        central_longitude=lon0, central_latitude=lat0,
    )),
    ("Orthographic (globe)", ccrs.Orthographic(
        central_longitude=lon0, central_latitude=lat0,
    )),
    ("Mercator", ccrs.Mercator(central_longitude=lon0)),
]

fig = plt.figure(figsize=(14, 12))

for i, (name, proj) in enumerate(projections, 1):
    ax = fig.add_subplot(2, 2, i, projection=proj)
    radar_geo["sweep_0"]["DBZ"].plot(
        x="x", y="y", cmap="ChaseSpectral",
        transform=ccrs.PlateCarree(), ax=ax,
        add_colorbar=False,
    )
    ax.coastlines()
    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
    ax.gridlines(draw_labels=isinstance(proj, (ccrs.PlateCarree, ccrs.Mercator)))
    ax.set_title(name)

fig.tight_layout()

/home/docs/checkouts/readthedocs.org/user_builds/xradar/conda/latest/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

../_images/7fb8ecfb9abce978adc73a74d8181fd78d7f119c0ec56888147e60ea3c7e0379.png

Summary#

target_crs accepts anything pyproj.CRS(...) accepts: int EPSG codes, "EPSG:xxxx" strings, WKT, or pyproj.CRS instances.
The spatial_ref / crs_wkt coordinate is updated to reflect the target CRS — sweep.xradar.get_crs() will return it.
Currently only x and y are reprojected; z stays as AEQD altitude above the radar site.
If you want separate lon / lat / alt coordinates (without overwriting x, y), see the Assign_GeoCoords notebook.