Chapter 2: Case Study - MRMS Cookbook

March 24-27, 2023 Tornado Outbreak¶

Next, title your notebook appropriately with a top-level Markdown header, # (see the very first cell above). Do not use this level header anywhere else in the notebook. Our book build process will use this title in the navbar, table of contents, etc. Keep it short, keep it descriptive.

Follow this with a --- cell to visually distinguish the transition to the prerequisites section.

Overview¶

The tornado outbreak of March 24–27, 2023 was a devastating multi-day severe weather event that swept across the Southern United States, particularly impacting Mississippi, Alabama, Tennessee, and Georgia. Triggered by a slow-moving upper-level trough interacting with moist, unstable air from the Gulf of Mexico, the outbreak produced 35 confirmed tornadoes, including a violent EF4 that tore through Rolling Fork, Midnight, and Silver City, Mississippi with peak winds of 195 mph. This EF4 tornado alone caused catastrophic damage and multiple fatalities and brought tornado emergencies ahead of widespread destruction.
Over the four-day span, the system also unleashed damaging straight-line winds, large hail, and flooding. In total, the outbreak resulted in 23 fatalities (plus two from non-tornadic causes), over 230 injuries, and an estimated $1.9 billion in damage. The event was notable not only for its intensity but also for its geographic breadth and the prolonged nature of the severe weather threat.

This chapter explores MRMS data from this tornado outbreak, specifically from the Rolling Fork–Silver City, MS tornado on March 24, 2023. The chapter investigates:

Reflectivity
Precipitation Rates
Rotation
Hail Swaths - Under Construction!
Storm Intensity (Vertically integrated liquid) - Under Construction!

To support this analysis, the chapter introduces a practical method for accessing MRMS data directly from an AWS server. It guides readers through defining a map and customizing its spatial extent to filter and process data relevant to the event. By showcasing multiple variations of each MRMS variable, the chapter highlights the accessibility and versatility of these datasets for visualizing high-impact weather events. A focused case study on the Rolling Fork tornado illustrates how several hours of MRMS data can be leveraged to gain insight into storm properties of a significant weather event.

Imports¶

This section brings in packages for handling AWS requests, manipulating MRMS files, and rendering plots on a map.

import sys
import s3fs
import urllib
import tempfile
import gzip
import xarray as xr
import xarray
import io
import numpy as np
import cartopy
import datetime
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

aws = s3fs.S3FileSystem(anon=True)

Build Map¶

This section uses Cartopy to build a blank map then define our extent for our case study.

# Set up the map projection
projection = ccrs.LambertConformal(central_longitude=-96, central_latitude=39)

# Create the figure and axes
fig, ax = plt.subplots(figsize=(12, 8), subplot_kw={'projection': projection})

# Set extent for CONUS (approximate)
ax.set_extent([-125, -66.5, 24, 50], crs=ccrs.PlateCarree())

# Add geographic features
ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))

# Optional: remove ticks
ax.set_xticks([])
ax.set_yticks([])

# Add title
plt.title("Blank Map", fontsize=18)

plt.show()

We will be looking specifically at the tornado outbreak that occurred in Dixie Alley, so let’s set our extents specifically to Dixie Alley.¶

lon_min, lon_max = -96, -80
lat_min, lat_max = 29, 38

# Set up the map projection
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create the figure and axes
fig, ax = plt.subplots(figsize=(12, 8), subplot_kw={'projection': projection})

# Set extent for CONUS (approximate)
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

# Add geographic features
ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))

# Optional: remove ticks
ax.set_xticks([])
ax.set_yticks([])

# Add title
plt.title("Blank Dixie Alley Map", fontsize=18)

plt.show()

Fetch Data¶

This section uses s3 to pull in the data from the AWS s3 server. Following data acquisition, the module uses xarray to filter the resulting DataArray for the size of the desired map.¶

For variable names, see link here.¶

# def fetch_mrms_data(variable: str, yyyymmdd: str, hh: str) -> xr.DataArray:
#     """
#     Downloads and loads MRMS data from NOAA PDS.

#     Parameters:
#         variable (str): MRMS product name (e.g., 'MergedReflectivityQC').
#         yyyymmdd (str): Date in YYYYMMDD format.
#         hh (str): Hour in HH format (00–23 UTC).

#     Returns:
#         xarray.DataArray: Decoded MRMS data array.
#     """
#     url = (
#         f"https://noaa-mrms-pds.s3.amazonaws.com/CONUS/{variable}/"
#         f"{yyyymmdd}/MRMS_{variable}_{yyyymmdd}-{hh}0000.grib2.gz"
#     )

