Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/history.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

## 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)
* ADD: ``open_cfradial2_datatree`` reader with grouped CfRadial2 compatibility normalization for common FM301/CfRadial2 naming differences ({issue}`93`, {issue}`287`) by [@syedhamidali](https://github.com/syedhamidali)
* ENH: Move station coordinates (``latitude``, ``longitude``, ``altitude``) to root node as coordinates for DataTree coordinate inheritance ({issue}`331`, {pull}`333`) by [@aladinor](https://github.com/aladinor)
* ENH: Add ``optional_groups`` parameter (default ``False``) to all ``open_*_datatree()`` functions to control inclusion of ``/radar_parameters``, ``/georeferencing_correction``, and ``/radar_calibration`` subgroups ({issue}`331`, {pull}`333`) by [@aladinor](https://github.com/aladinor)
Expand Down
290 changes: 290 additions & 0 deletions docs/notebooks/Georeference_TargetCRS.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,290 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "eb6eb2fb",
"metadata": {},
"source": [
"# Reproject Radar Coordinates via the Accessor\n",
"\n",
"This notebook demonstrates the ``target_crs`` option now exposed on\n",
"``.xradar.georeference()`` for ``xarray.Dataset`` and ``xarray.DataTree``.\n",
"\n",
"It lets you reproject radar ``x, y`` coordinates from the radar-native\n",
"Azimuthal Equidistant (AEQD) projection into any ``pyproj``-compatible CRS\n",
"in a single accessor call:\n",
"\n",
"```python\n",
"radar = radar.xradar.georeference(target_crs=4326)\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "8fe22ab0",
"metadata": {},
"source": [
"## Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ea52d5af",
"metadata": {},
"outputs": [],
"source": [
"import cmweather # noqa\n",
"import matplotlib.pyplot as plt\n",
"from open_radar_data import DATASETS\n",
"\n",
"import xradar as xd"
]
},
{
"cell_type": "markdown",
"id": "de72bd74",
"metadata": {},
"source": [
"## Read a sample radar file"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ddec4345",
"metadata": {},
"outputs": [],
"source": [
"filename = DATASETS.fetch(\"cfrad.20080604_002217_000_SPOL_v36_SUR.nc\")\n",
"radar = xd.io.open_cfradial1_datatree(filename, first_dim=\"auto\")\n",
"radar"
]
},
{
"cell_type": "markdown",
"id": "e87b83a0",
"metadata": {},
"source": [
"## 1. Default georeferencing (AEQD, meters)\n",
"\n",
"Without ``target_crs``, the accessor adds ``x``, ``y``, ``z`` in the\n",
"radar-native AEQD projection (distances in meters from the radar site)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "61077c76",
"metadata": {},
"outputs": [],
"source": [
"radar_aeqd = radar.xradar.georeference()\n",
"sweep = radar_aeqd[\"sweep_0\"].to_dataset(inherit=\"all_coords\")\n",
"crs = sweep.xradar.get_crs()\n",
"print(\n",
" \"projection:\",\n",
" crs.coordinate_operation.method_name if crs.coordinate_operation else crs.name,\n",
")\n",
"print(\"is_projected / is_geographic:\", crs.is_projected, \"/\", crs.is_geographic)\n",
"print(\"x units:\", sweep.x.attrs.get(\"units\"))\n",
"print(\"x range:\", float(sweep.x.min()), float(sweep.x.max()))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ce9c5d2",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(7, 6))\n",
"radar_aeqd[\"sweep_0\"][\"DBZ\"].plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", ax=ax)\n",
"ax.set_title(\"AEQD (default) — x, y in meters from radar\")\n",
"ax.set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"id": "103f0bf4",
"metadata": {},
"source": [
"## 2. Reproject to geographic lon/lat (EPSG:4326)\n",
"\n",
"Pass ``target_crs=4326`` to get ``x`` as longitude and ``y`` as latitude.\n",
"Attributes are updated automatically (``standard_name``, ``units``)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d2ea627c",
"metadata": {},
"outputs": [],
"source": [
"radar_geo = radar.xradar.georeference(target_crs=4326)\n",
"sweep = radar_geo[\"sweep_0\"].to_dataset(inherit=\"all_coords\")\n",
"crs = sweep.xradar.get_crs()\n",
"print(\"CRS name:\", crs.name, \"| EPSG:\", crs.to_epsg())\n",
"print(\"x attrs:\", dict(sweep.x.attrs))\n",
"print(\"y attrs:\", dict(sweep.y.attrs))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6ce55e7",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(8, 6))\n",
"radar_geo[\"sweep_0\"][\"DBZ\"].plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", ax=ax)\n",
"ax.set_xlabel(\"Longitude\")\n",
"ax.set_ylabel(\"Latitude\")\n",
"ax.set_title(\"EPSG:4326 — x=lon, y=lat\")"
]
},
{
"cell_type": "markdown",
"id": "9adc842c",
"metadata": {},
"source": [
"## 3. Reproject to a projected CRS (UTM, Web Mercator)\n",
"\n",
"Any ``pyproj``-compatible CRS input works: EPSG code, WKT string, or a\n",
"``pyproj.CRS`` instance."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b68b43c6",
"metadata": {},
"outputs": [],
"source": [
"# Web Mercator (meters)\n",
"radar_wm = radar.xradar.georeference(target_crs=\"EPSG:3857\")\n",
"sweep_wm = radar_wm[\"sweep_0\"].to_dataset(inherit=\"all_coords\")\n",
"crs = sweep_wm.xradar.get_crs()\n",
"print(\"CRS:\", crs.name, \"| EPSG:\", crs.to_epsg())"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fa85ef99",
"metadata": {},
"outputs": [],
"source": [
"import pyproj\n",
"\n",
"utm_crs = pyproj.CRS(\"EPSG:32633\") # UTM zone 33N (tweak for your radar site)\n",
"radar_utm = radar.xradar.georeference(target_crs=utm_crs)\n",
"sweep_utm = radar_utm[\"sweep_0\"].to_dataset(inherit=\"all_coords\")\n",
"crs = sweep_utm.xradar.get_crs()\n",
"print(\"CRS:\", crs.name, \"| EPSG:\", crs.to_epsg())"
]
},
{
"cell_type": "markdown",
"id": "342b18be",
"metadata": {},
"source": [
"## 4. Plot on a cartopy map\n",
"\n",
"Once data is in a known CRS, hand it off to cartopy for map plotting."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "316b50c8",
"metadata": {},
"outputs": [],
"source": [
"import cartopy.crs as ccrs\n",
"\n",
"radar_geo = radar.xradar.georeference(target_crs=4326)\n",
"\n",
"fig = plt.figure(figsize=(9, 8))\n",
"ax = fig.add_subplot(111, projection=ccrs.PlateCarree())\n",
"radar_geo[\"sweep_0\"][\"DBZ\"].plot(\n",
" x=\"x\",\n",
" y=\"y\",\n",
" cmap=\"ChaseSpectral\",\n",
" transform=ccrs.PlateCarree(),\n",
" ax=ax,\n",
" cbar_kwargs=dict(pad=0.05, shrink=0.7),\n",
")\n",
"ax.coastlines()\n",
"ax.gridlines(draw_labels=True)\n",
"ax.set_title(\"Georeferenced PPI on PlateCarree map\")"
]
},
{
"cell_type": "markdown",
"id": "d49ee5c1",
"metadata": {},
"source": [
"## 5. Works on a single Dataset too\n",
"\n",
"The accessor is available on ``Dataset`` and ``DataArray`` as well, not just\n",
"``DataTree``."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "014505ea",
"metadata": {},
"outputs": [],
"source": [
"ds = radar[\"sweep_0\"].to_dataset(inherit=\"all_coords\")\n",
"ds_geo = ds.xradar.georeference(target_crs=4326)\n",
"ds_geo[[\"DBZ\"]]"
]
},
{
"cell_type": "markdown",
"id": "7852a0d2",
"metadata": {},
"source": [
"## Notes\n",
"\n",
"* ``target_crs`` accepts anything ``pyproj.CRS(...)`` accepts: int\n",
" [EPSG codes](https://epsg.io/), ``\"EPSG:xxxx\"`` strings, WKT, or\n",
" ``pyproj.CRS`` instances.\n",
"* The ``spatial_ref`` / ``crs_wkt`` coordinate is updated to reflect the\n",
" target CRS — ``sweep.xradar.get_crs()`` will return it.\n",
"* Currently only ``x`` and ``y`` are reprojected; ``z`` stays as AEQD\n",
" altitude above the radar site.\n",
"* If you want **separate** ``lon`` / ``lat`` / ``alt`` coordinates (without\n",
" overwriting ``x, y``), see the [Assign_GeoCoords](Assign_GeoCoords.md) notebook."
]
}
],
"metadata": {
"jupytext": {
"main_language": "python"
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
1 change: 1 addition & 0 deletions docs/usage.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
21 changes: 21 additions & 0 deletions tests/test_accessors.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading
Loading