"""
ATMS 411/401 Homework 4 - Advanced Sounding Analysis
- Bypasses broken UWyo module with direct HTML fetcher.
- Falls back to NOAA IGRA if UWyo is down (with correct unit handling).
- Computes: Virtual CAPE, 50-mb ML Parcel, and Integrated Vapor Transport (IVT) using dz.
- Enforces strict monotonic pressure sorting to prevent massive negative CAPE bugs.
- Generates 5 strictly requested PNG files (Skew-Ts with Wind Barbs, Temp/Pres, Densities, Delta P).
"""
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
from siphon.simplewebservice.igra2 import IGRAUpperAir
def get_user_datetime():
print("--- Upper Air Sounding Fetcher ---")
print("Enter the date and time for the sounding (e.g., Year: 2025, Month: 9, Day: 16, Hour: 12)")
try:
year = int(input("Year (YYYY): "))
month = int(input("Month (MM): "))
day = int(input("Day (DD): "))
hour = int(input("Hour (00 or 12): "))
return datetime(year, month, day, hour)
except ValueError:
print("\nInvalid input. Defaulting to 2025-09-16 12Z.")
return datetime(2025, 9, 16, 12)
def direct_wyoming_fetch(dt, wy_id, region):
"""Directly scrapes the UWyo HTTP site, bypassing the broken Siphon module."""
year = dt.year
month = f"{dt.month:02d}"
dayhour = f"{dt.day:02d}{dt.hour:02d}"
url = f"http://weather.uwyo.edu/cgi-bin/sounding?region={region}&TYPE=TEXT%3ALIST&YEAR={year}&MONTH={month}&FROM={dayhour}&TO={dayhour}&STNM={wy_id}"
print(f" -> Fetching: {url}")
req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'})
try:
with urllib.request.urlopen(req, timeout=10) as response:
html = response.read().decode('utf-8')
except Exception as e:
raise RuntimeError(f"Connection failed: {e}")
if "Can't get" in html or "Sorry" in html:
raise RuntimeError("Data missing on UWyo.")
try:
raw_text = html.split('')[1].split('')[0] if '' in html else html.split('')[1].split('')[0]
except IndexError:
raise RuntimeError("Data table not found in the HTML.")
data_lines = [line for line in raw_text.split('\n') if line.strip() and not line.startswith('---') and 'PRES' not in line and 'hPa' not in line]
if not data_lines:
raise RuntimeError("Data table was empty.")
df = pd.read_fwf(io.StringIO('\n'.join(data_lines)), widths=[7]*11, header=None,
names=['pressure', 'height', 'temperature', 'dewpoint', 'relh', 'mixr', 'direction', 'speed', 'thta', 'thte', 'thtv'])
for col in df.columns:
df[col] = pd.to_numeric(df[col], errors='coerce')
direction = df['direction'].values * units.degrees
speed = df['speed'].values * units.knots
u, v = mpcalc.wind_components(speed, direction)
df['u_wind'] = u.m
df['v_wind'] = v.m
print(" [Success: Data retrieved from Wyoming]")
return df
def fetch_sounding(dt, wy_id, wy_region, igra_id):
try:
return direct_wyoming_fetch(dt, wy_id, wy_region)
except Exception as e:
print(f" [Wyoming failed: {e}]")
print(f" -> Falling back to NOAA IGRA server (ID: {igra_id})...")
try:
df, _ = IGRAUpperAir.request_data(dt, igra_id)
if 'dewpoint' not in df.columns and 'dewpoint_depression' in df.columns:
df['dewpoint'] = df['temperature'] - df['dewpoint_depression']
spd_col = 'speed' if 'speed' in df.columns else 'wind_speed'
dir_col = 'direction' if 'direction' in df.columns else 'wind_direction'
if spd_col in df.columns and dir_col in df.columns:
direction = df[dir_col].values * units.degrees
speed = df[spd_col].values * units('m/s')
u, v = mpcalc.wind_components(speed, direction)
df['u_wind'] = u.to('knots').m
df['v_wind'] = v.to('knots').m
else:
df['u_wind'], df['v_wind'] = np.nan, np.nan
print(" [Success: Data retrieved from NOAA IGRA]")
return df
except Exception as e_igra:
raise RuntimeError(f"All servers failed: {e_igra}")
def compute_metrics(df):
"""Calculates metrics using independent data pipelines to preserve maximum resolution for each variable."""
epsilon = 0.62198
Rd = 287.058 * units('J / (kg * K)')
# =========================================================================
# PIPELINE 1: Full-Resolution Thermodynamics (Skew-T & CAPE)
# =========================================================================
df_th = df.dropna(subset=['pressure', 'temperature', 'dewpoint']).copy()
df_th = df_th.sort_values('pressure', ascending=False).drop_duplicates('pressure')
if len(df_th) > 1:
valid_p_idx = np.concatenate(([True], np.diff(df_th['pressure'].values) < 0))
df_th = df_th[valid_p_idx]
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
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
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
# =========================================================================
# PIPELINE 2: Density Profile (Requires Height)
# =========================================================================
df_z = df.dropna(subset=['pressure', 'height', 'temperature', 'dewpoint']).copy()
df_z = df_z.sort_values('height', ascending=True).drop_duplicates('height')
if len(df_z) > 1:
valid_z_idx = np.concatenate(([True], np.diff(df_z['height'].values) > 0))
df_z = df_z[valid_z_idx]
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')
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')
# =========================================================================
# PIPELINE 3: IVT Integration (Requires Height AND Wind)
# =========================================================================
df_ivt = df.dropna(subset=['pressure', 'height', 'temperature', 'dewpoint', 'u_wind', 'v_wind']).copy()
df_ivt = df_ivt.sort_values('height', ascending=True).drop_duplicates('height')
if len(df_ivt) > 1:
valid_ivt_idx = np.concatenate(([True], np.diff(df_ivt['height'].values) > 0))
df_ivt = df_ivt[valid_ivt_idx]
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 (Requires Pressure and Wind)
# =========================================================================
df_wind = df.dropna(subset=['pressure', 'u_wind', 'v_wind']).copy()
df_wind = df_wind.sort_values('pressure', ascending=False).drop_duplicates('pressure')
if len(df_wind) > 1:
valid_p_wind_idx = np.concatenate(([True], np.diff(df_wind['pressure'].values) < 0))
df_wind = df_wind[valid_p_wind_idx]
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,
'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,
'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=(9, 9))
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')
# Plot Wind Barbs (Decimated by a factor of 4 to prevent overlapping ink blobs)
if len(metrics['p_wind']) > 0:
skew.plot_barbs(metrics['p_wind'][::3], metrics['u_wind'][::3], metrics['v_wind'][::3])
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
# Dynamically scale bottom Y-axis to guarantee surface data visibility (>1000mb if needed)
p_max = metrics['p_th'].m.max()
bottom_limit = max(1000, np.ceil(p_max / 10) * 10 + 10)
skew.ax.set_ylim(bottom_limit, 100)
skew.ax.set_xlim(-40, 60)
skew.ax.set_xlabel('Temperature (\xb0C)')
skew.ax.set_ylabel('Pressure (hPa)')
skew.ax.legend(loc='upper right')
skew.ax.set_title(f"Skew-T Log-P: {station_name}\n{dt.strftime('%Y-%m-%d %H:%M Z')}", fontsize=14)
filename = f"SkewT_{station_name.split(',')[0]}_{dt.strftime('%Y%m%d_%H')}Z.png"
plt.savefig(filename, dpi=150, bbox_inches='tight')
plt.close()
def plot_temp_pressure_height(mb, mr, dt):
if len(mb['p_z']) == 0 or len(mr['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(mb['p_z'].m, mb['z_z'].to('km').m, 'g-', label='Barrow P', linewidth=2)
ax1.plot(mr['p_z'].m, mr['z_z'].to('km').m, 'g--', label='Rochambeau P', linewidth=2)
ax2.plot(mb['T_z'].m, mb['z_z'].to('km').m, 'r-', label='Barrow T', linewidth=2)
ax2.plot(mr['T_z'].m, mr['z_z'].to('km').m, 'r--', label='Rochambeau T', linewidth=2)
fig.legend(loc='center right', bbox_to_anchor=(0.85, 0.75))
plt.title(f"Temperature and Pressure vs Height\n{dt.strftime('%Y-%m-%d %H:%M Z')}", y=1.08, fontsize=14)
filename = f"Temp_Pressure_vs_Height_{dt.strftime('%Y%m%d_%H')}Z.png"
plt.savefig(filename, dpi=150, bbox_inches='tight')
plt.close()
def plot_densities_height(mb, mr, dt):
if len(mb['rho_air']) == 0 or len(mr['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(mb['rho_wv'].m, mb['z_z'].to('km').m, 'b-', label='Barrow $\\rho_{wv}$', linewidth=2)
ax1.plot(mr['rho_wv'].m, mr['z_z'].to('km').m, 'b--', label='Rochambeau $\\rho_{wv}$', linewidth=2)
ax2.plot(mb['rho_air'].m, mb['z_z'].to('km').m, color='orange', linestyle='-', label='Barrow $\\rho_{air}$', linewidth=2)
ax2.plot(mr['rho_air'].m, mr['z_z'].to('km').m, color='orange', linestyle='--', label='Rochambeau $\\rho_{air}$', linewidth=2)
fig.legend(loc='center right', bbox_to_anchor=(0.85, 0.75))
plt.title(f"Densities vs Height\n{dt.strftime('%Y-%m-%d %H:%M Z')}", y=1.08, fontsize=14)
filename = f"Densities_vs_Height_{dt.strftime('%Y%m%d_%H')}Z.png"
plt.savefig(filename, dpi=150, bbox_inches='tight')
plt.close()
def plot_interpolated_pressure_diff(mb, mr, dt):
zb = mb['z_z'].m
zr = mr['z_z'].m
if len(zb) == 0 or len(zr) == 0: return
z_min = max(zb[0], zr[0])
z_max = min(zb[-1], zr[-1])
if z_max <= z_min: return
z_grid = np.arange(z_min, z_max, 11)
pb_interp = np.interp(z_grid, zb, mb['p_z'].m)
pr_interp = np.interp(z_grid, zr, mr['p_z'].m)
delta_p5 = 5 * (pr_interp - pb_interp)
fig, ax = plt.subplots(figsize=(8, 10))
ax.plot(delta_p5, z_grid / 1000.0, 'k-', linewidth=2)
ax.set_xlabel(r'$5 \times (P_{Rochambeau} - P_{Barrow})$ (hPa)', fontsize=12)
ax.set_ylabel('Height (km)', fontsize=12)
ax.set_title(f"Interpolated Pressure Difference vs Height\n{dt.strftime('%Y-%m-%d %H:%M Z')}", fontsize=14)
ax.grid(True, linestyle='--', alpha=0.7)
filename = f"Delta_Pressure5_{dt.strftime('%Y%m%d_%H')}Z.png"
plt.savefig(filename, dpi=150, bbox_inches='tight')
plt.close()
def main():
dt = get_user_datetime()
print(f"\nFetching data for {dt.strftime('%Y-%m-%d %H:%M Z')}...\n")
try:
print("Locating Barrow, AK data:")
df_b = fetch_sounding(dt, '70026', 'naconf', 'USM00070026')
print("\nLocating Rochambeau, FG data:")
df_r = fetch_sounding(dt, '81405', 'samer', 'FGM00081405')
except RuntimeError as e:
print(f"\nTERMINAL ERROR: {e}")
return
mb = compute_metrics(df_b)
mr = compute_metrics(df_r)
print("\n" + "="*50)
print(f"METRICS FOR {dt.strftime('%Y-%m-%d %H:%M Z')}")
print("="*50)
print("BARROW, AK:")
print(f" Precipitable Water : {mb['pw'].m:.2f} mm")
print(f" IVT u-component : {mb['ivt_u'].m:.2f} kg/m/s")
print(f" IVT v-component : {mb['ivt_v'].m:.2f} kg/m/s")
print(f" IVT Magnitude : {mb['ivt_mag'].m:.2f} kg/m/s")
print(f" IVT Direction : {mb['ivt_dir'].m:.1f} deg (Origin)")
print(f" 50mb ML Parcel Start : T={mb['T_ml']:.1f}, Td={mb['Td_ml']:.1f} (at {mb['p_ml']:.0f} hPa)")
print(f" ML CAPE (Virtual) : {mb['cape']:.1f} J/kg")
print(f" ML CIN : {mb['cin']:.1f} J/kg")
print("-" * 50)
print("ROCHAMBEAU, FG:")
print(f" Precipitable Water : {mr['pw'].m:.2f} mm")
print(f" IVT u-component : {mr['ivt_u'].m:.2f} kg/m/s")
print(f" IVT v-component : {mr['ivt_v'].m:.2f} kg/m/s")
print(f" IVT Magnitude : {mr['ivt_mag'].m:.2f} kg/m/s")
print(f" IVT Direction : {mr['ivt_dir'].m:.1f} deg (Origin)")
print(f" 50mb ML Parcel Start : T={mr['T_ml']:.1f}, Td={mr['Td_ml']:.1f} (at {mr['p_ml']:.0f} hPa)")
print(f" ML CAPE (Virtual) : {mr['cape']:.1f} J/kg")
print(f" ML CIN : {mr['cin']:.1f} J/kg")
print("="*50)
print("\nGenerating Skew-T diagrams with Wind Barbs...")
plot_skewt(mb, "Barrow, AK", dt)
plot_skewt(mr, "Rochambeau, FG", dt)
print("Generating Temperature & Pressure Overlay...")
plot_temp_pressure_height(mb, mr, dt)
print("Generating Densities Overlay...")
plot_densities_height(mb, mr, dt)
print("Generating 5x Pressure Difference Graph...")
plot_interpolated_pressure_diff(mb, mr, dt)
print("\nAll 5 plots saved successfully to your current working directory as PNG files.")
if __name__ == "__main__":
main()