#     response = urllib.request.urlopen(url)
#     compressed_file = response.read()

#     with tempfile.NamedTemporaryFile(suffix=".grib2") as f:
#         f.write(gzip.decompress(compressed_file))
#         data_in = xr.load_dataarray(f.name, engine='cfgrib', decode_timedelta=True)

#     return data_in

# def fetch_mrms_data(
#     variable: str,
#     yyyymmdd: str,
#     hh: str,
#     lon_min: float = None,
#     lat_min: float = None,
#     lon_max: float = None,
#     lat_max: float = None
# ) -> xr.DataArray:
#     """
#     Downloads and loads MRMS data from NOAA PDS, with optional spatial filtering.

#     Parameters:
#         variable (str): MRMS product name (e.g., 'MergedReflectivityQC').
#         yyyymmdd (str): Date in YYYYMMDD format.
#         hh (str): Hour in HH format (00–23 UTC).
#         lon_min, lat_min, lon_max, lat_max (float, optional): Bounding box for spatial subset. 

#     Returns:
#         xarray.DataArray: Decoded MRMS data array, optionally subset by lat/lon.

#     Example use: 
#         data = fetch_mrms_data('MergedReflectivityQC', '20230325', '02', lon_min=-96, lat_min=29, lon_max=-80, lat_max=38)
#     """
#     url = (
#         f"https://noaa-mrms-pds.s3.amazonaws.com/CONUS/{variable}/"
#         f"{yyyymmdd}/MRMS_{variable}_{yyyymmdd}-{hh}0000.grib2.gz"
#     )

#     response = urllib.request.urlopen(url)
#     compressed_file = response.read()

#     with tempfile.NamedTemporaryFile(suffix=".grib2") as f:
#         f.write(gzip.decompress(compressed_file))
#         data_in = xr.load_dataarray(f.name, engine='cfgrib', decode_timedelta=True)

#     # Optional spatial filtering
#     if all(v is not None for v in [lon_min, lat_min, lon_max, lat_max]):
#         data_in = data_in.sel(
#             latitude=slice(lat_max, lat_min),  # descending order
#             longitude=slice(360 - abs(lon_min), 360 - abs(lon_max))        
#         )

#     return data_in

def find_available_files(
    variable: str,
    yyyymmdd: str,
    hh: str
):
    
    files_list = []

    available_files = aws.ls(f'noaa-mrms-pds/CONUS/{variable}/{yyyymmdd}/', refresh=True)
    for file in available_files:
        file_hour = file[-15:-13]
        if file_hour == hh:
            files_list.append(file)

    if len(files_list) == 0:
        raise ValueError(f"No files found for {variable} on {yyyymmdd} at hour {hh}.")
    else:    
        return files_list

def fetch_mrms_data(    
    file: str,
    lon_min: float = None,
    lat_min: float = None,
    lon_max: float = None,
    lat_max: float = None
):
    url = (f"https://noaa-mrms-pds.s3.amazonaws.com/{file[14:]}")
    response = urllib.request.urlopen(url)
    compressed_file = response.read()        
    with tempfile.NamedTemporaryFile(suffix=".grib2") as f:
        f.write(gzip.decompress(compressed_file))
        data_in = xr.load_dataarray(f.name, engine='cfgrib', decode_timedelta=True)

    # Optional spatial filtering
    if all(v is not None for v in [lon_min, lat_min, lon_max, lat_max]):
        data_in = data_in.sel(
            latitude=slice(lat_max, lat_min),  # descending order
            longitude=slice(360 - abs(lon_min), 360 - abs(lon_max))        
        )

    return data_in

# response = urllib.request.urlopen("https://noaa-mrms-pds.s3.amazonaws.com/CONUS/CREF_1HR_MAX_00.50/20230325/MRMS_CREF_1HR_MAX_00.50_20230325-010000.grib2.gz")

# compressed_file = response.read()

# with tempfile.NamedTemporaryFile(suffix=".grib2") as f:
#             f.write(gzip.decompress(compressed_file))
#             data_in = xr.load_dataarray(f.name, engine='cfgrib', decode_timedelta=True)

Case Study - March 24, 2023¶

Rolling Fork - Silver City, MS Tornado¶

3/25/23 1z to 2z¶

# Lon mins and maxes for our projections:
lon_min, lon_max = -96, -80
lat_min, lat_max = 29, 38

