diff --git a/jobs/JLANDDA_PREP_OBS b/jobs/JLANDDA_PREP_OBS
index 3e6bb7ec..b456670c 100755
--- a/jobs/JLANDDA_PREP_OBS
+++ b/jobs/JLANDDA_PREP_OBS
@@ -73,6 +73,8 @@ fi
 mkdir -p ${COMOUT}
 export COMOUTobs=${COMOUTobs:-${COMOUT}/obs}
 mkdir -p ${COMOUTobs}
+export COMOUTplot="${COMOUTplot:-${COMOUT}/plot}"
+mkdir -p ${COMOUTplot}
 
 # Create a teomporary share directory
 export DATA_SHARE="${DATA_SHARE:-${DATAROOT}/DATA_SHARE/${PDY}}"
diff --git a/modulefiles/tasks/hera/task.prep_obs.lua b/modulefiles/tasks/hera/task.prep_obs.lua
index e7bcb54c..d57971da 100644
--- a/modulefiles/tasks/hera/task.prep_obs.lua
+++ b/modulefiles/tasks/hera/task.prep_obs.lua
@@ -10,7 +10,5 @@ load(pathJoin("py-cartopy", py_cartopy_ver))
 load(pathJoin("py-matplotlib", py_matplotlib_ver))
 load(pathJoin("py-netcdf4", py_netcdf4_ver))
 load(pathJoin("py-numpy", py_numpy_ver))
-load(pathJoin("py-pandas", py_pandas_ver))
 load(pathJoin("py-pyyaml", py_pyyaml_ver))
-load(pathJoin("py-xarray", py_xarray_ver))
 
diff --git a/scripts/exlandda_plot_stats.sh b/scripts/exlandda_plot_stats.sh
index 93ac7c14..69179a72 100755
--- a/scripts/exlandda_plot_stats.sh
+++ b/scripts/exlandda_plot_stats.sh
@@ -19,30 +19,28 @@ nHH=${NTIME:8:2}
 # Stats Plot
 ############################################################
 
-# Path to the directory containing the input file
-INPUTFP="${DATA_HOFX}"
 # Field variable
-FIELDVAR="OMA"
+field_var="OMA"
 # Field Range for scatter plot: [Low,High]
-FRLOW=-300
-FRHIGH=300
+field_range_low=-300
+field_range_high=300
 # Number of bins in histogram plot
-NBINS=100
+nbins=100
 # Plot type (scatter/histogram/both)
-PLOTTYPE="both"
+plottype="both"
 # Figure title
-FIGTITLE="GHCN Snow Depth (mm)::Obs-Ana::${PDY}"
+title_fig="GHCN Snow Depth (mm)::Obs-Ana::${PDY}"
 # Prefix of output file name
-PREOUTFN="hofx_oma_${PDY}"
+output_prefix="hofx_oma_${PDY}"
 
 cat > plot_hofx.yaml <<EOF
-hofx_files: '${INPUTFP}'
-field_var: '${FIELDVAR}'
-field_range: [${FRLOW},${FRHIGH}]
-nbins: ${NBINS}
-plottype: '${PLOTTYPE}'
-title_fig: '${FIGTITLE}'
-output_prefix: '${PREOUTFN}'
+hofx_files: '${DATA_HOFX}'
+field_var: '${field_var}'
+field_range: [${field_range_low},${field_range_high}]
+nbins: ${nbins}
+plottype: '${plottype}'
+title_fig: '${title_fig}'
+output_prefix: '${output_prefix}'
 machine: '${MACHINE}'
 EOF
 
@@ -52,35 +50,35 @@ if [ $? -ne 0 ]; then
 fi
 
 # Copy result files to COMOUT
-cp -p ${PREOUTFN}* ${COMOUTplot}
+cp -p ${output_prefix}* ${COMOUTplot}
 
 
 ############################################################
 # Time-history Plot
 ############################################################
 
