ENH: Expose target_crs in georeference accessor

docs/history.md

@@ -1,5 +1,9 @@
 # History
 
+## Development
+
+* ENH: Expose ``target_crs`` on ``.xradar.georeference()`` accessor for DataArray, Dataset and DataTree, enabling reprojection via e.g. ``radar.xradar.georeference(target_crs=4326)`` ({issue}`243`) by [@syedhamidali](https://github.com/syedhamidali)
+
 ## 0.12.0 (2026-04-21)
 
 * MNT: Unpin xarray, require ``xarray >= 2026.4.0`` in ``requirements.txt``, ``environment.yml``, ``ci/unittests.yml``, and ``ci/notebooktests.yml`` by [@aladinor](https://github.com/aladinor)

docs/notebooks/Georeference_TargetCRS.md

@@ -0,0 +1,171 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.19.1
+  main_language: python
+kernelspec:
+  display_name: Python 3
+  name: python3
+---
+
+# Reproject Radar Coordinates via the Accessor
+
+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:
+
+```python
+radar = radar.xradar.georeference(target_crs=4326)
+```
+
++++
+
+## Imports
+
+```{code-cell}
+import cmweather  # noqa
+import matplotlib.pyplot as plt
+from open_radar_data import DATASETS
+
+import xradar as xd
+```
+
+## Read a sample radar file
+
+```{code-cell}
+filename = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc")
+radar = xd.io.open_cfradial1_datatree(filename, first_dim="auto")
+radar
+```
+
++++
+
+## 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).
+
+```{code-cell}
+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()))
+```
+
+```{code-cell}
+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")
+```
+
++++
+
+## 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``).
+
+```{code-cell}
+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))
+```
+
+```{code-cell}
+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")
+```
+
++++
+
+## 3. Reproject to a projected CRS (UTM, Web Mercator)
+
+Any ``pyproj``-compatible CRS input works: EPSG code, WKT string, or a
+``pyproj.CRS`` instance.
+
+```{code-cell}
+# 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())
+```
+
+```{code-cell}
+import pyproj
+
+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())
+```
+
++++
+
+## 4. Plot on a cartopy map
+
+Once data is in a known CRS, hand it off to cartopy for map plotting.
+
+```{code-cell}
+import cartopy.crs as ccrs
+
+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")
+```
+
++++
+
+## 5. Works on a single Dataset too
+
+The accessor is available on ``Dataset`` and ``DataArray`` as well, not just
+``DataTree``.
+
+```{code-cell}
+ds = radar["sweep_0"].to_dataset(inherit="all_coords")
+ds_geo = ds.xradar.georeference(target_crs=4326)
+ds_geo[["DBZ"]]
+```
+
++++
+
+## Notes
+
+* ``target_crs`` accepts anything ``pyproj.CRS(...)`` accepts: int
+  [EPSG codes](https://epsg.io/), ``"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](Assign_GeoCoords.md) notebook.

docs/usage.md

@@ -44,6 +44,7 @@ notebooks/UF
 :caption: Examples
 
 notebooks/Assign_GeoCoords
+notebooks/Georeference_TargetCRS
 notebooks/CfRadial1_Export
 notebooks/Read-plot-Sigmet-data-from-AWS
 notebooks/plot-ppi

tests/test_accessors.py

@@ -51,6 +51,27 @@ def test_georeference_datatree():
     )
 
 
+def test_georeference_target_crs_dataset():
+    radar = xd.model.create_sweep_dataset()
+    geo = radar.xradar.georeference(target_crs=4326)
+    crs = geo.xradar.get_crs()
+    assert crs.is_geographic
+    assert geo.x.attrs["standard_name"] == "longitude"
+    assert geo.y.attrs["standard_name"] == "latitude"
+    assert_almost_equal(geo.x.values[0, 0], radar.longitude.values, decimal=3)
+    assert_almost_equal(geo.y.values[0, 0], radar.latitude.values, decimal=3)
+
+
+def test_georeference_target_crs_datatree():
+    radar = xd.model.create_sweep_dataset()
+    tree = xr.DataTree.from_dict({"sweep_0": radar})
+    geo = tree.xradar.georeference(target_crs=4326)["sweep_0"]
+    assert geo["x"].attrs["standard_name"] == "longitude"
+    assert geo["y"].attrs["standard_name"] == "latitude"
+    assert_almost_equal(geo["x"].values[0, 0], radar.longitude.values, decimal=3)
+    assert_almost_equal(geo["y"].values[0, 0], radar.latitude.values, decimal=3)
+
+
 def test_crs_dataarray():
     """Test the add_crs and get_crs methods on a DataArray."""
     # Create a sample DataArray with radar data

xradar/accessors.py

@@ -73,7 +73,7 @@ class XradarDataArrayAccessor(XradarAccessor):
     """Adds a number of xradar specific methods to xarray.DataArray objects."""
 
     def georeference(
-        self, earth_radius=None, effective_radius_fraction=None
+        self, earth_radius=None, effective_radius_fraction=None, target_crs=None
     ) -> xr.DataArray:
         """
         Parameters
@@ -83,6 +83,9 @@ def georeference(
             WGS84 ellipsoid.
         effective_radius_fraction: float
             Fraction of earth to use for the effective radius (default is 4/3).
+        target_crs: pyproj.CRS or int or str, optional
+            Coordinate reference system to reproject x, y to (e.g. ``4326`` for
+            geographic lon/lat). If not provided, radar-native AEQD is used.
         Returns
         -------
         da = xr.DataArray
@@ -94,6 +97,7 @@ def georeference(
             get_x_y_z,
             earth_radius=earth_radius,
             effective_radius_fraction=effective_radius_fraction,
+            target_crs=target_crs,
         )
 
     def add_crs(self) -> xr.DataArray:
@@ -123,7 +127,7 @@ class XradarDataSetAccessor(XradarAccessor):
     """Adds a number of xradar specific methods to xarray.Dataset objects."""
 
     def georeference(
-        self, earth_radius=None, effective_radius_fraction=None
+        self, earth_radius=None, effective_radius_fraction=None, target_crs=None
     ) -> xr.Dataset:
         """
         Add georeference information to an xarray dataset
@@ -134,6 +138,9 @@ def georeference(
             WGS84 ellipsoid.
         effective_radius_fraction: float
             Fraction of earth to use for the effective radius (default is 4/3).
+        target_crs: pyproj.CRS or int or str, optional
+            Coordinate reference system to reproject x, y to (e.g. ``4326`` for
+            geographic lon/lat). If not provided, radar-native AEQD is used.
         Returns
         -------
         da = xarray.Dataset
@@ -144,6 +151,7 @@ def georeference(
             get_x_y_z,
             earth_radius=earth_radius,
             effective_radius_fraction=effective_radius_fraction,
+            target_crs=target_crs,
         )
 
     def add_crs(self) -> xr.DataSet:
@@ -185,7 +193,7 @@ class XradarDataTreeAccessor(XradarAccessor):
     """Adds a number of xradar specific methods to xarray.DataTree objects."""
 
     def georeference(
-        self, earth_radius=None, effective_radius_fraction=None
+        self, earth_radius=None, effective_radius_fraction=None, target_crs=None
     ) -> xr.DataTree:
         """
         Add georeference information to an xradar datatree object
@@ -196,6 +204,9 @@ def georeference(
             WGS84 ellipsoid.
         effective_radius_fraction: float
             Fraction of earth to use for the effective radius (default is 4/3).
+        target_crs: pyproj.CRS or int or str, optional
+            Coordinate reference system to reproject x, y to (e.g. ``4326`` for
+            geographic lon/lat). If not provided, radar-native AEQD is used.
         Returns
         -------
         da = xarray.DataTree
@@ -206,6 +217,7 @@ def georeference(
             get_x_y_z_tree,
             earth_radius=earth_radius,
             effective_radius_fraction=effective_radius_fraction,
+            target_crs=target_crs,
         )
 
     def add_crs(self) -> xr.DataTree: