xjtulyc/metpy-atmosphere

MetPy Atmospheric Analysis Skill

This skill covers the full workflow for operational and research-grade atmospheric science in Python: programmatic data retrieval from UCAR THREDDS and NOAA NOMADS servers using Siphon, thermodynamic parcel analysis (CAPE, CIN, LCL, LFC, EL) using MetPy, skew-T log-P diagram generation, wind hodograph construction, and synoptic-scale composite analysis over gridded model output.

MetPy is the de-facto standard Python library for atmospheric calculations. It handles unit-aware arrays via Pint and provides a comprehensive set of thermodynamic, kinematic, and interpolation routines. Siphon acts as the data-access layer that queries THREDDS Data Server (TDS) catalogs and returns xarray datasets with CF-compliant metadata.

---

When to Use This Skill

• You need to retrieve real-time or archived NWP model data (GFS, NAM, RAP) programmatically without manual download.
• You are producing upper-air soundings or comparing observed radiosondes with model profiles.
• You need convective instability metrics (CAPE, CIN, LIFTED index, K-index) for a case study or climatological analysis.
• You want to create publication-quality synoptic maps (500 hPa heights, SLP, moisture flux).
• You are building an automated workflow that runs nightly and ingests the latest model cycle.

---

Background & Key Concepts

Skew-T Log-P Diagram

A thermodynamic diagram where temperature is plotted on a skewed x-axis and pressure decreases logarithmically upward. Dry adiabats, moist adiabats, and saturation mixing-ratio lines help visualize parcel ascent and atmospheric stability.

CAPE and CIN

• CAPE (Convective Available Potential Energy): the positive buoyancy energy a parcel gains ascending from the LFC to the EL. Units are J kg⁻¹. Values > 1000 J kg⁻¹ indicate significant severe-weather potential.
• CIN (Convective INhibition): the negative area below the LFC that must be overcome before free convection begins. Large CIN suppresses convection; small CIN allows it.

THREDDS / NOMADS

The UCAR THREDDS Data Server hosts model grids and observational datasets accessible via OPeNDAP, HTTP, and WMS. NOAA NOMADS is the National Operational Model Archive and Distribution System providing real-time GFS, NAM, and RAP output. Siphon wraps both services with a Pythonic catalog interface.

Hodograph

A polar diagram of the wind vector at successive altitude levels. The shape of the hodograph quantifies wind shear and storm-relative helicity (SRH), key ingredients for supercell thunderstorm development.

---

Environment Setup

Install Dependencies

# Create and activate a conda environment (recommended for Cartopy)
conda create -n atmos python=3.11
conda activate atmos
conda install -c conda-forge metpy siphon xarray cartopy matplotlib numpy pandas
pip install metpy>=1.5 siphon>=0.9

# Alternatively, pure pip (may require system GEOS/PROJ for Cartopy)
pip install metpy>=1.5 siphon>=0.9 xarray>=2023.6 matplotlib>=3.7 \
            numpy>=1.24 pandas>=2.0 cartopy>=0.22

Verify Installation

import metpy
import siphon
import xarray as xr
import cartopy

print(f"MetPy    {metpy.__version__}")
print(f"Siphon   {siphon.__version__}")
print(f"xarray   {xr.__version__}")
print(f"Cartopy  {cartopy.__version__}")

Optional: NOMADS API Access

NOMADS is public and requires no key. THREDDS catalogs at UCAR are also open. No authentication is needed for the examples below.

---

Core Workflow

Step 1 — Retrieve Sounding Data from THREDDS

"""
step1_siphon_sounding.py
------------------------
Retrieve the latest GFS 0-h analysis profile at a single point via the
UCAR THREDDS server using Siphon, then parse it into a tidy DataFrame.
"""

import os
import warnings
from datetime import datetime, timedelta, timezone

import numpy as np
import pandas as pd
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS

warnings.filterwarnings("ignore")


