Merge pull request #316 from mbruckner-work/develop_tropomifix

Changes to MOPITT pairing

Author: rschwant

docs/appendix/yaml.rst

@@ -89,6 +89,12 @@ or other interactive Python session,
 as it gives you some visual indication of the progress of multi-file data loading
 and some parts of the processing.
 
+**pairing_kwargs:** This is an optional argument. This dictionary allows for specifying keyword arguments for pairing methods.
+First level should be the observation type (e.g. "sat_grid_clm", "sat_swath_clm"). Then under the observation type label provide the specific pairing options for your application.
+   
+   * **apply_ak:** This is an optional argument used for pairing of satellite data. When no pairing keyword arguments are specified it will default to True. This should be set to True when application of satellite averaging kernels or apriori data to model observations is desired.
+   * **mod_to_overpass:** This is an optional argument used for pairing of satellite data. When set to True the model data will be pre-processed to the published local overpass time for the satellite. As of now, local overpass times are hard-wired.
+
 Models
 ------
 All input for each instance of the model class. First level should be the model 
@@ -120,9 +126,12 @@ data (e.g., surf_only: True).
 Typically this is set at the horizontal resolution of your model * 1.5. Setting 
 this to a smaller value will speed up the pairing process. 
 
-**apply_ak:** This is an optional argument used for pairing of satellite data. This
-should be set to True when application of satellite averaging kernels or apriori data 
-to model observations is desired. 
+**apply_ak:** Removed. Instead, specify ``pairing_kwargs`` in the analysis section.
+
+**is_global:** Optional boolean argument to specify if the model dataset is global or
+regional. Used in some satellite pairing methods to indicate if a longitude wrap should 
+be applied. Defaults to False when unspecified and xesmf-based satellite pairing methods
+will assume the model datset is regional.
 
 **mapping:** This is the mapping dictionary for all variables to be plotted. 
 For each observational dataset, add a mapping dictionary where the model 
@@ -204,6 +213,8 @@ See :doc:`../getting_started/downloading_obs` for more details.
 **obs_type:** The observation type. Options are: "pt_sfc" or point surface. Adding 
 options for Aircraft and Satellite observations are under development.
 
+**sat_type:** The satellite observation type. Options include: "mopitt_l3", "omps_l3", "omps_nm", "modis_l2", and "tropomi_l2_no2". Additional options are under development. 
+
 **data_proc:** This section stores all of the data processing information.
 
    * **filter_dict:** This is a dictionary used to filter the observation data 

docs/develop/development_team.rst

@@ -52,8 +52,8 @@ in regional coupled chemistry-meteorology models. My plans include testing and
 expanding the capability of MELODIES MONET for evaluating simulations with
 research and operational models of fire impacts on air quality and weather.
 
-**Margaret Bruckner:**
-I am a graduate student at the University of Wisconsin-Madison. My development plans
+**Maggie Bruckner:**
+I am an NRC postdoc associated with NOAA CSL. My development plans
 primarily focus on adding capabilities for comparison of satellite observations to model
 output and expanding data processing options. 
 

examples/jupyter_notebooks/UFS-AQM_mopitt_example.ipynb

