add lonboard example

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 ([email protected])\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
+}