add MOPITT CO L3 pairing and updates to satplot

examples/yaml/control_mopitt.yaml

@@ -5,19 +5,65 @@
 # Seting axis values - If set_axis = True in data_proc section of each plot_grp the yaxis for the plot will be set based on the values specified in the obs section for each variable. If set_axis is set to False, then defaults will be used. 'vmin_plot' and 'vmax_plot' are needed for 'timeseries', 'spatial_overlay', and 'boxplot'. 'vdiff_plot' is needed for 'spatial_bias' plots and'ty_scale' is needed for 'taylor' plots. 'nlevels' or the number of levels used in the contour plot can also optionally be provided for spatial_overlay plot. If set_axis = True and the proper limits are not provided in the obs section, a warning will print, and the plot will be created using the default limits.
 #-------------------------------------------------------------------
 analysis:
-  start_time: '2020-08-01'
-  end_time:   '2020-08-03'
-  output_dir: /glade/work/buchholz/melodies-monet/output
+  start_time: '2019-07-15'
+  end_time:   '2019-07-19'
+  output_dir: /ships19/aqda/mbruckner/Mopitt_mm_test
   debug: True
 
 #-------------------------------------------------------------------
-# model:
-
+model:
+  raqms:
+    files: /ships19/aqda/mbruckner/monet_plots/mopitt_example_raqms/uwhyb*nc
+    mod_type: 'raqms'
+    apply_ak: True # for satellite comparison, applies averaging kernels/apriori when true. Default to False
+    variables: #Opt
+      ico: # specifying to switch units to ppbv
+        need: True
+    mapping: #model species name : obs species name
+      mopitt_l3: 
+        ico: column #The mapping tables need to contain the same species for all models.
+    plot_kwargs: #Opt 
+      color: 'purple'
+      marker: '^'
+      linestyle: 'dotted'
 #-------------------------------------------------------------------
 obs:    
   mopitt_l3: # obs label
-    filename: /glade/work/buchholz/data/MOPITT/MOP03JM-2020*
+    filename: /ships19/aqda/mbruckner/Mopitt_mm_test/MOP03J-201907*he5
     obs_type: sat_grid_clm
-    variables: null
-
-
+    sat_type: mopitt_l3
+    variables: 
+      column:
+        ylabel_plot: '$10^{18} molec/cm^{2}$'    
+#-------------------------------------------------------------------
+plots:
+  plot_grp3:
+    type: 'gridded_spatial_bias' #'gridded_spatial_bias' # plot type
+    fig_kwargs: #For all spatial plots, specify map_kwargs here too.
+      states: True
+      figsize: [10, 5] # figure size 
+    text_kwargs: #Opt
+      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'] #List of domain names. If domain_type = all domain_name is used in plot title.
+    data: ['mopitt_l3_raqms'] # 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 vdiff_plot for each variable in obs.
+
+
+  plot_grp2:
+    type: 'taylor' # plot type
+    fig_kwargs: #Opt to define figure options
+      figsize: [8,8] # figure size if multiple plots
+    default_plot_kwargs: # Opt to define defaults for all plots. Model kwargs overwrite these.
+      linewidth: 2.0
+      markersize: 10.
+    text_kwargs: #Opt
+      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_l3_raqms'] # 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_omps_nm-raqms.yaml

@@ -12,7 +12,7 @@ analysis:
   read:
       paired:
           method: 'netcdf'
-          filenames: {'omps_nm_raqms':['firex_omps_201908*_20190*_omps_nm_raqms.nc4']}
+          filenames: {'omps_nm_raqms':['firex_omps_*_omps_nm_raqms.nc4']} 
 model:
   raqms: # model label
     files: /ships19/aqda/mbruckner/MELODIES-MONET-1/examples/jupyter_notebooks/raqms-files.txt
@@ -27,33 +27,35 @@ model:
         o3vmr: ozone_column #The mapping tables need to contain the same species for all models.
     plot_kwargs: #Opt 
       color: 'purple'
-      marker: '+'
+      marker: '^'
       linestyle: 'dotted'
 obs:
   omps_nm: # obs label
     filename: /ships19/aqda/mbruckner/OMPS-NPP/O3-daily/2019/nadir_mapper/OMPS-NPP_NMTO3-L2_v2.1_2019m08*t*.h5
     obs_type: sat_swath_clm
