Advanced: Solar Validation Against HALO-(AC)³ Aircraft Measurements

Contents

Advanced: Solar Validation Against HALO-(AC)³ Aircraft Measurements#

In this notebook, we’ll compare pyRadtran’s solar broadband simulations against measurements from multiple aircraft during the HALO-(AC)³ campaign. We validate downwelling solar irradiance against observations and compute statistical error metrics.

Setup#

import pyradtran
from pyradtran import load_config
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from pathlib import Path
import pandas as pd
import logging

# Configure logging for pyradtran
logging.getLogger('pyradtran').setLevel(logging.CRITICAL)

# ── Simulation parameters ─────────────────────────────────────────────────────
# Paths come from ~/.pyradtran/config.yaml (master config)
# radiosonde_base is overridden in the saved YAML to use local sounding data
cfg = load_config()

cfg.simulation_defaults.rte_solver           = "disort"
cfg.simulation_defaults.mol_abs_param        = "reptran medium"
cfg.simulation_defaults.source               = "solar"
cfg.simulation_defaults.wavelength_nm        = [290, 4500]
cfg.simulation_defaults.integrate_wavelength = True
cfg.simulation_defaults.h2o_source           = "radiosonde"
cfg.simulation_defaults.ozone_du             = None   # use radiosonde
cfg.simulation_defaults.h2o_mm               = None   # use radiosonde
cfg.simulation_defaults.surface_temperature_k = 253.15
cfg.simulation_defaults.output_altitudes_km  = [0.0]
cfg.simulation_defaults.output_columns       = ["zout", "lambda", "sza", "edir", "eglo", "edn", "eup", "enet", "albedo"]
cfg.execution.max_workers                    = 32
cfg.execution.cleanup_temp_files             = False
cfg.execution.timeout_seconds               = 3600

# Override radiosonde_base to use local sounding data bundled in the notebook
cfg.paths.radiosonde_base = str(Path("~/pyRadtran/book/notebooks/radiosonde/").expanduser())

config_path = Path("config/halo-ac3_bbr_all_aircraft.yaml")
cfg.to_yaml(config_path)
print(f"Config saved to {config_path}")

2026-04-19 23:44:35,870 - pyradtran.config - INFO - Configuration written to config/halo-ac3_bbr_all_aircraft.yaml

Config saved to config/halo-ac3_bbr_all_aircraft.yaml

Data Preparation#

# Load CSV and convert to xarray Dataset
df = pd.read_csv('data/HALO-AC3_HALO_P5_P6_aircraft_broadband_radiation_clear_sky_with_ocean_100s.csv', parse_dates=['time'])
df = df.set_index('time')

ds = xr.Dataset.from_dataframe(df).drop_duplicates('time').sortby('time')
ds_to_sim = ds.rename({'Lat': 'latitude', 'Lon': 'longitude', 'Alt': 'altitude'})

# Use uncorrected solar flux where available
ds_to_sim['F_down_solar'].loc[{'time': ds_to_sim.F_down_solar_uncorr.notnull()}] = ds_to_sim.F_down_solar_uncorr.loc[{'time': ds_to_sim.F_down_solar_uncorr.notnull()}]

# Compute albedo from upwelling/downwelling solar flux
albedo = (ds_to_sim['F_up_solar'] / ds_to_sim['F_down_solar'])
albedo = albedo.where(albedo < 1)  # Remove unphysical values > 1
ds_to_sim['albedo'] = ('time', albedo.astype(float).values)
ds_to_sim['albedo'] = ds_to_sim['albedo'].fillna(0.5)  # Default albedo for NaN values

# Convert altitude from meters to km
ds_to_sim['altitude'] = ds_to_sim['altitude'] / 1e3
ds_to_sim = ds_to_sim.assign_coords(altitude=ds_to_sim['altitude'])