@@ -0,0 +1,79 @@
+analysis:
+        start_time: '2023-06-27-00:00:00'
+        end_time: '2023-06-29-00:59:00'
+        output_dir: /scratch1/BMC/rcm2/mbruckner/mopitt/ufs_aqm_analysis/plots
+        debug: False
+        pairing_kwargs:
+                sat_grid_clm:
+                        apply_ak: True
+                        mod_to_overpass: True
+model:
+        ufsaqm: 
+                files: '/scratch1/BMC/rcm2/rhs/monet_example/AEROMMA/UFS-AQM/cmaq54_OriRave1/aqm.202306*/12/aqm.t12z.dyn.*.nc'
+                mod_type: 'rrfs'
+                is_global: False
+                mapping: 
+                        mopitt: 
+                                co: column
+                plot_kwargs: 
+                        color: 'purple'
+                        marker: '^'
+                        linestyle: 'dotted'
+obs:
+        mopitt:
+                filename: /scratch1/BMC/rcm2/mbruckner/mopitt/v009_daily/*he5 
+                obs_type: sat_grid_clm
+                sat_type: mopitt_l3
+                variables:
+                        column:
+                                ylabel_plot: '$10^{18} molec/cm^{2}$'
+                                vmin_plot: 0.0 # optMin for y-axis during plotting. To apply to a plot, change restrict_yaxis = True.
+                                vmax_plot: 16
+                                vdiff_plot: 15.0
+plots:
+        #plot_ts: 
+        #        type: 'timeseries'
+        #        fig_kwargs:
+        #                figsize: [12,6]
+        #        default_plot_kwargs: # Opt to define defaults for all plots. Model kwargs overwrite these.
+        #                linewidth: 2.0
+        #                markersize: 10.
+        #        text_kwargs: #Opt
+        #                fontsize: 18.
+        #        domain_type: ['all']
+        #        domain_name: ['CONUS']
+        #        data: ['mopitt_ufsaqm']
+        #        data_proc:
+        #                rem_obs_nan: True # True: Remove all points where model or obs variable is NaN. False: Remove only points where model variable is NaN.
+        #                ts_select_time: 'time' #Time used for avg and plotting: Options: 'time' for UTC or 'time_local'
+        #                ts_avg_window: 'h' # Options: None for no averaging or list pandas resample rule (e.g., 'h', 'D')
+        #                set_axis: False #If select True, add vmin_plot and vmax_plot for each variable in obs.
+        plot_grp1:
+                type: 'gridded_spatial_bias'
+                fig_kwargs:
+                        states: True
+                        figsize: [10,5]
+                text_kwargs:
+                        fontsize: 16.
+                domain_type: ['all']
+                domain_name: ['UFS-AQM']
+                data: ['mopitt_ufsaqm']
+                data_proc:
+                        rem_obs_nan: True
+                        set_axis: True
+      #plot_grp2:
+      #          type: 'taylor'
+      #          fig_kwargs:
+      #                  figsize: [8.8]
+      #          default_plot_kwargs:
+      #                  linewidth: 2.0
+      #                  markersize: 10.
+      #          text_kwargs:
+      #                  fontsize: 16.
+      #          domain_type: ['all'] #List of domain types: 'all' or any domain in obs file. (e.g., airnow: epa_region, state_name, siteid, etc.)
+      #          domain_name: ['Global'] # of domain names. If domain_type = all domain_name is used in plot title.
+      #          data: ['mopitt_ufsaqm'] # make this a list of pairs in obs_model where the obs is the obs label and model is the model_label
+      #          data_proc:
+      #                  rem_obs_nan: True # True: Remove all points where model or obs variable is NaN. False: Remove only points where model variable is NaN.
+      #                 set_axis: True #If select True, add ty_scale for each variable in obs.
+

examples/yaml/control_ufsaqm_mopitt.yaml

@@ -247,7 +247,7 @@ def open_sat_obs(self, time_interval=None, control_dict=None):
         None
         """
         from .util import time_interval_subset as tsub
-        
+        from glob import glob
         try:
             if self.sat_type == 'omps_l3':
                 print('Reading OMPS L3')
@@ -275,6 +275,13 @@ def open_sat_obs(self, time_interval=None, control_dict=None):
                 else: flst = self.file
                 self.obj = mio.sat._mopitt_l3_mm.open_dataset(flst, ['column','pressure_surf','apriori_col',
                                                                           'apriori_surf','apriori_prof','ak_col'])
+
+                # Determine if monthly or daily product and set as attribute
+                if any(mtype in glob(self.file)[0] for mtype in ('MOP03JM','MOP03NM','MOP03TM')):
+                    self.obj.attrs['monthly'] = True
+                else:
+                    self.obj.attrs['monthly'] = False
+                            
             elif self.sat_type == 'modis_l2':
                 # from monetio import modis_l2
                 print('Reading MODIS L2')
@@ -427,7 +434,7 @@ class model:
     def __init__(self):
         """Initialize a :class:`model` object."""
         self.model = None
-        self.apply_ak = False
+        self.is_global = False
         self.radius_of_influence = None
         self.mod_kwargs = {}
         self.file_str = None
@@ -716,6 +723,8 @@ def __init__(self):
         self.obs_gridded_dataset = None
         self.add_logo = True
         """bool, default=True : Add the MELODIES MONET logo to the plots."""
+        self.pairing_kwargs = {}
+
 
     def __repr__(self):
         return (
@@ -809,6 +818,10 @@ def read_control(self, control=None):
                 = [[time_stamps[n], time_stamps[n+1]]
                     for n in range(len(time_stamps)-1)]
 
+        # specific arguments for pairing options
+        if 'pairing_kwargs' in self.control_dict['analysis'].keys():
+            self.pairing_kwargs = self.control_dict['analysis']['pairing_kwargs']
+            
         # Enable Dask progress bars? (default: false)
         enable_dask_progress_bars = self.control_dict["analysis"].get(
             "enable_dask_progress_bars", False)
@@ -911,8 +924,8 @@ def open_models(self, time_interval=None,load_files=True):
                 # this is the model type (ie cmaq, rapchem, gsdchem etc)
                 m.model = self.control_dict['model'][mod]['mod_type']
                 # set the model label in the dictionary and model class instance
-                if "apply_ak" in self.control_dict['model'][mod].keys():
-                    m.apply_ak = self.control_dict['model'][mod]['apply_ak']
+                if "is_global" in self.control_dict['model'][mod].keys():
+                    m.is_global = self.control_dict['model'][mod]['is_global']
                 if 'radius_of_influence' in self.control_dict['model'][mod].keys():
                     m.radius_of_influence = self.control_dict['model'][mod]['radius_of_influence']
                 else:
@@ -1330,6 +1343,12 @@ def pair_data(self, time_interval=None):
                 # TODO: add other network types / data types where (ie flight, satellite etc)
                 # if sat_swath_clm (satellite l2 column products)
                 elif obs.obs_type.lower() == 'sat_swath_clm':
+                    # grab kwargs for pairing. Use default if not specified
+                    pairing_kws = {'apply_ak':True,'mod_to_overpass':False}
+                    for key in self.pairing_kwargs[obs.obs_type.lower()]:
+                        pairing_kws[key] = self.pairing_kwargs[obs.obs_type.lower()][key]
+                    if 'apply_ak' not in self.pairing_kwargs[obs.obs_type.lower()]:
+                        print('WARNING: The satellite pairing option apply_ak is being set to True because it was not specified in the YAML. Pairing will fail if there is no AK available.')
 
                     if obs.sat_type == 'omps_nm':
 
@@ -1340,7 +1359,7 @@ def pair_data(self, time_interval=None):
                         if 'time' in obs.obj.dims:
                             obs.obj = obs.obj.sel(time=slice(self.start_time,self.end_time))
                             obs.obj = obs.obj.swap_dims({'time':'x'})
-                        if mod.apply_ak == True:
+                        if pairing_kws['apply_ak'] == True:
                             model_obj = mod.obj[keys+['pres_pa_mid','surfpres_pa']]
 
                             paired_data = sutil.omps_nm_pairing_apriori(model_obj,obs.obj,keys)
@@ -1370,7 +1389,7 @@ def pair_data(self, time_interval=None):
                         model_obj = mutil.cal_model_no2columns(model_obj)
                         #obs_dat = obs.obj.sel(time=slice(self.start_time.date(),self.end_time.date())).copy()
 
-                        if mod.apply_ak == True:
+                        if pairing_kws['apply_ak'] == True:
                             paired_data = sutil.trp_interp_swatogrd_ak(obs.obj, model_obj)
                         else:
                             paired_data = sutil.trp_interp_swatogrd(obs.obj, model_obj)
@@ -1394,6 +1413,12 @@ def pair_data(self, time_interval=None):
 
                 # if sat_grid_clm (satellite l3 column products)
                 elif obs.obs_type.lower() == 'sat_grid_clm':
+                    # grab kwargs for pairing. Use default if not specified
+                    pairing_kws = {'apply_ak':True,'mod_to_overpass':False}
+                    for key in self.pairing_kwargs[obs.obs_type.lower()]:
+                        pairing_kws[key] = self.pairing_kwargs[obs.obs_type.lower()][key]
+                    if 'apply_ak' not in self.pairing_kwargs[obs.obs_type.lower()]:
+                        print('WARNING: The satellite pairing option apply_ak is being set to True because it was not specified in the YAML. Pairing will fail if there is no AK available.')
                     if len(keys) > 1: 
                         print('Caution: More than 1 variable is included in mapping keys.')
                         print('Pairing code is calculating a column for {}'.format(keys[0])) 
@@ -1416,13 +1441,22 @@ def pair_data(self, time_interval=None):
 
                     elif obs.sat_type == 'mopitt_l3':
                         from .util import satellite_utilities as sutil
-                        if mod.apply_ak: 
+                        
+                        if pairing_kws['apply_ak']: 
                             model_obj = mod.obj[keys+['pres_pa_mid']]
                             # trim to only data within analysis window, as averaging kernels can't be applied outside it
-                            obs_dat = obs.obj.sel(time=slice(self.start_time.date(),self.end_time.date()))#.copy()
-                            model_obj = model_obj.sel(time=slice(self.start_time.date(),self.end_time.date()))#.copy()
+                            obs_dat = obs.obj.sel(time=slice(self.start_time.date(),self.end_time.date()))
+                            model_obj = model_obj.sel(time=slice(self.start_time.date(),self.end_time.date()))
+                            
+                            # Sample model to observation overpass time
+                            if pairing_kws['mod_to_overpass']:
+                                print('sampling model to 10:30 local overpass time')
+                                overpass_datetime = pd.date_range(self.start_time.replace(hour=10,minute=30),
+                                                                  self.end_time.replace(hour=10,minute=30),freq='D')
+                                model_obj = sutil.mod_to_overpasstime(model_obj,overpass_datetime)
+                            
                             # interpolate model to observation, calculate column with averaging kernels applied
-                            paired = sutil.mopitt_l3_pairing(model_obj,obs_dat,keys[0])
+                            paired = sutil.mopitt_l3_pairing(model_obj,obs_dat,keys[0],global_model=mod.is_global)
                             p = pair()
                             p.type = obs.obs_type
                             p.obs = obs.label

melodies_monet/driver.py

@@ -291,18 +291,16 @@ def cal_amf_wrfchem(scatw, wrfpreslayer, tpreslev, troppres, wrfno2layer_molec,
     vertical_pres = []
     vertical_scatw = []
     vertical_wrfp = []
-
+    
     for llb in range(len(lb[0])):
         yy = lb[0][llb]
         xx = lb[1][llb]
         vertical_pres = tpreslev[:,yy,xx] # mli, update to new dimension
         vertical_scatw = scatw[yy,xx,:]
         vertical_wrfp = wrfpreslayer[yy,xx,:]
-
         f = interpolate.interp1d(np.log10(vertical_pres[:]),vertical_scatw[:], fill_value="extrapolate")# relationship between pressure to avk
         wrfavk[yy,xx,:] = f(np.log10(vertical_wrfp[:])) #wrf-chem averaging kernel
 
-
     for l in range(nz-1):
         # check if it's within tropopause
         preminus[:,:]         = wrfpreslayer[:,:,l] - troppres[:,:]