""" ATMS 411/401 Homework - Single-Site High-Resolution Sounding Analysis - Fetches NCEI National Weather Service Raw Upper Air Data (via UWyo BUFR API). - Computes: Virtual CAPE, 50-mb ML Parcel, IVT, LFC, and Equilibrium Level (EL). - Interpolates Y-axis pressure ticks to display corresponding physical elevation in meters. - Labels Water Vapor Mixing Ratio lines gracefully without cluttering the sounding. - Implements strict monotonic filtering to eliminate balloon pendulum/telemetry noise. - Generates 4 independent, beautiful single-site PNG files. """ import warnings warnings.filterwarnings("ignore") # Suppress MetPy NaN and integration warnings import urllib.request import io import pandas as pd import numpy as np from datetime import datetime import matplotlib.pyplot as plt import metpy.calc as mpcalc from metpy.plots import SkewT from metpy.units import units def get_user_inputs(): print("--- NCEI High-Resolution Raw Upper Air Fetcher ---") try: year = int(input("Year (YYYY) [e.g. 2025]: ")) month = int(input("Month (MM) [e.g. 9]: ")) day = int(input("Day (DD) [e.g. 16]: ")) hour = int(input("Hour (00 or 12): ")) site_id = input("Station Number (WMO ID, e.g., 72489 for Reno): ").strip() if not site_id: site_id = "72489" return datetime(year, month, day, hour), site_id except ValueError: print("\nInvalid input. Defaulting to Reno (72489) at 2025-09-16 12Z.") return datetime(2025, 9, 16, 12), "72489" def fetch_ncei_raw_bufr(dt, site_id): """Fetches high-resolution NWS Raw Upper Air Data via the modern WSGI API endpoint.""" url = f"https://weather.uwyo.edu/wsgi/sounding?datetime={dt.strftime('%Y-%m-%d')}%20{dt.strftime('%H')}:00:00&id={site_id}&src=BUFR&type=TEXT:CSV" print(f"\n -> Fetching NCEI Raw Upper Air Data: {url}") req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'}) try: with urllib.request.urlopen(req, timeout=15) as response: csv_data = response.read().decode('utf-8') except Exception as e: raise RuntimeError(f"Connection failed: {e}") if "error" in csv_data.lower() or len(csv_data.strip().split('\n')) < 5: raise RuntimeError("Data missing or server returned an error for this station/time.") df = pd.read_csv(io.StringIO(csv_data)) speed_col = [c for c in df.columns if 'speed' in c.lower()][0] is_knots = 'knot' in speed_col.lower() or 'kts' in speed_col.lower() col_map = {} for col in df.columns: c = col.lower() if 'pressure' in c and 'hpa' in c: col_map[col] = 'pressure' elif 'height' in c and 'm' in c: col_map[col] = 'height' elif 'dew' in c and 'c' in c: col_map[col] = 'dewpoint' elif 'temperature' in c and 'c' in c and 'dew' not in c and 'ice' not in c: col_map[col] = 'temperature' elif 'direction' in c: col_map[col] = 'direction' elif 'speed' in c: col_map[col] = 'speed' df = df.rename(columns=col_map) req_cols = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed'] for req_col in req_cols: if req_col not in df.columns: raise RuntimeError(f"Missing required column '{req_col}' in dataset.") direction = df['direction'].values * units.degrees speed = df['speed'].values * (units('knots') if is_knots else units('m/s')) u, v = mpcalc.wind_components(speed, direction) df['u_wind'] = u.to('knots').m df['v_wind'] = v.to('knots').m print(" [Success: High-resolution Raw BUFR data retrieved]") return df def get_clean_ascent(df, req_cols, check_height=True): """Filters out balloon swinging and sensor drops to enforce a perfectly smooth ascent profile.""" d = df.dropna(subset=req_cols).copy() d = d.sort_values('pressure', ascending=False).reset_index(drop=True) if len(d) == 0: return d valid_p = np.concatenate(([True], np.diff(d['pressure'].values) < 0)) d = d[valid_p].reset_index(drop=True) if check_height and 'height' in d.columns: h_max = -np.inf keep = [] for h in d['height'].values: if h > h_max: keep.append(True) h_max = h else: keep.append(False) d = d[keep] return d def compute_metrics(df): """Calculates metrics using independent data pipelines to preserve maximum resolution.""" epsilon = 0.62198 Rd = 287.058 * units('J / (kg * K)') # ========================================================================= # PIPELINE 1: Full-Resolution Thermodynamics (Skew-T, CAPE, LFC, EL) # ========================================================================= df_th = get_clean_ascent(df, ['pressure', 'temperature', 'dewpoint'], check_height=False) if len(df_th) > 1: p_th = df_th['pressure'].values * units.hPa T_th = df_th['temperature'].values * units.degC Td_th = df_th['dewpoint'].values * units.degC pw = mpcalc.precipitable_water(p_th, Td_th).to('mm') try: p_ml, T_ml, Td_ml = mpcalc.mixed_parcel(p_th, T_th, Td_th, depth=50 * units.hPa) parcel_prof = mpcalc.parcel_profile(p_th, T_ml, Td_ml).to('degC') cape_ml, cin_ml = mpcalc.cape_cin(p_th, T_th, Td_th, parcel_prof) cape_val = max(0.0, cape_ml.m) cin_val = cin_ml.m # Exact Surface-Based LFC and EL lfc_sb_val = mpcalc.lfc(p_th, T_th, Td_th) sb_lfc_p = np.nanmean(lfc_sb_val[0].m) if lfc_sb_val[0] is not None else np.nan el_sb_val = mpcalc.el(p_th, T_th, Td_th) sb_el_p = np.nanmean(el_sb_val[0].m) if el_sb_val[0] is not None else np.nan # Exact Mixed-Layer LFC and EL (Override surface point to force ML LCL calculation) T_env_ml = np.copy(T_th.m) Td_env_ml = np.copy(Td_th.m) T_env_ml[0] = T_ml.m Td_env_ml[0] = Td_ml.m lfc_ml_val = mpcalc.lfc(p_th, T_env_ml * units.degC, Td_env_ml * units.degC, parcel_temperature_profile=parcel_prof) ml_lfc_p = np.nanmean(lfc_ml_val[0].m) if lfc_ml_val[0] is not None else np.nan el_ml_val = mpcalc.el(p_th, T_env_ml * units.degC, Td_env_ml * units.degC, parcel_temperature_profile=parcel_prof) ml_el_p = np.nanmean(el_ml_val[0].m) if el_ml_val[0] is not None else np.nan except Exception: p_ml, T_ml, Td_ml = [np.nan * units.hPa, np.nan * units.degC, np.nan * units.degC] parcel_prof = np.full(len(p_th), np.nan) * units.degC cape_val, cin_val = 0.0, 0.0 sb_lfc_p, sb_el_p, ml_lfc_p, ml_el_p = np.nan, np.nan, np.nan, np.nan else: p_th, T_th, Td_th = [np.array([]) * u for u in (units.hPa, units.degC, units.degC)] pw, cape_val, cin_val = 0.0 * units.mm, 0.0, 0.0 p_ml, T_ml, Td_ml = [np.nan * units.hPa, np.nan * units.degC, np.nan * units.degC] parcel_prof = np.array([]) * units.degC sb_lfc_p, sb_el_p, ml_lfc_p, ml_el_p = np.nan, np.nan, np.nan, np.nan # ========================================================================= # PIPELINE 2: Density Profile & LFC Height Interpolation # ========================================================================= df_z = get_clean_ascent(df, ['pressure', 'height', 'temperature', 'dewpoint'], check_height=True) if len(df_z) > 1: p_z = df_z['pressure'].values * units.hPa z_z = df_z['height'].values * units.meters T_z = df_z['temperature'].values * units.degC Td_z = df_z['dewpoint'].values * units.degC e_vap = mpcalc.saturation_vapor_pressure(Td_z).to('Pa') p_dry = p_z.to('Pa') - e_vap rho_dry = (p_dry / (Rd * T_z.to('kelvin'))).to('kg/m**3') w = epsilon * e_vap / p_dry Tv = mpcalc.virtual_temperature(T_z, w) rho_air = (p_z.to('Pa') / (Rd * Tv.to('kelvin'))).to('kg/m**3') rho_wv = (rho_dry * w).to('kg/m**3') theta = mpcalc.potential_temperature(p_z, T_z) theta_e = mpcalc.equivalent_potential_temperature(p_z, T_z, Td_z) # Interpolate LFC and EL heights using logarithm of pressure (Hypsometric equation) p_log_asc = np.log(p_z.m)[::-1] z_asc = z_z.m[::-1] sb_lfc_z = np.interp(np.log(sb_lfc_p), p_log_asc, z_asc) if not np.isnan(sb_lfc_p) else np.nan sb_el_z = np.interp(np.log(sb_el_p), p_log_asc, z_asc) if not np.isnan(sb_el_p) else np.nan ml_lfc_z = np.interp(np.log(ml_lfc_p), p_log_asc, z_asc) if not np.isnan(ml_lfc_p) else np.nan ml_el_z = np.interp(np.log(ml_el_p), p_log_asc, z_asc) if not np.isnan(ml_el_p) else np.nan else: p_z, z_z, T_z, Td_z = [np.array([]) * u for u in (units.hPa, units.meters, units.degC, units.degC)] rho_air, rho_wv = np.array([]) * units('kg/m**3'), np.array([]) * units('kg/m**3') theta, theta_e = np.array([]) * units.kelvin, np.array([]) * units.kelvin sb_lfc_z, sb_el_z, ml_lfc_z, ml_el_z = np.nan, np.nan, np.nan, np.nan # ========================================================================= # PIPELINE 3: IVT Integration # ========================================================================= df_ivt = get_clean_ascent(df, ['pressure', 'height', 'temperature', 'dewpoint', 'u_wind', 'v_wind'], check_height=True) if len(df_ivt) > 1: p_ivt = df_ivt['pressure'].values * units.hPa z_ivt = df_ivt['height'].values * units.meters T_ivt = df_ivt['temperature'].values * units.degC Td_ivt = df_ivt['dewpoint'].values * units.degC u_ivt = df_ivt['u_wind'].values * units.knots v_ivt = df_ivt['v_wind'].values * units.knots e_vap_ivt = mpcalc.saturation_vapor_pressure(Td_ivt).to('Pa') p_dry_ivt = p_ivt.to('Pa') - e_vap_ivt rho_dry_ivt = (p_dry_ivt / (Rd * T_ivt.to('kelvin'))).to('kg/m**3') w_ivt = epsilon * e_vap_ivt / p_dry_ivt rho_wv_ivt = (rho_dry_ivt * w_ivt).to('kg/m**3') ivt_u_val = np.trapz(rho_wv_ivt.m * u_ivt.m_as('m/s'), z_ivt.m) ivt_v_val = np.trapz(rho_wv_ivt.m * v_ivt.m_as('m/s'), z_ivt.m) ivt_u = ivt_u_val * units('kg / (m * s)') ivt_v = ivt_v_val * units('kg / (m * s)') ivt_mag = np.sqrt(ivt_u**2 + ivt_v**2) ivt_dir = mpcalc.wind_direction(ivt_u.m * units('m/s'), ivt_v.m * units('m/s')) else: ivt_mag, ivt_u, ivt_v, ivt_dir = [0 * units('kg / (m * s)')] * 3 + [0 * units.degrees] # ========================================================================= # PIPELINE 4: Wind Barbs # ========================================================================= df_wind = get_clean_ascent(df, ['pressure', 'u_wind', 'v_wind'], check_height=False) if len(df_wind) > 1: p_wind = df_wind['pressure'].values * units.hPa u_wind = df_wind['u_wind'].values * units.knots v_wind = df_wind['v_wind'].values * units.knots else: p_wind, u_wind, v_wind = [np.array([]) * u for u in (units.hPa, units.knots, units.knots)] return { 'pw': pw, 'ivt_u': ivt_u, 'ivt_v': ivt_v, 'ivt_mag': ivt_mag, 'ivt_dir': ivt_dir, 'p_ml': p_ml, 'T_ml': T_ml, 'Td_ml': Td_ml, 'cape': cape_val, 'cin': cin_val, 'sb_lfc_p': sb_lfc_p, 'sb_lfc_z': sb_lfc_z, 'sb_el_p': sb_el_p, 'sb_el_z': sb_el_z, 'ml_lfc_p': ml_lfc_p, 'ml_lfc_z': ml_lfc_z, 'ml_el_p': ml_el_p, 'ml_el_z': ml_el_z, 'p_th': p_th, 'T_th': T_th, 'Td_th': Td_th, 'parcel_prof': parcel_prof, 'p_z': p_z, 'z_z': z_z, 'T_z': T_z, 'rho_air': rho_air, 'rho_wv': rho_wv, 'theta': theta, 'theta_e': theta_e, 'p_wind': p_wind, 'u_wind': u_wind, 'v_wind': v_wind } def plot_skewt(metrics, station_name, dt): if len(metrics['p_th']) == 0: return fig = plt.figure(figsize=(10, 10)) skew = SkewT(fig, rotation=45) skew.plot(metrics['p_th'], metrics['T_th'], 'r', label='Temp', linewidth=1.5) skew.plot(metrics['p_th'], metrics['Td_th'], 'g', label='Dewpoint', linewidth=1.5) skew.plot(metrics['p_th'], metrics['parcel_prof'], 'k:', linewidth=2, label='50mb ML Parcel') if len(metrics['p_wind']) > 0: step = max(1, len(metrics['p_wind']) // 45) skew.plot_barbs(metrics['p_wind'][::step], metrics['u_wind'][::step], metrics['v_wind'][::step]) skew.plot_dry_adiabats() skew.plot_moist_adiabats() skew.plot_mixing_lines() # --------------------------------------------------------- # Add values to the water vapor mixing ratio lines (green dashes) # MetPy defaults: [0.0004, 0.001, 0.002, 0.004, 0.007, 0.01, 0.016, 0.024, 0.032] w_vals = np.array([0.0004, 0.001, 0.002, 0.004, 0.007, 0.01, 0.016, 0.024, 0.032]) * units('dimensionless') w_labels = [0.4, 1, 2, 4, 7, 10, 16, 24, 32] # Place labels safely at 925 hPa and 625 hPa so they don't overlap the main gridlines for p_label in [925 * units.hPa, 625 * units.hPa]: e_label = mpcalc.vapor_pressure(p_label, w_vals) t_label = mpcalc.dewpoint(e_label) for w_num, t in zip(w_labels, t_label): if -38 <= t.m <= 43: # Keep comfortably inside the X-axis bounds skew.ax.text(t.m, p_label.m, f'{w_num:g}', fontsize=8, color='forestgreen', weight='bold', alpha=0.8, ha='center', va='center', bbox=dict(facecolor='white', edgecolor='none', alpha=0.6, pad=1.2)) # --------------------------------------------------------- # Configure limits (-40C to 45C horizontally) p_max_th = metrics['p_th'].m.max() bottom_limit = max(1000, np.ceil(p_max_th / 10) * 10 + 10) skew.ax.set_ylim(bottom_limit, 100) skew.ax.set_xlim(-40, 45) # Inject physical elevation into the Y-axis pressure labels yticks = np.arange(100, 1051, 100) valid_ticks = [t for t in yticks if bottom_limit >= t >= 100] if len(metrics['p_z']) > 0: p_z_max = metrics['p_z'].m.max() new_labels = [] for t in valid_ticks: if t > p_z_max + 2: # Ignore subterranean pressures (e.g. 1000hPa in Reno) new_labels.append(f"{t}") else: z_interp = np.interp(np.log(t), np.log(metrics['p_z'].m)[::-1], metrics['z_z'].m[::-1]) new_labels.append(f"{t}\n({int(z_interp)}m)") skew.ax.set_yticks(valid_ticks) skew.ax.set_yticklabels(new_labels) skew.ax.set_xlabel('Temperature (\xb0C)') skew.ax.set_ylabel('Pressure (hPa) & Elevation (m)') skew.ax.legend(loc='upper right') skew.ax.set_title(f"Skew-T Log-P (Raw High-Res): {station_name}\n{dt.strftime('%Y-%m-%d %H:%M Z')}", fontsize=14) filename = f"SkewT_{station_name}_{dt.strftime('%Y%m%d_%H')}Z.png" plt.savefig(filename, dpi=150, bbox_inches='tight') plt.close() def plot_temp_pressure_height(metrics, station_name, dt): if len(metrics['p_z']) == 0: return fig, ax1 = plt.subplots(figsize=(8, 10)) ax2 = ax1.twiny() ax1.set_xlabel('Pressure (hPa)', fontsize=12, color='darkgreen') ax1.tick_params(axis='x', labelcolor='darkgreen') ax1.set_ylabel('Height (km)', fontsize=12) ax2.set_xlabel('Temperature (\xb0C)', fontsize=12, color='darkred') ax2.tick_params(axis='x', labelcolor='darkred') ax1.plot(metrics['p_z'].m, metrics['z_z'].to('km').m, 'g-', label='Pressure', linewidth=2) ax2.plot(metrics['T_z'].m, metrics['z_z'].to('km').m, 'r-', label='Temperature', linewidth=2) lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right') plt.title(f"Temperature and Pressure vs Height\nSite: {station_name} | {dt.strftime('%Y-%m-%d %H:%M Z')}", y=1.08, fontsize=14) filename = f"Temp_Pressure_vs_Height_{station_name}_{dt.strftime('%Y%m%d_%H')}Z.png" plt.savefig(filename, dpi=150, bbox_inches='tight') plt.close() def plot_densities_height(metrics, station_name, dt): if len(metrics['rho_air']) == 0: return fig, ax1 = plt.subplots(figsize=(8, 10)) ax2 = ax1.twiny() ax1.set_xlabel('Water Vapor Density (kg/m$^3$)', fontsize=12, color='blue') ax1.tick_params(axis='x', labelcolor='blue') ax1.set_ylabel('Height (km)', fontsize=12) ax2.set_xlabel('Air Density using Virtual Temp (kg/m$^3$)', fontsize=12, color='orange') ax2.tick_params(axis='x', labelcolor='orange') ax1.plot(metrics['rho_wv'].m, metrics['z_z'].to('km').m, 'b-', label='$\\rho_{wv}$', linewidth=2) ax2.plot(metrics['rho_air'].m, metrics['z_z'].to('km').m, color='orange', linestyle='-', label='$\\rho_{air}$', linewidth=2) lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right') plt.title(f"Densities vs Height\nSite: {station_name} | {dt.strftime('%Y-%m-%d %H:%M Z')}", y=1.08, fontsize=14) filename = f"Densities_vs_Height_{station_name}_{dt.strftime('%Y%m%d_%H')}Z.png" plt.savefig(filename, dpi=150, bbox_inches='tight') plt.close() def plot_theta_profile(metrics, station_name, dt): """Plot Potential Temp and Equivalent Potential Temp strictly up to 12 km.""" mask = metrics['z_z'].m <= 12000 if not np.any(mask): return z = metrics['z_z'][mask].to('km') theta = metrics['theta'][mask] theta_e = metrics['theta_e'][mask] fig, ax = plt.subplots(figsize=(8, 10)) ax.plot(theta.m, z.m, 'b-', linewidth=2, label=r'Potential Temp ($\theta$)') ax.plot(theta_e.m, z.m, 'r-', linewidth=2, label=r'Equivalent Potential Temp ($\theta_e$)') ax.set_ylabel('Height (km)', fontsize=12) ax.set_xlabel('Temperature (Kelvin)', fontsize=12) ax.set_title(f"Thermodynamic Profile (Surface to 12 km)\nSite: {station_name} | {dt.strftime('%Y-%m-%d %H:%M Z')}", fontsize=14) ax.grid(True, linestyle='--', alpha=0.7) ax.legend(loc='upper right') ax.set_ylim(z.m.min(), 12) filename = f"Theta_Profile_{station_name}_{dt.strftime('%Y%m%d_%H')}Z.png" plt.savefig(filename, dpi=150, bbox_inches='tight') plt.close() def main(): dt, site_id = get_user_inputs() try: df = fetch_ncei_raw_bufr(dt, site_id) except RuntimeError as e: print(f"\nTERMINAL ERROR: {e}") return m = compute_metrics(df) print("\n" + "="*50) print(f"METRICS FOR SITE {site_id} | {dt.strftime('%Y-%m-%d %H:%M Z')}") print("="*50) print(f" Precipitable Water : {m['pw'].m:.2f} mm") print(f" IVT u-component : {m['ivt_u'].m:.2f} kg/m/s") print(f" IVT v-component : {m['ivt_v'].m:.2f} kg/m/s") print(f" IVT Magnitude : {m['ivt_mag'].m:.2f} kg/m/s") print(f" IVT Direction : {m['ivt_dir'].m:.1f} deg (Origin)") print(f" 50mb ML Parcel Start : T={m['T_ml']:.1f}, Td={m['Td_ml']:.1f} (at {m['p_ml']:.0f} hPa)") print(f" ML CAPE (Virtual) : {m['cape']:.1f} J/kg") print(f" ML CIN : {m['cin']:.1f} J/kg") print("-" * 50) print(f" Surface-Based LFC : {m['sb_lfc_z']:.0f} m (at {m['sb_lfc_p']:.0f} hPa)" if not np.isnan(m['sb_lfc_z']) else " Surface-Based LFC : None") print(f" Surface-Based EL : {m['sb_el_z']:.0f} m (at {m['sb_el_p']:.0f} hPa)" if not np.isnan(m['sb_el_z']) else " Surface-Based EL : None") print(f" Mixed-Layer LFC : {m['ml_lfc_z']:.0f} m (at {m['ml_lfc_p']:.0f} hPa)" if not np.isnan(m['ml_lfc_z']) else " Mixed-Layer LFC : None") print(f" Mixed-Layer EL : {m['ml_el_z']:.0f} m (at {m['ml_el_p']:.0f} hPa)" if not np.isnan(m['ml_el_z']) else " Mixed-Layer EL : None") print("="*50) print("\nGenerating Skew-T diagram with Wind Barbs and Elevation Ticks...") plot_skewt(m, site_id, dt) print("Generating Temperature & Pressure vs Height...") plot_temp_pressure_height(m, site_id, dt) print("Generating Densities vs Height...") plot_densities_height(m, site_id, dt) print("Generating Theta & Theta_e vs Height (to 12km)...") plot_theta_profile(m, site_id, dt) print("\nAll 4 plots saved successfully to your current working directory as PNG files.") if __name__ == "__main__": main()