{ "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": [ "
Important Note:
\n", "\n",
" This radar data is from a mobile research radar called Doppler on Wheels (DOW), and its site coordinates (latitude, longitude) \n",
" often vary slightly during operation, as can be seen from the shape ((731)
) of the data in the above cell, while we expect it to \n",
" be of unity shapes or empty, i.e., (1)
or ()
. As a result, multiple site coordinate values can exist, creating a challenge for assigning \n",
" consistent x, y, z
or lat, lon, and alt
coordinates using the current georeferencing system in xradar
. \n",
" To address this, a custom function like fix_sitecoords
(defined above) can be created, leveraging the map_over_sweeps
\n",
" function to standardize the site coordinates.\n",
"