def fetch_gfs_sounding(
    lat: float,
    lon: float,
    date: datetime | None = None,
    variables: list[str] | None = None,
) -> pd.DataFrame:
    """
    Fetch a vertical profile from the GFS 0.25-degree analysis via THREDDS/NCSS.

    Parameters
    ----------
    lat, lon : float
        Station latitude and longitude in decimal degrees.
    date : datetime, optional
        UTC time of the desired model run. Defaults to the most recent 00Z cycle.
    variables : list of str, optional
        List of NCSS variable names. Defaults to temperature, dewpoint, and wind.

    Returns
    -------
    pd.DataFrame
        Columns: pressure [hPa], temperature [degC], dewpoint [degC],
        u_wind [m/s], v_wind [m/s], height [m].
    """
    if date is None:
        now = datetime.now(tz=timezone.utc)
        # Roll back to the most recent 00Z or 12Z cycle
        hour = (now.hour // 12) * 12
        date = now.replace(hour=hour, minute=0, second=0, microsecond=0)

    if variables is None:
        variables = [
            "Temperature_isobaric",
            "Dewpoint_temperature_isobaric",
            "u-component_of_wind_isobaric",
            "v-component_of_wind_isobaric",
            "Geopotential_height_isobaric",
        ]

    catalog_url = (
        "https://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/"
        "Global_0p25deg/catalog.xml"
    )

    print(f"Connecting to THREDDS catalog: {catalog_url}")
    cat = TDSCatalog(catalog_url)
    # Select the latest dataset
    ds_name = sorted(cat.datasets.keys())[-1]
    print(f"Using dataset: {ds_name}")

    ncss = NCSS(cat.datasets[ds_name].access_urls["NetcdfSubset"])
    query = ncss.query()
    query.lonlat_point(lon, lat)
    query.variables(*variables)
    query.vertical_level_range(100, 1050)  # hPa
    query.accept("netcdf4")

    raw = ncss.get_data(query)

    # Parse pressure levels and profile arrays
    pressure = raw.variables["pressure"][:]  # hPa
    temp_k = raw.variables["Temperature_isobaric"][0, :, 0, 0]
    td_k = raw.variables["Dewpoint_temperature_isobaric"][0, :, 0, 0]
    u = raw.variables["u-component_of_wind_isobaric"][0, :, 0, 0]
    v = raw.variables["v-component_of_wind_isobaric"][0, :, 0, 0]
    z = raw.variables["Geopotential_height_isobaric"][0, :, 0, 0]

    df = pd.DataFrame({
        "pressure": pressure,
        "temperature": temp_k - 273.15,
        "dewpoint": td_k - 273.15,
        "u_wind": u,
        "v_wind": v,
        "height": z,
    }).sort_values("pressure", ascending=False).reset_index(drop=True)

    print(f"Retrieved {len(df)} pressure levels for ({lat:.2f}N, {lon:.2f}E)")
    return df


if __name__ == "__main__":
    # Example: profile over central Oklahoma (thunderstorm alley)
    sounding_df = fetch_gfs_sounding(lat=35.5, lon=-97.5)
    print(sounding_df.head(10).to_string(index=False))
    sounding_df.to_csv("gfs_sounding_okc.csv", index=False)
    print("Saved: gfs_sounding_okc.csv")

Step 2 — Skew-T Log-P Diagram and Parcel Analysis

"""
step2_skewt_parcel.py
---------------------
Build a skew-T log-P diagram from a sounding DataFrame and overlay
the parcel ascent path.  Compute CAPE, CIN, LCL, LFC, and EL.
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units


def load_sounding(csv_path: str) -> tuple:
    """
    Load a sounding CSV (columns: pressure, temperature, dewpoint, u_wind, v_wind)
    and attach MetPy units.

    Returns
    -------
    tuple of pint Quantity arrays: pressure, temperature, dewpoint, u, v
    """
    df = pd.read_csv(csv_path)

    # Drop levels with missing temperature or dewpoint
    df = df.dropna(subset=["temperature", "dewpoint"]).copy()
    df = df.sort_values("pressure", ascending=False).reset_index(drop=True)

    pressure = df["pressure"].values * units.hPa
    temperature = df["temperature"].values * units.degC
    dewpoint = df["dewpoint"].values * units.degC
    u = df["u_wind"].values * units("m/s")
    v = df["v_wind"].values * units("m/s")

    return pressure, temperature, dewpoint, u, v


def plot_skewt(
    pressure, temperature, dewpoint, u, v,
    title: str = "Skew-T Log-P Diagram",
    output_path: str = "skewt.png",
) -> plt.Figure:
    """
    Produce a complete skew-T log-P diagram with parcel trace and indices.

    The figure includes:
    - Observed temperature and dewpoint profiles
    - Surface-based parcel ascent path
    - CAPE (green shading) and CIN (red shading) areas
    - Wind barbs on the right side
    - Text box with convective indices

    Parameters
    ----------
    pressure, temperature, dewpoint : pint Quantity arrays
        Mandatory profile data.
    u, v : pint Quantity arrays
        Wind components in m/s.
    title : str
        Plot title.
    output_path : str
        File path to save the PNG figure.

    Returns
    -------
    matplotlib.figure.Figure
    """
    # ----- Parcel analysis -----
    parcel_prof = mpcalc.parcel_profile(pressure, temperature[0], dewpoint[0])

    cape, cin = mpcalc.cape_cin(pressure, temperature, dewpoint, parcel_prof)
    lcl_p, lcl_t = mpcalc.lcl(pressure[0], temperature[0], dewpoint[0])
    lfc_p, lfc_t = mpcalc.lfc(pressure, temperature, dewpoint)
    el_p, el_t = mpcalc.el(pressure, temperature, dewpoint)

    # ----- Figure layout -----
    fig = plt.figure(figsize=(10, 12))
    skew = SkewT(fig, rotation=45)

    # Background lines
    skew.plot_dry_adiabats(alpha=0.25, colors="orange")
    skew.plot_moist_adiabats(alpha=0.25, colors="green")
    skew.plot_mixing_lines(alpha=0.25, colors="blue")

    # Observed profiles
    skew.plot(pressure, temperature, "r", linewidth=2, label="Temperature")
    skew.plot(pressure, dewpoint, "g", linewidth=2, label="Dewpoint")

    # Parcel trace
    skew.plot(pressure, parcel_prof, "k--", linewidth=1.5, label="Parcel")

    # CAPE / CIN shading
    skew.shade_cape(pressure, temperature, parcel_prof)
    skew.shade_cin(pressure, temperature, parcel_prof, dewpoint)

    # Key levels
    skew.ax.axhline(lcl_p.m, color="purple", linestyle=":", linewidth=1.2, label="LCL")
    if not np.isnan(lfc_p.m):
        skew.ax.axhline(lfc_p.m, color="darkorange", linestyle="-.", linewidth=1.2, label="LFC")
    if not np.isnan(el_p.m):
        skew.ax.axhline(el_p.m, color="cyan", linestyle="-.", linewidth=1.2, label="EL")

    # Wind barbs
    wind_slice = slice(None, None, 5)
    skew.plot_barbs(
        pressure[wind_slice], u[wind_slice], v[wind_slice],
        x_clip_radius=0.12, y_clip_radius=0.12,
    )

    # Indices text box
    indices_text = (
        f"CAPE: {cape.m:.0f} J kg⁻¹\n"
        f"CIN:  {cin.m:.0f} J kg⁻¹\n"
        f"LCL:  {lcl_p.m:.0f} hPa / {lcl_t.m:.1f} °C\n"
        f"LFC:  {lfc_p.m:.0f} hPa\n"
        f"EL:   {el_p.m:.0f} hPa"
    )
    skew.ax.text(
        0.02, 0.98, indices_text,
        transform=skew.ax.transAxes,
        fontsize=9, verticalalignment="top",
        bbox=dict(boxstyle="round,pad=0.4", facecolor="wheat", alpha=0.8),
    )

    skew.ax.set_xlim(-40, 50)
    skew.ax.set_ylim(1050, 100)
    skew.ax.set_xlabel("Temperature (°C)")
    skew.ax.set_ylabel("Pressure (hPa)")
    skew.ax.set_title(title, fontsize=13)
    skew.ax.legend(loc="upper right", fontsize=8)

    fig.tight_layout()
    fig.savefig(output_path, dpi=150, bbox_inches="tight")
    print(f"Saved: {output_path}")
    print(f"\nConvective indices:")
    print(f"  CAPE = {cape.m:.1f} J/kg")
    print(f"  CIN  = {cin.m:.1f} J/kg")
    print(f"  LCL  = {lcl_p.m:.1f} hPa, {lcl_t.m:.1f} °C")

    return fig


if __name__ == "__main__":
    p, t, td, u, v = load_sounding("gfs_sounding_okc.csv")
    fig = plot_skewt(
        p, t, td, u, v,
        title="GFS Analysis – Central Oklahoma\nSkew-T Log-P",
        output_path="skewt_okc.png",
    )
    plt.show()

Step 3 — Wind Hodograph and Storm-Relative Helicity

"""
step3_hodograph.py
------------------
Construct a wind hodograph and compute storm-relative helicity (SRH)
for the 0–1 km and 0–3 km layers.
"""

import numpy as np
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import Hodograph
from metpy.units import units
import pandas as pd


def plot_hodograph(
    pressure, u, v, height,
    title: str = "Wind Hodograph",
    output_path: str = "hodograph.png",
) -> plt.Figure:
    """
    Plot a wind hodograph colored by altitude layer and annotate SRH values.

    Parameters
    ----------
    pressure : pint Quantity (hPa)
    u, v : pint Quantity (m/s)
    height : pint Quantity (m AGL)
    title : str
    output_path : str

    Returns
    -------
    matplotlib.figure.Figure
    """
    # Compute SRH for key layers
    srh_1km = mpcalc.storm_relative_helicity(
        height, u, v, depth=1000 * units.m
    )
    srh_3km = mpcalc.storm_relative_helicity(
        height, u, v, depth=3000 * units.m
    )

    # Bulk wind shear 0–6 km
    shear_u, shear_v = mpcalc.bulk_shear(pressure, u, v, depth=6000 * units.m)
    bws = np.sqrt(shear_u**2 + shear_v**2).to("kt")

    fig, ax = plt.subplots(1, 1, figsize=(8, 8))
    h = Hodograph(ax, component_range=80)
    h.add_grid(increment=10)

    # Color-coded by height layer
    intervals = np.array([0, 1000, 3000, 6000, 12000]) * units.m
    colors = ["purple", "red", "darkorange", "gold"]
    h.plot_colormapped(u, v, height, intervals=intervals, colors=colors)

    # Annotate wind speeds at key levels (surface, 850, 700, 500 hPa)
    key_levels = [0, 850, 700, 500]
    p_vals = pressure.m
    for target_p in key_levels:
        idx = np.argmin(np.abs(p_vals - target_p))
        ax.annotate(
            f"{target_p} hPa" if target_p > 0 else "Sfc",
            xy=(u[idx].m, v[idx].m),
            fontsize=7, color="black",
            xytext=(3, 3), textcoords="offset points",
        )

    # Annotations
    info = (
        f"0–1 km SRH: {srh_1km[0].m:.0f} m²/s²\n"
        f"0–3 km SRH: {srh_3km[0].m:.0f} m²/s²\n"
        f"0–6 km BWS: {bws.m:.1f} kt"
    )
    ax.text(
        0.02, 0.98, info,
        transform=ax.transAxes, fontsize=9,
        verticalalignment="top",
        bbox=dict(boxstyle="round,pad=0.4", facecolor="lightyellow", alpha=0.9),
    )

    ax.set_title(title, fontsize=12)
    fig.tight_layout()
    fig.savefig(output_path, dpi=150, bbox_inches="tight")
    print(f"Saved: {output_path}")
    return fig


if __name__ == "__main__":
    df = pd.read_csv("gfs_sounding_okc.csv").dropna()
    df = df.sort_values("pressure", ascending=False).reset_index(drop=True)

    p = df["pressure"].values * units.hPa
    u = df["u_wind"].values * units("m/s")
    v = df["v_wind"].values * units("m/s")
    z = df["height"].values * units.m

    fig = plot_hodograph(p, u, v, z,
                         title="GFS – OKC Hodograph",
                         output_path="hodograph_okc.png")
    plt.show()

Step 4 — Synoptic Composite Map (500 hPa Heights + Vorticity)

"""
step4_synoptic_map.py
---------------------
Download a GFS 500 hPa analysis grid from NOMADS and plot geopotential
height contours plus absolute vorticity with a Cartopy basemap.
"""

import os
import warnings
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy.calc as mpcalc
from metpy.units import units

warnings.filterwarnings("ignore")

# NOMADS real-time GFS 0.25-deg via OPeNDAP — no authentication required
# Adjust the cycle date/time as needed; this URL pattern works for recent runs.
NOMADS_BASE = (
    "https://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs{date}/gfs_0p25_00z"
)


def fetch_500hpa_grid(date_str: str) -> xr.Dataset:
    """
    Open GFS 500 hPa geopotential height and wind via OPeNDAP from NOMADS.

    Parameters
    ----------
    date_str : str
        Date in YYYYMMDD format, e.g. '20240601'.

    Returns
    -------
    xr.Dataset
        Dataset with hgtprs (geopotential height) and ugrdprs/vgrdprs (winds).
    """
    url = NOMADS_BASE.format(date=date_str)
    print(f"Opening: {url}")
    ds = xr.open_dataset(url, engine="netcdf4")
    # Select 500 hPa level (lev coordinate in hPa) and first forecast hour
    ds_500 = ds.sel(lev=500, method="nearest").isel(time=0)
    return ds_500


def plot_500hpa_map(
    ds: xr.Dataset,
    extent: list[float] = [-130, -60, 20, 60],
    output_path: str = "500hpa_analysis.png",
) -> plt.Figure:
    """
    Plot 500 hPa geopotential height contours and absolute vorticity fill.

    Parameters
    ----------
    ds : xr.Dataset
        Dataset with hgtprs, ugrdprs, vgrdprs, lat, lon variables.
    extent : list of float
        Map extent [lon_min, lon_max, lat_min, lat_max].
    output_path : str
        Output file path.

    Returns
    -------
    matplotlib.figure.Figure
    """
    lat = ds["lat"].values
    lon = ds["lon"].values
    lon2d, lat2d = np.meshgrid(lon, lat)

    hgt = ds["hgtprs"].values          # geopotential height in meters
    u = ds["ugrdprs"].values * units("m/s")
    v = ds["vgrdprs"].values * units("m/s")

    # Compute absolute vorticity
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
    avor = mpcalc.absolute_vorticity(u, v, dx=dx, dy=dy,
                                      latitude=lat2d * units.degrees)
    avor_scaled = avor.m * 1e5  # units of 10^-5 s^-1

    # ---- Figure ----
    proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35)
    fig, ax = plt.subplots(1, 1, figsize=(14, 9), subplot_kw={"projection": proj})
    ax.set_extent(extent, crs=ccrs.PlateCarree())

    ax.add_feature(cfeature.COASTLINE, linewidth=0.6)
    ax.add_feature(cfeature.BORDERS, linewidth=0.4)
    ax.add_feature(cfeature.STATES, linewidth=0.3, alpha=0.5)

    # Vorticity fill
    vort_levels = np.arange(0, 50, 4)
    cf = ax.contourf(
        lon2d, lat2d, avor_scaled,
        levels=vort_levels, cmap="YlOrRd",
        transform=ccrs.PlateCarree(), extend="max", alpha=0.7,
    )
    plt.colorbar(cf, ax=ax, orientation="horizontal", pad=0.04,
                 label="Absolute Vorticity (×10⁻⁵ s⁻¹)", shrink=0.7)

    # Height contours
    hgt_levels = np.arange(4800, 6000, 60)
    cs = ax.contour(
        lon2d, lat2d, hgt,
        levels=hgt_levels, colors="black", linewidths=1.0,
        transform=ccrs.PlateCarree(),
    )
    ax.clabel(cs, fmt="%d", fontsize=8, inline=True)

    gl = ax.gridlines(draw_labels=True, linewidth=0.4, color="gray", alpha=0.5)
    gl.top_labels = False
    gl.right_labels = False

    ax.set_title("GFS 500 hPa Geopotential Height (dam) and Absolute Vorticity",
                 fontsize=13, pad=10)

    fig.tight_layout()
    fig.savefig(output_path, dpi=150, bbox_inches="tight")
    print(f"Saved: {output_path}")
    return fig

---

Advanced Usage

Composite Soundings: CAPE Climatology Over a Season

"""
advanced_cape_climatology.py
-----------------------------
Iterate over a list of dates, download GFS soundings, compute CAPE/CIN,
and build a monthly climatology DataFrame.
"""

import warnings
import numpy as np
import pandas as pd
from datetime import datetime, timedelta
import metpy.calc as mpcalc
from metpy.units import units

warnings.filterwarnings("ignore")


def compute_cape_cin_from_df(df: pd.DataFrame) -> dict:
    """
    Compute CAPE and CIN from a sounding DataFrame.

    Parameters
    ----------
    df : pd.DataFrame
        Columns: pressure [hPa], temperature [°C], dewpoint [°C].

    Returns
    -------
    dict with keys 'cape', 'cin', 'lcl_p', 'lifted_index'.
    """
    df = df.dropna(subset=["temperature", "dewpoint"])
    df = df.sort_values("pressure", ascending=False).reset_index(drop=True)

    p = df["pressure"].values * units.hPa
    t = df["temperature"].values * units.degC
    td = df["dewpoint"].values * units.degC

    parcel = mpcalc.parcel_profile(p, t[0], td[0])
    cape, cin = mpcalc.cape_cin(p, t, td, parcel)
    lcl_p, _ = mpcalc.lcl(p[0], t[0], td[0])
    li = mpcalc.lifted_index(p, t, parcel)

    return {
        "cape": float(cape.m),
        "cin": float(cin.m),
        "lcl_p": float(lcl_p.m),
        "lifted_index": float(li.m),
    }


def seasonal_cape_summary(records: list[dict]) -> pd.DataFrame:
    """
    Build a summary table from a list of per-sounding dicts.

    Each dict must have at least 'date', 'cape', 'cin', 'lifted_index'.
    """
    df = pd.DataFrame(records)
    df["date"] = pd.to_datetime(df["date"])
    df["month"] = df["date"].dt.month

    monthly = (
        df.groupby("month")[["cape", "cin", "lifted_index"]]
        .agg(["mean", "std", "max"])
        .round(1)
    )
    return monthly


# --- Demonstration with synthetic data (replace with real sounding fetches) ---
if __name__ == "__main__":
    rng = np.random.default_rng(42)
    n = 180  # six months of daily soundings
    start = datetime(2023, 4, 1)
    records = []
    for i in range(n):
        # Synthetic CAPE: higher in summer, lower in spring/autumn
        day = start + timedelta(days=i)
        base_cape = 800 + 600 * np.sin(np.pi * i / 180)
        records.append({
            "date": day,
            "cape": max(0, base_cape + rng.normal(0, 200)),
            "cin": abs(rng.normal(80, 40)),
            "lifted_index": rng.normal(-3 - base_cape / 600, 1),
        })

    summary = seasonal_cape_summary(records)
    print("Monthly CAPE/CIN/LI summary (Apr – Sep 2023):")
    print(summary.to_string())

Custom Parcel: Most-Unstable CAPE

"""
advanced_mucape.py
------------------
Compute Most-Unstable CAPE (MUCAPE) by lifting the parcel from the
most buoyant 100-hPa mixed layer in the lowest 300 hPa.
"""

import numpy as np
import metpy.calc as mpcalc
from metpy.units import units


def most_unstable_cape(pressure, temperature, dewpoint) -> tuple:
    """
    Find the most unstable parcel in the lowest 300 hPa and return its CAPE/CIN.

    Parameters
    ----------
    pressure, temperature, dewpoint : pint Quantity arrays
        Full sounding arrays sorted surface-first (decreasing altitude).

    Returns
    -------
    tuple: (mucape, mucin, parcel_pressure)
        All pint Quantities.
    """
    # Restrict search to lowest 300 hPa
    sfc_p = pressure[0]
    mask = pressure >= (sfc_p - 300 * units.hPa)

    p_layer = pressure[mask]
    t_layer = temperature[mask]
    td_layer = dewpoint[mask]

    best_cape = 0.0 * units("J/kg")
    best_cin = 0.0 * units("J/kg")
    best_idx = 0

    for i in range(len(p_layer)):
        try:
            parcel = mpcalc.parcel_profile(pressure, t_layer[i], td_layer[i])
            c, ci = mpcalc.cape_cin(pressure, temperature, dewpoint, parcel)
            if c > best_cape:
                best_cape = c
                best_cin = ci
                best_idx = i
        except Exception:
            continue

    return best_cape, best_cin, p_layer[best_idx]


if __name__ == "__main__":
    # Example: load the previously saved sounding
    import pandas as pd

    df = pd.read_csv("gfs_sounding_okc.csv").dropna()
    df = df.sort_values("pressure", ascending=False).reset_index(drop=True)

    p = df["pressure"].values * units.hPa
    t = df["temperature"].values * units.degC
    td = df["dewpoint"].values * units.degC

    mucape, mucin, mu_p = most_unstable_cape(p, t, td)
    print(f"MUCAPE = {mucape.m:.1f} J/kg")
    print(f"MUCIN  = {mucin.m:.1f} J/kg")
    print(f"Most unstable parcel lifted from {mu_p.m:.0f} hPa")

---

Troubleshooting

Connection Errors to THREDDS / NOMADS

The UCAR THREDDS server and NOMADS occasionally time out under heavy load or during maintenance windows. Strategies:

import time
from siphon.catalog import TDSCatalog

def robust_catalog(url: str, retries: int = 3, wait: float = 5.0) -> TDSCatalog:
    """Retry a TDS catalog fetch up to `retries` times with a delay."""
    for attempt in range(1, retries + 1):
        try:
            return TDSCatalog(url)
        except Exception as exc:
            print(f"Attempt {attempt}/{retries} failed: {exc}")
            if attempt < retries:
                time.sleep(wait)
    raise RuntimeError(f"Could not connect to {url} after {retries} attempts.")

Units Mismatch Errors

MetPy raises pint.errors.DimensionalityError when incompatible units are combined. Always attach units immediately after loading data:

# Wrong — will raise error later:
temp_raw = df["temperature"].values        # plain numpy array

# Correct:
temp = df["temperature"].values * units.degC

Parcel Analysis Fails Near the Surface

If mpcalc.parcel_profile raises or returns NaN arrays, check that the surface level (lowest-pressure index) has valid temperature and dewpoint:

print(f"Surface: p={p[0]:.1f}, T={t[0]:.1f}, Td={td[0]:.1f}")
assert not np.isnan(t[0].m), "Surface temperature is NaN"
assert not np.isnan(td[0].m), "Surface dewpoint is NaN"
assert td[0] <= t[0], "Dewpoint exceeds temperature at surface"

Cartopy Installation Issues on Windows

Cartopy requires GEOS and PROJ. On Windows, use conda:

conda install -c conda-forge cartopy

Large Grid Memory Usage

For global GFS grids (1440 × 721), load only the variables and region you need:

ds = xr.open_dataset(url, engine="netcdf4")
ds_sub = ds.sel(lat=slice(20, 60), lon=slice(230, 300))

---

External Resources

• MetPy documentation: https://unidata.github.io/MetPy/latest/
• MetPy example gallery: https://unidata.github.io/MetPy/latest/examples/index.html
• Siphon documentation: https://unidata.github.io/siphon/latest/
• UCAR THREDDS catalog browser: https://thredds.ucar.edu/thredds/catalog.html
• NOAA NOMADS server: https://nomads.ncep.noaa.gov/
• AMS Glossary of Meteorology: https://glossary.ametsoc.org/wiki/Main_Page
• Bluestein (1992) Synoptic-Dynamic Meteorology in Midlatitudes, Oxford University Press
• Doswell & Rasmussen (1994) The effect of hodograph smoothing on derived thermodynamic parameters. Wea. Forecasting, 9, 576–593.

---

Examples

Example 1: Full Sounding Workflow — OKC Tornado Outbreak Day

This example combines Steps 1–3 to reproduce a textbook severe-weather sounding from a historical tornado outbreak day by querying archived GFS analyses.

"""
example1_tornado_day_sounding.py
---------------------------------
Retrieve, process, and visualize a sounding profile for a tornado-outbreak
scenario: 27 April 2011 (Super Outbreak).  Demonstrates the full MetPy
sounding workflow from data retrieval to skew-T and hodograph figures.
"""

import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT, Hodograph
from metpy.units import units

warnings.filterwarnings("ignore")

# Synthetic high-CAPE sounding representative of 27 Apr 2011 OKC environment
# Replace with real sounding fetch via step1_siphon_sounding.py for live data.
PRESSURE_HPA = np.array([
    1003, 982, 964, 950, 925, 900, 875, 850, 825, 800, 775, 750,
    725, 700, 675, 650, 625, 600, 575, 550, 525, 500, 475, 450,
    425, 400, 375, 350, 325, 300, 275, 250, 225, 200, 175, 150, 125, 100,
])
TEMPERATURE_C = np.array([
    28.0, 25.5, 23.0, 21.8, 19.5, 17.2, 14.8, 12.5, 10.2, 7.8, 5.4, 3.2,
    1.0, -1.8, -4.5, -7.0, -9.8, -12.5, -15.0, -18.0, -21.5, -25.0, -28.5, -32.0,
    -36.0, -40.0, -44.5, -49.0, -53.5, -58.0, -61.0, -55.0, -51.0, -50.0, -53.0, -56.0, -60.0, -66.0,
])
DEWPOINT_C = np.array([
    22.0, 21.5, 21.0, 20.8, 19.0, 17.0, 14.5, 12.0, 9.0, 5.0, 1.5, -1.0,
    -3.5, -7.0, -11.0, -15.0, -19.5, -23.0, -27.0, -31.0, -35.5, -40.0, -44.0, -48.0,
    -52.0, -56.0, -60.0, -65.0, -68.0, -70.0, -72.0, -70.0, -68.0, -66.0, -66.0, -68.0, -70.0, -75.0,
])
U_WIND_MS = np.array([
    -4, -5, -6, -6, -7, -8, -9, -10, -11, -12, -13, -14,
    -14, -15, -16, -17, -18, -20, -22, -24, -26, -28, -30, -32,
    -33, -34, -36, -38, -40, -42, -43, -44, -45, -44, -42, -40, -38, -35,
], dtype=float)
V_WIND_MS = np.array([
    2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 12,
    13, 14, 15, 15, 16, 16, 17, 17, 18, 18, 18, 17,
    16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3,
], dtype=float)
HEIGHT_M = np.array([
    87, 298, 505, 650, 902, 1160, 1424, 1695, 1973, 2260, 2555, 2859,
    3172, 3493, 3826, 4170, 4527, 4895, 5278, 5675, 6086, 6515, 6960, 7425,
    7910, 8415, 8945, 9500, 10082, 10696, 11343, 12028, 12756, 13534, 14368, 15272, 16267, 17476,
])

p = PRESSURE_HPA * units.hPa
t = TEMPERATURE_C * units.degC
td = DEWPOINT_C * units.degC
u = U_WIND_MS * units("m/s")
v = V_WIND_MS * units("m/s")
z = HEIGHT_M * units.m

# ---- Parcel analysis ----
parcel = mpcalc.parcel_profile(p, t[0], td[0])
cape, cin = mpcalc.cape_cin(p, t, td, parcel)
lcl_p, lcl_t = mpcalc.lcl(p[0], t[0], td[0])
lfc_p, _ = mpcalc.lfc(p, t, td)
el_p, _ = mpcalc.el(p, t, td)
srh_1 = mpcalc.storm_relative_helicity(z, u, v, depth=1000 * units.m)
srh_3 = mpcalc.storm_relative_helicity(z, u, v, depth=3000 * units.m)

print(f"CAPE={cape.m:.0f} J/kg | CIN={cin.m:.0f} J/kg")
print(f"LCL={lcl_p.m:.0f} hPa | LFC={lfc_p.m:.0f} hPa | EL={el_p.m:.0f} hPa")
print(f"0-1 km SRH={srh_1[0].m:.0f} m²/s² | 0-3 km SRH={srh_3[0].m:.0f} m²/s²")

# ---- Combined skew-T + hodograph figure ----
fig = plt.figure(figsize=(16, 10))
skew = SkewT(fig, rotation=45, rect=(0.05, 0.05, 0.55, 0.90))

skew.plot_dry_adiabats(alpha=0.25, colors="orange")
skew.plot_moist_adiabats(alpha=0.25, colors="green")
skew.plot_mixing_lines(alpha=0.25, colors="blue")
skew.plot(p, t, "r", linewidth=2, label="Temperature")
skew.plot(p, td, "g", linewidth=2, label="Dewpoint")
skew.plot(p, parcel, "k--", linewidth=1.5, label="Parcel")
skew.shade_cape(p, t, parcel)
skew.shade_cin(p, t, parcel, td)
skew.ax.axhline(lcl_p.m, color="purple", linestyle=":", linewidth=1.2)
skew.ax.axhline(lfc_p.m, color="darkorange", linestyle="-.", linewidth=1.2)
skew.ax.axhline(el_p.m, color="cyan", linestyle="-.", linewidth=1.2)
skew.plot_barbs(p[::4], u[::4], v[::4])
skew.ax.set_xlim(-40, 50)
skew.ax.set_ylim(1050, 100)
skew.ax.set_title("27 Apr 2011 Super Outbreak – OKC Sounding\nSkew-T Log-P", fontsize=12)
skew.ax.legend(loc="upper right", fontsize=8)

# Hodograph inset
ax_hod = fig.add_axes([0.63, 0.45, 0.33, 0.48])
h = Hodograph(ax_hod, component_range=60)
h.add_grid(increment=10)
intervals = np.array([0, 1000, 3000, 6000]) * units.m
colors = ["purple", "red", "darkorange"]
h.plot_colormapped(u, v, z, intervals=intervals, colors=colors)
ax_hod.set_title("Hodograph (m/s)", fontsize=10)

# Summary text
summary = (
    f"CAPE: {cape.m:.0f} J/kg\n"
    f"CIN:  {cin.m:.0f} J/kg\n"
    f"0-1 SRH: {srh_1[0].m:.0f} m²/s²\n"
    f"0-3 SRH: {srh_3[0].m:.0f} m²/s²"
)
fig.text(0.63, 0.30, summary, fontsize=10,
         bbox=dict(boxstyle="round,pad=0.4", facecolor="wheat", alpha=0.85))

fig.savefig("example1_tornado_sounding.png", dpi=150, bbox_inches="tight")
print("Saved: example1_tornado_sounding.png")
plt.show()

Example 2: Multi-Station CAPE Comparison Along a Dryline

"""
example2_dryline_cape_map.py
-----------------------------
Simulate a spatial survey of CAPE values across a dryline gradient and
plot them on a Cartopy map to illustrate spatial gradients in instability.
Replace the synthetic data section with real Siphon sounding fetches.
"""

import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy.calc as mpcalc
from metpy.units import units

warnings.filterwarnings("ignore")

# Representative station locations along a north-south dryline transect
STATIONS = {
    "AMA": (35.2, -101.7),   # Amarillo, TX (dry side)
    "OKC": (35.4, -97.6),    # Oklahoma City (moist side)
    "LBB": (33.7, -101.8),   # Lubbock, TX (dry)
    "MKC": (39.1, -94.6),    # Kansas City (moist)
    "DDC": (37.8, -99.9),    # Dodge City, KS (transition zone)
    "SGF": (37.2, -93.4),    # Springfield, MO (moist)
    "SHV": (32.4, -93.8),    # Shreveport, LA (very moist)
    "OUN": (35.2, -97.5),    # Norman, OK (moist)
}

# Synthetic CAPE values representing a classic plains dryline scenario
# West of dryline (~101°W): low CAPE; east of dryline: high CAPE
CAPE_VALUES = {
    "AMA": 150,
    "OKC": 2800,
    "LBB": 200,
    "MKC": 2200,
    "DDC": 900,
    "SGF": 3100,
    "SHV": 3500,
    "OUN": 2950,
}
CIN_VALUES = {k: np.random.default_rng(i).integers(10, 120)
              for i, k in enumerate(STATIONS)}


def plot_dryline_cape(
    stations: dict,
    cape_values: dict,
    cin_values: dict,
    output_path: str = "example2_dryline_cape.png",
) -> plt.Figure:
    """
    Plot CAPE values at sounding stations on a plains regional map,
    with circle size proportional to CAPE and color indicating CIN.
    """
    proj = ccrs.LambertConformal(central_longitude=-98, central_latitude=37)
    fig, ax = plt.subplots(figsize=(12, 9), subplot_kw={"projection": proj})
    ax.set_extent([-108, -88, 28, 46], crs=ccrs.PlateCarree())

    ax.add_feature(cfeature.LAND, facecolor="#f0f0e8")
    ax.add_feature(cfeature.OCEAN, facecolor="#cce5ff")
    ax.add_feature(cfeature.COASTLINE, linewidth=0.6)
    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
    ax.add_feature(cfeature.STATES, linewidth=0.8, edgecolor="gray")

    # Approximate dryline position
    dryline_lons = [-101.5, -101.0, -100.5, -100.8, -101.2]
    dryline_lats = [28, 32, 35, 38, 43]
    ax.plot(dryline_lons, dryline_lats,
            transform=ccrs.PlateCarree(),
            color="brown", linewidth=2.5, linestyle="--",
            label="Approximate Dryline", zorder=5)

    # Scatter stations
    lons = [v[1] for v in stations.values()]
    lats = [v[0] for v in stations.values()]
    capes = [cape_values[k] for k in stations]
    cins = [cin_values[k] for k in stations]

    sc = ax.scatter(
        lons, lats,
        s=[c / 15 for c in capes],
        c=cins, cmap="Reds_r", vmin=0, vmax=150,
        transform=ccrs.PlateCarree(),
        zorder=6, edgecolors="black", linewidths=0.7,
        label="Stations (size ∝ CAPE)",
    )
    plt.colorbar(sc, ax=ax, orientation="vertical", label="CIN (J/kg)", shrink=0.6)

    for name, (lat, lon) in stations.items():
        cape = cape_values[name]
        ax.annotate(
            f"{name}\n{cape} J/kg",
            xy=(lon, lat),
            xytext=(5, 5), textcoords="offset points",
            fontsize=7, transform=ccrs.PlateCarree(),
            zorder=7,
        )

    ax.set_title(
        "Synthetic Dryline CAPE Survey — Southern Great Plains\n"
        "Circle size ∝ CAPE; color = CIN (J/kg)",
        fontsize=12,
    )
    ax.legend(loc="lower right", fontsize=9)
    gl = ax.gridlines(draw_labels=True, linewidth=0.3, color="gray", alpha=0.5)
    gl.top_labels = False
    gl.right_labels = False

    fig.tight_layout()
    fig.savefig(output_path, dpi=150, bbox_inches="tight")
    print(f"Saved: {output_path}")
    return fig


if __name__ == "__main__":
    fig = plot_dryline_cape(STATIONS, CAPE_VALUES, CIN_VALUES)
    plt.show()