Reflectivity_Script.py

# Reflectivity_Script.py
# Michael Wasserstein
# 12/28/2025
#
# Retrieves the most recent KMTX (WSR-88D) level-2 scan from the NOAA NEXRAD
# AWS archive, computes the composite reflectivity, produces a plan-view map,
# and generates range-height cross sections along the azimuths toward two
# mountain sites (Snowbasin / SNI and Ben Lomond Peak / BLP).
#
# External dependency: map_script_2.py (not included here) must be importable
# from the same directory or your PYTHONPATH.  It supplies the topography grid
# (lons, lats, topo) and the Great Salt Lake geometry (lakes_gdf).
#
# AWS access uses the legacy boto library (boto.s3.connection).  If you prefer
# the newer boto3 SDK, replace the S3Connection block with a boto3 equivalent.

import os
import datetime
import tempfile
import numpy as np
import pandas as pd

import pyart
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

from boto.s3.connection import S3Connection

# map_script_2 provides: lons, lats, topo (topography grid) and lakes_gdf
from map_script_2 import *


# --------------------------------------------------------------------------- #
# Configuration
# --------------------------------------------------------------------------- #
radar_id = 'KMTX'

# Output directory for saved figures
YYYY     = datetime.datetime.utcnow().strftime('%Y')
mm       = datetime.datetime.utcnow().strftime('%m')
dd       = datetime.datetime.utcnow().strftime('%d')
FigDir   = f'/path/to/your/output/{YYYY}{mm}{dd}/{radar_id}/'

os.makedirs(FigDir, exist_ok=True)

# Azimuth, distance-to-site (km), site lat/lon, cross-section end lat/lon
site_dict = {
    'SNI': [97.9,  49.8,  41.1995,  -111.8590, 41.13031, -111.25546],  # Snowbasin middle bowl
    'BLP': [73.1,  43.96, 41.37603, -111.94405, 41.51578, -111.32309], # Ben Lomond Peak SNOTEL
}

# Matplotlib style
props  = dict(boxstyle='square', facecolor='white', alpha=0.8,  ec="gray")
props2 = dict(boxstyle='square', facecolor='white', alpha=0.9,  ec="none")

topo_levels = np.arange(1200, 3000, 500)  # contour levels for topography (m MSL)

# --------------------------------------------------------------------------- #
# Retrieve most recent KMTX scan from AWS
# --------------------------------------------------------------------------- #
conn   = S3Connection(anon=True)
bucket = conn.get_bucket('unidata-nexrad-level2')

curr_time = datetime.datetime.utcnow()

aws_directory = '{:04d}/{:02d}/{:02d}/'.format(
    curr_time.year, curr_time.month, curr_time.day) + radar_id + '/'
file_list = bucket.list(aws_directory, "/")

# Step back one day at a time if no files exist for today (e.g., just after 00 UTC)
while len(list(file_list)) == 0:
    print('No radar files for {}/{}/{}, stepping back one day'.format(
        curr_time.year, curr_time.month, curr_time.day))
    curr_time = curr_time - datetime.timedelta(days=1)
    aws_directory = '{:04d}/{:02d}/{:02d}/'.format(
        curr_time.year, curr_time.month, curr_time.day) + radar_id + '/'
    file_list = bucket.list(aws_directory, "/")

radfile = list(file_list)[-1]
if "MDM" in radfile.key:   # skip the MDM supplemental file
    radfile = list(file_list)[-2]

localfile = tempfile.NamedTemporaryFile()
radfile.get_contents_to_filename(localfile.name)
radar = pyart.io.read_nexrad_archive(localfile.name)
localfile.close()

# Extract metadata
date           = radar.time['units'][14:].split('T')[0]
time           = radar.time['units'][14:].split('T')[1]
valid_time     = date + ' ' + time
save_time      = pd.to_datetime(valid_time).strftime('%Y%m%d%H%M%S')
radar_location = radar.metadata['instrument_name']
type_of_data   = radar.metadata['original_container']
radar_lat      = radar.latitude['data'][0]
radar_lon      = radar.longitude['data'][0]
rad_altitude   = radar.altitude['data'][0]

# --------------------------------------------------------------------------- #
# Composite reflectivity map
# --------------------------------------------------------------------------- #
composite     = pyart.retrieve.composite_reflectivity(radar)
lat_composite = composite.get_gate_lat_lon_alt(sweep=0)[0]
lon_composite = composite.get_gate_lat_lon_alt(sweep=0)[1]
composite_ref = composite.get_field(field_name='composite_reflectivity', sweep=0)

cmap = plt.get_cmap('pyart_LangRainbow12')
cmap.set_under('none')

