Updates to ASIA-AQ Vertical Profile Plots

CHANGELOG.md

@@ -18,7 +18,7 @@
 - Rearranged and renamed some of the ARCSIX codes
 - refactored ARCSIX mission post-processing code for suitability for archive upload
 - refactored INSPYRE trajectory preparation and plotting scripts
-
+- Expanded capability for ASIA-AQ vertical profile plots
 - updated monthly_mod.py to reduce the time it takes to run
 
 ## [v2.3.0] - 2025-07-25
@@ -30,6 +30,7 @@
 ### Changed
 - removed all NNR scripts from GAAS_App and relocated them in AeroML repository
 
+
 ## [v2.2.0] - 2025-07-15
 
 ### Added

src/Components/missions/ASIA-AQ/profile_plotting/plotaerosolprofiles.py

@@ -1,84 +1,183 @@
 #!/usr/bin/env python
 
+import sys
 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']]
+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))])
+# 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']
 
-# 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
+# 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()

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'

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'
   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})'
+  units: '($\mu$g m$^{-3}$)'
   fullname: 'Black Carbon'
-  unitconversion: 1000
+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'