+    sat_type: omps_nm
     variables: #Opt 
-
+      ozone_column:
+        ylabel_plot: 'DU'
 
 plots:
-  plot_grp1:
-    type: 'timeseries' # plot type
-    fig_kwargs: #Opt to define figure options
-      figsize: [12,6] # figure size if multiple plots
-    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'] #List of domain types: 'all' or any domain in obs file. (e.g., airnow: epa_region, state_name, siteid, etc.) 
-    domain_name: ['Global'] #List of domain names. If domain_type = all domain_name is used in plot title.
-    data: ['omps_nm_raqms'] # 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.
-      ts_select_time: 'time' #Time used for avg and plotting: Options: 'time' for UTC or 'time_local'
-      ts_avg_window: 'min' # 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: 'timeseries' # plot type
+#    fig_kwargs: #Opt to define figure options
+#      figsize: [12,6] # figure size if multiple plots
+#    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'] #List of domain types: 'all' or any domain in obs file. (e.g., airnow: epa_region, state_name, siteid, etc.) 
+#    domain_name: ['Global'] #List of domain names. If domain_type = all domain_name is used in plot title.
+#    data: ['omps_nm_raqms_grid'] # 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.
+#      ts_select_time: 'time' #Time used for avg and plotting: Options: 'time' for UTC or 'time_local'
+#      ts_avg_window: 'min' # 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_grp2:
     type: 'taylor' # plot type
     fig_kwargs: #Opt to define figure options
@@ -76,10 +78,9 @@ plots:
       figsize: [10, 5] # figure size 
     text_kwargs: #Opt
       fontsize: 16.
-      #label: '$\\Delta$ DU'
     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'] #List of domain names. If domain_type = all domain_name is used in plot title.
     data: ['omps_nm_raqms'] # 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 vdiff_plot for each variable in obs.
+      set_axis: True #If select True, add vdiff_plot for each variable in obs.

melodies_monet/driver.py

@@ -217,7 +217,7 @@ def open_sat_obs(self,time_interval=None):
         try:
             if self.sat_type == 'omps_l3':
                 print('Reading OMPS L3')
-                self.obj = mio.sat._omps_l3_mm.read_OMPS_l3(self.file)
+                self.obj = mio.sat._omps_l3_mm.open_dataset(self.file)
             elif self.sat_type == 'omps_nm':
                 print('Reading OMPS_NM')
                 if time_interval is not None:
@@ -236,7 +236,8 @@ def open_sat_obs(self,time_interval=None):
 
             elif self.sat_type == 'mopitt_l3':
                 print('Reading MOPITT')
-                self.obj = mio.sat._mopitt_l3_mm.read_mopittdataset(self.file, 'column')
+                self.obj = mio.sat._mopitt_l3_mm.open_dataset(self.file, ['column','pressure_surf','apriori_col',
+                                                                          'apriori_surf','apriori_prof','ak_col'])
             elif self.sat_type == 'modis_l2':
                 from monetio import modis_l2
                 print('Reading MODIS L2')
@@ -393,7 +394,7 @@ def glob_files(self):
             self.files = [tutorial.fetch_example(example_id)]
         else:
             self.files = sort(glob(self.file_str))
-
+            
         # add option to read list of files from text file
         _, extension = os.path.splitext(self.file_str)
         if extension.lower() == '.txt':
@@ -493,7 +494,7 @@ def open_model_files(self, time_interval=None, control_dict=None):
                 self.mod_kwargs.update({'scrip_file' : self.scrip_file})            
                 self.obj = mio.models._cesm_se_mm.open_mfdataset(self.files,**self.mod_kwargs)
                 #self.obj, self.obj_scrip = read_cesm_se.open_mfdataset(self.files,**self.mod_kwargs)
-                #self.obj.monet.scrip = self.obj_scrip
+                #self.obj.monet.scrip = self.obj_scrip      
             elif 'raqms' in self.model.lower():
                 if len(self.files) > 1:
                     self.obj = mio.raqms.open_mfdataset(self.files,**self.mod_kwargs)
