From f155b546fa1318eab2d9697b40a8ceef7226f6fa Mon Sep 17 00:00:00 2001 From: acollow Date: Fri, 27 Feb 2026 14:01:08 -0500 Subject: [PATCH 1/5] expand variables for optics, allow for plotting of multiple experiments --- .../profile_plotting/plotaerosolprofiles.py | 207 +++++++++++++----- .../ASIA-AQ/profile_plotting/plotconfig.yaml | 28 ++- .../ASIA-AQ/profile_plotting/variablemap.yaml | 57 ++++- 3 files changed, 219 insertions(+), 73 deletions(-) diff --git a/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py b/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py index 12fe972d..b714d975 100755 --- a/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py +++ b/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py @@ -1,84 +1,185 @@ #!/usr/bin/env python +import sys +sys.path.insert(0, '/discover/nobackup/acollow/GMAOpyobs/src') + import numpy as np import xarray as xr import matplotlib.pyplot as plt import yaml from pyobs.icartt import ICARTT -import os, sys +import os from glob import glob import netCDF4 as nc +import re -#This code will generate a flight average vertical profile for a single ASIA-AQ flight. The syntax for running -#is ./plotaerosolprofiles.py VAR yyyymmdd where VAR must be a configuration listed in variablemap.yaml and the -#date must be a date in which an ASIA-AQ DC8 flight occurred as listed in plotconfig.yaml. Prior to running the -#python path needs to be set using "export PYTHONPATH=../AeroApps/install/lib/Python/" (or the full path to your -#AeroApps build). - -yamlkey_var=sys.argv[1] -flightdate=sys.argv[2] +yamlkey_var = sys.argv[1] +flightdate = sys.argv[2] with open('variablemap.yaml') as f: - var = yaml.safe_load(f) + var = yaml.safe_load(f) with open('plotconfig.yaml') as f: - config = yaml.safe_load(f) - -sampledfile=glob(config['geospath'] + '/asiaaq-mrg60_dc8_*' + flightdate + '*' + var[yamlkey_var]['sampling'] + '.nc4') -ictfile=glob(config['ictpath'] + '/asiaaq-mrg60_dc8_*' + flightdate + '*ict') - -# Load data -ds = xr.open_dataset(sampledfile[0], engine='netcdf4') -alt_geos = ds['H'].values -vardata_geos= ds[var[yamlkey_var]['geosvar']].values -if yamlkey_var=='OC': - vardata_geos=vardata_geos+ds['BR'].values -ds.close -vardata_geos=vardata_geos*var[yamlkey_var]['unitconversion'] + config = yaml.safe_load(f) +def evaluate_geosvar_expression(expr, ds): + """ + Evaluate a simple expression like 'OC + BR' or 'SCA / EXT' + Supports +, -, *, / operators + ds can be a single dataset or a dictionary of datasets + """ + # If it's just a simple variable name (no operators), return it directly + if not any(op in expr for op in ['+', '-', '*', '/']): + var_name = expr.strip() + if isinstance(ds, dict): + # For multi-file variables, extract the base variable name + return ds[var_name].values + else: + return ds[var_name].values + + # Find all variable names in the expression + var_names = re.findall(r'[A-Za-z_][A-Za-z0-9_]*', expr) + + # Create a namespace with the variables from the dataset(s) + if isinstance(ds, dict): + # Multiple datasets case + namespace = {name: ds[name].values for name in var_names if name in ds} + else: + # Single dataset case + namespace = {name: ds[name].values for name in var_names if name in ds} + + # Evaluate the expression + return eval(expr, {"__builtins__": {}}, namespace) + +# Load observational data once +ictfile = glob(config['ictpath'] + '/asiaaq-mrg60_dc8_*' + flightdate + '*ict') m = ICARTT(ictfile) #print([key for key in m.__dict__.keys()]) -ALT,obs_ts=m.Altitude_AGL_m_DIGANGI, m.__dict__[var[yamlkey_var]['obsvar']] - -# Get the indices of minimum distances for all columns at once -indices = np.array([np.argmin(np.abs(alt_geos[n,:] - ALT[n])) for n in range(len(ALT))]) +ALT, obs_ts = m.Altitude_AGL_m_DIGANGI, m.__dict__[var[yamlkey_var]['obsvar']] -# Use advanced indexing to get the corresponding values from GEOS and filter for available obs -geos_ts = np.array([vardata_geos[n,indices[n]] for n in range(len(ALT))]) -geos_ts[np.isnan(obs_ts)] = np.nan +# Apply unit conversion if specified in YAML +if 'obs_unit_conversion' in var[yamlkey_var]: + obs_ts = obs_ts * var[yamlkey_var]['obs_unit_conversion'] -#Find layer stats +# Define bin edges for vertical profiles bin_edges = np.arange(0, 4001, 250) heightbins = np.digitize(ALT, bin_edges) +height = np.arange(0.125, 4.1, step=0.25) + +# Calculate observational statistics once obs = np.array([np.nanmean(obs_ts[heightbins == h]) if np.any(heightbins == h) - else np.nan for h in range(1, 17)]) -geos = np.array([np.nanmean(geos_ts[heightbins == h]) if np.any(heightbins == h) - else np.nan for h in range(1, 17)]) -obs25 = np.array([np.nanpercentile(obs_ts[heightbins == h],25) if np.any(heightbins == h) - else np.nan for h in range(1, 17)]) -geos25 = np.array([np.nanpercentile(geos_ts[heightbins == h],25) if np.any(heightbins == h) - else np.nan for h in range(1, 17)]) -obs75 = np.array([np.nanpercentile(obs_ts[heightbins == h],75) if np.any(heightbins == h) - else np.nan for h in range(1, 17)]) -geos75 = np.array([np.nanpercentile(geos_ts[heightbins == h],75) if np.any(heightbins == h) - else np.nan for h in range(1, 17)]) - - -#Plot + else np.nan for h in range(1, 17)]) +obs25 = np.array([np.nanpercentile(obs_ts[heightbins == h], 25) if np.any(heightbins == h) + else np.nan for h in range(1, 17)]) +obs75 = np.array([np.nanpercentile(obs_ts[heightbins == h], 75) if np.any(heightbins == h) + else np.nan for h in range(1, 17)]) + +# Store model data for each experiment +model_data = {} +colors = ['red', 'blue', 'green', 'orange', 'purple', 'brown'] # Add more colors as needed + +# Loop through experiments +for exp_idx, exp in enumerate(config['experiments']): + exp_name = exp['name'] + exp_path = exp['path'] + + # Check if sampling is a dictionary (multi-file case) or a string (single file) + sampling_config = var[yamlkey_var]['sampling'] + + if isinstance(sampling_config, dict): + # Multi-file case: load multiple files and create a combined dataset dictionary + ds_dict = {} + for var_alias, sampling_name in sampling_config.items(): + sampledfile = glob(exp_path + '/asiaaq-mrg60_dc8_*' + flightdate + '*' + + sampling_name + '.nc4') + + if not sampledfile: + print(f"Warning: No sampled file found for {var_alias} in experiment {exp_name}") + continue + + # Load the dataset + ds_temp = xr.open_dataset(sampledfile[0], engine='netcdf4') + + # Store altitude from first file + if 'alt_geos' not in locals(): + alt_geos = ds_temp['H'].values + + # Extract the actual variable name from the alias (e.g., 'SCA_rh8' -> 'SCA') + actual_var = var_alias.split('_')[0] + ds_dict[var_alias] = ds_temp[actual_var] + ds_temp.close() + + # Evaluate expression using the dictionary of datasets + vardata_geos = evaluate_geosvar_expression(var[yamlkey_var]['geosvar'], ds_dict) + + else: + # Single file case (original behavior) + sampledfile = glob(exp_path + '/asiaaq-mrg60_dc8_*' + flightdate + '*' + + sampling_config + '.nc4') + + if not sampledfile: + print(f"Warning: No sampled file found for experiment {exp_name}") + continue + + # Load data + ds = xr.open_dataset(sampledfile[0], engine='netcdf4') + alt_geos = ds['H'].values + + # Evaluate the geosvar expression (handles both simple variables and expressions) + vardata_geos = evaluate_geosvar_expression(var[yamlkey_var]['geosvar'], ds) + + ds.close() + + # Get the indices of minimum distances for all columns at once + indices = np.array([np.argmin(np.abs(alt_geos[n, :] - ALT[n])) for n in range(len(ALT))]) + + # Use advanced indexing to get the corresponding values from GEOS and filter for available obs + geos_ts = np.array([vardata_geos[n, indices[n]] for n in range(len(ALT))]) + geos_ts[np.isnan(obs_ts)] = np.nan + + # Find layer stats for this experiment + geos = np.array([np.nanmean(geos_ts[heightbins == h]) if np.any(heightbins == h) + else np.nan for h in range(1, 17)]) + geos25 = np.array([np.nanpercentile(geos_ts[heightbins == h], 25) if np.any(heightbins == h) + else np.nan for h in range(1, 17)]) + geos75 = np.array([np.nanpercentile(geos_ts[heightbins == h], 75) if np.any(heightbins == h) + else np.nan for h in range(1, 17)]) + + # Store results + model_data[exp_name] = { + 'mean': geos, + 'p25': geos25, + 'p75': geos75, + 'color': exp.get('color', colors[exp_idx % len(colors)]) + } + +# Plot fig = plt.figure(figsize=(10, 8)) fig.tight_layout(pad=1.0) plt.rcParams['font.size'] = '14' -ax2 = plt.subplot(111) -height=np.arange(0.125, 4.1, step=0.25) -plt.plot(obs,height,'black') -plt.plot(geos,height,'red') -ax2.fill_betweenx(height,obs25,obs75,facecolor='dimgray') -ax2.fill_betweenx(height,np.squeeze(geos25),np.squeeze(geos75),facecolor='red',alpha=0.3) +ax2 = plt.subplot(111) + +# Plot observations +plt.plot(obs, height, 'black', linewidth=2, label='Obs') +ax2.fill_betweenx(height, obs25, obs75, facecolor='dimgray', alpha=0.3, + label='Obs 25th-75th') + +# Plot each model experiment +for exp_name, data in model_data.items(): + plt.plot(data['mean'], height, color=data['color'], linewidth=2, label=exp_name) + ax2.fill_betweenx(height, np.squeeze(data['p25']), np.squeeze(data['p75']), + facecolor=data['color'], alpha=0.2, + label=f"{exp_name} 25th-75th") + plt.xlabel(var[yamlkey_var]['fullname'] + ' ' + var[yamlkey_var]['units']) plt.ylabel('Height (km)') -ax2.legend(['Obs',config['expname'],'25th-75th Percentile'],fontsize=12) +ax2.legend(fontsize=10, loc='best') + country = config['flightdate'][int(flightdate)]['country'] datestr = config['flightdate'][int(flightdate)]['datestring'] plt.title(f"{country}, {datestr}") -plt.savefig(config['expname'] + '_' + var[yamlkey_var]['geosvar'] + flightdate + '.png') +# Save with all experiment names in filename +exp_names = '_'.join([exp['name'] for exp in config['experiments']]) +plt.savefig(f"{exp_names}_{yamlkey_var}_{flightdate}.png", + dpi=300, bbox_inches='tight') plt.show() diff --git a/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml b/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml index 8bb24416..c7b720df 100644 --- a/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml +++ b/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml @@ -1,6 +1,13 @@ -geospath: '/discover/nobackup/acollow/AeroApps/c360R_v11.6.1_asiaaqbaseline/sampled' -ictpath: '/discover/nobackup/projects/gmao/iesa/aerosol/campaigns/ASIA-AQ/archive/merges/v20250423' -expname: 'c360R_v11.6.1_asiaaqbaseline' +ictpath: '/discover/nobackup/projects/gmao/iesa/aerosol/campaigns/ASIA-AQ/archive/merges/v20260218' + +experiments: + - name: 'Baseline' + path: '/discover/nobackup/projects/gocart/hbian/proposal/2022-AsiaAQ/aaqoxgmi/sampled' + color: 'red' + - name: 'Updated Optics' + path: '/discover/nobackup/projects/gocart/hbian/proposal/2022-AsiaAQ/aaqoxh24/sampled' + color: 'blue' + flightdate: 20240206: country: 'Philippines' @@ -13,7 +20,7 @@ flightdate: datestring: '2/11/24' 20240213: country: 'Philippines' - datestring: '2/13/24' + datestring: '2/13/24' 20240215: country: 'Taiwan Overpass' datestring: '2/15/24' @@ -22,16 +29,16 @@ flightdate: datestring: '2/17/24' 20240226: country: 'South Korea' - datestring: '2/26/24' + datestring: '2/26/24' 20240308: country: 'South Korea' - datestring: '3/8/24' + datestring: '3/8/24' 20240310: country: 'South Korea' - datestring: '3/10/24' + datestring: '3/10/24' 20240311: country: 'South Korea' - datestring: '3/11/24' + datestring: '3/11/24' 20240313: country: 'Taiwan Overpass' datestring: '3/13/24' @@ -40,7 +47,7 @@ flightdate: datestring: '3/16/24' 20240318: country: 'Thailand' - datestring: '3/18/24' + datestring: '3/18/24' 20240321: country: 'Thailand' datestring: '3/21/24' @@ -50,3 +57,6 @@ flightdate: 20240327: country: 'Taiwan Overpass' datestring: '3/27/24' + + + diff --git a/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml b/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml index a665a11e..b577527a 100644 --- a/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml +++ b/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml @@ -4,46 +4,81 @@ SU: sampling: 'ams' units: '($\mu$g m$^{-3}$)' fullname: 'Sulfate' - unitconversion: 1 NH4: geosvar: 'NH4' obsvar: 'Ammonium_PM15_AMS_60s_JIMENEZ' sampling: 'ams' units: '($\mu$g m$^{-3}$)' fullname: 'Ammonium' - unitconversion: 1 NI: geosvar: 'NI' obsvar: 'Nitrate_PM15_AMS_60s_JIMENEZ' - sampling: 'ams' + sampling: 'ams' units: '($\mu$g m$^{-3}$)' fullname: 'Nitrate' - unitconversion: 1 OC: geosvar: 'OC' + geosvar_combined: ['OC', 'BR'] # List of variables to sum obsvar: 'OA_PM15_AMS_60s_JIMENEZ' sampling: 'ams' units: '($\mu$g m$^{-3}$)' fullname: 'Organic Carbon' - unitconversion: 1 SS: geosvar: 'SS' obsvar: 'Seasalt_PM15_AMS_60s_JIMENEZ' sampling: 'ams' units: '($\mu$g m$^{-3}$)' fullname: 'Sea Salt' - unitconversion: 1 BC: geosvar: 'BC' obsvar: 'BC_mass_MOORE' sampling: 'bc_sp2' - units: '(ng m$^{-3})' - fullname: 'Black Carbon' - unitconversion: 1000 + units: '($\mu$g m$^{-3}$)' + fullname: 'Black Carbon' +SO2: + geosvar: 'SO2' + obsvar: 'SO2_CIT_CROUNSE' + sampling: 'inst3_aer_Nv' + units: 'ppt' + fullname: 'Sulfur Dioxide' RH: geosvar: 'RH' obsvar: 'RHw_DLH_DISKIN' sampling: 'inst3_aer_Nv' - units: '(%)' + units: '%' fullname: 'Relative Humidity' - unitconversion: 100 +Ext532dry: + geosvar: 'EXT' + obsvar: 'Ext532_submicron_dry_MOORE' + obs_unit_conversion: 0.001 # Convert Mm^-1 to km^-1 + sampling: '532_submicron_aop_STP_rh2' + units: '(km$^{-1}$)' + fullname: '532 nm Dry Extinction' +fRH550: + geosvar: 'SCA_rh8 / SCA_rh2' + obsvar: 'fRH550_RH20to80_MOORE' + sampling: + SCA_rh8: '550_submicron_aop_STP_rh8' + SCA_rh2: '550_submicron_aop_STP_rh2' + units: '' + fullname: 'fRH at 550 nm' +Ext532amb: + geosvar: 'EXT' + obsvar: 'Ext532_submicron_amb_MOORE' + obs_unit_conversion: 0.001 # Convert Mm^-1 to km^-1 + sampling: '532_submicron_aop_STP' + units: '(km$^{-1}$)' + fullname: '532 nm Ambient Extinction' +Scat550dry: + geosvar: 'SCA' + obsvar: 'Sc550_submicron_MOORE' + obs_unit_conversion: 0.001 # Convert Mm^-1 to km^-1 + sampling: '550_submicron_aop_STP_rh2' + units: '(km$^{-1}$)' + fullname: '550 nm Dry Scattering' +SSA550dry: + geosvar: 'SCA / EXT' + obsvar: 'SSA_dry_550nm_MOORE' + sampling: '550_submicron_aop_STP_rh2' + units: '' + fullname: '550 nm Dry SSA' From f4daa01cb5e415248afc2df8ca52e497c63ca074 Mon Sep 17 00:00:00 2001 From: Allie Collow <31133739+acollow@users.noreply.github.com> Date: Wed, 27 May 2026 15:16:16 -0400 Subject: [PATCH 2/5] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index eee019e3..2b3417ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ - Extract TROPOMI pixels over MPL sites - ASIA-AQ mission scripts - capability to produce flight averaged vertical profile plots for ASIA-AQ +- Expanded capability for ASIA-AQ vertical profile plots ### Fixed ### Changed From 208f81623111a75c8ecfd15e0b9343a375850b9e Mon Sep 17 00:00:00 2001 From: Allie Collow <31133739+acollow@users.noreply.github.com> Date: Wed, 27 May 2026 15:40:53 -0400 Subject: [PATCH 3/5] removing training spaces --- .../missions/ASIA-AQ/profile_plotting/variablemap.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml b/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml index b577527a..f076c45c 100644 --- a/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml +++ b/src/Components/missions/ASIA-AQ/profile_plotting/variablemap.yaml @@ -13,7 +13,7 @@ NH4: NI: geosvar: 'NI' obsvar: 'Nitrate_PM15_AMS_60s_JIMENEZ' - sampling: 'ams' + sampling: 'ams' units: '($\mu$g m$^{-3}$)' fullname: 'Nitrate' OC: @@ -34,7 +34,7 @@ BC: obsvar: 'BC_mass_MOORE' sampling: 'bc_sp2' units: '($\mu$g m$^{-3}$)' - fullname: 'Black Carbon' + fullname: 'Black Carbon' SO2: geosvar: 'SO2' obsvar: 'SO2_CIT_CROUNSE' From e4b1dfb233d259a7d8aaa18e5c241cbab61137d9 Mon Sep 17 00:00:00 2001 From: Allie Collow <31133739+acollow@users.noreply.github.com> Date: Wed, 27 May 2026 15:42:50 -0400 Subject: [PATCH 4/5] remove trailing spaces --- .../ASIA-AQ/profile_plotting/plotconfig.yaml | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml b/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml index c7b720df..9aef559d 100644 --- a/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml +++ b/src/Components/missions/ASIA-AQ/profile_plotting/plotconfig.yaml @@ -20,7 +20,7 @@ flightdate: datestring: '2/11/24' 20240213: country: 'Philippines' - datestring: '2/13/24' + datestring: '2/13/24' 20240215: country: 'Taiwan Overpass' datestring: '2/15/24' @@ -29,16 +29,16 @@ flightdate: datestring: '2/17/24' 20240226: country: 'South Korea' - datestring: '2/26/24' + datestring: '2/26/24' 20240308: country: 'South Korea' - datestring: '3/8/24' + datestring: '3/8/24' 20240310: country: 'South Korea' - datestring: '3/10/24' + datestring: '3/10/24' 20240311: country: 'South Korea' - datestring: '3/11/24' + datestring: '3/11/24' 20240313: country: 'Taiwan Overpass' datestring: '3/13/24' @@ -47,7 +47,7 @@ flightdate: datestring: '3/16/24' 20240318: country: 'Thailand' - datestring: '3/18/24' + datestring: '3/18/24' 20240321: country: 'Thailand' datestring: '3/21/24' @@ -57,6 +57,3 @@ flightdate: 20240327: country: 'Taiwan Overpass' datestring: '3/27/24' - - - From 890391418ed50bf6b408c311a403ff8371f4b820 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <69310046+patricia-nasa@users.noreply.github.com> Date: Fri, 29 May 2026 12:41:41 -0400 Subject: [PATCH 5/5] Update plotaerosolprofiles.py --- .../missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py b/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py index b714d975..19d3f09f 100755 --- a/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py +++ b/src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py @@ -1,8 +1,6 @@ #!/usr/bin/env python import sys -sys.path.insert(0, '/discover/nobackup/acollow/GMAOpyobs/src') - import numpy as np import xarray as xr import matplotlib.pyplot as plt