-FN_DATA_ANAL_PREFIX="analysis_"
-FN_DATA_ANAL_SUFFIX=".log"
-FN_DATA_FCST_PREFIX="forecast_"
-FN_DATA_FCST_SUFFIX=".log"
-OUT_TITLE_ANAL_BASE="Land-DA::Analysis::QC SnowDepthGHCN::"
-OUT_FN_ANAL_BASE="landda_timehistory_"
-OUT_TITLE_TIME="Land-DA::Wall-clock time"
-OUT_FN_TIME="landda_timehistory_wtime"
+fn_data_anal_prefix="analysis_"
+fn_data_anal_suffix=".log"
+fn_data_fcst_prefix="forecast_"
+fn_data_fcst_suffix=".log"
+out_title_anal_base="Land-DA::Analysis::QC SnowDepthGHCN::"
+out_fn_anal_base="landda_timehistory_"
+out_title_time="Land-DA::Wall-clock time"
+out_fn_time="landda_timehistory_wtime"
 
 cat > plot_timehistory.yaml <<EOF
 path_data: '${LOGDIR}'
 work_dir: '${DATA}'
-fn_data_anal_prefix: '${FN_DATA_ANAL_PREFIX}'
-fn_data_anal_suffix: '${FN_DATA_ANAL_SUFFIX}'
-fn_data_fcst_prefix: '${FN_DATA_FCST_PREFIX}'
-fn_data_fcst_suffix: '${FN_DATA_FCST_SUFFIX}'
+fn_data_anal_prefix: '${fn_data_anal_prefix}'
+fn_data_anal_suffix: '${fn_data_anal_suffix}'
+fn_data_fcst_prefix: '${fn_data_fcst_prefix}'
+fn_data_fcst_suffix: '${fn_data_fcst_suffix}'
 nprocs_anal: '${NPROCS_ANALYSIS}'
 nprocs_fcst: '${NPROCS_FORECAST}'
-out_title_anal_base: '${OUT_TITLE_ANAL_BASE}'
-out_fn_anal_base: '${OUT_FN_ANAL_BASE}'
-out_title_time: '${OUT_TITLE_TIME}'
-out_fn_time: '${OUT_FN_TIME}'
+out_title_anal_base: '${out_title_anal_base}'
+out_fn_anal_base: '${out_fn_anal_base}'
+out_title_time: '${out_title_time}'
+out_fn_time: '${out_fn_time}'
 EOF
 
 ${USHlandda}/plot_analysis_timehistory.py
@@ -89,28 +87,28 @@ if [ $? -ne 0 ]; then
 fi
 
 # Copy result files to COMOUT
-cp -p ${OUT_FN_ANAL_BASE}* ${COMOUTplot}
-cp -p ${OUT_FN_TIME}* ${COMOUTplot}
+cp -p ${out_fn_anal_base}* ${COMOUTplot}
+cp -p ${out_fn_time}* ${COMOUTplot}
 
 
 ############################################################
 # Restart Plot
 ############################################################
 
-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}_"
+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}_"
 
 cat > plot_restart.yaml <<EOF
 path_data: '${COMIN}/RESTART'
 work_dir: '${DATA}'
-fn_data_base: '${FN_DATA_BASE}'
-fn_data_ext: '${FN_DATA_EXT}'
-soil_lvl_number: '${SOIL_LEVEL_NUMBER}'
-out_title_base: '${OUT_TITLE_BASE}'
-out_fn_base: '${OUT_FN_BASE}'
+fn_data_base: '${fn_data_base}'
+fn_data_ext: '${fn_data_ext}'
+soil_lvl_number: '${soil_level_number}'
+out_title_base: '${out_title_base}'
+out_fn_base: '${out_fn_base}'
 machine: '${MACHINE}'
 EOF
 
@@ -120,5 +118,5 @@ if [ $? -ne 0 ]; then
 fi
 
 # Copy result files to COMOUT
-cp -p ${OUT_FN_BASE}* ${COMOUTplot}
+cp -p ${out_fn_base}* ${COMOUTplot}
 
diff --git a/scripts/exlandda_prep_obs.sh b/scripts/exlandda_prep_obs.sh
index 1bc770a5..a8272e9c 100755
--- a/scripts/exlandda_prep_obs.sh
+++ b/scripts/exlandda_prep_obs.sh
@@ -28,6 +28,7 @@ if [ "${OBS_GHCN}" = "YES" ]; then
   # check obs is available
   if [ -f "${obs_fp}" ]; then
     echo "GHCN observation file: ${obs_fp}"
+    cp -p "${obs_fp}" .
     cp -p "${obs_fp}" "${COMOUTobs}/${out_fn}"
   else
     input_ghcn_file="${DATA_GHCN_RAW}/${YYYP}.csv"
@@ -36,14 +37,36 @@ if [ "${OBS_GHCN}" = "YES" ]; then
       echo "GHCN raw data file: ${YYYP}.csv"
       err_exit "GHCN raw data file does not exist in designated path !!!"
     fi
-    output_ioda_file="${obs_fn}"
     ghcn_station_file="${DATA_GHCN_RAW}/ghcnd-stations.txt"
 
-    ${USHlandda}/ghcn_snod2ioda.py -i ${input_ghcn_file} -o ${output_ioda_file} -f ${ghcn_station_file} -d ${YYYP}${MP}${DP}${HP} -m maskout
+    ${USHlandda}/ghcn_snod2ioda.py -i ${input_ghcn_file} -o ${obs_fn} -f ${ghcn_station_file} -d ${YYYP}${MP}${DP}${HP} -m maskout
     if [ $? -ne 0 ]; then
       err_exit "Generation of GHCN obs file failed !!!"
     fi
-    cp -p "${output_ioda_file}" "${COMOUTobs}/${out_fn}"
+    cp -p "${obs_fn}" "${COMOUTobs}/${out_fn}"
+  fi
+
+  ############################################################
+  # Observation File Plot
+  ############################################################
+
+  out_title_base="Land-DA::Obs::GHCN::${PDY}::"
+  out_fn_base="landda_obs_ghcn_${PDY}_"
 
+  cat > plot_obs_ghcn.yaml <<EOF
+work_dir: '${DATA}'
+fn_input: '${obs_fn}'
+out_title_base: '${out_title_base}'
+out_fn_base: '${out_fn_base}'
+machine: '${MACHINE}'
+EOF
+
+  ${USHlandda}/plot_obs_ghcn.py
+  if [ $? -ne 0 ]; then
+    err_exit "Observation file plot failed"
   fi
+
+  # Copy result file to COMOUT
+  cp -p ${out_fn_base}* ${COMOUTplot}
+
 fi
diff --git a/ush/plot_obs_ghcn.py b/ush/plot_obs_ghcn.py
new file mode 100755
index 00000000..f7587d6c
--- /dev/null
+++ b/ush/plot_obs_ghcn.py
@@ -0,0 +1,212 @@
+#!/usr/bin/env python3
+
+###################################################################### CHJ #####
+## Name		  : plot_obs_ghcn.py
+## Usage	  : Plot GHCN observation ioda file of land-DA workflow
+## Input file : ghcn_snwd_ioda.nc
+## NOAA/EPIC
+## History ===============================
+## V000: 2024/12/03: Chan-Hoo Jeon : Preliminary version
+###################################################################### CHJ #####
+
+import os, sys
+import yaml
+import numpy as np
+import netCDF4 as nc
+import cartopy
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+
+# Main part (will be called at the end) ============================= CHJ =====
+def main():
+# =================================================================== CHJ =====
+
+    yaml_file="plot_obs_ghcn.yaml"
+    with open(yaml_file, 'r') as f:
+        yaml_data=yaml.load(f, Loader=yaml.FullLoader)
+    f.close()
+    print("YAML_DATA:",yaml_data)
+
+    work_dir=yaml_data['work_dir']
+    fn_input=yaml_data['fn_input']
+    out_title_base=yaml_data['out_title_base']
+    out_fn_base=yaml_data['out_fn_base']
+    machine=yaml_data['machine']
+
+    if yaml_data['machine']=='hera':
+        cartopy.config['data_dir']='/scratch2/NAGAPE/epic/UFS_Land-DA_Dev/inputs/NaturalEarth'
+    elif yaml_data['machine']=='orion' or yaml_data['machine']=='hercules':
+        cartopy.config['data_dir']='/work/noaa/epic/UFS_Land-DA_Dev/inputs/NaturalEarth'
+    elif yaml_data['machine']=='singularity':
+        cartopy.config['data_dir']='SINGULARITY_WORKING_DIR/land-DA_workflow/fix/NaturalEarth'
+
+    print(' ===== INPUT: '+fn_input+' ================================')
+    # open the data file
+    fpath=os.path.join(work_dir,fn_input)
+    try: mdat=nc.Dataset(fpath)
+    except: raise Exception('Could NOT find the file',fpath)
+    print(mdat)
+    print(mdat.groups['MetaData'])
+    print(mdat.groups['ObsError'])
+    print(mdat.groups['ObsValue'])
+    print(mdat.groups['PreQC'])
+
+    longitude=mdat.groups['MetaData'].variables['longitude'][:]
+    latitude=mdat.groups['MetaData'].variables['latitude'][:]
+    datetime=mdat.groups['MetaData'].variables['dateTime'][:]
+    stationElevation=mdat.groups['MetaData'].variables['stationElevation'][:]
+#    height=mdat.groups['MetaData'].variables['height'][:]
+    stationID=mdat.groups['MetaData'].variables['stationIdentification'][:]
+
+    # Longitude 0:360 => -180:180
+#    lon_max=np.max(lon)
+#    if lon_max>180:
+#        lon=(lon+180)%360-180
+
+    lon=longitude
+    lat=latitude
+
+    # Highest and lowest longitudes and latitudes for plot extent
+    lon_min=np.min(lon)
+    lon_max=np.max(lon)
+    lat_min=np.min(lat)
+    lat_max=np.max(lat)
+#    extent=[lon_min,lon_max,lat_min,lat_max]
+    # for CONUS
+#    extent=[-125,-66,23,53]
+    # for Northern Hemisphere
+    extent=[-179,179,0,82.5]
+    print(extent)
+
+#    c_lon=np.mean(extent[:2])
+    c_lon=-77.0369 # D.C.
+    print(' c_lon=',c_lon)
+
+    # Variables
+#    vars_out=["ObsValue","ObsError","PreQC"]
+    vars_out=["ObsValue"]
+    for svar in vars_out:
+        svar_plot(svar,mdat,lon,lat,extent,c_lon,out_title_base,out_fn_base,work_dir)
+    
+
+# Variable plot =============================================== CHJ =====
+def svar_plot(svar,mdat,lon,lat,extent,c_lon,out_title_base,out_fn_base,work_dir):
+# ============================================================= CHJ =====
+
+    print(' ===== '+svar+' === Total Snow Depth =====================')
+    # Extract data array
+    sfld=mdat.groups[svar].variables['totalSnowDepth'][:]
+
+    svar="SnowDepth"
+    out_title_fld=out_title_base+svar
+    out_fn=out_fn_base+svar
+
+    cs_cmap='gist_ncar_r'
+    lb_ext='neither'
+    tick_ln=1.5
+    tick_wd=0.45
+    tlb_sz=3
+    scat_sz=1.0
+    n_rnd=2
+    cmap_range='fixed'
+
+    print(' svar name=',svar)
+
+    # Max and Min of the field
+    fmax=np.max(sfld)
+    fmin=np.min(sfld)
+    print(' fld_max=',fmax)
+    print(' flx_min=',fmin)
+
+    # Make the colormap range symmetry
+    print(' cmap range=',cmap_range)
+    if cmap_range=='symmetry':
+        tmp_cmp=max(abs(fmax),abs(fmin))
+        cs_min=round(-tmp_cmp,n_rnd)
+        cs_max=round(tmp_cmp,n_rnd)
+    elif cmap_range=='round':
+        cs_min=round(fmin,n_rnd)
+        cs_max=round(fmax,n_rnd)
+    elif cmap_range=='real':
+        cs_min=fmin
+        cs_max=fmax
+    elif cmap_range=='fixed':
+        cs_min=0
+        cs_max=1000.0
+    else:
+        sys.exit('ERROR: wrong colormap-range flag !!!')
+
+    print(' cs_max=',cs_max)
+    print(' cs_min=',cs_min)
+    print(' extent=',extent)
+
+    # Plot field
+    fig,ax=plt.subplots(1,1,subplot_kw=dict(projection=ccrs.Robinson(c_lon)))
+    ax.set_extent(extent, ccrs.PlateCarree())
+    # Call background plot
+    back_plot(ax)
+    ax.set_title(out_title_fld,fontsize=9)
+    cs=ax.scatter(lon,lat,transform=ccrs.PlateCarree(),c=sfld,cmap=cs_cmap,
+                  vmin=cs_min,vmax=cs_max,s=scat_sz)
+    divider=make_axes_locatable(ax)
+    ax_cb=divider.new_horizontal(size="3%",pad=0.1,axes_class=plt.Axes)
+    fig.add_axes(ax_cb)
+    cbar=plt.colorbar(cs,cax=ax_cb,extend=lb_ext)
+    cbar.ax.tick_params(labelsize=8)
+    cbar.set_label(svar,fontsize=8)
+
+    # Output figure
+    ndpi=300
+    out_file(work_dir,out_fn,ndpi)
+
+
+# Background plot ==================================================== CHJ =====
+def back_plot(ax):
+# ==================================================================== CHJ =====
+
+    # Resolution of background natural earth data ('50m' or '110m')
+    back_res='50m'
+
+    fline_wd=0.5  # line width
+    falpha=0.7 # transparency
+
+    # natural_earth
+    land=cfeature.NaturalEarthFeature('physical','land',back_res,
+                      edgecolor='face',facecolor=cfeature.COLORS['land'],
+                      alpha=falpha)
+    lakes=cfeature.NaturalEarthFeature('physical','lakes',back_res,
+                      edgecolor='blue',facecolor='none',
+                      linewidth=fline_wd,alpha=falpha)
+    coastline=cfeature.NaturalEarthFeature('physical','coastline',
+                      back_res,edgecolor='black',facecolor='none',
+                      linewidth=fline_wd,alpha=falpha)
+    states=cfeature.NaturalEarthFeature('cultural','admin_1_states_provinces',
+                      back_res,edgecolor='green',facecolor='none',
+                      linewidth=fline_wd,linestyle=':',alpha=falpha)
+    borders=cfeature.NaturalEarthFeature('cultural','admin_0_countries',
+                      back_res,edgecolor='red',facecolor='none',
+                      linewidth=fline_wd,alpha=falpha)
+
+#    ax.add_feature(land)
+#    ax.add_feature(lakes)
+#    ax.add_feature(states)
+#    ax.add_feature(borders)
+    ax.add_feature(coastline)
+
+
+# Output file ======================================================= CHJ =====
+def out_file(work_dir,out_file,ndpi):
+# =================================================================== CHJ =====
+    # Output figure
+    fp_out=os.path.join(work_dir,out_file)
+    plt.savefig(fp_out+'.png',dpi=ndpi,bbox_inches='tight')
+    plt.close('all')
+
+
+# Main call ========================================================= CHJ =====
+if __name__=='__main__':
+    main()
+