add python scripts for computing stats of ObsFcstAna output (via intermediate monthly files)

GEOSldas_App/util/postproc/ObsFcstAna_stats/helper/compute_monthly_stats.py

@@ -0,0 +1,87 @@
+import numpy as np
+from datetime import datetime, timedelta
+from dateutil.relativedelta import relativedelta
+from helper.read_GEOSldas import read_ObsFcstAna, read_tilecoord, read_obs_param
+
+def compute_monthly_stats(expdir,expid,domain,this_month,tc,obs_param,var_list):
+
+    n_tile = tc['N_tile']
+    n_spec = len(obs_param)
+
+    start_time = this_month.replace(day=1,hour=3) 
+    end_time = start_time + relativedelta(months=1)
+
+    data_sum = {}
+    data2_sum = {}
+
+    N_data = np.zeros((n_tile, n_spec))
+    oxf_sum = np.zeros((n_tile, n_spec))
+    oxa_sum = np.zeros((n_tile, n_spec))
+    fxa_sum = np.zeros((n_tile, n_spec))
+
+    for var in var_list:
+        data_sum[var] = np.zeros((n_tile, n_spec))
+        data2_sum[var] = np.zeros((n_tile, n_spec))
+
+    date_time = start_time
+    while date_time < end_time:
+
+        fname =  fname = expdir+expid+'/output/'+domain+'/ana/ens_avg/Y'+ \
+                          date_time.strftime('%Y') + '/M' + \
+                          date_time.strftime('%m') + '/'  + \
+                          expid+'.ens_avg.ldas_ObsFcstAna.' + \
+                          date_time.strftime('%Y%m%d_%H%M') +'z.bin'
+
+        OFA = read_ObsFcstAna(fname)
+
+        if len(OFA['obs_tilenum'] > 0):
+            # Initialize full size variable to keep values
+            data_tile={}
+            for var in var_list:
+                data_tile[var] = np.zeros((n_tile, n_spec)) +np.nan
+
+            for ispec in np.arange(n_spec):
+                # check species overall "assim" flag for masking
+                this_species = int(obs_param[ispec]['species'])
+                masked_data = {}
+                if obs_param[ispec]['assim'] == 'T':
+                    masked_tilenum = OFA['obs_tilenum'][np.logical_and(OFA['obs_species'] == this_species, OFA['obs_assim']==1)]
+                    for var in var_list:
+                        masked_data[var] = OFA[var][np.logical_and(OFA['obs_species'] == this_species, OFA['obs_assim']==1)]
+                else:
+                    masked_tilenum = OFA['obs_tilenum'][OFA['obs_species'] == this_species]
+                    for var in var_list:
+                        masked_data[var] = OFA[var][OFA['obs_species'] == this_species]          
+
+                tile_idx = np.where(np.isin(tc['tile_id'], masked_tilenum))[0]
+
+                for var in var_list:
+                    data_tile[var][tile_idx, ispec] = masked_data[var]         
+
+            is_valid = ~np.isnan(data_tile['obs_obs'])
+            N_data[is_valid] += 1
+            oxf_sum[is_valid] += data_tile['obs_obs'][is_valid] * data_tile['obs_fcst'][is_valid]
+            oxa_sum[is_valid] += data_tile['obs_obs'][is_valid] * data_tile['obs_ana'][is_valid]
+            fxa_sum[is_valid] += data_tile['obs_fcst'][is_valid] * data_tile['obs_ana'][is_valid]
+            for var in var_list:
+                data_sum[var][is_valid] += data_tile[var][is_valid]
+                data2_sum[var][is_valid] += data_tile[var][is_valid] **2
+
+        date_time = date_time + timedelta(seconds=10800)
+
+    return N_data, data_sum, data2_sum, oxf_sum, oxa_sum, fxa_sum
+
+if __name__ == '__main__':
+    date_time = datetime(2015,5,1)
+    expdir = '/gpfsm/dnb05/projects/p51/SMAP_Nature/'
+    expid = 'SPL4SM_Vv8010'
+    domain = 'SMAP_EASEv2_M09_GLOBAL'
+    var_list = ['obs_obs', 'obs_obsvar', 'obs_fcst', 'obs_fcstvar', 'obs_ana', 'obs_anavar']
+    ftc = expdir+expid+'/output/'+domain+'/rc_out/'+expid+'.ldas_tilecoord.bin'
+    tc = read_tilecoord(ftc)
+
+    fop = expdir+expid+'/output/'+domain+'/rc_out/Y2015/M04/'+expid+'.ldas_obsparam.20150401_0000z.txt'
+    obs_param = read_obs_param(fop)
+
+    N_data, data_sum, data2_sum, oxf_sum, oxa_sum, fxa_sum = \
+           compute_monthly_stats(expdir,expid,domain,date_time,tc,obs_param,var_list)

GEOSldas_App/util/postproc/ObsFcstAna_stats/helper/plot.py

@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits import basemap
+
+def plotMap(
+    data,
+    *,
+    ax=None,
+    lat=None,
+    lon=None,
+    title=None,
+    cRange=None,
+    figsize=(8, 4),
+    clbar=True,
+    cRangeint=False,
+    cmap=plt.cm.jet,
+    bounding=None,
+    prj="cyl",
+):
+
+    if cRange is not None:
+        vmin = cRange[0]
+        vmax = cRange[1]
+    else:
+        temp = flatData(data)
+        vmin = np.percentile(temp, 5)
+        vmax = np.percentile(temp, 95)
+        if cRangeint is True:
+            vmin = int(round(vmin))
+            vmax = int(round(vmax))
+    if ax is None:
+        fig = plt.figure(figsize=figsize)
+        ax = fig.subplots()
+
+    if bounding is None:
+        bounding = [
+            np.min(lat) - 0.5,
+            np.max(lat) + 0.5,
+            np.min(lon) - 0.5,
+            np.max(lon) + 0.5,
+        ]
+
+    mm = basemap.Basemap(
+        llcrnrlat=bounding[0],
+        urcrnrlat=bounding[1],
+        llcrnrlon=bounding[2],
+        urcrnrlon=bounding[3],
+        projection=prj,
+        resolution="c",
+        ax=ax,
+    )
+    mm.drawcoastlines()
+    #mm.drawstates(linestyle="dashed")
+    #mm.drawcountries(linewidth=1.0, linestyle="-.")
+    cs = mm.pcolormesh(lon, lat, data, cmap=cmap, vmin=vmin, vmax=vmax)
+    if clbar is True:
+        cb = mm.colorbar(cs, pad="5%", location="bottom")
+        if 'normalized' in title:
+            cb.set_ticks(np.linspace(vmin,vmax,6))
+            cb.set_ticklabels([f'{10**x:.2f}' for x in np.linspace(vmin, vmax, 6)])
+
+    if title is not None:
+        ax.set_title(title)
+    if ax is None:
+        return fig, ax, mm
+    else:
+        return mm, cs
+

GEOSldas_App/util/postproc/ObsFcstAna_stats/helper/read_GEOSldas.py

@@ -0,0 +1,224 @@
+import numpy as np
+import struct
+import os
+import numpy as np
+
+def read_obs_param(fname):
+    print(f"Reading {fname}")
+
+    with open(fname, 'r') as fid:
+        N_obs_param = int(fid.readline().strip())
+
+        obs_param = []
+        for _ in range(N_obs_param):
+            param = {}
+            param['descr'] = fid.readline().strip().strip('"')
+            param['species'] = float(fid.readline().strip())
+            param['orbit'] = float(fid.readline().strip())
+            param['pol'] = float(fid.readline().strip())
+            param['N_ang'] = int(float(fid.readline().strip()))
+
+            param['ang'] = np.array([float(x) for x in fid.readline().split()])
+
+            param['freq'] = float(fid.readline().strip())
+            param['FOV'] = float(fid.readline().strip())
+            param['FOV_units'] = fid.readline().strip().strip('"')
+            param['assim'] = fid.readline().strip()
+            param['scale'] = fid.readline().strip()
+            param['getinnov'] = fid.readline().strip()
+            param['RTM_ID'] = float(fid.readline().strip())
+            param['bias_Npar'] = float(fid.readline().strip())
+            param['bias_trel'] = float(fid.readline().strip())
+            param['bias_tcut'] = float(fid.readline().strip())
+            param['nodata'] = float(fid.readline().strip())
+            param['varname'] = fid.readline().strip().strip('"')
+            param['units'] = fid.readline().strip().strip('"')
+            param['path'] = fid.readline().strip().strip('"')
+            param['name'] = fid.readline().strip().strip('"')
+
+            if 'Vv80' in fname or 'OL80' in fname:
+                fid.readline()  # Skip two lines
+                fid.readline()
+
+            param['scalepath'] = fid.readline().strip().strip('"')
+            param['scalename'] = fid.readline().strip().strip('"')
+            param['flistpath'] = fid.readline().strip().strip('"')
+            param['flistname'] = fid.readline().strip().strip('"')
+            param['errstd'] = float(fid.readline().strip())
+            param['std_normal_max'] = float(fid.readline().strip())
+            param['zeromean'] = fid.readline().strip()
+            param['coarsen_pert'] = fid.readline().strip()
+            param['xcorr'] = float(fid.readline().strip())
+            param['ycorr'] = float(fid.readline().strip())
+            param['adapt'] = float(fid.readline().strip())
+
+            obs_param.append(param)
+
+    print(f"Done reading obs_param for {N_obs_param} species")
+
+    return obs_param
+
+def read_tilecoord(fname):
+    int_precision = 'i'
+    float_precision = 'f'
+
+    # Determine endianness
+    machfmt = '<'  # '>' for big-endian, '<' for little-endian
+
+    print(f"reading from {fname}")
+
+    tile_coord = {}
+
+    with open(fname, 'rb') as ifp:
+        fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+        tile_coord['N_tile'] = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+        fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+        Nt = tile_coord['N_tile']
+
+        fields = ['tile_id', 'typ', 'pfaf', 'com_lon', 'com_lat', 'min_lon', 'max_lon',
+                      'min_lat', 'max_lat', 'i_indg', 'j_indg', 'frac_cell', 'frac_pfaf',
+                      'area', 'elev']
+
+        for field in fields:
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            dtype = int_precision if field in ['tile_id', 'typ', 'pfaf', 'i_indg', 'j_indg'] else float_precision
+            tile_coord[field] = np.frombuffer(ifp.read(Nt * 4), dtype=f'{machfmt}{dtype}')
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+    print("done reading file")
+    return tile_coord
+
+def read_ObsFcstAna(fname, isLDASsa=False):
+
+    # Initialize outputs
+    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 = []
+
+    # Determine endianness
+    machfmt = '>' if isLDASsa else '<'  # '>' for big-endian, '<' for little-endian
+
+    if os.path.exists(fname):
+        print(f"reading from {fname}")
+
+        with open(fname, 'rb') as ifp:
+            # Read N_obs and time stamp entry
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            N_obs = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            year, month, day, hour, minute, second, dofyr, pentad = struct.unpack(f'{machfmt}8i', ifp.read(32))
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            date_time = {
+                'year': year,
+                'month': month,
+                'day': day,
+                'hour': hour,
+                'min': minute,
+                'sec': second,
+                'dofyr': dofyr,
+                'pentad': pentad
+            }
+
+            # Read observation assim flag
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            tmp_data = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}i').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            obs_assim = np.zeros(N_obs, dtype=int)
+            obs_assim[tmp_data != 0] = 1
+
+            # Read species information
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_species = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}i').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read tile number information
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_tilenum = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}i').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read longitude
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_lon = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read latitude
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_lat = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read observation value
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_obs = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read observation variance
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_obsvar = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read observation-space model forecast value
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_fcst = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read observation-space model forecast variance
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_fcstvar = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read observation-space analysis value
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_ana = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # Read observation-space analysis variance
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+            obs_anavar = np.frombuffer(ifp.read(N_obs * 4), dtype=f'{machfmt}f').copy()
+            fortran_tag = struct.unpack(f'{machfmt}i', ifp.read(4))[0]
+
+            # No-data check
+            obs_obsvar[obs_obsvar == nodata] = np.nan
+            obs_fcst[obs_fcst == nodata] = np.nan
+            obs_fcstvar[obs_fcstvar == nodata] = np.nan
+            obs_ana[obs_ana == nodata] = np.nan
+            obs_anavar[obs_anavar == nodata] = np.nan
+
+    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}
+