diff --git a/docs/history.md b/docs/history.md index d48cb738..536cf982 100644 --- a/docs/history.md +++ b/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) diff --git a/docs/notebooks/Georeference_TargetCRS.md b/docs/notebooks/Georeference_TargetCRS.md new file mode 100644 index 00000000..935f7da1 --- /dev/null +++ b/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. diff --git a/docs/usage.md b/docs/usage.md index ecd81e9e..e986d7d8 100644 --- a/docs/usage.md +++ b/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 diff --git a/tests/test_accessors.py b/tests/test_accessors.py index b8120888..37059720 100644 --- a/tests/test_accessors.py +++ b/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 diff --git a/xradar/accessors.py b/xradar/accessors.py index c8a18140..ac9a75cf 100644 --- a/xradar/accessors.py +++ b/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: