Add NetCDF4 ObsFcstAna postprocessing support

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added ObsFcstAna postprocessing support for NetCDF4 diagnostics, with fallback to legacy binary ObsFcstAna files.
 - Added support for river routing.
 - Added optional NetCDF4 output mode for ObsFcstAna, including NetCDF metadata and runtime context. Changed namelist variable "out_ObsFcstAna" from logical to integer.
 - Added support for ensemble simulations for routing; for now, landice hardwired to NensLandice=1.
@@ -26,6 +27,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- Fixed `read_obs_param()` parsing for the current obsparam format by reading forecast variable names and units.
 - Fixed crashes in debug mode; including adding an extra mask to work around a compiler bug in "where(elemental)".
 - Fixed string matching for EASE tile file to accommodate new "EASE*-Pfafstetter" tile file for runoff routing purposes.
 - Fixed GEOSlandpert build when MKL is unavailable by enabling MKL-specific code paths only when MKL is detected.
@@ -183,4 +185,3 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Inaugural version.  0-diff vs. GEOSldas v18.0.0.
 
 -----------------------------
-

GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py

@@ -14,7 +14,7 @@
 from datetime                    import datetime, timedelta
 from dateutil.relativedelta      import relativedelta
 from netCDF4                     import Dataset, date2num
-from read_GEOSldas               import read_ObsFcstAna, read_tilecoord, read_obs_param
+from read_GEOSldas               import read_ObsFcstAna, read_ObsFcstAna_nc4, read_tilecoord, read_obs_param
 
 from helper.write_nc4            import write_sums_nc4, write_stats_nc4
 
@@ -94,6 +94,7 @@ def compute_monthly_sums(self, date_time):
         n_spec = len(obsparam_list[0])
 
         date_time = date_time.replace(hour=int(self.da_t0), minute=int(np.mod(self.da_t0,1)*60))
+        month_string = date_time.strftime('%Y%m')
         stop_time = date_time + relativedelta(months=1) 
 
         data_sum  = {}
@@ -108,19 +109,28 @@ def compute_monthly_sums(self, date_time):
             data_sum[ var] = np.zeros((n_tile, n_spec))
             data2_sum[var] = np.zeros((n_tile, n_spec))
 
+        n_files_read = 0
+
         while date_time < stop_time:
 
             # read the list of experiments at each time step (OFA="ObsFcstAna")
             OFA_list = []
             for i in range(len(expdir_list)):
-                fname = expdir_list[i]+expid_list[i]+'/output/'+self.domain+'/ana/ens_avg/Y'+ \
-                                  date_time.strftime('%Y') + '/M' + \
-                                  date_time.strftime('%m') + '/'  + \
-                                  expid_list[i]+'.ens_avg.ldas_ObsFcstAna.' + \
-                                  date_time.strftime('%Y%m%d_%H%M') +'z.bin'
-                if os.path.isfile(fname):
-                    print('read '+fname)
-                    OFA_list.append(read_ObsFcstAna(fname))
+                fname_base = expdir_list[i]+expid_list[i]+'/output/'+self.domain+'/ana/ens_avg/Y'+ \
+                                      date_time.strftime('%Y') + '/M' + \
+                                      date_time.strftime('%m') + '/'  + \
+                                      expid_list[i]+'.ens_avg.ldas_ObsFcstAna.' + \
+                                      date_time.strftime('%Y%m%d_%H%M') +'z'
+                fname_nc4 = fname_base + '.nc4'
+                fname_bin = fname_base + '.bin'
+                if os.path.isfile(fname_nc4):
+                    print('read '+fname_nc4)
+                    OFA_list.append(read_ObsFcstAna_nc4(fname_nc4))
+                    n_files_read += 1
+                elif os.path.isfile(fname_bin):
+                    print('read '+fname_bin)
+                    OFA_list.append(read_ObsFcstAna(fname_bin))
+                    n_files_read += 1
 
             data_all=[]
             for OFA, obs_param in zip(OFA_list,obsparam_list):
@@ -181,6 +191,9 @@ def compute_monthly_sums(self, date_time):
 
             date_time = date_time + timedelta(seconds=da_dt)
 
+        if n_files_read == 0:
+            raise FileNotFoundError('No ObsFcstAna .nc4 or .bin files found for ' + month_string)
+
         return N_data, data_sum, data2_sum, oxf_sum, oxa_sum, fxa_sum
 
     # ----------------------------------------------------------------------------------------------------------
@@ -212,8 +225,13 @@ def save_monthly_sums(self):
             if  not os.path.isfile(fout):
                 print('computing monthly sums...')
                 # compute monthly sums
