Visualizing Severe Weather in Wisconsin#

Motivation#

The Outbreak#

Last Friday (March 31, 2023), a severe weather outbreak impacted a large portion of the Central United States, including the Upper Midwest. There were a total of:

  • 142 tornado reports

  • 368 wind reports, with 7 of those being high wind reports, in excess of 65 knots (~75 mph)

  • 184 hail reports, with 14 of those bing large hail reports, in excess of 2 inches in diameter

Here is a map of the storm reports provided by the National Weather Service:

from IPython.display import IFrame

IFrame(
    src="https://www.spc.noaa.gov/climo/gmf.php?rpt=230331_rpts_filtered",
    width=700,
    height=500,
)

Focusing on Southern Wisconsin#

In this blog post, we explore a single area - Walworth County, Wisconsin, which is located in Southeastern Wisconsin. This county is home to Lake Geneva, Wisconsin, a popular tourist destination, especially for those living in the Chicagoland area.

3 confirmed tornadoes moved through Walworth county, with those confirmations from the local National Weather Service (NWS Milwaukee). They also found a considerable amount of wind and hail damage around the area.

NWS Milwaukee assembled an event summary (Link to NWS Event Summary) which includes the following radar imagery, which was produced using data from their radar and GRLevel software: NWS Milwaukee Radar Image

Reproducing Event Images#

What if we wanted to reproduce these plots, using free, open-source software, using open data? We can do this using Py-ART. We encourage people check out the Radar Cookbook, which walks through how to get started with the software stack, and details how to work with NEXRAD data, as we are in this case.

The focus here is to build upon the workflow described there, and automate this as much as possible, to enable local NWS offices to easily create post-event evaluation images.

Visualizing the Event#

Imports#

We import a variety of python packages, which can be installed using

conda install -c conda-forge arm_pyart metpy cartopy pandas numpy shapely
from datetime import datetime, timedelta
import warnings

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature
import fsspec
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
import numpy as np
import pandas as pd
import pyart

warnings.filterwarnings("ignore")
## 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 as part of the Atmospheric
## Radiation Measurement (ARM) Climate Research Facility, an Office of
## Science 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
/Users/mgrover/miniforge3/envs/pyart-dev/lib/python3.10/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Access the Data#

We start first with some extra map features - in this case, we add interstates to the map (credit to the Census database for the data, and Randy Chase for the code)

Download the Road Data#

I am running this on a Mac, so I used curl, but if you are on a Linux machine you would want to use wget here.

!curl https://www2.census.gov/geo/tiger/TIGER2016/PRIMARYROADS/tl_2016_us_primaryroads.zip --output tl_2016_us_primaryroads.zip
!unzip tl_2016_us_primaryroads.zip
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 25.7M    0 25.7M    0     0  21.0M      0 --:--:--  0:00:01 --:--:-- 21.0M
Archive:  tl_2016_us_primaryroads.zip
 extracting: tl_2016_us_primaryroads.cpg  
  inflating: tl_2016_us_primaryroads.dbf  
  inflating: tl_2016_us_primaryroads.prj  
  inflating: tl_2016_us_primaryroads.shp  
  inflating: tl_2016_us_primaryroads.shp.ea.iso.xml  
  inflating: tl_2016_us_primaryroads.shp.iso.xml  
  inflating: tl_2016_us_primaryroads.shp.xml  
  inflating: tl_2016_us_primaryroads.shx  

Add the Interstate Date as a feature#

We can add this data as a feature to be used on the maps using Shapely

# Read shape file
reader = shpreader.Reader("./tl_2016_us_primaryroads.shp")
names = []
geoms = []

# Loop through and find the interstates, looking for an I in the full name
for rec in reader.records():
    if rec.attributes["FULLNAME"][0] == "I":
        names.append(rec)
        geoms.append(rec.geometry)

# make interstate feature
interstate_feature = ShapelyFeature(
    geoms, ccrs.PlateCarree(), edgecolor="darkblue", lw=0.5, facecolor="none"
)

Investigate the NEXRAD Data Available#

The NEXRAD data is available through Amazon Web Services (AWS), and we are using the Level2 data, which has been quality controlled to an extent, but has not been fully cleaned (ex. velocity field may still be aliased).

fs = fsspec.filesystem("s3", anon=True)
site = "KMKX"
year = "2023"
month = "04"
day = "01"
hour = "01"
files = sorted(
    fs.glob(
        f"s3://noaa-nexrad-level2/{year}/{month}/{day}/{site}/{site}{year}{month}{day}_{hour}*06"
    )
)
files
['noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_010517_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_011057_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_011638_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_012219_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_012759_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_013341_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_013922_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_014502_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_015042_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_015623_V06',
 'noaa-nexrad-level2/2023/04/01/KMKX/KMKX20230401_015928_V06']

Setup our Plotting Functions#

Now that we have our data, we setup a function to process data, which:

  • Reads the data into a Py-ART radar object

  • Subset the lowest elevations

  • Dealias the radial velocity

  • Plot the velocity

  • Plot the reflectivity

  • Save the plot to a file