@@ -935,6 +936,7 @@ def pair_data(self, time_interval=None):
                         print('After pairing: ', paired_data)
                     # this outputs as a pandas dataframe.  Convert this to xarray obj
                     p = pair()
+                    print('saving pair')
                     p.obs = obs.label
                     p.model = mod.label
                     p.model_vars = keys
@@ -949,13 +951,15 @@ def pair_data(self, time_interval=None):
                 # if sat_swath_clm (satellite l2 column products)
                 elif obs.obs_type.lower() == 'sat_swath_clm':
 
-                    if obs.label == 'omps_nm':
+                    if obs.sat_type == 'omps_nm':
 
                         from .util import satellite_utilities as sutil
 
-                        #necessary observation index things 
-                        #the along track coordinate dim sometimes needs to be time and other times an unassigned 'x'
-                        obs.obj = obs.obj.swap_dims({'time':'x'})
+                        # necessary observation index things 
+                        ## the along track coordinate dim sometimes needs to be time and other times an unassigned 'x'
+                        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:
                             model_obj = mod.obj[keys+['pres_pa_mid','surfpres_pa']]
 
@@ -976,22 +980,46 @@ def pair_data(self, time_interval=None):
                         self.paired[label] = p
                 # if sat_grid_clm (satellite l3 column products)
                 elif obs.obs_type.lower() == 'sat_grid_clm':
-                    if obs.label == 'omps_l3':
+                    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])) 
+                    if obs.sat_type == 'omps_l3':
                         from .util import satellite_utilities as sutil
                         # trim obs array to only data within analysis window
-                        obs_dat = obs.obj.sel(time=slice(self.start_time.date(),self.end_time.date())).copy()
-                        model_obsgrid = sutil.omps_l3_daily_o3_pairing(mod.obj,obs_dat,'o3vmr')
-                        # combine model and observations into paired dataset
-                        obs_dat['o3vmr'] = (['time','x','y'],model_obsgrid.sel(time=slice(self.start_time.date(),self.end_time.date())).data)
+                        obs_dat = obs.obj.sel(time=slice(self.start_time.date(),self.end_time.date()))#.copy()
+                        mod_dat = mod.obj.sel(time=slice(self.start_time.date(),self.end_time.date()))
+                        paired_obsgrid = sutil.omps_l3_daily_o3_pairing(mod_dat,obs_dat,keys[0])
+
                         p = pair()
                         p.type = obs.obs_type
                         p.obs = obs.label
                         p.model = mod.label
                         p.model_vars = keys
                         p.obs_vars = obs_vars
-                        p.obj = obs_dat
+                        p.obj = paired_obsgrid
                         label = '{}_{}'.format(p.obs,p.model)
                         self.paired[label] = p
+                    elif obs.sat_type == 'mopitt_l3':
+                        from .util import satellite_utilities as sutil
+                        if mod.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()
+                            # interpolate model to observation, calculate column with averaging kernels applied
+                            paired = sutil.mopitt_l3_pairing(model_obj,obs_dat,keys[0])
+                            p = pair()
+                            p.type = obs.obs_type
+                            p.obs = obs.label
+                            p.model = mod.label
+                            p.model_vars = keys
+                            p.model_vars[0] += '_column_model'
+                            p.obs_vars = obs_vars
+                            p.obj = paired
+                            label ='{}_{}'.format(p.obs,p.model)
+                            self.paired[label] = p
+                        else:
+                            print("Pairing without averaging kernel has not been enabled for this dataset")
     def concat_pairs(self):
         """Read and concatenate all observation and model time interval pair data,
         populating the :attr:`paired` dict.
@@ -1001,7 +1029,7 @@ def concat_pairs(self):
         None
         """
         pass
-    
+
     ### TODO: Create the plotting driver (most complicated one)
     # def plotting(self):
     def plotting(self):
@@ -1024,8 +1052,10 @@ def plotting(self):
         pair_keys = list(self.paired.keys())
         if self.paired[pair_keys[0]].type.lower() == 'pt_sfc':
             from .plots import surfplots as splots,savefig
+            if self.debug: print('using surfplots')
         else:
             from .plots import satplots as splots,savefig
+            if self.debug: print('using satplots')
 
         # Disable figure count warning
         initial_max_fig = plt.rcParams["figure.max_open_warning"]
