Create a lonboard map of NEXRAD reflectivity

An example which follows the “Create a plot of NEXRAD reflectivity”, but plots with the lonboard deck.gl python interface.

Author: Sam Gardner (samuel.gardner@ttu.edu) License: BSD 3 clause

import cmweather
import numpy as np
from arro3.core import Array, Field, Schema, Table, fixed_size_list_array, list_array
from lonboard import Map, SolidPolygonLayer
from lonboard.colormap import apply_continuous_cmap

import pyart
from pyart.testing import get_test_data

## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar data. Py-ART is partly supported
## by the U.S. Department of Energy Office of Science as part of
## the Atmospheric Radiation Measurement (ARM) User Facility.
##
## If you use this software to prepare a publication, please cite:
##
##     JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119

Download and read the test data file

We use Py-ART’s built-in test data fetcher to grab a sample NEXRAD Level 2 archive file, then read it in and extract the first sweep.

filename = get_test_data("Level2_KATX_20130717_1950.ar2v")
radar = pyart.io.read_nexrad_archive(filename)

# extract the first sweep of data
first_sweep = radar.extract_sweeps([0])

Compute gate corner coordinates

Extract the range, azimuth, and elevation of gates in the first sweep, then convert the gate center coordinates into latitude/longitude edges (i.e. the corners of each radar gate) so we can later build polygons for them.

range_center = first_sweep.range["data"]
az_center = first_sweep.azimuth["data"]
elevation = first_sweep.elevation["data"]

x_edge, y_edge, _ = pyart.core.antenna_vectors_to_cartesian(
    range_center, az_center, elevation, edges=True
)


radar_lat = radar.latitude["data"][0]
radar_lon = radar.longitude["data"][0]
projparams = {"proj": "pyart_aeqd", "lon_0": radar_lon, "lat_0": radar_lat}
lons, lats = pyart.core.cartesian_to_geographic(x_edge, y_edge, projparams)

Normalize reflectivity for colormapping

Lonboard colormaps expect data to be scaled to the range 0-1, so we normalize the reflectivity data using a vmin/vmax pair. The ChaseSpectral colormap aligns its lightness breaks with meteorologically relevant reflectivity values when used with vmin = -10 dBZ and vmax = 80 dBZ.

Extract reflectivity data from the first sweep, fill masked values with NaN, and scale it to 0-1. Rendering is much faster if data below vmin and NaN values are excluded before plotting, so we also build a mask of points to keep.

reflectivity = first_sweep.fields["reflectivity"]["data"].filled(np.nan)

chasespectral_vmin = -10
chasespectral_vmax = 80


reflectivity_normalized = np.ravel(
    (reflectivity - chasespectral_vmin) / (chasespectral_vmax - chasespectral_vmin)
)


reflectivity_normalized[reflectivity_normalized < 0] = np.nan
points_to_keep = ~np.isnan(reflectivity_normalized)
reflectivity_normalized = reflectivity_normalized[points_to_keep]

Build geoarrow polygons from gate corners

The function below takes the corners of each radar gate and returns them as a geoarrow table. Geoarrow is a way to represent geospatial data in Apache Arrow, which is what lonboard uses under the hood for rendering geospatial data.

This function was originally written by Kyle Barron — see this lonboard discussion for more details.

def lon_lat_to_arrow(lon_edge, lat_edge, data_mask):
    # bottom-left: (i, j)
    bl_lon = lon_edge[:-1, :-1]
    bl_lat = lat_edge[:-1, :-1]
    # bottom-right: (i+1, j)
    br_lon = lon_edge[1:, :-1]
    br_lat = lat_edge[1:, :-1]
    # top-right: (i+1, j+1)
    tr_lon = lon_edge[1:, 1:]
    tr_lat = lat_edge[1:, 1:]
    # top-left: (i, j+1)
    tl_lon = lon_edge[:-1, 1:]
    tl_lat = lat_edge[:-1, 1:]

    lons_to_poly = np.array(
        [bl_lon.flatten(), br_lon.flatten(), tr_lon.flatten(), tl_lon.flatten()]
    )
    lats_to_poly = np.array(
        [bl_lat.flatten(), br_lat.flatten(), tr_lat.flatten(), tl_lat.flatten()]
    )

    poly_coords = np.stack([lons_to_poly, lats_to_poly], axis=-1).transpose(1, 0, 2)
    poly_coords = poly_coords[data_mask, :, :]
    poly_coords_closed = np.concat((poly_coords, poly_coords[:, 0:1, :]), axis=1)

    poly_coords_arrow = fixed_size_list_array(
        Array.from_numpy(np.ravel(poly_coords_closed)), 2
    )

    ring_offsets = Array.from_numpy(
        np.arange(0, (poly_coords_closed.shape[0] + 1) * 5, 5, dtype=np.int32)
    )
    arrow_rings = list_array(ring_offsets, poly_coords_arrow)

    geom_offsets = Array.from_numpy(
        np.arange(poly_coords_closed.shape[0] + 1, dtype=np.int32)
    )
    arrow_geoms = list_array(geom_offsets, arrow_rings)

    extension_metadata = {"ARROW:extension:name": "geoarrow.polygon"}
    field = Field(
        "geometry", arrow_geoms.type, nullable=True, metadata=extension_metadata
    )
    schema = Schema([field])
    table = Table.from_arrays([arrow_geoms], schema=schema)
    return table

Build the map

Apply the lon_lat_to_arrow function to our gate corners and mask, apply the ChaseSpectral colormap to the normalized reflectivity values, and render everything as a SolidPolygonLayer on a lonboard Map.

table = lon_lat_to_arrow(lons, lats, points_to_keep)
colors = apply_continuous_cmap(
    reflectivity_normalized, cmweather.cm_colorblind.ChaseSpectral
)
layer = SolidPolygonLayer(table=table, get_fill_color=colors)
m = Map(layer)
m