plot_configuration = {
    "reflectivity": {
        "vmin": -20,
        "vmax": 70,
        "cmap": "pyart_ChaseSpectral",
        "label": "Horizontal Reflectivity (dBZ)",
    },
    "velocity": {
        "vmin": -50,
        "vmax": 50,
        "cmap": "twilight_shifted",
        "label": "Radial Velocity (m/s)",
    },
    "cross_correlation_ratio": {
        "vmin": 0.7,
        "vmax": 1.0,
        "cmap": "pyart_ChaseSpectral",
        "label": "Cross Correlation Ratio",
    },
}
def plot_radar_ppi(
    file,
    left_field="velocity",
    right_field="reflectivity",
    plot_config=plot_configuration,
):
    # Read the radar data
    radar = pyart.io.read_nexrad_archive(f"s3://{file}")

    # Extract the lowest two sweeps - to help with speeding up the data cleaning
    radar = radar.extract_sweeps([0, 1])

    if "velocity" in [left_field, right_field]:
        # Dealias the radial velocities, using a reflectivity threshold to only look at precip regions
        nyq = radar.instrument_parameters["nyquist_velocity"]["data"][0]
        gatefilter = pyart.filters.GateFilter(radar)
        gatefilter.exclude_below("reflectivity", 0)
        vel_dealias = pyart.correct.dealias_region_based(
            radar, nyquist_vel=nyq, centered=True, gatefilter=gatefilter
        )
        radar.add_field("corrected_velocity", vel_dealias, replace_existing=True)

    # Add a subplot for the velocity
    fig = plt.figure(figsize=(14, 8))
    graph = pyart.graph.RadarMapDisplay(radar)

    # Create a plot number to iterate over
    plot_number = 0
    for field in (left_field, right_field):
        # Use the configuration for the field
        config = plot_config[field]

        # Assume 0 sweep (lowest elevation)
        sweep = 0

        # User a new name for the corrected velocity, and use the next sweep
        if (field == "velocity") & (vel_dealias is not None):
            field = "corrected_velocity"
            sweep = 1

        # Plot the first field
        ax = plt.subplot(121 + plot_number, projection=ccrs.PlateCarree())

        # Prepare the title
        title = graph.generate_title(field, 0)
        time = pd.to_datetime(radar.time["units"][-20:])
        time_steps = radar.time["data"]
        index = radar.get_start_end(0)[-1]
        time += timedelta(seconds=time_steps[index])
        central_time = str(time.astimezone("US/Central"))[:-13]
        new_title = f"{title[0:13]} \n {central_time} CDT"
        graph.plot_ppi_map(
            field,
            cmap=config["cmap"],
            sweep=sweep,
            vmin=config["vmin"],
            vmax=config["vmax"],
            ax=ax,
            colorbar_orient="horizontal",
            colorbar_label=config["label"],
            projection=ccrs.PlateCarree(),
        )

        # Add a point near the start of the tornado in Sharon
        ax.scatter(
            -88.7, 42.53, label="Sharon, Wisconsin (Near Start of Tornado)", color="k"
        )

        # Add a point for Delavan Lake, Wisconsin
        ax.scatter(
            -88.6105,
            42.60666,
            label="Delavan Lake, Wisconsin (Near End of Tornado)",
            color="lightgrey",
        )

        # Add the interstates
        ax.add_feature(
            interstate_feature,
            edgecolor="lavenderblush",
            facecolor="None",
            linewidth=2,
            label="Interstate 43",
        )

        # Add the counties
        ax.add_feature(USCOUNTIES, linewidth=0.5)

        # Add gridlines
        gl = ax.gridlines(
            crs=ccrs.PlateCarree(),
            draw_labels=True,
            linewidth=1,
            color="gray",
            alpha=0.3,
            linestyle="--",
        )
        gl.right_labels = False
        gl.top_labels = False

        # Add a label for the interstates
        plt.plot(0, 0, color="lavenderblush", label="Interstate 43")
        plt.xlim(-88.8, -88.2)
        plt.ylim(42.3, 42.9)
        plt.legend(loc="upper right")

        # Add one to the plot number
        plot_number += 1

    plt.title(new_title)
    plt.savefig(f"ppi_velocity_reflectivity_{time}.png", dpi=300)
    plt.show()
    plt.close()

Run our Plotting Functions#

Now that we have our plotting functions, we can apply these to the first volume of our file list, when the NWS identified the Sharon to Delavan Lake tornado NWS Tornado Damage Survey

plot_radar_ppi(files[0])
../../_images/50d142cd84cb9f3b4f2a224163b90f9904ab4f4bba2cfc62d01e3b9874806d67.png
plot_radar_ppi(files[0], right_field="cross_correlation_ratio")
../../_images/7d730a3d5f6e619091f430c4d673102bce40dd02112fe4b73dd96734dacb3650.png

Conclusion#

Within this post, we examined how to use data from the National Weather Service to identify tornado signatures, and create a flexible plotting function to visualize its path. We even added interstates for reference!

We hope that this will be applicable to NWS offices interested in making colorblind-friendly maps of severe weather events, and automate the data selection and visualization process (where possible).