ADD: Example for GeoCoords (#251)

openradar · Nov 28, 2024 · 0ea8d58 · 0ea8d58
1 parent 836a4c5
commit 0ea8d58
Show file tree

Hide file tree

Showing 3 changed files with 363 additions and 0 deletions.
diff --git a/docs/history.md b/docs/history.md
@@ -7,6 +7,7 @@
 * FIX: Keeping attributes at each variable when using `open_nexradlevel2_datatree`. ({issue}`241`) ({pull}`242`) by [@aladinor](https://github.com/aladinor)
 * FIX: Correctly read transition rays in RHI scans ({issue}`247`) ({pull}`250`) by [@rcjackson](https://github.com/rcjackson)
 * FIX: Correctly open NEXRAD files when split cut mode is enable ({issue} `245`) ({pull}`246`) by [@aladinor](https://github.com/aladinor)
+* ADD: Example Notebook for assigning geocoords. ({issue}`243`) and ({pull}`251`) by [@syedhamidali](https://github.com/syedhamidali)
 
 ## 0.8.0 (2024-11-04)
 

diff --git a/docs/usage.md b/docs/usage.md
@@ -41,6 +41,7 @@ notebooks/NexradLevel2
 :maxdepth: 2
 :caption: Examples
 
+notebooks/Assign_GeoCoords
 notebooks/CfRadial1_Export
 notebooks/Read-plot-Sigmet-data-from-AWS
 notebooks/plot-ppi

diff --git a/examples/notebooks/Assign_GeoCoords.ipynb b/examples/notebooks/Assign_GeoCoords.ipynb
@@ -0,0 +1,361 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# Assign GeoCoords"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1",
+   "metadata": {},
+   "source": [
+    "This notebook demonstrates how to work with geocoordinates (longitude and latitude) instead of radar-centric x and y coordinates."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cmweather  # noqa: F401\n",
+    "import matplotlib.pyplot as plt\n",
+    "import xarray as xr\n",
+    "from open_radar_data import DATASETS\n",
+    "\n",
+    "import xradar as xd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3",
+   "metadata": {},
+   "source": [
+    "## Define Functions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_geocoords(ds):\n",
+    "    \"\"\"\n",
+    "    Converts Cartesian coordinates (x, y, z) in a radar dataset to geographic\n",
+    "    coordinates (longitude, latitude, altitude) using CRS transformation.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    ds : xarray.Dataset\n",
+    "        Radar dataset with Cartesian coordinates.\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    xarray.Dataset\n",
+    "        Dataset with added 'lon', 'lat', and 'alt' coordinates and their attributes.\n",
+    "    \"\"\"\n",
+    "    from pyproj import CRS, Transformer\n",
+    "\n",
+    "    # Convert the dataset to georeferenced coordinates\n",
+    "    ds = ds.xradar.georeference()\n",
+    "    # Define source and target coordinate reference systems (CRS)\n",
+    "    src_crs = ds.xradar.get_crs()\n",
+    "    trg_crs = CRS.from_user_input(4326)  # EPSG:4326 (WGS 84)\n",
+    "    # Create a transformer for coordinate conversion\n",
+    "    transformer = Transformer.from_crs(src_crs, trg_crs)\n",
+    "    # Transform x, y, z coordinates to latitude, longitude, and altitude\n",
+    "    trg_y, trg_x, trg_z = transformer.transform(ds.x, ds.y, ds.z)\n",
+    "    # Assign new coordinates with appropriate attributes\n",
+    "    ds = ds.assign_coords(\n",
+    "        {\n",
+    "            \"lon\": (ds.x.dims, trg_x, xd.model.get_longitude_attrs()),\n",
+    "            \"lat\": (ds.y.dims, trg_y, xd.model.get_latitude_attrs()),\n",
+    "            \"alt\": (ds.z.dims, trg_z, xd.model.get_altitude_attrs()),\n",
+    "        }\n",
+    "    )\n",
+    "    return ds\n",
+    "\n",
+    "\n",
+    "def fix_sitecoords(ds):\n",
+    "    coords = [\"longitude\", \"latitude\", \"altitude\", \"altitude_agl\"]\n",
+    "    for coord in coords:\n",
+    "        # Compute median excluding NaN\n",
+    "        data = ds[coord].median(skipna=True).item()\n",
+    "        attrs = ds[coord].attrs if coord in ds else {}\n",
+    "        ds = ds.assign_coords({coord: xr.DataArray(data=data, attrs=attrs)})\n",
+    "    return ds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5",
+   "metadata": {},
+   "source": [
+    "## Load Data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "file1 = DATASETS.fetch(\"cfrad.20080604_002217_000_SPOL_v36_SUR.nc\")\n",
+    "file2 = DATASETS.fetch(\"cfrad.20211011_201557.188_to_20211011_201617.720_DOW8_PPI.nc\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "## Example #1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtree1 = xd.io.open_cfradial1_datatree(file1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtree1 = dtree1.xradar.georeference()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "display(dtree1[\"sweep_0\"].ds)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11",
+   "metadata": {},
+   "source": [
+    "## Assign lat, lon, and alt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtree1 = dtree1.xradar.map_over_sweeps(get_geocoords)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "display(dtree1[\"sweep_0\"].ds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = dtree1[\"sweep_0\"].to_dataset()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(1, 2, figsize=(11, 4))\n",
+    "ds[\"DBZ\"].plot(x=\"x\", y=\"y\", vmin=-10, vmax=75, cmap=\"HomeyerRainbow\", ax=ax[0])\n",
+    "\n",
+    "ds[\"DBZ\"].plot(x=\"lon\", y=\"lat\", vmin=-10, vmax=75, cmap=\"HomeyerRainbow\", ax=ax[1])\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "16",
+   "metadata": {},
+   "source": [
+    "## Example #2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtree2 = xd.io.open_cfradial1_datatree(file2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "18",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    dtree2 = dtree2.xradar.georeference()\n",
+    "except Exception:\n",
+    "    print(\"Georeferencing failed!\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Longitude\", dtree2[\"sweep_0\"][\"longitude\"].shape)\n",
+    "print(\"Latitude\", dtree2[\"sweep_0\"][\"latitude\"].shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20",
+   "metadata": {},
+   "source": [
+    "<div class=\"alert alert-info\">\n",
+    "    <p style=\"font-weight:bold; margin:0;\">Important Note:</p>\n",
+    "    <p>\n",
+    "        This radar data is from a mobile research radar called <b>Doppler on Wheels (DOW)</b>, and its site coordinates (latitude, longitude) \n",
+    "        often vary slightly during operation, as can be seen from the shape (<code>(731)</code>) of the data in the above cell, while we expect it to \n",
+    "        be of unity shapes or empty, i.e., <code>(1)</code> or <code>()</code>. As a result, multiple site coordinate values can exist, creating a challenge for assigning \n",
+    "        consistent <code>x, y, z</code> or <code>lat, lon, and alt</code> coordinates using the current georeferencing system in <code>xradar</code>. \n",
+    "        To address this, a custom function like <code>fix_sitecoords</code> (defined above) can be created, leveraging the <code>map_over_sweeps</code> \n",
+    "        function to standardize the site coordinates.\n",
+    "    </p>\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "21",
+   "metadata": {},
+   "source": [
+    "## Fix Coords"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtree2 = dtree2.xradar.map_over_sweeps(fix_sitecoords)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "23",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dtree2 = dtree2.xradar.map_over_sweeps(get_geocoords)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "24",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "display(dtree2[\"sweep_0\"].ds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "25",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = dtree2[\"sweep_0\"].to_dataset()\n",
+    "ref = ds.where(ds.DBZHC >= 5)[\"DBZHC\"]\n",
+    "\n",
+    "fig = plt.figure(figsize=(6, 5))\n",
+    "ax = plt.axes()\n",
+    "pl = ref.plot(\n",
+    "    x=\"lon\",\n",
+    "    y=\"lat\",\n",
+    "    vmin=-10,\n",
+    "    vmax=70,\n",
+    "    cmap=\"HomeyerRainbow\",\n",
+    "    ax=ax,\n",
+    ")\n",
+    "\n",
+    "ax.minorticks_on()\n",
+    "ax.grid(ls=\"--\", lw=0.5)\n",
+    "ax.set_aspect(\"auto\")\n",
+    "title = (\n",
+    "    dtree2.attrs[\"instrument_name\"]\n",
+    "    + \" \"\n",
+    "    + str(ds.time.mean().dt.strftime(\"%Y-%m-%d %H:%M:%S\").values)\n",
+    ")\n",
+    "ax.set_title(title)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "26",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "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.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}