diff --git a/examples/plotting/lonboard_nexrad_reflectivity.ipynb b/examples/plotting/lonboard_nexrad_reflectivity.ipynb new file mode 100644 index 0000000000..6f25917cf0 --- /dev/null +++ b/examples/plotting/lonboard_nexrad_reflectivity.ipynb @@ -0,0 +1,183 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "====================================\n", + "Create a lonboard map of NEXRAD reflectivity\n", + "====================================\n", + "\n", + "An example which follows the \"Create a plot of NEXRAD reflectivity\", but plots with the\n", + "lonboard deck.gl python interface.\n", + "\n", + "\"\"\"\n", + "\n", + "print(__doc__)\n", + "\n", + "# Author: Sam Gardner (samuel.gardner@ttu.edu)\n", + "# License: BSD 3 clause\n", + "\n", + "import cmweather\n", + "import numpy as np\n", + "from arro3.core import Array, Field, Schema, Table, fixed_size_list_array, list_array\n", + "from lonboard import Map, SolidPolygonLayer\n", + "from lonboard.colormap import apply_continuous_cmap\n", + "\n", + "import pyart\n", + "from pyart.testing import get_test_data\n", + "\n", + "# download and open test data file\n", + "filename = get_test_data(\"Level2_KATX_20130717_1950.ar2v\")\n", + "radar = pyart.io.read_nexrad_archive(filename)\n", + "# extract the first sweep of data\n", + "first_sweep = radar.extract_sweeps([0])" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# extract the range, azimuth, and elevation of gates in the first sweep\n", + "range_center = first_sweep.range[\"data\"]\n", + "az_center = first_sweep.azimuth[\"data\"]\n", + "elevation = first_sweep.elevation[\"data\"]\n", + "# Convert the gate center coordinates to latitude and longitude edges\n", + "x_edge, y_edge, _ = pyart.core.antenna_vectors_to_cartesian(\n", + " range_center, az_center, elevation, edges=True\n", + ")\n", + "radar_lat = radar.latitude[\"data\"][0]\n", + "radar_lon = radar.longitude[\"data\"][0]\n", + "projparams = {\"proj\": \"pyart_aeqd\", \"lon_0\": radar_lon, \"lat_0\": radar_lat}\n", + "lons, lats = pyart.core.cartesian_to_geographic(x_edge, y_edge, projparams)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Extract reflectivity data from the first sweep and fill masked values with NaN\n", + "reflectivity = first_sweep.fields[\"reflectivity\"][\"data\"].filled(np.nan)\n", + "# lonboard colormaps data over a range of 0-1, so we normalize the reflectivity data by our vmin and vmax\n", + "# The ChaseSpectral colormap aligns lightness breaks with relevant meteorological phenomena\n", + "# when used with a vmin of -10 dBZ and a vmax of 80 dBZ\n", + "chasespectral_vmin = -10\n", + "chasespectral_vmax = 80\n", + "# Scale the reflectivity data to the range 0-1\n", + "reflectivity_normalized = np.ravel(\n", + " (reflectivity - chasespectral_vmin) / (chasespectral_vmax - chasespectral_vmin)\n", + ")\n", + "# Plotting takes much less time if data below the vmin and NaN values are not plotted\n", + "reflectivity_normalized[reflectivity_normalized < 0] = np.nan\n", + "points_to_keep = ~np.isnan(reflectivity_normalized)\n", + "reflectivity_normalized = reflectivity_normalized[points_to_keep]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# This function gets the corners of each radar gate and returns them as a geoarrow table\n", + "# Geoarrow tables are a way to represent geospatial data in Apache Arrow, which is used by lonboard\n", + "# for rendering geospatial data. This function originally written by Kyle Barron\n", + "# See https://github.com/developmentseed/lonboard/discussions/643\n", + "def lon_lat_to_arrow(lon_edge, lat_edge, data_mask):\n", + " # bottom-left: (i, j)\n", + " bl_lon = lon_edge[:-1, :-1]\n", + " bl_lat = lat_edge[:-1, :-1]\n", + "\n", + " # bottom-right: (i+1, j)\n", + " br_lon = lon_edge[1:, :-1]\n", + " br_lat = lat_edge[1:, :-1]\n", + "\n", + " # top-right: (i+1, j+1)\n", + " tr_lon = lon_edge[1:, 1:]\n", + " tr_lat = lat_edge[1:, 1:]\n", + "\n", + " # top-left: (i, j+1)\n", + " tl_lon = lon_edge[:-1, 1:]\n", + " tl_lat = lat_edge[:-1, 1:]\n", + "\n", + " lons_to_poly = np.array(\n", + " [bl_lon.flatten(), br_lon.flatten(), tr_lon.flatten(), tl_lon.flatten()]\n", + " )\n", + " lats_to_poly = np.array(\n", + " [bl_lat.flatten(), br_lat.flatten(), tr_lat.flatten(), tl_lat.flatten()]\n", + " )\n", + " poly_coords = np.stack([lons_to_poly, lats_to_poly], axis=-1).transpose(1, 0, 2)\n", + " poly_coords = poly_coords[data_mask, :, :]\n", + "\n", + " poly_coords_closed = np.concat((poly_coords, poly_coords[:, 0:1, :]), axis=1)\n", + "\n", + " poly_coords_arrow = fixed_size_list_array(\n", + " Array.from_numpy(np.ravel(poly_coords_closed)), 2\n", + " )\n", + " ring_offsets = Array.from_numpy(\n", + " np.arange(0, (poly_coords_closed.shape[0] + 1) * 5, 5, dtype=np.int32)\n", + " )\n", + " arrow_rings = list_array(ring_offsets, poly_coords_arrow)\n", + " geom_offsets = Array.from_numpy(\n", + " np.arange(poly_coords_closed.shape[0] + 1, dtype=np.int32)\n", + " )\n", + " arrow_geoms = list_array(geom_offsets, arrow_rings)\n", + " extension_metadata = {\"ARROW:extension:name\": \"geoarrow.polygon\"}\n", + " field = Field(\n", + " \"geometry\", arrow_geoms.type, nullable=True, metadata=extension_metadata\n", + " )\n", + " schema = Schema([field])\n", + " table = Table.from_arrays([arrow_geoms], schema=schema)\n", + " return table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the geoarrow table from the lon and lat data\n", + "arrow_table = lon_lat_to_arrow(lons, lats, points_to_keep)\n", + "# Create a solid polygon layer with the geoarrow table and the reflectivity data\n", + "layer = SolidPolygonLayer(\n", + " table=arrow_table,\n", + " get_fill_color=apply_continuous_cmap(\n", + " reflectivity_normalized, cmweather.cm_colorblind.ChaseSpectral\n", + " ),\n", + ")\n", + "# Plot the data\n", + "m = Map(layers=[layer])\n", + "m" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}