if not os.path.exists(FigDir + 'composite_reflectivity_' + save_time + '.png'):

    fig, ax1 = plt.subplots(
        1, 1, figsize=(10, 8),
        subplot_kw={'projection': ccrs.PlateCarree()},
        facecolor='white', edgecolor='k'
    )

    plot = ax1.pcolormesh(lon_composite, lat_composite, composite_ref,
                          cmap=cmap, vmin=0, vmax=40, zorder=15)

    ax1.add_geometries(lakes_gdf.geometry, ccrs.PlateCarree(),
                       zorder=150, facecolor='none', edgecolor='black', linewidth=2.5)

    ax1.contour(lons, lats, topo, zorder=100, levels=topo_levels,
                cmap='binary', linewidths=1.3)

    ax1.set_title('Valid: ' + valid_time, loc='right')
    ax1.set_title(radar_location + ' ' + type_of_data + '\nComposite Reflectivity', loc='left')

    for site in ('SNI', 'BLP'):
        ax1.plot([radar_lon, site_dict[site][5]], [radar_lat, site_dict[site][4]],
                 zorder=101, color='#676e69', linewidth=4,
                 solid_capstyle='round', linestyle='-')
        ax1.text(site_dict[site][3], site_dict[site][2], site,
                 bbox=props, fontsize=17, zorder=5000,
                 horizontalalignment='center', verticalalignment='center',
                 transform=ccrs.PlateCarree())

    ax1.set_xlim(-113, -111)
    ax1.set_ylim(40.3, 41.8)

    cax  = plt.axes([0.92, 0.132, 0.02, 0.73])
    cbar = plt.colorbar(plot, cax=cax, orientation='vertical')
    cbar.ax.set_ylabel('dBZ', rotation=0, fontsize=16, labelpad=15)
    cbar.ax.tick_params(labelsize=14, size=8)

    plt.savefig(FigDir + 'composite_reflectivity_' + save_time + '.png',
                dpi=300, bbox_inches='tight')
    plt.close()


# --------------------------------------------------------------------------- #
# Range-height cross sections
# --------------------------------------------------------------------------- #
def produce_cross_sections(site, radar):
    """
    Plot a range-height cross section of reflectivity along a fixed azimuth.

    The terrain profile along the azimuth is loaded from pre-computed .npy files
    (BLP_distance_along_azimuth.npy / BLP_elevation_along_azimuth.npy, etc.)
    that are included with this repository.

    Inputs:
    site  - Site key ('SNI' or 'BLP')
    radar - Py-ART radar object for the current scan
    """
    target_az      = site_dict[site][0]
    dist_to_point  = site_dict[site][1]
    site_lat       = site_dict[site][2]
    site_lon       = site_dict[site][3]

    rad_graph = pyart.graph.RadarMapDisplay(radar)

    distance_along_azimuth  = np.load(f'data/{site}_distance_along_azimuth.npy')
    elevation_along_azimuth = np.load(f'data/{site}_elevation_along_azimuth.npy')

    cmap          = plt.get_cmap('pyart_LangRainbow12')
    yticks_labels = [-1000, 0, 1000, 2000, 3000, 4000, 5000, 6000]
    yticks        = np.arange(-1, 6.5, 1)
    vmin, vmax    = 0, 40

    fig, ax1 = plt.subplots(1, 1, facecolor='white', edgecolor='k', figsize=(16, 14))

    plot = rad_graph.plot_azimuth_to_rhi(
        field='reflectivity', target_azimuth=target_az, ax=ax1,
        vmin=vmin, vmax=vmax, cmap=cmap, title='',
        colorbar_flag=False, colorbar_label='dBZ'
    )

    # Terrain profile (subtract radar altitude because y-axis is distance above radar)
    ax1.plot(distance_along_azimuth,
             (elevation_along_azimuth - rad_altitude) / 1000,
             color='black', linewidth=2)

    ax1.set_title('Valid: ' + valid_time, loc='right', fontsize=24)
    ax1.set_title(
        radar_location + ' ' + type_of_data +
        '\nCross Section along ' + str(target_az) + '° Azimuth',
        loc='left', fontsize=24
    )

    ax1.set_xlim(0, 80)
    ax1.set_ylim(-1, 6)
    ax1.set_yticks(yticks, labels=yticks_labels)
    ax1.tick_params(axis='both', size=8, labelsize=20)
    ax1.set_xlabel('Distance from Radar (km)', fontsize=24)
    ax1.set_ylabel('Distance above Radar (m)', fontsize=24)

    # Second y-axis showing height MSL
    ax2 = ax1.twinx()
    ax2.tick_params(axis='both', size=8, labelsize=20)
    ax2.set_ylim(rad_altitude - 1000, rad_altitude + 6000)
    ax2.tick_params(size=8)
    ax2.set_ylabel('Height (m MSL)', fontsize=24, rotation=270, labelpad=30)

    cax  = plt.axes([0.99, 0.11, 0.02, 0.77])
    cbar = fig.colorbar(plot, cmap=cmap, cax=cax)
    cbar.ax.set_ylabel('dBZ', rotation=0, fontsize=24, labelpad=30)
    cbar.ax.tick_params(labelsize=14, size=8)
    cbar.ax.set_yticks(np.linspace(0, 1.0, 9),
                       labels=np.arange(vmin, vmax + 0.01, 5).astype(int),
                       fontsize=24)

    # Dashed vertical line marking the site location
    ax1.vlines(dist_to_point, -5, 20, color='black', linestyle='--')
    ax1.text(dist_to_point, 2.9, site,
             bbox=props, fontsize=24, zorder=20, horizontalalignment='center')

    plt.savefig(FigDir + f'cross_section_{site}_' + save_time + '.png',
                dpi=300, bbox_inches='tight')
    plt.close()


for site in ('BLP', 'SNI'):
    if not os.path.exists(FigDir + f'cross_section_{site}_' + save_time + '.png'):
        produce_cross_sections(site=site, radar=radar)