Enhance post-processing plots

jobs/JLANDDA_PLOT_STATS

@@ -67,12 +67,16 @@ export COMIN="${COMIN:-${COMROOT}/${NET}/${model_ver}/${RUN}.${PDY}}"
 export COMOUT="${COMOUT:-${COMROOT}/${NET}/${model_ver}/${RUN}.${PDY}}"
 
 mkdir -p ${COMOUT}
+export COMOUThofx="${COMOUThofx:-${COMOUT}/hofx}"
+mkdir -p ${COMOUThofx}
 export COMOUTplot="${COMOUTplot:-${COMOUT}/plot}"
 mkdir -p ${COMOUTplot}
 
 # Create a teomporary share directory
 export DATA_HOFX="${DATA_HOFX:-${DATAROOT}/DATA_SHARE/hofx}"
 mkdir -p ${DATA_HOFX}
+export DATA_HOFX_OMA="${DATA_HOFX_OMA:-${DATAROOT}/DATA_SHARE/hofx_oma}"
+mkdir -p ${DATA_HOFX_OMA}
 #
 #-----------------------------------------------------------------------
 #

scripts/exlandda_analysis.sh

@@ -202,7 +202,7 @@ done
 
 if [ -d diags ]; then
   cp -p diags/* ${COMOUThofx}
-  ln -nsf ${COMOUThofx}/* ${DATA_HOFX}
+  ln -nsf ${COMOUThofx}/*.nc ${DATA_HOFX}
 fi
 
 

scripts/exlandda_plot_stats.sh

@@ -48,6 +48,8 @@ plottype: '${plottype}'
 title_fig: '${title_fig}'
 output_prefix: '${output_prefix}'
 cartopy_ne_path: '${FIXlandda}/NaturalEarth'
+hofx_data_path: '${DATA_HOFX_OMA}'
+cdate: '${YYYY}-${MM}-${DD}-${HH}'
 EOF
 
   ${USHlandda}/hofx_analysis_stats.py
@@ -57,6 +59,7 @@ EOF
 
   # Copy result files to COMOUT
   cp -p ${output_prefix}* ${COMOUTplot}
+  cp -p "${DATA_HOFX_OMA}/hofx_oma_timehis"* ${COMOUThofx}
 fi
 
 
@@ -82,6 +85,7 @@ nprocs_anal: '${NPROCS_ANALYSIS}'
 nprocs_fcst: '${nprocs_forecast}'
 obs_type: '${OBS_TYPE}'
 out_fn_base: '${out_fn_base}'
+hofx_data_path: '${DATA_HOFX_OMA}'
 EOF
 
   ${USHlandda}/plot_analysis_timehistory.py
@@ -95,14 +99,15 @@ fi
 
 
 ###########################################################
-# Restart Plot
+# Plot restart tiles
 ###########################################################
 if [ "${DO_PLOT_RESTART}" = "YES" ]; then
   fn_data_base="ufs_land_restart.${nYYYY}-${nMM}-${nDD}_${nHH}-00-00.tile"
   fn_data_ext=".nc"
   soil_level_number="1"
   out_title_base="Land-DA::restart::${nYYYY}-${nMM}-${nDD}_${nHH}::"
   out_fn_base="landda_out_restart_${nYYYY}-${nMM}-${nDD}_${nHH}_"
+  plot_cs_cmap="gist_ncar_r"
 
   cat > plot_restart.yaml <<EOF
 path_data: '${COMIN}/RESTART'
@@ -113,6 +118,7 @@ soil_lvl_number: '${soil_level_number}'
 out_title_base: '${out_title_base}'
 out_fn_base: '${out_fn_base}'
 cartopy_ne_path: '${FIXlandda}/NaturalEarth'
+plot_cs_cmap: '${plot_cs_cmap}'
 EOF
 
   ${USHlandda}/plot_forecast_restart.py
@@ -134,6 +140,14 @@ if [ "${DO_PLOT_COMBINE_TILES}" = "YES" ]; then
   soil_level_number="1"
   out_title_base="Land-DA::${nYYYY}-${nMM}-${nDD}_${nHH}::"
   out_fn_base="landda_out_combined_${nYYYY}-${nMM}-${nDD}_${nHH}_"
+  # Number of mesh (grid) points in longitudinal direction
+  nlon_plot=400
+  # Number of mesh (grid) points in latitudinal direction
+  nlat_plot=200
+  # SciPy griddata methods: linear, nearest, cubic
+  griddata_method="nearest"
+  # matplolib pcolormesh shading options: flat, nearest, auto, gouraud
+  shading_option="auto"
 
   cat > plot_combine_tiles.yaml <<EOF
 path_data: '${COMIN}/RESTART'
@@ -144,6 +158,10 @@ soil_lvl_number: '${soil_level_number}'
 out_title_base: '${out_title_base}'
 out_fn_base: '${out_fn_base}'
 cartopy_ne_path: '${FIXlandda}/NaturalEarth'
+nlon_plot: ${nlon_plot}
+nlat_plot: ${nlat_plot}
+griddata_method: '${griddata_method}'
+shading_option: '${shading_option}'
 EOF
 
   ${USHlandda}/plot_combine_tiles.py

ush/hofx_analysis_stats.py

@@ -66,7 +66,7 @@ def get_obs_stats(fdir, plottype, jedi_exe):
 
     return total_oma,total_omb,total_lat,total_lon
 
-def plot_scatter():
+def plot_scatter(hofx_data_path,cdate):
     print("===== PLOT: SCATTER =====")
 
     # Set the path to Natural Earth dataset
@@ -79,7 +79,23 @@ def plot_scatter():
     print("Mean |OMA|=",field_mean)
     print("STDV |OMA|=",field_std)
     print("Max |OMA|=",field_max)
-    print("Min |OMA|=",field_min)    
+    print("Min |OMA|=",field_min)
+
+    # Print out OMA values to file
+    hofx_data_fp=os.path.join(hofx_data_path,"hofx_oma_timehis_abs.txt")
+    if os.path.exists(hofx_data_fp):
+        # Remove line for same date
+        with open(hofx_data_fp, 'r') as f:
+            lines = f.readlines()
+        with open(hofx_data_fp, 'w') as f:
+            for line in lines:
+                columns = line.strip().split(' ')
+                if columns and columns[0].strip() != cdate:
+                    f.write(line)
+
+    with open(hofx_data_fp, 'a') as f:
+        print(cdate,field_mean,field_std,field_max,field_min, file=f)
+
     crs=ccrs.PlateCarree()
     fig=plt.figure(figsize=(8,5))
     ax=plt.subplot(111, projection=crs)
@@ -100,7 +116,7 @@ def plot_scatter():
     plt.savefig(output_fn,dpi=200,bbox_inches='tight')
     plt.close('all')
 
-def plot_histogram():
+def plot_histogram(hofx_data_path,cdate):
     print("===== PLOT: HISTOGRAM =====")    
     field_mean=float("{:.2f}".format(np.mean(field)))
     field_std=float("{:.2f}".format(np.std(field)))
@@ -110,8 +126,23 @@ def plot_histogram():
     print("STDV OMA=",field_std)
     print("Max OMA=",field_max)
     print("Min OMA=",field_min)
-    nbins=yaml_data['nbins']
 
+    # Print out OMA values to file
+    hofx_data_fp=os.path.join(hofx_data_path,"hofx_oma_timehis.txt")
+    if os.path.exists(hofx_data_fp):
+        # Remove line for same date
+        with open(hofx_data_fp, 'r') as f:
+            lines = f.readlines()
+        with open(hofx_data_fp, 'w') as f:
+            for line in lines:
+                columns = line.strip().split(' ')
+                if columns and columns[0].strip() != cdate:
+                    f.write(line)
+
+    with open(hofx_data_fp, 'a') as f:
+        print(cdate,field_mean,field_std,field_max,field_min, file=f)
+
+    nbins=yaml_data['nbins']
     opt_xlimit='auto'
     if opt_xlimit=='auto':
         fld_min=int(field_min)
@@ -138,6 +169,8 @@ def plot_histogram():
         yaml_data=yaml.load(f, Loader=yaml.FullLoader)
     f.close()
     print("YAML_DATA:",yaml_data)
+    hofx_data_path=yaml_data['hofx_data_path']
+    cdate=yaml_data['cdate']
 
     oma,omb,lat,lon=get_obs_stats(yaml_data['hofx_files'],yaml_data['plottype'],yaml_data['jedi_exe'])
     if yaml_data['field_var']=='OMA':
@@ -146,6 +179,6 @@ def plot_histogram():
         field=omb    
 
     if yaml_data['plottype']=='scatter' or yaml_data['plottype']=='both': 
-        plot_scatter()
+        plot_scatter(hofx_data_path,cdate)
     if yaml_data['plottype']=='histogram' or yaml_data['plottype']=='both':
-        plot_histogram()
+        plot_histogram(hofx_data_path,cdate)

ush/plot_analysis_timehistory.py

@@ -9,6 +9,7 @@
 ## V000: 2024/10/14: Chan-Hoo Jeon : Preliminary version
 ## V001: 2024/10/15: Chan-Hoo Jeon : Add wall-clock time plot
 ## V002: 2024/10/31: Chan-Hoo Jeon : Fix input log file name issue
+## V003: 2025/02/26: Chan-Hoo Jeon : Add h(x) Obs-ana plot
 ###################################################################### CHJ #####
 
 import os, sys
@@ -43,6 +44,7 @@ def main():
     nprocs_fcst = yaml_data['nprocs_fcst']
     obs_type = yaml_data['obs_type']
     out_fn_base = yaml_data['out_fn_base']
+    hofx_data_path = yaml_data['hofx_data_path']
 
     var_list = ["totalSnowDepth"]
     nprocs_anal = int(nprocs_anal)
@@ -53,6 +55,7 @@ def main():
         var_dict_anal = get_data_analysis(path_data,fn_data_anal_prefix,fn_data_anal_suffix,jedi_exe,nprocs_anal,var_nm)
         var_dict_fcst = get_data_forecast(path_data,fn_data_fcst_prefix,fn_data_fcst_suffix,nprocs_fcst)
         plot_data(var_dict_anal,var_dict_fcst,jedi_exe,obs_type,out_fn_base,work_dir,var_nm)
+        plot_his_oma(var_dict_anal,out_fn_base,work_dir,var_nm,hofx_data_path)
 
 
 # Get data from files =============================================== CHJ =====
@@ -342,6 +345,61 @@ def plot_his_qc(dfa,min_var,max_var,rms_var,out_title_qc,out_fn_qc,work_dir,qc_t
     out_file(work_dir,out_fn_qc,ndpi)
 
 
+# Plot time-history of H(x) OMA data ================================ CHJ =====
+def plot_his_oma(var_dict_anal,out_fn_base,work_dir,var_nm,hofx_data_path):
+# =================================================================== CHJ =====
+
+    dfa = pd.DataFrame(var_dict_anal)
+
+    oma_fp = os.path.join(hofx_data_path,"hofx_oma_timehis.txt")
+    with open(oma_fp, 'r') as f:
+        lines = f.readlines()
+    column_data = [line.strip().split(' ') for line in lines]
+    num_columns = len(column_data[0]) if column_data else 0
+    columns = [[] for _ in range(num_columns)]
+    for row in column_data:
+        for i, value in enumerate(row):
+            if i>0:
+                value = float(value)
+            columns[i].append(value)
+
+    print(columns[1])
+
+    out_title_oma = f'''Land-DA::OMA (observation-analysis)::{var_nm}'''
+    out_fn_oma = f'''{out_fn_base}_oma_{var_nm}'''
+
+    # figsize=(width,height) in inches
+    fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(6,6))
+    fig.suptitle(out_title_oma,fontsize=txt_fnt+1,y=0.97)
+
+    axes[0].plot(dfa['Date'],columns[1],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz,label='Mean')
+    axes[0].set_ylabel('OMA: Mean', fontsize=txt_fnt-1)
+    axes[0].tick_params(axis="y",labelsize=txt_fnt-2)
+#    axes[0].legend(fontsize=txt_fnt-1, loc='center')
+    axes[0].grid(linewidth=0.2)
+
+    axes[1].plot(dfa['Date'],columns[2],'s-.',color='red',mfc='none',linewidth=ln_wdth,markersize=mk_sz,label='STD')
+    axes[1].set_ylabel('OMA: STDV', fontsize=txt_fnt-1)
+    axes[1].tick_params(axis="y",labelsize=txt_fnt-2)
+#    axes[1].legend(fontsize=txt_fnt-1, loc='center')
+    axes[1].grid(linewidth=0.2)
+
+    axes[2].plot(dfa['Date'],dfa['nobs_in'],'o-',color='blue',linewidth=ln_wdth,markersize=mk_sz,label='N_obs:raw')
+    axes[2].plot(dfa['Date'],dfa['nobs_QC'],'s-.',color='red',mfc='none',linewidth=ln_wdth,markersize=mk_sz,label='N_obs:QC')
+    axes[2].set_xlabel('Date', fontsize=txt_fnt-1)
+    axes[2].set_ylabel('Number of observations', fontsize=txt_fnt-1)
+    axes[2].tick_params(axis="x",labelsize=txt_fnt-2)
+    axes[2].tick_params(axis="y",labelsize=txt_fnt-2)
+    axes[2].legend(fontsize=txt_fnt-1, loc='center right')
+    axes[2].grid(linewidth=0.2)
+
+    plt.xticks(rotation=30, ha='right')
+    plt.tight_layout()
+    # Output figure
+    ndpi = 300
+    out_file(work_dir,out_fn_oma,ndpi)
+
+
 # Output file ======================================================= CHJ =====
 def out_file(work_dir,out_file,ndpi):
 # =================================================================== CHJ =====

ush/plot_combine_tiles.py

@@ -45,6 +45,12 @@ def main():
     out_title_base=yaml_data['out_title_base']
     out_fn_base=yaml_data['out_fn_base']
     cartopy_ne_path=yaml_data['cartopy_ne_path']
+    nlon_plot=yaml_data['nlon_plot']
+    nlat_plot=yaml_data['nlat_plot']
+    nlon_plot=int(nlon_plot)
+    nlat_plot=int(nlat_plot)
+    griddata_method=yaml_data['griddata_method']
+    shading_option=yaml_data['shading_option']
 
     # Set the path to Natural Earth dataset
     cartopy.config['data_dir']=cartopy_ne_path
@@ -59,7 +65,8 @@ def main():
     # plot restart file
     for var_nm in var_list:
         sfc_data = get_data(path_data,fn_data_base,fn_data_ext,var_nm,soil_lvl_num)
-        plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir)
+        plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir,
+                  nlon_plot,nlat_plot,griddata_method,shading_option)
 
 
 # geo lon/lat from orography ======================================== CHJ =====
@@ -149,7 +156,8 @@ def get_data(path_data,fn_data_base,fn_data_ext,var_nm,soil_lvl_num):
 
 
 # Plot data ========================================================= CHJ =====
-def plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir):
+def plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir,
+              nlon_plot,nlat_plot,griddata_method,shading_option):
 # =================================================================== CHJ =====
 
     print(' ===== Plotting data ================================')
@@ -159,8 +167,8 @@ def plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir):
     print('Fill value:',sfc_fill_value)
 
     # Define and interpolate the grid/data for mesh plot
-    num_glon_mesh=200
-    num_glat_mesh=100
+    num_glon_mesh=nlon_plot
+    num_glat_mesh=nlat_plot
     lon_min=np.min(glon)
     lon_max=np.max(glon)
     lat_min=round(np.min(glat))
@@ -177,7 +185,7 @@ def plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir):
     print('Min glat mesh:', np.min(glat_m))
     print(glon_m.shape)
     print(glat_m.shape)
-    sfc_data_m = griddata((glon,glat),sfc_data,(glon_m,glat_m),method='nearest',fill_value=sfc_fill_value)
+    sfc_data_m = griddata((glon,glat),sfc_data,(glon_m,glat_m),method=griddata_method,fill_value=sfc_fill_value)
     sfc_data_m_masked = np.ma.masked_where(sfc_data_m == sfc_fill_value, sfc_data_m)
     print(sfc_data_m_masked.shape)
 
@@ -201,7 +209,7 @@ def plot_data(glon,glat,sfc_data,var_nm,out_title_base,out_fn_base,work_dir):
     back_plot(ax)
 
     cs=ax.pcolormesh(glon_m,glat_m,sfc_data_m_masked,cmap=cs_cmap,rasterized=True,
-       vmin=cs_min,vmax=cs_max,transform=ccrs.PlateCarree())
+       shading=shading_option,vmin=cs_min,vmax=cs_max,transform=ccrs.PlateCarree())
 
     divider=make_axes_locatable(ax)
     ax_cb=divider.new_horizontal(size="3%",pad=0.1,axes_class=plt.Axes)