{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\"wradlib" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# wradlib - clutter and beamblockage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview\n", "\n", "Within this notebook, we will cover:\n", "\n", "1. Reading data using xradar\n", "1. Clutter detection\n", "1. Beam Blockage calculation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "| Concepts | Importance | Notes |\n", "| --- | --- | --- |\n", "| [Xarray Basics](https://tutorial.xarray.dev/intro.html) | Helpful | Basic Dataset/DataArray |\n", "| [Matplotlib Basics](https://foundations.projectpythia.org/core/matplotlib/matplotlib-basics.html) | Helpful | Basic Plotting |\n", "| [Intro to Cartopy](https://foundations.projectpythia.org/core/cartopy/cartopy.html) | Helpful | Projections |\n", "\n", "- **Time to learn**: 10 minutes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import wradlib as wrl\n", "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "import cartopy\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "import xradar as xd\n", "import hvplot\n", "import hvplot.xarray" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## retrieve data from s3 bucket\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import urllib.request\n", "from pathlib import Path\n", "\n", "# Set the URL for the cloud\n", "URL = \"https://js2.jetstream-cloud.org:8001/\"\n", "path = \"pythia/radar/erad2024\"\n", "!mkdir -p data\n", "files = [\n", " \"20240522_MeteoSwiss_ARPA_Lombardia/Data/Xband/DES_VOL_RAW_20240522_1600.nc\",\n", " \"wradlib/desio_dem.tif\",\n", "]\n", "for file in files:\n", " file_remote = os.path.join(path, file)\n", " file_local = os.path.join(\"data\", Path(file).name)\n", " if not os.path.exists(file_local):\n", " print(f\"downloading, {file_local}\")\n", " urllib.request.urlretrieve(f\"{URL}{file_remote}\", file_local)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Open CfRadial1 Volume\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "reindex = dict(angle_res=1, direction=1, start_angle=0, stop_angle=360)\n", "dtree = xd.io.open_cfradial1_datatree(\"data/DES_VOL_RAW_20240522_1600.nc\")\n", "display(dtree.load())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Get first sweep" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "swp = (\n", " dtree[\"sweep_0\"]\n", " .to_dataset()\n", " .wrl.georef.georeference(crs=wrl.georef.get_earth_projection())\n", " .set_coords(\"sweep_mode\")\n", ")\n", "swp.x.attrs = xd.model.get_longitude_attrs()\n", "swp.y.attrs = xd.model.get_latitude_attrs()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "display(swp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get Digital Elevation Map (DEM)\n", "\n", "If we have access to the NASA EarthData GESDISC, we can use the BearerToken to retrieve SRTM data corresponding to the actual radar domain. Or we can choose the precompiled GeoTiff." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extent = [swp.x.min().values, swp.x.max().values, swp.y.min().values, swp.y.max().values]\n", "# import os\n", "# os.environ[\"WRADLIB_EARTHDATA_BEARER_TOKEN\"] = \"\"\n", "# dem = wrl.io.get_srtm(extent)\n", "# wrl.io.write_raster_dataset(\"desio_dem.tif\", dem)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dem = (\n", " xr.open_dataset(\"data/desio_dem.tif\", engine=\"rasterio\")\n", " .isel(band=0)\n", " .rename(band_data=\"DEM\")\n", " .reset_coords(\"band\", drop=True)\n", ")\n", "display(dem)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Extract radar parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "radar_parameters = dtree[\"radar_parameters\"]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bw = radar_parameters[\"radar_beam_width_h\"]\n", "bw" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prepare DEM for Polar Processing\n", "\n", "Here the power of [xr.apply_ufunc](https://docs.xarray.dev/en/stable/generated/xarray.apply_ufunc.html) is shown, a wrapper to xarray-ify numpy functions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def interpolate_dem(obj, dem, **kwargs):\n", " dim0 = obj.wrl.util.dim0()\n", "\n", " def wrapper(sx, sy, dx, dy, dem, *args, **kwargs):\n", " y, x = np.meshgrid(dy, dx)\n", " rastercoords = np.dstack([x, y])\n", " polcoords = np.dstack([sx, sy])\n", " return wrl.ipol.cart_to_irregular_spline(rastercoords, dem, polcoords, **kwargs)\n", "\n", " out = xr.apply_ufunc(\n", " wrapper,\n", " obj.x,\n", " obj.y,\n", " dem.x,\n", " dem.y,\n", " dem,\n", " input_core_dims=[[dim0, \"range\"], [dim0, \"range\"], [\"x\"], [\"y\"], [\"y\", \"x\"]],\n", " output_core_dims=[[dim0, \"range\"]],\n", " dask=\"parallelized\",\n", " vectorize=True,\n", " kwargs=kwargs,\n", " dask_gufunc_kwargs=dict(allow_rechunk=True),\n", " )\n", " out.name = \"DEM\"\n", " return obj.assign(DEM=out)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "swp = interpolate_dem(swp, dem.DEM, order=3, prefilter=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "swp.DEM.wrl.vis.plot(cmap=\"terrain\", vmin=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot scan strategy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nrays = swp.azimuth.size\n", "nbins = swp.range.size\n", "range_res = 300.0\n", "ranges = np.arange(nbins) * range_res\n", "elevs = dtree.root.sweep_fixed_angle.values\n", "sitecoords = (\n", " dtree.root.longitude.values.item(),\n", " dtree.root.latitude.values.item(),\n", " dtree.root.altitude.values.item(),\n", ")\n", "\n", "ax = wrl.vis.plot_scan_strategy(\n", " ranges,\n", " elevs,\n", " sitecoords,\n", " beamwidth=radar_parameters[\"radar_beam_width_h\"].values,\n", " terrain=None,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use `terrain=swp.DEM.sel(azimuth=0, method=\"nearest\")` to get some arbitrary ray." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ax = wrl.vis.plot_scan_strategy(\n", " ranges,\n", " elevs,\n", " sitecoords,\n", " beamwidth=radar_parameters[\"radar_beam_width_h\"].values,\n", " terrain=swp.DEM.sel(azimuth=0, method=\"nearest\"),\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate clutter map" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "clmap = swp.DBZ_TOT.wrl.classify.filter_gabella(\n", " wsize=5,\n", " thrsnorain=0.0,\n", " tr1=21.0, # 21.,\n", " n_p=23.0, # 23,\n", " tr2=1.3,\n", " rm_nans=False,\n", ")\n", "swp = swp.assign({\"CMAP\": clmap})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Reflectivities, Clutter and Cluttermap" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(15, 12))\n", "ax1 = fig.add_subplot(221)\n", "from osgeo import osr\n", "\n", "wgs84 = osr.SpatialReference()\n", "wgs84.ImportFromEPSG(4326)\n", "# swp = swp.sel(range=slice(0, 100000)).set_coords(\"sweep_mode\").wrl.georef.georeference(crs=wgs84)\n", "swp.DBZ_TOT.plot(x=\"x\", y=\"y\", ax=ax1, vmin=0, vmax=60)\n", "ax1.set_title(\"Reflectivity raw\")\n", "ax2 = fig.add_subplot(222)\n", "swp.CMAP.plot(x=\"x\", y=\"y\", ax=ax2)\n", "ax2.set_title(\"Cluttermap\")\n", "ax3 = fig.add_subplot(223)\n", "swp.DBZ_TOT.where(swp.CMAP == 1).plot(x=\"x\", y=\"y\", ax=ax3, vmin=0, vmax=60)\n", "ax3.set_title(\"Clutter\")\n", "ax4 = fig.add_subplot(224)\n", "swp.DBZ_TOT.where(swp.CMAP < 1).plot(x=\"x\", y=\"y\", ax=ax4, vmin=0, vmax=60)\n", "ax4.set_title(\"Reflectivity clutter removed\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Compare with corrected reflectivity from signal processor \n", "\n", "plus additional simple RHOHV filter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(15, 6))\n", "ax1 = fig.add_subplot(121)\n", "swp.DBZ.plot(x=\"x\", y=\"y\", ax=ax1, vmin=0, vmax=60)\n", "ax1.set_title(\"Reflectivity corr\")\n", "ax2 = fig.add_subplot(122)\n", "swp.DBZ_TOT.where((swp.CMAP < 1) & (swp.RHOHV >= 0.8)).plot(\n", " x=\"x\", y=\"y\", ax=ax2, vmin=0, vmax=60\n", ")\n", "ax2.set_title(\"Reflectivity clutter removed\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(8, 5))\n", "ax1 = fig.add_subplot(111)\n", "swp.CMAP.where(swp.CMAP == 1).plot(x=\"x\", y=\"y\", vmin=0, vmax=1, cmap=\"turbo\")\n", "ax1.set_title(\"Reflectivity corr\")\n", "dem.DEM.plot(ax=ax1, zorder=-2, cmap=\"terrain\", vmin=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Use hvplot for zooming and panning" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to rechunk the coordinates as hvplot needs chunked variables and coords." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cl = (\n", " swp.CMAP.where(swp.CMAP == 1)\n", " .chunk(chunks={})\n", " .hvplot.quadmesh(\n", " x=\"x\", y=\"y\", cmap=\"turbo\", width=600, height=500, clim=(0, 1), rasterize=True\n", " )\n", ")\n", "dm = dem.DEM.chunk(chunks={}).hvplot(\n", " x=\"x\", y=\"y\", cmap=\"terrain\", width=600, height=500, rasterize=True\n", ")\n", "dm * cl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## BeamBlockage Calculation\n", "\n", "Can you xarray-ify the following, too?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "beamradius = wrl.util.half_power_radius(swp.range, bw)\n", "PBB = wrl.qual.beam_block_frac(swp.DEM.values, swp.z.values, beamradius)\n", "PBB = np.ma.masked_invalid(PBB)\n", "CBB = wrl.qual.cum_beam_block_frac(PBB)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "swp = swp.assign(\n", " CBB=([\"azimuth\", \"range\"], CBB),\n", " PBB=([\"azimuth\", \"range\"], PBB),\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# just a little helper function to style x and y axes of our maps\n", "def annotate_map(ax, cm=None, title=\"\"):\n", " # ticks = (ax.get_xticks() / 1000).astype(int)\n", " # ax.set_xticklabels(ticks)\n", " # ticks = (ax.get_yticks() / 1000).astype(int)\n", " # ax.set_yticklabels(ticks)\n", " # ax.set_xlabel(\"Kilometers\")\n", " # ax.set_ylabel(\"Kilometers\")\n", " if not cm is None:\n", " plt.colorbar(cm, ax=ax)\n", " if not title == \"\":\n", " ax.set_title(title)\n", " ax.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib as mpl\n", "\n", "sitecoords = (swp.longitude.values, swp.latitude.values, swp.altitude.values)\n", "r = swp.range.values\n", "az = swp.azimuth.values\n", "\n", "alt = swp.z.values\n", "fig = plt.figure(figsize=(15, 12))\n", "\n", "# create subplots\n", "ax1 = plt.subplot2grid((2, 2), (0, 0))\n", "ax2 = plt.subplot2grid((2, 2), (0, 1))\n", "ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)\n", "\n", "# azimuth angle\n", "angle = 0\n", "\n", "# Plot terrain (on ax1)\n", "dem_pm = swp.DEM.wrl.vis.plot(ax=ax1, cmap=mpl.cm.terrain, vmin=0.0, add_colorbar=False)\n", "swp.sel(azimuth=0, method=\"nearest\").plot.scatter(\n", " x=\"x\", y=\"y\", marker=\".\", color=\"r\", alpha=0.2, ax=ax1\n", ")\n", "ax1.plot(swp.longitude.values, swp.latitude.values, \"ro\")\n", "annotate_map(\n", " ax1,\n", " dem_pm,\n", " \"Terrain within {0} km range\".format(np.max(swp.range.values / 1000.0) + 0.1),\n", ")\n", "\n", "# Plot CBB (on ax2)\n", "cbb = swp.CBB.wrl.vis.plot(ax=ax2, cmap=mpl.cm.PuRd, vmin=0, vmax=1, add_colorbar=False)\n", "annotate_map(ax2, cbb, \"Beam-Blockage Fraction\")\n", "\n", "# Plot single ray terrain profile on ax3\n", "(bc,) = ax3.plot(\n", " swp.range / 1000.0, swp.z[angle, :], \"-b\", linewidth=3, label=\"Beam Center\"\n", ")\n", "(b3db,) = ax3.plot(\n", " swp.range.values / 1000.0,\n", " (swp.z[angle, :] + beamradius),\n", " \":b\",\n", " linewidth=1.5,\n", " label=\"3 dB Beam width\",\n", ")\n", "ax3.plot(swp.range / 1000.0, (swp.z[angle, :] - beamradius), \":b\")\n", "ax3.fill_between(swp.range / 1000.0, 0.0, swp.DEM[angle, :], color=\"0.75\")\n", "ax3.set_xlim(0.0, np.max(swp.range / 1000.0) + 0.1)\n", "ax3.set_ylim(0.0, 3000)\n", "ax3.set_xlabel(\"Range (km)\")\n", "ax3.set_ylabel(\"Altitude (m)\")\n", "ax3.grid()\n", "\n", "axb = ax3.twinx()\n", "(bbf,) = axb.plot(swp.range / 1000.0, swp.CBB[angle, :], \"-k\", label=\"BBF\")\n", "axb.set_ylabel(\"Beam-blockage fraction\")\n", "axb.set_ylim(0.0, 1.0)\n", "axb.set_xlim(0.0, np.max(swp.range / 1000.0) + 0.1)\n", "\n", "\n", "legend = ax3.legend(\n", " (bc, b3db, bbf),\n", " (\"Beam Center\", \"3 dB Beam width\", \"BBF\"),\n", " loc=\"upper left\",\n", " fontsize=10,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some EyeCandy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def height_formatter(x, pos):\n", " x = (x - 6370000) / 1000\n", " fmt_str = \"{:g}\".format(x)\n", " return fmt_str\n", "\n", "\n", "def range_formatter(x, pos):\n", " x = x / 1000.0\n", " fmt_str = \"{:g}\".format(x)\n", " return fmt_str" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 6))\n", "\n", "cgax, caax, paax = wrl.vis.create_cg(fig=fig, rot=0, scale=1)\n", "\n", "# azimuth angle\n", "angle = 0\n", "\n", "# fix grid_helper\n", "er = 6370000\n", "gh = cgax.get_grid_helper()\n", "gh.grid_finder.grid_locator2._nbins = 80\n", "gh.grid_finder.grid_locator2._steps = [1, 2, 4, 5, 10]\n", "\n", "# calculate beam_height and arc_distance for ke=1\n", "# means line of sight\n", "bhe = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=1.0)\n", "ade = wrl.georef.bin_distance(r, 0, sitecoords[2], re=er, ke=1.0)\n", "nn0 = np.zeros_like(r)\n", "# for nice plotting we assume earth_radius = 6370000 m\n", "ecp = nn0 + er\n", "# theta (arc_distance sector angle)\n", "thetap = -np.degrees(ade / er) + 90.0\n", "\n", "# zero degree elevation with standard refraction\n", "bh0 = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er)\n", "\n", "# plot (ecp is earth surface normal null)\n", "(bes,) = paax.plot(thetap, ecp, \"-k\", linewidth=3, label=\"Earth Surface NN\")\n", "(bc,) = paax.plot(thetap, ecp + alt[angle, :], \"-b\", linewidth=3, label=\"Beam Center\")\n", "(bc0r,) = paax.plot(thetap, ecp + bh0, \"-g\", label=\"0 deg Refraction\")\n", "(bc0n,) = paax.plot(thetap, ecp + bhe, \"-r\", label=\"0 deg line of sight\")\n", "(b3db,) = paax.plot(\n", " thetap, ecp + alt[angle, :] + beamradius, \":b\", label=\"+3 dB Beam width\"\n", ")\n", "paax.plot(thetap, ecp + alt[angle, :] - beamradius, \":b\", label=\"-3 dB Beam width\")\n", "\n", "# orography\n", "paax.fill_between(thetap, ecp, ecp + swp.DEM[angle, :], color=\"0.75\")\n", "\n", "# shape axes\n", "cgax.set_xlim(0, np.max(ade))\n", "cgax.set_ylim([ecp.min() - 1000, ecp.max() + 2500])\n", "caax.grid(True, axis=\"x\")\n", "cgax.grid(True, axis=\"y\")\n", "cgax.axis[\"top\"].toggle(all=False)\n", "caax.yaxis.set_major_locator(\n", " mpl.ticker.MaxNLocator(steps=[1, 2, 4, 5, 10], nbins=20, prune=\"both\")\n", ")\n", "caax.xaxis.set_major_locator(mpl.ticker.MaxNLocator())\n", "caax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(height_formatter))\n", "caax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(range_formatter))\n", "\n", "caax.set_xlabel(\"Range (km)\")\n", "caax.set_ylabel(\"Altitude (km)\")\n", "\n", "legend = paax.legend(\n", " (bes, bc0n, bc0r, bc, b3db),\n", " (\n", " \"Earth Surface NN\",\n", " \"0 deg line of sight\",\n", " \"0 deg std refraction\",\n", " \"Beam Center\",\n", " \"3 dB Beam width\",\n", " ),\n", " loc=\"lower left\",\n", " fontsize=10,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "We've just learned how to use $\\omega radlib$'s Gabella clutter detection for single sweeps. We've looked into digital elevation maps and beam blockage calculations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Resources and references\n", "\n", "- [xarray](https://docs.xarray.dev)\n", "- [dask](https://docs.dask.org/)\n", "- [GDAL](https://gdal.org)\n", "- [xradar backends](https://docs.openradarscience.org/projects/xradar/en/stable/importers.html)\n", "- [CfRadial1](https://ncar.github.io/CfRadial/)" ] } ], "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.4" }, "nbdime-conflicts": { "local_diff": [ { "diff": [ { "diff": [ { "key": 0, "op": "addrange", "valuelist": [ "Python 3" ] }, { "key": 0, "length": 1, "op": "removerange" } ], "key": "display_name", "op": "patch" } ], "key": "kernelspec", "op": "patch" } ], "remote_diff": [ { "diff": [ { "diff": [ { "key": 0, "op": "addrange", "valuelist": [ "Python3" ] }, { "key": 0, "length": 1, "op": "removerange" } ], "key": "display_name", "op": "patch" } ], "key": "kernelspec", "op": "patch" } ] }, "toc-autonumbering": false }, "nbformat": 4, "nbformat_minor": 4 }