{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Create a lonboard map of NEXRAD reflectivity\n", "\n", "An example which follows the \"Create a plot of NEXRAD reflectivity\", but\n", "plots with the lonboard deck.gl python interface.\n", "\n", "Author: Sam Gardner (samuel.gardner@ttu.edu)\n", "License: BSD 3 clause\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "## You are using the Python ARM Radar Toolkit (Py-ART), an open source\n", "## library for working with weather radar data. Py-ART is partly supported\n", "## by the U.S. Department of Energy Office of Science as part of\n", "## the Atmospheric Radiation Measurement (ARM) User Facility.\n", "##\n", "## If you use this software to prepare a publication, please cite:\n", "##\n", "## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119\n", "\n" ] } ], "source": [ "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Download and read the test data file\n", "\n", "We use Py-ART's built-in test data fetcher to grab a sample NEXRAD Level 2\n", "archive file, then read it in and extract the first sweep." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "filename = get_test_data(\"Level2_KATX_20130717_1950.ar2v\")\n", "radar = pyart.io.read_nexrad_archive(filename)\n", "\n", "# extract the first sweep of data\n", "first_sweep = radar.extract_sweeps([0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute gate corner coordinates\n", "\n", "Extract the range, azimuth, and elevation of gates in the first sweep, then\n", "convert the gate center coordinates into latitude/longitude edges (i.e.\n", "the corners of each radar gate) so we can later build polygons for them." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "range_center = first_sweep.range[\"data\"]\n", "az_center = first_sweep.azimuth[\"data\"]\n", "elevation = first_sweep.elevation[\"data\"]\n", "\n", "x_edge, y_edge, _ = pyart.core.antenna_vectors_to_cartesian(\n", " range_center, az_center, elevation, edges=True\n", ")\n", "\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": "markdown", "metadata": {}, "source": [ "## Normalize reflectivity for colormapping\n", "\n", "Lonboard colormaps expect data to be scaled to the range 0-1, so we\n", "normalize the reflectivity data using a vmin/vmax pair. The\n", "ChaseSpectral colormap aligns its lightness breaks with meteorologically\n", "relevant reflectivity values when used with vmin = -10 dBZ and\n", "vmax = 80 dBZ.\n", "\n", "\n", "Extract reflectivity data from the first sweep, fill masked values with\n", "NaN, and scale it to 0-1. Rendering is much faster if data below vmin\n", "and NaN values are excluded before plotting, so we also build a mask of\n", "points to keep." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "reflectivity = first_sweep.fields[\"reflectivity\"][\"data\"].filled(np.nan)\n", "\n", "chasespectral_vmin = -10\n", "chasespectral_vmax = 80\n", "\n", "\n", "reflectivity_normalized = np.ravel(\n", " (reflectivity - chasespectral_vmin) / (chasespectral_vmax - chasespectral_vmin)\n", ")\n", "\n", "\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": "markdown", "metadata": {}, "source": [ "## Build geoarrow polygons from gate corners\n", "\n", "The function below takes the corners of each radar gate and returns them as\n", "a geoarrow table. Geoarrow is a way to represent geospatial data in\n", "Apache Arrow, which is what lonboard uses under the hood for rendering\n", "geospatial data.\n", "\n", "\n", "This function was originally written by Kyle Barron — see this\n", "lonboard discussion\n", "for more details." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "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", " # bottom-right: (i+1, j)\n", " br_lon = lon_edge[1:, :-1]\n", " br_lat = lat_edge[1:, :-1]\n", " # top-right: (i+1, j+1)\n", " tr_lon = lon_edge[1:, 1:]\n", " tr_lat = lat_edge[1:, 1:]\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", "\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", " 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", "\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", "\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", "\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": "markdown", "metadata": {}, "source": [ "## Build the map\n", "\n", "Apply the lon_lat_to_arrow function to our gate corners and mask, apply\n", "the ChaseSpectral colormap to the normalized reflectivity values, and\n", "render everything as a SolidPolygonLayer on a lonboard Map." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "73b0e071a53f4286b7fa8aa18a4b168f", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(, VBox(children=(ErrorOutput(), ErrorOutput(), Erro…" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "table = lon_lat_to_arrow(lons, lats, points_to_keep)\n", "colors = apply_continuous_cmap(\n", " reflectivity_normalized, cmweather.cm_colorblind.ChaseSpectral\n", ")\n", "layer = SolidPolygonLayer(table=table, get_fill_color=colors)\n", "m = Map(layer)\n", "m" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.14.8" } }, "nbformat": 4, "nbformat_minor": 4 }