# Define output altitude grid (0–12 km in 1 km steps)
ds_to_sim['altitude'] = ('altitude', np.arange(0, 13, 1))

Simulation Results#

Spectral radiative transfer calculations completed using DISORT solver with multi-level output altitude grid.

Comparing pyRadtran with Published Results#

Finally, we compare our pyRadtran simulations against the aircraft observations and the published simulation results from the dataset (Becker et al., 2023). This three-way comparison validates both the measurement pipeline and our simulation approach:

X = ds_sim_alt['eglo']
Y = ds_to_sim['F_down_solar']

fig, ax_scatter = plt.subplots(1, 1, figsize=(6, 6))

# --- Scatter: pyRadtran simulated vs observed irradiance ---
cmap = plt.get_cmap('viridis', n_colors)
scatter = ax_scatter.scatter(X / 1e3, Y, c=airplane_codes, cmap=cmap, s=25, alpha=0.7, edgecolor='k')
ax_scatter.set_xlabel('Simulated Irradiance (pyRadtran) (W/m²)')
ax_scatter.set_ylabel('Observed Irradiance (W/m²)')

n_colors = len(airplane_types)

# 1:1 reference line and axis limits
ax_scatter.set_xlim(0, 600)
ax_scatter.set_ylim(0, 600)
ax_scatter.plot(np.arange(0, 600, 100), np.arange(0, 600, 100), color='gray', linestyle='--', linewidth=0.5)

# Inset colorbar for aircraft categories
cax = inset_axes(ax_scatter, width="30%", height="5%", loc='upper right')
norm = mcolors.BoundaryNorm(np.arange(n_colors + 1) - 0.5, n_colors)
cb = fig.colorbar(scatter, cax=cax, orientation='horizontal', ticks=np.arange(n_colors), norm=norm)
cb.set_ticklabels(airplane_types, fontweight='bold')
cb.ax.tick_params(labelsize=8)
cax.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=True, labeltop=False)
cax.spines[['top', 'right', 'left', 'bottom']].set_visible(False)
for a in [ax_scatter]:
    a.grid(True, linestyle='--', alpha=0.5)
    a.spines[['top', 'right']].set_visible(False)

# --- Error metrics ---
nan_mask = ~np.isnan(X) & ~np.isnan(Y)
X = X[nan_mask].values / 1000
Y = Y[nan_mask].values

rmse = np.sqrt(mean_squared_error(X, Y))
bias = np.mean(X - Y)
r2 = r2_score(X, Y)
mae = mean_absolute_error(X, Y)

stats_text = f'RMSE: {rmse:.2f} W/m²\nMAE:  {mae:.2f} W/m²\nBias:  {bias:.2f} W/m²\nR²: {r2:.2f}'
ax_scatter.text(0.05, 0.95, stats_text, transform=ax_scatter.transAxes, fontsize=10, verticalalignment='top')

Text(0.05, 0.95, 'RMSE: 18.29 W/m²\nMAE:  13.94 W/m²\nBias:  -1.18 W/m²\nR²: 0.98')

../_images/0956d13106dd9e3e7a459c0f1160990c11d865d2c6d84068b365f910ec2ae16a.png

df_to_compare = ds[['F_down_solar', 'F_down_solar_sim']]
df_to_compare['F_down_solar_sim_pyRadtran'] = ('time', ds_sim_alt.eglo.values / 1000)  # Convert to W/m²

df_to_compare = df_to_compare.dropna('time', how='any').to_dataframe()
df_to_compare

### make this cool pandas scatter comparison plot where you can compare multiple columns 
import seaborn as sns
import matplotlib.pyplot as plt
#sns.set(style="whitegrid")
g = sns.pairplot(df_to_compare, kind='scatter', diag_kind='kde', markers=["o", "x"], height=3, aspect=1.2, corner=True)

../_images/40696ebf1ecc12dc4a7c2ebea6152ce6845412cca735da07c507bf9d18a4bf4f.png