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. That tornado alone caused catastrophic damage and multiple fatalities, prompting tornado emergencies and 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, using:

Reflectivity
Rotation
Hail Swaths
Local Storm Reports

Imports¶

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

Build Map¶

This section uses Cartopy to build a map.¶

# 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()

Great, now we can highlight our case study!¶

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

# 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
cref1z = fetch_mrms_data('CREF_1HR_MAX_00.50', '20230325', '01')
cref2z = fetch_mrms_data('CREF_1HR_MAX_00.50', '20230325', '02')
# 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)	- 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¶

precip_1z = fetch_mrms_data('PrecipRate_00.00', '20230325', '01') # Precip Rate
precip_2z = fetch_mrms_data('PrecipRate_00.00', '20230325', '02')
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
QPE_1z = fetch_mrms_data('MultiSensor_QPE_01H_Pass1_00.00', '20230325', '01') # QPE: Quantified Precip Estimation - Offered hourly.
QPE_2z = fetch_mrms_data('MultiSensor_QPE_01H_Pass1_00.00', '20230325', '02') 
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
QPE_1z_p2 = fetch_mrms_data('MultiSensor_QPE_01H_Pass2_00.00', '20230325', '01') # QPE: Quantified Precip Estimation - last hour, 2nd pass
QPE_2z_p2 = fetch_mrms_data('MultiSensor_QPE_01H_Pass2_00.00', '20230325', '02') 
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¶

# RadarOnly_QPE_01H_00.00
RQPE_1z = fetch_mrms_data('RadarOnly_QPE_01H_00.00', '20230325', '01') # RadarOnly_QPE: Radar Only Quantified Precip Estimation - last hour
RQPE_2z = fetch_mrms_data('RadarOnly_QPE_01H_00.00', '20230325', '02') 
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 and turbulence.¶

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 Az Shear 3-6km AGL	Mid-level azimuthal shear (3–6 km AGL); captures elevated storm rotation.	Instantaneous	MergedAzShear_3-6kmAGL_00.50
Rotation Track 30min	30-min accumulation of low-level rotation; useful for short-term tracking.	30 minutes	RotationTrack30min_00.50
Rotation Track 60min	60-min accumulation of low-level rotation; highlights sustained activity.	60 minutes	RotationTrack60min_00.50
Rotation Track ML 30min	ML-enhanced 30-min rotation track; filters noise, boosts confidence.	30 minutes	RotationTrackML30min_00.50
Rotation Track ML 60min	ML-enhanced 60-min rotation track; detects short-lived intense rotation.	60 minutes	RotationTrackML60min_00.50

# Naming convention is completely different! Would have to completely rewrite my method to get this - get plots similar to those above

# MergedAzShear_0-2kmAGL_00.50

# MergedAzShear_3-6kmAGL_00.50 -  for the sake of time, hold on this

# RotationTrack30min_00.50 -  for the sake of time, hold on this
# RotationTrack60min_00.50 -  for the sake of time, hold on this
# 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 - Maximum Expected Size of Hail - for the sake of time, hold on this
# SHI_00.50 - Severe Hail Index
# VII_00.50 - Vergically Integrated Ice
# POSH_00.50

Storm Intensity - Vertically Integrated Liquid¶

# VIL_00.50		
# VIL_Density_00.50			
# VIL_Max_120min_00.50

Summary¶

Add one final --- marking the end of your body of content, and then conclude with a brief single paragraph summarizing at a high level the key pieces that were learned and how they tied to your objectives. Look to reiterate what the most important takeaways were.

What’s next?¶

Let Jupyter book tie this to the next (sequential) piece of content that people could move on to down below and in the sidebar. However, if this page uniquely enables your reader to tackle other nonsequential concepts throughout this book, or even external content, link to it here!

Resources and references¶

Finally, be rigorous in your citations and references as necessary. Give credit where credit is due. Also, feel free to link to relevant external material, further reading, documentation, etc. Then you’re done! Give yourself a quick review, a high five, and send us a pull request. A few final notes:

Kernel > Restart Kernel and Run All Cells... to confirm that your notebook will cleanly run from start to finish
Kernel > Restart Kernel and Clear All Outputs... before committing your notebook, our machines will do the heavy lifting
Take credit! Provide author contact information if you’d like; if so, consider adding information here at the bottom of your notebook
Give credit! Attribute appropriate authorship for referenced code, information, images, etc.
Only include what you’re legally allowed: no copyright infringement or plagiarism

Thank you for your contribution!

MRMS Cookbook

Project Pythia Notebook Template