-                mN_data, mdata_sum, mdata2_sum, moxf_sum, moxa_sum, mfxa_sum = \
-                    self.compute_monthly_sums(date_time)
+                try:
+                    mN_data, mdata_sum, mdata2_sum, moxf_sum, moxa_sum, mfxa_sum = \
+                        self.compute_monthly_sums(date_time)
+                except FileNotFoundError as e:
+                    print(f'WARNING: skipping {date_time.strftime("%Y%m")} - {e}')
+                    date_time = date_time + relativedelta(months=1)
+                    continue
 
                 # save monthly sums in nc4 file
                 write_sums_nc4(fout, mN_data,mdata_sum, mdata2_sum, moxf_sum, moxa_sum, mfxa_sum, obsparam_list[0])

GEOSldas_App/util/shared/python/read_GEOSldas.py

@@ -40,6 +40,8 @@ def read_obs_param(fname):
             param['nodata']         = float(fid.readline().strip())
             param['varname']        =       fid.readline().strip().strip('"')
             param['units']          =       fid.readline().strip().strip('"')
+            param['fcstvarname']    =       fid.readline().strip().strip('"')
+            param['fcstunits']      =       fid.readline().strip().strip('"')
             param['path']           =       fid.readline().strip().strip('"')
             param['name']           =       fid.readline().strip().strip('"')
             param['maskpath']       =       fid.readline().strip().strip('"')
@@ -292,4 +294,109 @@ def read_ObsFcstAna(fname, isLDASsa=False):
             'obs_ana'    : obs_ana, 
             'obs_anavar' : obs_anavar}
 
+# ----------------------------------------------------------------------------
+#
+# reader for GEOSldas ObsFcstAna file (NetCDF-4)
+
+def read_ObsFcstAna_nc4(fname):
+
+    from datetime import datetime
+    from netCDF4 import Dataset
+
+    nodata = -9999
+
+    date_time = {
+        'year'  : nodata,
+        'month' : nodata,
+        'day'   : nodata,
+        'hour'  : nodata,
+        'min'   : nodata,
+        'sec'   : nodata,
+        'dofyr' : nodata,
+        'pentad': nodata
+    }
+
+    obs_assim   = []
+    obs_species = []
+    obs_tilenum = []
+    obs_lon     = []
+    obs_lat     = []
+    obs_obs     = []
+    obs_obsvar  = []
+    obs_fcst    = []
+    obs_fcstvar = []
+    obs_ana     = []
+    obs_anavar  = []
+
+    if os.path.exists(fname):
+        print(f"reading from {fname}")
+
+        date_string = os.path.basename(fname).split('.')[-2].rstrip('z')
+        dt = datetime.strptime(date_string, '%Y%m%d_%H%M')
+        dofyr = int(dt.strftime('%j'))
+        is_leap_year = (dt.year % 4 == 0 and (dt.year % 100 != 0 or dt.year % 400 == 0))
+        pentad = (dofyr-1)//5 + 1
+        if is_leap_year and dofyr >= 59:
+            pentad = (dofyr-2)//5 + 1
+
+        date_time = {
+            'year'  : dt.year,
+            'month' : dt.month,
+            'day'   : dt.day,
+            'hour'  : dt.hour,
+            'min'   : dt.minute,
+            'sec'   : dt.second,
+            'dofyr' : dofyr,
+            'pentad': pentad
+        }
+
+        with Dataset(fname, 'r') as ncid:
+            def get_int(name):
+                return np.array(ncid.variables[name][:], dtype=int)
+
+            def get_float(name, mask_nodata=True):
+                var = ncid.variables[name]
+                data = var[:]
+                if np.ma.isMaskedArray(data):
+                    if mask_nodata:
+                        data = data.filled(np.nan)
+                    else:
+                        data = data.filled(getattr(var, '_FillValue', np.nan))
+                data = np.array(data, dtype=np.float32)
+
+                if mask_nodata:
+                    fill_value = getattr(var, '_FillValue', None)
+                    if fill_value is not None:
+                        data[data == fill_value] = np.nan
+
+                return data
+
+            obs_assim   = get_int('assim_flag')
+            obs_species = get_int('species')
+            obs_tilenum = get_int('tilenum')
+            obs_lon     = get_float('lon', mask_nodata=False)
+            obs_lat     = get_float('lat', mask_nodata=False)
+            obs_obs     = get_float('obs', mask_nodata=False)
+            obs_obsvar  = get_float('obsvar')
+            obs_fcst    = get_float('fcst')
+            obs_fcstvar = get_float('fcstvar')
+            obs_ana     = get_float('ana')
+            obs_anavar  = get_float('anavar')
+
+    else:
+        print(f"file does not exist: {fname}")
+
+    return {'date_time'  : date_time,
+            'obs_assim'  : obs_assim,
+            'obs_species': obs_species,
+            'obs_tilenum': obs_tilenum,
+            'obs_lon'    : obs_lon,
+            'obs_lat'    : obs_lat,
+            'obs_obs'    : obs_obs,
+            'obs_obsvar' : obs_obsvar,
+            'obs_fcst'   : obs_fcst,
+            'obs_fcstvar': obs_fcstvar,
+            'obs_ana'    : obs_ana,
+            'obs_anavar' : obs_anavar}
+
 # ================ EOF =================================================