@@ -1093,9 +1123,10 @@ def plotting(self):
 
                                 pairdf_all = p.obj.swap_dims({'x':'time'})
                             # squash lat/lon dimensions into single dimension
-                            elif obs_type == 'sat_grid_clm':
-                                pairdf_all = p.obj.stack(ll=['x','y'])
-                                pairdf_all = pairdf_all.rename_dims({'ll':'y'})
+                            ## 2024-03 MEB rechecking necessity of this.
+                            #elif obs_type == 'sat_grid_clm':
+                            #    pairdf_all = p.obj.stack(ll=['x','y'])
+                            #    pairdf_all = pairdf_all.rename_dims({'ll':'y'})
                             else:
                                 pairdf_all = p.obj
                             # Select only the analysis time window.
@@ -1283,8 +1314,10 @@ def plotting(self):
                             else:
                                 vmin = None
                                 vmax = None
-                            # Select time to use as index.
-                            pairdf = pairdf.set_index(grp_dict['data_proc']['ts_select_time'])
+                            # Select time to use as index. 
+                            # 2024-03-01 MEB needs to only apply if pandas. fails for xarray
+                            if obs_type == 'pt_sfc':
+                                pairdf = pairdf.set_index(grp_dict['data_proc']['ts_select_time'])
                             # Specify ts_avg_window if noted in yaml file. 
                             if 'ts_avg_window' in grp_dict['data_proc'].keys():
                                 a_w = grp_dict['data_proc']['ts_avg_window']

melodies_monet/plots/satplots.py

@@ -668,8 +668,10 @@ def make_spatial_bias_gridded(df, column_o=None, label_o=None, column_m=None,
                       domain_type=None, domain_name=None, fig_dict=None, 
                       text_dict=None,debug=False):
 
-    """Creates difference plot for satellite and model data. Needs to be altered for cases where more than 1 overpass for a location,
-    eg. more than 1 day of data."""
+    """Creates difference plot for satellite and model data.
+        For data in swath format, overplots all differences
+        For data on regular grid, mean difference.
+    """
     if debug == False:
         plt.ioff()
 
@@ -691,8 +693,10 @@ def make_spatial_bias_gridded(df, column_o=None, label_o=None, column_m=None,
         ylabel = column_o
 
     #Take the difference for the model output - the sat output
-    diff_mod_min_obs = (df[column_o] - df[column_m]).squeeze()
-
+    diff_mod_min_obs = (df[column_m] - df[column_o]).squeeze()
+    #Take mean over time, 
+    if len(diff_mod_min_obs.dims) == 3:
+        diff_mod_min_obs = diff_mod_min_obs.mean('time')
 
     #Determine the domain
     if domain_type == 'all' and domain_name == 'CONUS':
@@ -739,7 +743,7 @@ def make_spatial_bias_gridded(df, column_o=None, label_o=None, column_m=None,
     c = ax.axes.scatter(df.longitude,df.latitude,c=diff_mod_min_obs,cmap=cmap,s=2,norm=norm)
     plt.gcf().canvas.draw() 
     plt.tight_layout(pad=0)
-    plt.title(title_add + label_o + ' - ' + label_m,fontweight='bold',**text_kwargs)
+    plt.title(title_add + label_m + ' - ' + label_o,fontweight='bold',**text_kwargs)
     ax.axes.set_extent(map_kwargs['extent'],crs=ccrs.PlateCarree())    
 
     #Uncomment these lines if you update above just to verify colorbars are identical.
@@ -757,7 +761,7 @@ def make_spatial_bias_gridded(df, column_o=None, label_o=None, column_m=None,
     position_m = model_ax.get_position()
     position_c = cax.get_position()
     cax.set_position([position_c.x0, position_m.y0, position_c.x1 - position_c.x0, (position_m.y1-position_m.y0)*1.1])
-    cax.set_ylabel(ylabel,fontweight='bold',**text_kwargs)
+    cax.set_ylabel('$\Delta$'+ylabel,fontweight='bold',**text_kwargs)
     cax.tick_params(labelsize=text_kwargs['fontsize']*0.8,length=10.0,width=2.0,grid_linewidth=2.0)    
 
     #plt.tight_layout(pad=0)