#!/usr/bin/env python """ HRRR sounding via Herbie: - User chooses lat/lon and UTC date/hour - Gets surface elevation automatically (USGS EPQS) or asks user - Downloads HRRR pressure-level fields (prs product) - Computes Td, mixing ratio, theta, theta-e, theta-v - Saves UW-style sounding file and Skew-T PNG - Defines parcel at the actual surface pressure (site elevation) and lifts it to 100 hPa - Computes CAPE, CIN, LCL, LFC, EL - Only outputs/plots data from site elevation upward (geopotential height >= surface elevation) and down to 100 hPa """ import os import sys import requests import numpy as np import pandas as pd import matplotlib.pyplot as plt from herbie import Herbie import xarray as xr import metpy.calc as mpcalc from metpy.plots import SkewT from metpy.units import units def get_surface_elevation(lat, lon): """Try USGS EPQS service; fallback to user input.""" try: url = f"https://epqs.nationalmap.gov/v1/json?x={lon}&y={lat}&units=Meters&wkid=4326" r = requests.get(url, timeout=10) r.raise_for_status() elev = r.json()["value"]["elevation"] print(f"Retrieved surface elevation: {elev:.1f} m") return float(elev) except Exception as e: print("Could not retrieve elevation automatically:", e) elev = float(input("Enter surface elevation in meters: ").strip()) return elev def main(): # ----------------------------- # User input # ----------------------------- print("Enter location and UTC time for HRRR sounding:") try: lat = float(input("Latitude (e.g., 39.0968): ").strip()) lon = float(input("Longitude (e.g., -120.0324): ").strip()) date_str = input("Date (YYYY-MM-DD, e.g., 2025-06-21): ").strip() hour_str = input("Hour UTC (00–23, e.g., 19): ").strip() timestamp = f"{date_str} {int(hour_str):02d}:00" dt = pd.to_datetime(timestamp) except Exception as e: print("Input error:", e) sys.exit(1) # Surface elevation surface_elev = get_surface_elevation(lat, lon) # ----------------------------- # Herbie setup # ----------------------------- save_dir = os.path.join(os.path.expanduser("~"), "data", "hrrr") os.makedirs(save_dir, exist_ok=True) H = Herbie(timestamp, model="hrrr", product="prs", fxx=0, save_dir=save_dir) # Helper to read an isobaric field via cfgrib filter_by_keys def get_isobaric(short_name, rename=None): ds = H.xarray( backend_kwargs={ "filter_by_keys": {"shortName": short_name, "typeOfLevel": "isobaricInhPa"} } ) if isinstance(ds, list): ds = xr.merge(ds) if rename: var_names = list(ds.data_vars) if len(var_names) == 1 and var_names[0] != rename: ds = ds.rename({var_names[0]: rename}) elif rename in var_names: pass else: ds = ds.rename({var_names[0]: rename}) return ds # Fetch essential fields ds_T = get_isobaric("t", "T") # Temperature (K) ds_u = get_isobaric("u", "u") # u wind (m/s) ds_v = get_isobaric("v", "v") # v wind (m/s) if ds_T is None or ds_u is None or ds_v is None: raise RuntimeError("Missing essential fields (T, u, v). Try a different cycle/time or verify installation.") # Optional humidity and height ds_RH = get_isobaric("r", "RH") # Relative humidity (%) ds_gh = get_isobaric("gh", "gh") # Geopotential height (m) ds_q = None if ds_RH is not None else get_isobaric("q", "q") # Specific humidity (kg/kg) fallback # Merge datasets ds_all = ds_T for d in [ds_u, ds_v, ds_RH, ds_gh, ds_q]: if d is not None: ds_all = xr.merge([ds_all, d], compat="override", join="inner") # ----------------------------- # Nearest gridpoint on 2D lat/lon # ----------------------------- # Find names lat_name = "latitude" if "latitude" in ds_all else ("lat" if "lat" in ds_all else None) lon_name = "longitude" if "longitude" in ds_all else ("lon" if "lon" in ds_all else None) y_name = "y" if "y" in ds_all.dims else ("Y" if "Y" in ds_all.dims else None) x_name = "x" if "x" in ds_all.dims else ("X" if "X" in ds_all.dims else None) if any(n is None for n in [lat_name, lon_name, y_name, x_name]): raise KeyError("Could not locate 2D latitude/longitude or grid dims (y/x) in dataset.") lats = ds_all[lat_name].values lons = ds_all[lon_name].values target_lon = lon % 360 if np.nanmax(lons) > 180 else lon dist2 = (lats - lat) ** 2 + (lons - target_lon) ** 2 iy, ix = np.unravel_index(np.nanargmin(dist2), dist2.shape) pt = ds_all.isel(**{y_name: iy, x_name: ix}) # ----------------------------- # Extract base vertical profile # ----------------------------- p_name = "isobaricInhPa" if p_name not in pt.coords: raise RuntimeError(f"Pressure coordinate '{p_name}' not found in dataset.") p = pt[p_name].values * units.hPa T = pt["T"].values * units.kelvin u = pt["u"].values * units("m/s") v = pt["v"].values * units("m/s") # Humidity: RH preferred; fallback to q -> RH if "RH" in pt: RH = (pt["RH"].values / 100.0) * units.dimensionless elif "q" in pt: q = pt["q"].values * units("kg/kg") w_from_q = (q / (1.0 * units.dimensionless - q)).to("kg/kg") p_pa = p.to("Pa") e = (w_from_q * p_pa) / (0.622 * units.dimensionless + w_from_q) es = mpcalc.saturation_vapor_pressure(T).to("Pa") RH = (e / es).to("dimensionless") else: raise RuntimeError("Neither RH nor q available; cannot compute humidity.") # Constrain RH and compute derived quantities RH = np.clip(RH.magnitude, 0.0, 1.0) * units.dimensionless Td = mpcalc.dewpoint_from_relative_humidity(T, RH) w = mpcalc.mixing_ratio_from_relative_humidity(p, T, RH).to("g/kg") # Geopotential height: prefer model; else hypsometric if "gh" in pt: z = pt["gh"].values # meters else: Tv_for_z = T * (1.0 + 0.61 * w.to("kg/kg")) z = np.zeros_like(p.magnitude) Rd = 287.0 g0 = 9.80665 for i in range(1, len(p)): Tv_mean = 0.5 * (Tv_for_z[i - 1] + Tv_for_z[i]) z[i] = z[i - 1] + (Rd * Tv_mean / g0) * np.log(p[i - 1].magnitude / p[i].magnitude) # ----------------------------- # Clip profile: actual elevation upward, down to 100 hPa # ----------------------------- mask_clip = (z >= surface_elev) & (p.magnitude >= 100.0) if not np.any(mask_clip): raise RuntimeError("No HRRR pressure levels found above the site elevation with pressure >= 100 hPa.") p = p[mask_clip] T = T[mask_clip] Td = Td[mask_clip] u = u[mask_clip] v = v[mask_clip] RH = RH[mask_clip] w = w[mask_clip] z = z[mask_clip] # ----------------------------- # Thermodynamic potentials and wind diagnostics # ----------------------------- theta = mpcalc.potential_temperature(p, T) theta_e = mpcalc.equivalent_potential_temperature(p, T, Td) Tv = T * (1.0 + 0.61 * w.to("kg/kg")) theta_v = mpcalc.potential_temperature(p, Tv) drct = (np.degrees(np.arctan2(-u.magnitude, -v.magnitude)) + 360.0) % 360.0 sped = np.sqrt(u**2 + v**2).to("m/s") # Build dataframe, sorted surface -> top df = pd.DataFrame({ "PRES": p.magnitude, "HGHT": z, "TEMP": T.to("degC").magnitude, "DWPT": Td.to("degC").magnitude, "RELH": (RH.magnitude * 100.0), "MIXR": w.magnitude, # g/kg "DRCT": drct, "SPED": sped.magnitude, # m/s "THTA": theta.magnitude, "THTE": theta_e.magnitude, "THTV": theta_v.magnitude, }).sort_values(by="PRES", ascending=False).reset_index(drop=True) # ----------------------------- # Parcel: start at actual surface pressure (site elevation) # ----------------------------- p_prof = (df["PRES"].values * units.hPa) T_env_prof = (df["TEMP"].values * units.degC) Td_env_prof = (df["DWPT"].values * units.degC) T_sfc = df["TEMP"].iloc[0] * units.degC Td_sfc = df["DWPT"].iloc[0] * units.degC parcel_prof = mpcalc.parcel_profile(p_prof, T_sfc, Td_sfc) # CAPE/CIN and parcel levels cape, cin = mpcalc.cape_cin(p_prof, T_env_prof, Td_env_prof, parcel_prof) p_lcl, T_lcl = mpcalc.lcl(p_prof[0], T_sfc, Td_sfc) lfc_p, _ = mpcalc.lfc(p_prof, T_env_prof, Td_env_prof, parcel_prof) el_p, _ = mpcalc.el(p_prof, T_env_prof, Td_env_prof, parcel_prof) # Map pressures to heights via interpolation on clipped profile pres_full = df["PRES"].values hght_full = df["HGHT"].values pres_asc = pres_full[::-1] hght_asc = hght_full[::-1] def p_to_z(p_hpa): try: val = p_hpa.magnitude except Exception: return np.nan * units.meter if np.isnan(val): return np.nan * units.meter return float(np.interp(val, pres_asc, hght_asc)) * units.meter lcl_z = p_to_z(p_lcl) lfc_z = p_to_z(lfc_p) el_z = p_to_z(el_p) # Helper to format quantities safely def fmt_qty(q, fmt=".1f"): try: val = q.magnitude if np.isnan(val): return "NaN" return format(val, fmt) except Exception: return "NaN" # ----------------------------- # Save UW-format text file (only clipped levels) # ----------------------------- txt_name = f"HRRR_sounding_{dt.strftime('%Y%m%d_%H')}Z_{lat:.3f}_{lon:.3f}.txt" with open(txt_name, "w", encoding="utf-8") as f: f.write(f"Observations,for,Model,HRRR,at,{dt.strftime('%H')},UTC,{dt.day},{dt.strftime('%b')},{dt.year},,,\n") f.write("\n") f.write(f"Model derived sounding at lat: {lat:.4f}, lon: {lon:.4f}\n") f.write(f"Surface elevation (m): {surface_elev:.1f}\n") f.write("\n") f.write("----------------------------------------------------------------------------- \n") f.write(",PRES,HGHT,TEMP,DWPT,RELH,MIXR,DRCT,SPED,THTA,THTE,THTV\n") f.write(",hPa,m,C,C,%,g/kg,deg,m/s,K,K,K\n") f.write("----------------------------------------------------------------------------- \n") for _, row in df.iterrows(): f.write( f",{row['PRES']:.1f},{row['HGHT']:.1f},{row['TEMP']:.1f},{row['DWPT']:.1f}," f"{row['RELH']:.0f},{row['MIXR']:.2f},{row['DRCT']:.0f},{row['SPED']:.1f}," f"{row['THTA']:.1f},{row['THTE']:.0f},{row['THTV']:.1f}\n" ) f.write("----------------------------------------------------------------------------- \n") f.write("Parcel diagnostics:\n") f.write(f"- Parcel origin: Actual surface (site elevation)\n") f.write(f"- Surface (hPa,m): {df.loc[0, 'PRES']:.1f}, {df.loc[0, 'HGHT']:.1f}\n") f.write(f"- CAPE (J/kg): {fmt_qty(cape)}\n") f.write(f"- CIN (J/kg): {fmt_qty(cin)}\n") f.write(f"- LCL (hPa,m): {fmt_qty(p_lcl)}, {fmt_qty(lcl_z)}\n") f.write(f"- LFC (hPa,m): {fmt_qty(lfc_p)}, {fmt_qty(lfc_z)}\n") f.write(f"- EL (hPa,m): {fmt_qty(el_p)}, {fmt_qty(el_z)}\n") print(f"UW-style sounding saved: {os.path.abspath(txt_name)}") # ----------------------------- # Skew-T plot (only clipped range), styled # ----------------------------- # Wind barbs from df speed/direction dir_rad = np.deg2rad(df["DRCT"].values) u_barb = (-df["SPED"].values * np.sin(dir_rad)) * units("m/s") v_barb = (-df["SPED"].values * np.cos(dir_rad)) * units("m/s") fig = plt.figure(figsize=(8, 10)) skew = SkewT(fig) # Environmental T and Td: black and thicker skew.plot(df["PRES"].values * units.hPa, df["TEMP"].values * units.degC, color="k", linewidth=2.5, label="Temperature") skew.plot(df["PRES"].values * units.hPa, df["DWPT"].values * units.degC, color="k", linewidth=2.5, linestyle="--", label="Dewpoint") # Parcel trajectory: purple and thicker (from actual surface pressure to 100 hPa) skew.plot(p_prof, parcel_prof, color="purple", linewidth=3.0, label="Parcel") # Wind barbs skew.plot_barbs(df["PRES"].values * units.hPa, u_barb, v_barb) # Thin thermodynamic lines (version-safe) if hasattr(skew, "plot_dry_adiabats"): skew.plot_dry_adiabats(linewidth=0.5) elif hasattr(skew, "add_dry_adiabats"): skew.add_dry_adiabats(linewidth=0.5) if hasattr(skew, "plot_moist_adiabats"): skew.plot_moist_adiabats(linewidth=0.5) elif hasattr(skew, "add_moist_adiabats"): skew.add_moist_adiabats(linewidth=0.5) if hasattr(skew, "plot_mixing_lines"): skew.plot_mixing_lines(linewidth=0.5) elif hasattr(skew, "add_mixing_lines"): skew.add_mixing_lines(linewidth=0.5) # Y-limits: actual surface pressure to 100 hPa skew.ax.set_ylim(df["PRES"].max(), 100) skew.ax.set_xlim(-40, 50) plt.title( f"HRRR Skew-T at {lat:.3f}, {lon:.3f} on {dt.strftime('%Y-%m-%d %H')}Z\n" f"Surface elevation {surface_elev:.0f} m " f"CAPE={fmt_qty(cape, '.0f')} J/kg CIN={fmt_qty(cin, '.0f')} J/kg", fontsize=12 ) plt.legend(loc="best") plt.tight_layout() png_name = f"HRRR_skewt_{dt.strftime('%Y%m%d_%H')}Z_{lat:.3f}_{lon:.3f}.png" plt.savefig(png_name, dpi=300) plt.show() print(f"Skew-T PNG saved: {os.path.abspath(png_name)}") # Console diagnostics print("\nDiagnostics (parcel from actual surface):") print(f"CAPE: {fmt_qty(cape)} J/kg") print(f"CIN: {fmt_qty(cin)} J/kg") print(f"LCL: {fmt_qty(p_lcl)} hPa, {fmt_qty(lcl_z)} m") print(f"LFC: {fmt_qty(lfc_p)} hPa, {fmt_qty(lfc_z)} m") print(f"EL: {fmt_qty(el_p)} hPa, {fmt_qty(el_z)} m") if __name__ == "__main__": main()