Nice combo example #136

from datetime import datetime

import cartopy.crs as ccrs
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
from metpy.io import parse_metar_file
import metpy.plots as mpplots
from metpy.remote import GOESArchive
import numpy as np
from siphon.catalog import TDSCatalog
import xarray as xr

# Parse the mesoscale discussion bounds
md_text = ('37190506 37460487 37670470 37810422 37860215 37720171 '
           '36970146 36160161 35430215 35170275 35700516 36000550 '
           '36730529 37190506')
md_time = datetime(2023, 5, 17, 18, 29)
md_pts = np.array([(int('-1' + item[-4:]) / 100, int(item[:4]) / 100) for item in md_text.split()])

# Grab the corresponding visible satellite image
prod = GOESArchive(satellite=16).get_product('ABI-L1b-RadC', md_time, channel=2)
vis_img = prod.parse().metpy.parse_cf('Rad')

# Get some METAR data
metar_cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/noaaport/text/metar/catalog.xml')
sfc_obs = parse_metar_file(metar_cat.datasets.filter_time_nearest(md_time).remote_open(mode='t'))

# Keep only a single ob per station
sfc_obs = sfc_obs.set_index('date_time').groupby('station_id').first()

# Get MRMS data--need to work around some xarray + TDS incompatibilities
mrms_cat = TDSCatalog('https://thredds-test.unidata.ucar.edu/thredds/catalog/grib/NCEP/MRMS/CONUS/BaseRef/MRMS_CONUS_BaseReflectivity_20230517_1800.grib2/catalog.xml')
opendap_url = mrms_cat.datasets[0].access_urls['OPENDAP']
nc = xr.open_dataset(opendap_url, drop_variables='reftime')
ref = nc.metpy.parse_cf('MergedBaseReflectivityQC_altitude_above_msl')

# Select the data that matches our MD time. Also only plot > 25 dBz
ref = ref.sel(reftime4=md_time, method='nearest').squeeze()
ref.data[ref < 25] = np.nan

# Generate a base figure
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-102, central_latitude=37))
ax.set_extent((-107, -100, 34, 39))

# Plot satellite and radar as images
ax.imshow(vis_img, extent=(vis_img.metpy.x[0], vis_img.metpy.x[-1], vis_img.metpy.y[0], vis_img.metpy.y[-1]),
          cmap='gray', transform=vis_img.metpy.cartopy_crs, origin='lower')
ref_norm, ref_cmap = mpplots.colortables.get_with_steps('NWSStormClearReflectivity', -20, 0.5)
ax.imshow(ref,
          extent=(ref.metpy.x[0] - 360, ref.metpy.x[-1] - 360, ref.metpy.y[0], ref.metpy.y[-1]),
          transform=ref.metpy.cartopy_crs, cmap=ref_cmap, norm=ref_norm, origin='lower')

# Add base map features
ax.add_feature(mpplots.USCOUNTIES, edgecolor='lightgrey', linewidth=0.5)
ax.add_feature(mpplots.USSTATES)

# Plot the surface data
sp = mpplots.StationPlot(ax, sfc_obs.longitude, sfc_obs.latitude, transform=ccrs.PlateCarree(), clip_on=True)
sp.plot_parameter('NW', sfc_obs.air_temperature, color='red')
sp.plot_parameter('SW', sfc_obs.dew_point_temperature, color='green')
sp.plot_symbol('C', sfc_obs.cloud_coverage, mpplots.sky_cover)
sp.plot_barb(sfc_obs.eastward_wind, sfc_obs.northward_wind)

# Plot the MD
md = mpatches.Polygon(md_pts, transform=ccrs.PlateCarree(), path_effects=[mpplots.ScallopedStroke()],
                      edgecolor='brown', facecolor='none', linewidth=2)
ax.add_artist(md)

# Add a bit more decoration
ax.set_title(f'SPC Mesoscale Discussion Valid: {md_time:%Y-%m-%d %H:%M}Z')
mpplots.add_metpy_logo(fig, 25, 20, size='small')