Maximum 1-Hour Composite Reflectivity¶

The MRMS Max 1-Hour Composite Reflectivity product represents the highest reflectivity value observed within the past hour across all radar scans, providing a time-integrated view of storm intensity. It helps forecasters identify areas of persistent or intense convection, especially useful for tracking severe weather like hail or heavy rainfall. This product is derived from a seamless mosaic of multiple radars, quality-controlled to remove non-meteorological artifacts.¶

#### March 24, 2023 - Rolling Fork - Silver City, MS Tornado -- EF4, 71 minutes long, est winds 195 mph
## 3/25/23 1z to 2z, so we'll grab two hours of data shortly

#Grab 2 hours of data for plotting
cref1files = find_available_files('CREF_1HR_MAX_00.50', '20230325', '01')
cref1z = fetch_mrms_data(cref1files[0])    

cref2files = find_available_files('CREF_1HR_MAX_00.50', '20230325', '02')
cref2z = fetch_mrms_data(cref2files[0])

# Mask fill values for both datasets
masked1 = np.ma.masked_where(cref1z == -99.0, cref1z)
masked2 = np.ma.masked_where(cref2z == -99.0, cref2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        cref2z.longitude, cref2z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Reflectivity (dBZ)')
plt.suptitle('Max 1HR Composite Reflectivity:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')

plt.show()

Surface Precip Rate¶

To describe Surface Precip Rate, there are three variables that can be used:

Surface Precipitation Rate Products¶

Variable Name	Description	Temporal Resolution	Filename Pattern
Instantaneous PrecipRate	- Estimates current rainfall intensity - Derived from dual-pol radar	- Every 2 minutes	`PrecipRate_00.00`
MultiSensor QPE (Pass 1 & Pass 2)	- Combines radar and precip gauge data - Available in 1-pass and 2-pass versions - Used for hourly accumulation	- Hourly (Pass 1 and Pass 2)	`MRMS_QPE_01H_Pass1_00.00`
RadarOnly QPE	- Estimates surface rainfall rate using dual-polarization radar reflectivity. - Captures rapid changes in precipitation intensity at high temporal resolution.	- Every 2 minutes - Available in 15 minute as well as (1, 3, 6, 12, 24, 48) hour intervals - QPE since 12z also available	`RadarOnly_QPE_01H_00.00`

Instantaneous Precip Rate¶

precip1files= find_available_files('PrecipRate_00.00', '20230325', '01') # Precip Rate
precip_1z = fetch_mrms_data(precip1files[0]) 
precip2files = find_available_files('PrecipRate_00.00', '20230325', '02')
precip_2z = fetch_mrms_data(precip2files[0])
masked1 = np.ma.masked_where(precip_1z <= 0, precip_1z)
masked2 = np.ma.masked_where(precip_2z <= 0, precip_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        precip_1z.longitude, precip_2z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Precipitation (mm)')
plt.suptitle('Instantaneous Precipitation', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')

plt.show()

MultiSensorQPE - 1 Hour - Pass 1¶

# MultiSensor_QPE_01H_Pass1_00.00
QPE1files = find_available_files('MultiSensor_QPE_01H_Pass1_00.00', '20230325', '01') # QPE: Quantified Precip Estimation - Offered hourly.
QPE_1z = fetch_mrms_data(QPE1files[0]) 
QPE2files = find_available_files('MultiSensor_QPE_01H_Pass1_00.00', '20230325', '02')
QPE_2z = fetch_mrms_data(QPE2files[0])
masked1 = np.ma.masked_where(QPE_1z <= 0, QPE_1z)
masked2 = np.ma.masked_where(QPE_2z <= 0, QPE_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        QPE_1z.longitude, QPE_2z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Precipitation (mm)')
plt.suptitle('Multi-Sensor Quantified Precipitation Estimate (Pass 1):', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')

plt.show()

MultiSensorQPE - 1 Hour - Pass 2¶

# MultiSensor_QPE_01H_Pass2_00.00
QPE1zp2files = find_available_files('MultiSensor_QPE_01H_Pass2_00.00', '20230325', '01') # QPE: Quantified Precip Estimation - last hour, 2nd pass
QPE_1z_p2 = fetch_mrms_data(QPE1zp2files[0]) 
QPE2files = find_available_files('MultiSensor_QPE_01H_Pass2_00.00', '20230325', '02')
QPE_2z_p2 = fetch_mrms_data(QPE2files[0])
masked1 = np.ma.masked_where(QPE_1z_p2 <= 0, QPE_1z_p2)
masked2 = np.ma.masked_where(QPE_2z_p2 <= 0, QPE_2z_p2)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        QPE_1z_p2.longitude, QPE_2z_p2.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Precipitation (mm)')
plt.suptitle('Multi-Sensor Quantified Precipitation Estimate (Pass 2):', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')

plt.show()

Radar Only QPE - Last Hour¶

RQPE1files = find_available_files('RadarOnly_QPE_01H_00.00', '20230325', '01') # RadarOnly_QPE: Radar Only Quantified Precip Estimation - last hour
RQPE_1z = fetch_mrms_data(RQPE1files[0])
RQPE2files = find_available_files('RadarOnly_QPE_01H_00.00', '20230325', '02')
RQPE_2z = fetch_mrms_data(RQPE2files[0])
masked1 = np.ma.masked_where(RQPE_1z <= 0, RQPE_1z)
masked2 = np.ma.masked_where(RQPE_2z <= 0, RQPE_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        RQPE_1z.longitude, RQPE_2z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Precipitation (mm)')
plt.suptitle('Radar Only Quantified Precipitation Estimate:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')

plt.show()

Rotation¶

These rotation products can be combined to assess both the intensity and persistence of storm-scale rotation across multiple atmospheric layers and time scales. By layering instantaneous azimuthal shear with rotation tracks—especially ML-enhanced versions—forecasters and researchers can better identify evolving mesocyclones, discriminate between transient and sustained rotation, and refine environmental risk assessments for severe weather.¶

Variable	Description	Temporal Resolution	Filename Pattern
Merged AzShear (0-2km AGL)	Low-level azimuthal shear (0–2 km AGL); highlights near-surface rotation.	Instantaneous	MergedAzShear_0-2kmAGL_00.50
Merged AzShear (3-6km AGL)	Mid-level azimuthal shear (3–6 km AGL); captures elevated storm rotation.	Instantaneous	MergedAzShear_3-6kmAGL_00.50
Rotation Track (30-min)	30-min accumulation of low-level rotation; useful for short-term tracking.	30 minutes	RotationTrack30min_00.50
Rotation Track (60-min)	60-min accumulation of low-level rotation; highlights sustained activity.	60 minutes	RotationTrack60min_00.50
Rotation Track ML (30-min)	ML-enhanced 30-min rotation track; filters noise, boosts confidence.	30 minutes	RotationTrackML30min_00.50
Rotation Track ML (60-min)	ML-enhanced 60-min rotation track; detects short-lived intense rotation.	60 minutes	RotationTrackML60min_00.50

Merged AzShear (0-2km AGL)¶

# MergedAzShear_0-2kmAGL_00.50
Azshrfiles1_2km = find_available_files('MergedAzShear_0-2kmAGL_00.50', '20230325', '01') # Merged AzShr (0-2km AGL)
Azshr_1z = fetch_mrms_data(Azshrfiles1_2km[0])
Azshrfiles2_2km = find_available_files('MergedAzShear_0-2kmAGL_00.50', '20230325', '02')
Azshr_2z = fetch_mrms_data(Azshrfiles2_2km[0])
masked1 = np.ma.masked_where(Azshr_1z <= 0, Azshr_1z)
masked2 = np.ma.masked_where(Azshr_2z <= 0, Azshr_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        Azshr_1z.longitude, Azshr_1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Azimuthal Shear (s^-1)')
plt.suptitle('Merged AzShear (0-2km AGL):', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')

plt.show()

Merged AzShear (3-6 km AGL)¶

# MergedAzShear_3-6kmAGL_00.50
Azshrfiles1_3km = find_available_files('MergedAzShear_3-6kmAGL_00.50', '20230325', '01') # Merged AzShr (3-6km AGL)
Azshr3_1z = fetch_mrms_data(Azshrfiles1_3km[0])
Azshrfiles2_3km = find_available_files('MergedAzShear_3-6kmAGL_00.50', '20230325', '02')
Azshr3_2z = fetch_mrms_data(Azshrfiles2_3km[0])
masked1 = np.ma.masked_where(Azshr3_1z <= 0, Azshr3_1z)
masked2 = np.ma.masked_where(Azshr3_2z <= 0, Azshr3_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        Azshr3_1z.longitude, Azshr3_1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Azimuthal Shear (s^-1)')
plt.suptitle('Merged AzShear (3-6km AGL):', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Rotation Track (30-minute)¶

# RotationTrack30min_00.50 
RoTrack_30mFiles1 = find_available_files('RotationTrack30min_00.50', '20230325', '01') # Rotation Tracks (30-min)
RoTrack_30m_1z = fetch_mrms_data(RoTrack_30mFiles1[0])
RoTrack_30mFiles2 = find_available_files('RotationTrack30min_00.50', '20230325', '02')
RoTrack_30m_2z = fetch_mrms_data(RoTrack_30mFiles2[0])
masked1 = np.ma.masked_where(RoTrack_30m_1z <= 0, RoTrack_30m_1z)
masked2 = np.ma.masked_where(RoTrack_30m_2z <= 0, RoTrack_30m_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        RoTrack_30m_1z.longitude, RoTrack_30m_1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Rotation Track (30-minutes):', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Rotation Track (60 minutes)¶

# RotationTrack60min_00.50
RoTrack_60mFiles1 = find_available_files('RotationTrack60min_00.50', '20230325', '01') # Rotation Tracks (60-min)
RoTrack_60m_1z = fetch_mrms_data(RoTrack_60mFiles1[0])
RoTrack_60mFiles2 = find_available_files('RotationTrack30min_00.50', '20230325', '02')
RoTrack_60m_2z = fetch_mrms_data(RoTrack_60mFiles2[0])
masked1 = np.ma.masked_where(RoTrack_60m_1z <= 0, RoTrack_60m_1z)
masked2 = np.ma.masked_where(RoTrack_60m_2z <= 0, RoTrack_60m_2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        RoTrack_60m_1z.longitude, RoTrack_60m_1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Rotation Track (Last 60-minutes):', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

#Coming soon:
# RotationTrackML60min_00.50 -  for the sake of time, hold on this
# RotationTrackML30min_00.50 -  for the sake of time, hold on this

Hail Swaths¶

# MESH_00.50 - Maximum Expected Size of Hail
MeshFiles1 = find_available_files('MESH_00.50', '20230325', '01') # MESH
Mesh1z = fetch_mrms_data(MeshFiles1[0])
MeshFiles2 = find_available_files('MESH_00.50', '20230325', '02')
Mesh2z = fetch_mrms_data(MeshFiles2[0])
masked1 = np.ma.masked_where(Mesh1z <= 0, Mesh1z)
masked2 = np.ma.masked_where(Mesh2z <= 0, Mesh2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        Mesh1z.longitude, Mesh1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Maximum Expected Size of Hail:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Severe Hail Index¶

# SHI_00.50 - Severe Hail Index
SHIFiles1 = find_available_files('SHI_00.50', '20230325', '01') # SHI
SHI1z = fetch_mrms_data(SHIFiles1[0])
SHIFiles2 = find_available_files('SHI_00.50', '20230325', '02')
SHI2z = fetch_mrms_data(SHIFiles2[0])
masked1 = np.ma.masked_where(SHI1z <= 0, SHI1z)
masked2 = np.ma.masked_where(SHI2z <= 0, SHI2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        SHI1z.longitude, SHI1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Severe Hail Index:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Vertically Integrated Ice¶

# VII_00.50 - Vertically Integrated Ice
VIIFiles1 = find_available_files('VII_00.50', '20230325', '01')
VII1z = fetch_mrms_data(VIIFiles1[0])
VIIFiles2 = find_available_files('VII_00.50', '20230325', '02')
VII2z = fetch_mrms_data(VIIFiles2[0])
masked1 = np.ma.masked_where(VII1z <= 0, VII1z)
masked2 = np.ma.masked_where(VII2z <= 0, VII2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        VII1z.longitude, VII1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Vertically Integrated Ice:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Probability of Severe Hail¶

# POSH_00.50
POSHFiles1 = find_available_files('POSH_00.50', '20230325', '01')
POSH1z = fetch_mrms_data(POSHFiles1[0])
POSHFiles2 = find_available_files('POSH_00.50', '20230325', '02')
POSH2z = fetch_mrms_data(POSHFiles2[0])
masked1 = np.ma.masked_where(POSH1z <= 0, POSH1z)
masked2 = np.ma.masked_where(POSH2z <= 0, POSH2z)

projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        POSH1z.longitude, POSH1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Probability of Severe Hail', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Storm Intensity - Vertically Integrated Liquid¶

# VIL_00.50
VILFiles1 = find_available_files('VIL_00.50', '20230325', '01')
VIL1z = fetch_mrms_data(VILFiles1[0])
VILFiles2 = find_available_files('VIL_00.50', '20230325', '02')
VIL2z = fetch_mrms_data(VILFiles2[0])
masked1 = np.ma.masked_where(VIL1z <= 0, VIL1z)
masked2 = np.ma.masked_where(VIL2z <= 0, VIL2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        VIL1z.longitude, VIL1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Vertically Integrated Liquid:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Vertically Integrated Liquid Density¶

# VIL_Density_00.50
VILDFiles1 = find_available_files('VIL_Density_00.50', '20230325', '01')
VILD1z = fetch_mrms_data(VILDFiles1[0])
VILDFiles2 = find_available_files('VIL_Density_00.50', '20230325', '02')
VILD2z = fetch_mrms_data(VILDFiles2[0])
masked1 = np.ma.masked_where(VILD1z <= 0, VILD1z)
masked2 = np.ma.masked_where(VILD2z <= 0, VILD2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        VILD1z.longitude, VILD1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('Vertically Integrated Liquid Density:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

120-Minute Vertically Integrated Liquid Maximum¶

# VIL_Max_120min_00.50
VILMFiles1 = find_available_files('VIL_Max_120min_00.50', '20230325', '01')
VILM1z = fetch_mrms_data(VILMFiles1[0])
VILMFiles2 = find_available_files('VIL_Max_120min_00.50', '20230325', '02')
VILM2z = fetch_mrms_data(VILMFiles2[0])
masked1 = np.ma.masked_where(VILM1z <= 0, VILM1z)
masked2 = np.ma.masked_where(VILM2z <= 0, VILM2z)

# Define bounds for Dixie Alley
projection = ccrs.LambertConformal(central_longitude=-88, central_latitude=34)

# Create side-by-side subplots
fig, axes = plt.subplots(
    1, 2, figsize=(16, 8),
    subplot_kw={'projection': projection},
    gridspec_kw={'bottom': 0.2}  # leave room for shared colorbar
)

meshes = []
for ax, masked, title in zip(axes, [masked1, masked2], ["(a) 3/25/2023 @ 01z", "(b) 3/25/2023 @ 02z"]):
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='gray')
    ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', edgecolor='black')
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    
    mesh = ax.pcolormesh(
        VILM1z.longitude, VILM1z.latitude, masked,
        cmap='turbo', transform=ccrs.PlateCarree(), shading='auto'
    )
    ax.set_title(title, fontsize=15)
    meshes.append(mesh)

# Add shared colorbar beneath both plots
cbar_ax = fig.add_axes([0.25, 0.25, 0.5, 0.02])  # [left, bottom, width, height]
cbar = fig.colorbar(meshes[0], cax=cbar_ax, orientation='horizontal')
cbar.set_label('Units')
plt.suptitle('120-Minute Vertically Integrated Liqud Maximum:', fontsize='20', x=0.5, y=0.85, horizontalalignment='center', verticalalignment='top')
plt.show()

Conclusion¶

This chapter demonstrated how MRMS data can be efficiently accessed, processed, and visualized to analyze severe weather events, using the Rolling Fork–Silver City tornado as a case study. By leveraging AWS-hosted datasets and customizing spatial extents, users can explore multiple dimensions of storm structure—from reflectivity and rotation to hail swaths and precipitation rates.

The workflow presented here not only underscores the value of MRMS products for post-event analysis but also equips researchers and practitioners with tools to rapidly assess and interpret high-impact weather. As severe weather threats continue to evolve, accessible and scalable data pipelines like this one are essential for advancing situational awareness and scientific understanding.

Next Steps¶

In the next chapter, we shift focus from tornadic activity to flood impacts by examining the Texas Floods of early July 2025. Using the same MRMS data pipeline introduced here, we explore how radar reflectivity can be used to assess flood severity and spatial extent.

This upcoming case study highlights the adaptability of MRMS datasets for analyzing diverse weather hazards—demonstrating how a unified data approach can support both convective and flood event investigations.

Resources and references¶

NOAA/NSSL, 2025: Multi-Radar/Multi-Sensor System (MRMS). National Severe Storms Laboratory. https://www.nssl.noaa.gov/projects/mrms/

MRMS Cookbook

Chapter 1: Multi Radar/ Multi Sensor (MRMS) System: Overview, Case Studies, and More!

MRMS Cookbook

Chapter 3: July 2025 Central TX Floods