cleaned up exptag[_list] and outid; added/clarified comments in ObsFcstAna post-proc package

gmao-rreichle · gmao-rreichle · commit 990bce2a94dd · 2025-06-11T15:54:51.000-04:00
diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_maps.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_maps.py
@@ -2,6 +2,7 @@
 
 """
 Sample script for plotting maps of long-term data assimilation diagnostics.
+Computes Nobs-weighted avg of each metric across all species.
 Requires saved files with monthly sums (see Get_ObsFcstAna_stat.py).
 Stats of *normalized* OmFs are approximated!
 """
@@ -29,7 +30,7 @@ def plot_OmF_maps(postproc_obj, stats, fig_path='./'):
      
     Nmin = 20
 
-    # Then computer metrics of O-F, O-A, etc. based on above computed
+    # Then compute metrics of O-F, O-A, etc. based on above computed
     N_data = stats['N_data']
     O_mean = stats['obs_mean']
     # mean(x-y) = E[x] - E[y]   
@@ -66,7 +67,10 @@ def plot_OmF_maps(postproc_obj, stats, fig_path='./'):
     OmA_stdv[     N_data < Nmin] = np.nan
     N_data[       N_data < Nmin] = 0
 
-    # Combine metrics of individual species using weighted averaging
+    # Compute Nobs-weighted avg of each metric across all species.
+    # Typically used for SMAP Tb_h/h from asc and desc overpasses,
+    # or ASCAT soil moisture from Metop-A/B/C.
+    # DOES NOT MAKE SENSE IF, SAY, SPECIES HAVE DIFFERENT UNITS!
     Nobs_data     = np.nansum(              N_data, axis=1)
     OmF_mean      = np.nansum(OmF_mean     *N_data, axis=1)/Nobs_data
     OmF_stdv      = np.nansum(OmF_stdv     *N_data, axis=1)/Nobs_data
@@ -76,7 +80,7 @@ def plot_OmF_maps(postproc_obj, stats, fig_path='./'):
     OmA_stdv      = np.nansum(OmA_stdv     *N_data, axis=1)/Nobs_data
 
     # Plotting
-    expid = postproc_obj.outid
+    exptag = postproc_obj.exptag
 
     fig, axes = plt.subplots(2,2, figsize=(18,10))
     plt.rcParams.update({'font.size':14})
@@ -89,23 +93,23 @@ def plot_OmF_maps(postproc_obj, stats, fig_path='./'):
                 # crange is [cmin, cmax]
                 crange =[0, np.ceil((end_time-start_time).days/150)*300]
                 colormap = plt.get_cmap('jet',20)
-                title_txt = expid + ' Tb Nobs '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m')
+                title_txt = exptag + ' Tb Nobs '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m')
                 units = '[-]'
             if i == 0 and j ==1:
                 tile_data = OmF_mean
                 crange =[-3, 3]
                 colormap = plt.get_cmap('bwr', 15) 
-                title_txt = expid + ' Tb O-F mean '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m')
+                title_txt = exptag + ' Tb O-F mean '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m')
             if i == 1 and j == 0:
                 tile_data = OmF_stdv
                 crange =[0, 15]
                 colormap = plt.get_cmap ('jet',15)
-                title_txt = expid + ' Tb O-F stdv '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m')
+                title_txt = exptag + ' Tb O-F stdv '+ start_time.strftime('%Y%m')+'_'+end_time.strftime('%Y%m')
             if i == 1 and j == 1:
                 tile_data = OmF_norm_stdv
                 crange =[0, 15]
                 colormap = plt.get_cmap ('jet',15)
-                title_txt = expid + ' Tb normalized O-F stdv (approx!) '+ start_time.strftime('%Y%m%d')+'_'+end_time.strftime('%Y%m%d')
+                title_txt = exptag + ' Tb normalized O-F stdv (approx!) '+ start_time.strftime('%Y%m%d')+'_'+end_time.strftime('%Y%m%d')
 
             colormap.set_bad(color='0.9') # light grey, 0-black, 1-white
 
