diff --git a/CHANGELOG.md b/CHANGELOG.md index 2b78a700..ffae5bb9 100644 --- a/CHANGELOG.md +++ b/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. ----------------------------- - diff --git a/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py b/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py index 00058e06..c713d80c 100644 --- a/GEOSldas_App/util/postproc/ObsFcstAna_stats/postproc_ObsFcstAna.py +++ b/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]) diff --git a/GEOSldas_App/util/shared/python/read_GEOSldas.py b/GEOSldas_App/util/shared/python/read_GEOSldas.py index ccc1aa18..81ac4cf1 100644 --- a/GEOSldas_App/util/shared/python/read_GEOSldas.py +++ b/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 =================================================