@@ -159,7 +163,7 @@ def plot_OmF_maps(postproc_obj, stats, fig_path='./'):
 
     plt.tight_layout()
     # Save figure to file
-    fig.savefig(fig_path+'Map_OmF_'+ expid +'_'+start_time.strftime('%Y%m')+'_'+\
+    fig.savefig(fig_path+'Map_OmF_'+ exptag +'_'+start_time.strftime('%Y%m')+'_'+\
                         (end_time+timedelta(days=-1)).strftime('%Y%m')+'.png')
     #plt.show()
     plt.close(fig)
diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_timeseries.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/Plot_stats_timeseries.py
@@ -2,6 +2,7 @@
 
 """
 Sample script for plotting monthly time series of spatially averaged data assimilation diagnostics.
+Computes Nobs-weighted avg of each metric across all species.
 Requires saved files with monthly sums (see Get_ObsFcstAna_stat.py).
 """
 
@@ -17,11 +18,11 @@
     
 def Plot_monthly_OmF_bars(postproc_obj, fig_path='./'):
     import pickle
-    expid = postproc_obj.outid
+    exptag      = postproc_obj.exptag
     start_month = postproc_obj.start_time
-    end_month = postproc_obj.end_time
+    end_month   = postproc_obj.end_time
 
-    stats_file = fig_path + 'spatial_stats_'+expid+'_'+start_month.strftime('%Y%m')+ \
+    stats_file = fig_path + 'spatial_stats_'+exptag+'_'+start_month.strftime('%Y%m')+ \
             '_'+(end_month+timedelta(days=-1)).strftime('%Y%m')+'.pkl'
     
     if os.path.isfile(stats_file):
@@ -47,7 +48,10 @@ def Plot_monthly_OmF_bars(postproc_obj, fig_path='./'):
 
             OmFm,OmFs,OmAm,OmAs,Nobsm = postproc_obj.calc_spatial_stats_from_sums(current_month)
 
-            # Average individual species into a single value
+            # Compute Nobs-weighted avg of each metric across all species.
+            # Typically used for SMAP Tb_h/h from asc and desc overpasses,
+            # or ASCAT soil moisture from Metop-A/B/C.
+            # DOES NOT MAKE SENSE IF, SAY, SPECIES HAVE DIFFERENT UNITS!
             Nobsm = np.nansum(     Nobsm)
             OmFm  = np.nansum(OmFm*Nobsm)/Nobsm
             OmFs  = np.nansum(OmFs*Nobsm)/Nobsm
@@ -81,13 +85,13 @@ def Plot_monthly_OmF_bars(postproc_obj, fig_path='./'):
     plt.grid(True, linestyle='--', alpha=0.5)
 
     plt.xticks(ticks=date_vec[::2], labels=date_vec[::2])
-    plt.title(expid+ 'monthly '+plot_var)
+    plt.title(exptag+ 'monthly '+plot_var)
     plt.xlim(-1, len(date_vec)+1)
     plt.ylim(-.1, 2.)
 
     plt.tight_layout()
     #plt.show()
-    plt.savefig(fig_path+'Bars_'+plot_var+expid+date_vec[0]+'_'+\
+    plt.savefig(fig_path+'Bars_'+plot_var+exptag+date_vec[0]+'_'+\
                     date_vec[-1]+'.png')
     plt.close()
 
diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py
@@ -20,11 +20,11 @@
 
 class postproc_ObsFcstAna:
     
-    def __init__(self, exp_list, start_time, end_time, sum_path='./', outid=None):
+    def __init__(self, exp_list, start_time, end_time, sum_path='./'):
         self.exp_list      = exp_list
         self.expdir_list   = [item['expdir'] for item in exp_list]
         self.expid_list    = [item['expid']  for item in exp_list]
-        self.exptag_list   = [item['exptag'] for item in exp_list]
+        self.exptag        = exp_list[0]['exptag']
         self.domain        = exp_list[0]['domain']
         self.start_time    = start_time 
         self.end_time      = end_time
@@ -45,11 +45,7 @@ def __init__(self, exp_list, start_time, end_time, sum_path='./', outid=None):
                 else:
                     self.obs_from = exp_idx
         if self.obs_from < 0: self.obs_from = 0       # by default, obs data are from exp_list[0]
-        print(f"obs data are from {self.exptag_list[self.obs_from]}")   
-
-        # Tag of the output sums files. If not provided, use main exptag 
-        if outid is None:
-            self.outid =  self.exptag_list[0]
+        print(f"obs data are from {exp_list[self.obs_from]['expid']}")   
 
         # Verify the configuration every time when current class is initialized
         # to avoid saving sums with different configs in the same directory
@@ -63,7 +59,7 @@ def __init__(self, exp_list, start_time, end_time, sum_path='./', outid=None):
             config_list.append({var:exp[var] for var in config_verify if var in exp})
 
         # File of configuration for verification
-        f_config = self.sum_path + '/' + self.outid + '_config.yaml'
+        f_config = self.sum_path + '/' + self.exptag + '_config.yaml'
 
         # Save a new file or compare current configuration with previously saved 
         if not os.path.exists(f_config):
@@ -194,7 +190,6 @@ def compute_monthly_sums(self, date_time):
     def save_monthly_sums(self):
         expdir_list    = self.expdir_list
         expid_list     = self.expid_list
-        exptag_list    = self.exptag_list
         domain         = self.domain
         var_list       = self.var_list
         tc             = self.tilecoord
@@ -208,7 +203,7 @@ def save_monthly_sums(self):
             mo_path = self.sum_path + '/Y'+ date_time.strftime('%Y') + '/M' + date_time.strftime('%m') + '/'
             os.makedirs(mo_path, exist_ok=True)
     
-            fout = self.outid + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4'
+            fout = self.exptag + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4'
 
             fout = mo_path + fout
             
@@ -270,7 +265,7 @@ def calc_temporal_stats_from_sums(self, write_to_nc=True, fout_stats='./stats.nc
                 
                 fpath = self.sum_path + '/Y'+ date_time.strftime('%Y') + '/M' + date_time.strftime('%m') + '/'
                     
-                fout  = self.outid + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4'
+                fout  = self.exptag + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4'
 
                 fout  = fpath + fout
                                 
@@ -355,7 +350,7 @@ def calc_spatial_stats_from_sums(self, date_time):
         var_list = ['obs_obs', 'obs_fcst','obs_ana']
         
         mo_path = self.sum_path + '/Y'+ date_time.strftime('%Y') + '/M' + date_time.strftime('%m') + '/'            
-        fnc4_sums = mo_path + self.outid + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4'
+        fnc4_sums = mo_path + self.exptag + '.ens_avg.ldas_ObsFcstAna_sums.' + date_time.strftime('%Y%m') +'.nc4'
         
         mdata_sum = {}
         mdata2_sum = {}
@@ -405,7 +400,7 @@ def calc_spatial_stats_from_sums(self, date_time):
             # var(x) = E[x2] - (E[x])^2
             data_var[var] = data2_mean[var] - data_mean[var]**2
 
-        # Computer metrics of O-F, O-A, etc. based on above stats
+        # Compute metrics of O-F, O-A, etc. based on above stats
         O_mean = data_mean['obs_obs']
         F_mean = data_mean['obs_fcst']
         A_mean = data_mean['obs_ana']
diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/user_config.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/user_config.py
@@ -26,21 +26,23 @@ def get_config():
 
     exp_main = {
         'expdir'        : '/discover/nobackup/projects/gmao/merra/iau/merra_land/SMAP_runs/SMAP_Nature_v11/',
-        'expid'         : 'DAv8_SMOSSMAP',
-        'exptag'        : 'DAMulti_SMAP', 
-        'domain'        : 'SMAP_EASEv2_M36_GLOBAL',
-        'da_t0'         : 3,                              # (fractional) UTC hour of first land analysis 
-        'da_dt'         : 10800,                          # ObsFcstAna file interval in seconds
-        'species_list'  : [5,6,7,8],                      # indices of species to be processed 
-        'obsparam_time' : "20150401_0000"                 # time stamp of obsparam file (YYYYMMDD_HHMM)
+        'expid'         : 'DAv8_SMOSSMAP',             # GEOSldas exp ID of simulation
+        'exptag'        : 'DAMulti_SMAP',              # string used in output file names for sums & stats,
+                                                       #   can be same as expid or different -- e.g., reflect info
+                                                       #   about subset of included species or about cross-masking
+        'domain'        : 'SMAP_EASEv2_M36_GLOBAL',  
+        'da_t0'         : 3,                           # (fractional) UTC hour of first land analysis 
+        'da_dt'         : 10800,                       # ObsFcstAna file interval in seconds
+        'species_list'  : [5,6,7,8],                   # indices of species to be processed 
+        'obsparam_time' : "20150401_0000"              # time stamp of obsparam file (YYYYMMDD_HHMM)
     }
 
     # Optional experiment(s) can be added for cross-masking or extracting obs from a different experiment.
     #
     # All optional experiments and the main experiment must have identical tilecoords.
     # If the default "species" number/order do not match, set "species_list" accordingly to force a match.
     # Output will be cross-masked between all specified experiments.    
-
+    #
     # Forecasts and analyses are always from the main experiment.
     # Observations are from the experiment with 'use_obs' set to True (default is exp_main).  The most
     #   likely use case for reading obs from a supplemental experiment is when computing OmF etc diagnostics
@@ -50,21 +52,25 @@ def get_config():
     exp_sup1 = {
         'expdir'        : '/discover/nobackup/projects/gmao/merra/iau/merra_land/SMAP_runs/SMAP_Nature_v11/',
         'expid'         : 'DAv8_M36',
-        'exptag'        : 'DASMAP_SMAP', 
         'domain'        : 'SMAP_EASEv2_M36_GLOBAL',
-        'use_obs'       : True,                           # if True, use obs data from this exp
-        'species_list'  : [1,2,3,4],                      # indices of species to be processed 
-        'obsparam_time' : "20150401_0000"                 # time stamp of obsparam file (YYYYMMDD_HHMM)
+        'use_obs'       : True,                        # if True, use obs data from this exp
+        'species_list'  : [1,2,3,4],                   # indices of species to be processed,
+                                                       #   must identify same species as selected in main exp
+        'obsparam_time' : "20150401_0000"              # time stamp of obsparam file (YYYYMMDD_HHMM)
     }
 
-    # Convert experiments input to a list; first entry must be exp_main 
+    # Convert experiments input to a list; first entry must be exp_main: 
 
-    exp_list = [exp_main, exp_sup1]
+    #exp_list = [exp_main]               # no cross-masking
+    exp_list = [exp_main, exp_sup1]      # cross-mask exp_main with exp_sup1
     
-    # Top level directory to store monthly sum files; can use the experiment directory or a different path;
-    # /Yyyyy/Mmm/ is added automatically for individual months
+    # Top level directory for all output from this package:
+
+    out_path = '/discover/nobackup/[USERNAME]/SMAP_test/'
 
-    out_path = '/discover/nobackup/qliu/SMAP_test/'
+    # Directory for monthly sum files:
+    # - Can use the experiment directory or a different path.
+    # - Automatically appends /Yyyyy/Mmm/ for individual months.
     
     sum_path = out_path+'/'+exp_main['expid']+'/output/'+exp_main['domain']+'/ana/ens_avg/'
 
