Enable the manual operation of the marine verification tool

scripts/exgdas_global_marine_analysis_vrfy.py → scripts/old/exgdas_global_marine_analysis_vrfy.py

@@ -1,4 +1,3 @@
-#!/usr/bin/env python3
 ################################################################################
 #  UNIX Script Documentation Block
 #                      .                                             .

ush/eva/marine_eva_post.py

@@ -12,7 +12,9 @@
 vminmax = {'seaSurfaceTemperature': {'vmin': -2.0, 'vmax': 2.0},
            'seaIceFraction': {'vmin': -0.2, 'vmax': 0.2},
            'seaSurfaceSalinity': {'vmin': -0.2, 'vmax': 0.2},  # TODO: this should be changed
-           'absoluteDynamicTopography': {'vmin': -0.2, 'vmax': 0.2}}
+           'absoluteDynamicTopography': {'vmin': -0.2, 'vmax': 0.2},
+           'waterTemperature': {'vmin': -2.0, 'vmax': 2.0},
+           'salinity': {'vmin': -0.2, 'vmax': 0.2}}
 
 
 def marine_eva_post(inputyaml, outputdir, diagdir):

ush/eva/marine_gdas_plots.yaml

@@ -73,7 +73,7 @@ graphics:
         data variable: experiment::OmBQC::${variable}
     figure:
       layout: [1,1]
-      figure size: [11,5]
+      figure size: [20,10]
       title: 'OmB post QC | @NAME@ @CYCLE@ | ${variable_title}'
       output name: map_plots/@NAME@/${variable}/@CHANNELVAR@/@NAME@_${variable}@[email protected]
       tight_layout: true
@@ -94,11 +94,11 @@ graphics:
           data:
             variable: experiment::OmBQC::${variable}
             @CHANNELKEY@
-          markersize: 1
+          markersize: 0.01
           label: '$(variable)'
           colorbar: true
           # below may need to be edited/removed
-          cmap: ${dynamic_cmap}
+          cmap: 'seismic'
           vmin: ${dynamic_vmin}
           vmax: ${dynamic_vmax}
 

ush/soca/soca_vrfy.py

@@ -99,17 +99,17 @@ def plotHorizontalSlice(config):
 
     fig, ax = plt.subplots(figsize=(8, 5), subplot_kw={'projection': projs[config['proj']]})
 
-    # Plot the filled contours
-    contourf_plot = ax.contourf(np.squeeze(grid.lon),
+    # Use pcolor to plot the data
+    pcolor_plot = ax.pcolormesh(np.squeeze(grid.lon),
                                 np.squeeze(grid.lat),
                                 slice_data,
-                                levels=100,
                                 vmin=bounds[0], vmax=bounds[1],
                                 transform=ccrs.PlateCarree(),
-                                cmap=config['colormap'])
+                                cmap=config['colormap'],
+                                zorder=0)
 
     # Add colorbar for filled contours
-    cbar = fig.colorbar(contourf_plot, ax=ax, shrink=0.75, orientation='horizontal')
+    cbar = fig.colorbar(pcolor_plot, ax=ax, shrink=0.75, orientation='horizontal')
     cbar.set_label(label_colorbar)
 
     # Add contour lines with specified linewidths
@@ -120,16 +120,20 @@ def plotHorizontalSlice(config):
                levels=contour_levels,
                colors='black',
                linewidths=0.1,
-               transform=ccrs.PlateCarree())
+               transform=ccrs.PlateCarree(),
+               zorder=2)
 
-    ax.coastlines()  # TODO: make this work on hpc
+    try:
+        ax.coastlines()  # TODO: make this work on hpc
+    except Exception as e:
+        print(f"Warning: could not add coastlines. {e}")
     ax.set_title(title)
     if config['proj'] == 'South':
         ax.set_extent([-180, 180, -90, -50], ccrs.PlateCarree())
     if config['proj'] == 'North':
         ax.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
     # ax.add_feature(cartopy.feature.LAND)  # TODO: make this work on hpc
-    plt.savefig(figname, bbox_inches='tight', dpi=600)
+    plt.savefig(figname, bbox_inches='tight', dpi=300)
     plt.close(fig)
 
 
@@ -184,7 +188,7 @@ def plotZonalSlice(config):
     os.makedirs(dirname, exist_ok=True)
     figname = os.path.join(dirname, config['variable'] +
                            'zonal_lat_' + str(int(lat)) + '_' + str(int(config['max depth'])) + 'm')
-    plt.savefig(figname, bbox_inches='tight', dpi=600)
+    plt.savefig(figname, bbox_inches='tight', dpi=300)
     plt.close(fig)
 
 
@@ -239,7 +243,7 @@ def plotMeridionalSlice(config):
     os.makedirs(dirname, exist_ok=True)
     figname = os.path.join(dirname, config['variable'] +
                            'meridional_lon_' + str(int(lon)) + '_' + str(int(config['max depth'])) + 'm')
-    plt.savefig(figname, bbox_inches='tight', dpi=600)
+    plt.savefig(figname, bbox_inches='tight', dpi=300)
     plt.close(fig)
 
 

utils/soca/fig_gallery/exgdas_global_marine_analysis_vrfy_manual.py

@@ -0,0 +1,210 @@
+import os
+import numpy as np
+import gen_eva_obs_yaml
+import marine_eva_post
+import diag_statistics
+from multiprocessing import Process
+from soca_vrfy import statePlotter, plotConfig
+import subprocess
+
+comout = os.getenv('COM_OCEAN_ANALYSIS')
+com_ice_history = os.getenv('COM_ICE_HISTORY_PREV')
+com_ocean_history = os.getenv('COM_OCEAN_HISTORY_PREV')
+cyc = os.getenv('cyc')
+RUN = os.getenv('RUN')
+
+bcyc = str((int(cyc) - 3) % 24).zfill(2)
+gcyc = str((int(cyc) - 6) % 24).zfill(2)
+grid_file = os.path.join(comout, f'{RUN}.t'+bcyc+'z.ocngrid.nc')
+layer_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocninc.nc')
+
+# for eva
+diagdir = os.path.join(comout, 'diags')
+HOMEgfs = os.getenv('HOMEgfs')
+
+# Get flags from environment variables (set in the bash driver)
+run_ensemble_analysis = os.getenv('RUN_ENSENBLE_ANALYSIS', 'OFF').upper() == 'ON'
+run_bkgerr_analysis = os.getenv('RUN_BACKGROUND_ERROR_ANALYSIS', 'OFF').upper() == 'ON'
+run_bkg_analysis = os.getenv('RUN_BACKGROUND_ANALYSIS', 'OFF').upper() == 'ON'
+run_increment_analysis = os.getenv('RUN_INCREMENT_ANLYSIS', 'OFF').upper() == 'ON'
+
+# Initialize an empty list for the main config
+configs = [plotConfig(grid_file=grid_file,
+                      data_file=os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc'),
+                      variables_horiz={'ave_ssh': [-1.8, 1.3],
+                                       'Temp': [-1.8, 34.0],
+                                       'Salt': [32, 40]},
+                      colormap='nipy_spectral',
+                      comout=os.path.join(comout, 'vrfy', 'ana')),   # ocean surface analysis
+           plotConfig(grid_file=grid_file,
+                      data_file=os.path.join(comout, f'{RUN}.t'+cyc+'z.iceana.nc'),
+                      variables_horiz={'aice_h': [0.0, 1.0],
+                                       'hi_h': [0.0, 4.0],
+                                       'hs_h': [0.0, 0.5]},
+                      colormap='jet',
+                      projs=['North', 'South', 'Global'],
+                      comout=os.path.join(comout, 'vrfy', 'ana'))]   # sea ice analysis
+
+# Define each config and add to main_config if its flag is True
+if run_ensemble_analysis:
+    config_ens = [plotConfig(grid_file=grid_file,
+                             data_file=os.path.join(comout, f'{RUN}.t{cyc}z.ocn.recentering_error.nc'),
+                             variables_horiz={'ave_ssh': [-1, 1]},
+                             colormap='seismic',
+                             comout=os.path.join(comout, 'vrfy', 'recentering_error')),   # recentering error
+                  plotConfig(grid_file=grid_file,
+                             data_file=os.path.join(comout, f'{RUN}.t{cyc}z.ocn.ssh_steric_stddev.nc'),
+                             variables_horiz={'ave_ssh': [0, 0.8]},
+                             colormap='gist_ncar',
+                             comout=os.path.join(comout, 'vrfy', 'bkgerr', 'ssh_steric_stddev')),   # ssh steric stddev
+                  plotConfig(grid_file=grid_file,
+                             data_file=os.path.join(comout, f'{RUN}.t{cyc}z.ocn.ssh_unbal_stddev.nc'),
+                             variables_horiz={'ave_ssh': [0, 0.8]},
+                             colormap='gist_ncar',
+                             comout=os.path.join(comout, 'vrfy', 'bkgerr', 'ssh_unbal_stddev')),   # ssh unbal stddev
+                  plotConfig(grid_file=grid_file,
+                             data_file=os.path.join(comout, f'{RUN}.t{cyc}z.ocn.ssh_total_stddev.nc'),
+                             variables_horiz={'ave_ssh': [0, 0.8]},
+                             colormap='gist_ncar',
+                             comout=os.path.join(comout, 'vrfy', 'bkgerr', 'ssh_total_stddev')),   # ssh total stddev
+                  plotConfig(grid_file=grid_file,
+                             data_file=os.path.join(comout, f'{RUN}.t{cyc}z.ocn.steric_explained_variance.nc'),
+                             variables_horiz={'ave_ssh': [0, 1]},
+                             colormap='seismic',
+                             comout=os.path.join(comout, 'vrfy', 'bkgerr', 'steric_explained_variance'))]   # steric explained variance
+    configs.extend(config_ens)
+
+if run_bkgerr_analysis:
+    config_bkgerr = [plotConfig(grid_file=grid_file,
+                                layer_file=layer_file,
+                                data_file=os.path.join(comout, os.path.pardir, os.path.pardir,
+                                                       'bmatrix', 'ocean', f'{RUN}.t'+cyc+'z.ocean.bkgerr_stddev.nc'),
+                                lats=np.arange(-60, 60, 10),
+                                lons=np.arange(-280, 80, 30),
+                                variables_zonal={'Temp': [0, 2],
+                                                 'Salt': [0, 0.2],
+                                                 'u': [0, 0.2],
+                                                 'v': [0, 0.2]},
+                                variables_meridional={'Temp': [0, 2],
+                                                      'Salt': [0, 0.2],
+                                                      'u': [0, 0.2],
+                                                      'v': [0, 0.2]},
+                                variables_horiz={'Temp': [0, 2],
+                                                 'Salt': [0, 0.2],
+                                                 'u': [0, 0.2],
+                                                 'v': [0, 0.2],
+                                                 'ave_ssh': [0, 0.1]},
+                                colormap='jet',
+                                comout=os.path.join(comout, 'vrfy', 'bkgerr'))]   # ocn bkgerr stddev
+    configs.extend(config_bkgerr)
+
+if run_bkg_analysis:
+    config_bkg = [plotConfig(grid_file=grid_file,
+                             data_file=os.path.join(com_ice_history, f'{RUN}.ice.t{gcyc}z.inst.f006.nc'),
+                             variables_horiz={'aice_h': [0.0, 1.0],
+                                              'hi_h': [0.0, 4.0],
+                                              'hs_h': [0.0, 0.5]},
+                             colormap='jet',
+                             projs=['North', 'South', 'Global'],
+                             comout=os.path.join(comout, 'vrfy', 'bkg')),   # sea ice background
+                  plotConfig(grid_file=grid_file,
+                             layer_file=layer_file,
+                             data_file=os.path.join(com_ocean_history, f'{RUN}.ocean.t{gcyc}z.inst.f006.nc'),
+                             lats=np.arange(-60, 60, 10),
+                             lons=np.arange(-280, 80, 30),
+                             variables_zonal={'Temp': [-1.8, 34.0],
+                                              'Salt': [32, 40]},
+                             variables_meridional={'Temp': [-1.8, 34.0],
+                                                   'Salt': [32, 40]},
+                             variables_horiz={'ave_ssh': [-1.8, 1.3],
+                                              'Temp': [-1.8, 34.0],
+                                              'Salt': [32, 40]},
+                             colormap='nipy_spectral',
+                             comout=os.path.join(comout, 'vrfy', 'bkg'))]
+    configs.extend(config_bkg)
+
+if run_increment_analysis:
+    config_incr = [plotConfig(grid_file=grid_file,
+                              layer_file=layer_file,
+                              data_file=os.path.join(comout, f'{RUN}.t'+cyc+'z.ocninc.nc'),
+                              lats=np.arange(-60, 60, 10),
+                              lons=np.arange(-280, 80, 30),
+                              variables_zonal={'Temp': [-0.5, 0.5],
+                                               'Salt': [-0.1, 0.1]},
+                              variables_horiz={'Temp': [-0.5, 0.5],
+                                               'Salt': [-0.1, 0.1],
+                                               'ave_ssh': [-0.1, 0.1]},
+                              variables_meridional={'Temp': [-0.5, 0.5],
+                                                    'Salt': [-0.1, 0.1]},
+                              colormap='seismic',
+                              comout=os.path.join(comout, 'vrfy', 'incr')),   # ocean increment
+                   plotConfig(grid_file=grid_file,
+                              data_file=os.path.join(comout, f'{RUN}.t'+cyc+'z.ice.incr.nc'),
+                              lats=np.arange(-60, 60, 10),
+                              variables_horiz={'aice_h': [-0.2, 0.2],
+                                               'hi_h': [-0.5, 0.5],
+                                               'hs_h': [-0.1, 0.1]},
+                              colormap='seismic',
+                              projs=['North', 'South'],
+                              comout=os.path.join(comout, 'vrfy', 'incr'))]   # sea ice increment
+    configs.extend(config_incr)
+
+
+# plot marine analysis vrfy
+
+def plot_marine_vrfy(config):
+    ocnvrfyPlotter = statePlotter(config)
+    ocnvrfyPlotter.plot()
+
+
+# Number of processes
+num_processes = len(configs)
+
+# Create a list to store the processes
+processes = []
+
+# Iterate over configs
+for config in configs[:num_processes]:
+    process = Process(target=plot_marine_vrfy, args=(config,))
+    process.start()
+    processes.append(process)
+
+# Wait for all processes to finish
+for process in processes:
+    process.join()
+
+#######################################
+# eva plots
+#######################################
+
+evadir = os.path.join(HOMEgfs, 'sorc', f'{RUN}.cd', 'ush', 'eva')
+marinetemplate = os.path.join(evadir, 'marine_gdas_plots.yaml')
+varyaml = os.path.join(comout, 'yaml', 'var_original.yaml')
+
+# it would be better to refrence the dirs explicitly with the comout path
+# but eva doesn't allow for specifying output directories
+os.chdir(os.path.join(comout, 'vrfy'))
+if not os.path.exists('preevayamls'):
+    os.makedirs('preevayamls')
+if not os.path.exists('evayamls'):
+    os.makedirs('evayamls')
+
+gen_eva_obs_yaml.gen_eva_obs_yaml(varyaml, marinetemplate, 'preevayamls')
+
+files = os.listdir('preevayamls')
+for file in files:
+    infile = os.path.join('preevayamls', file)
+    marine_eva_post.marine_eva_post(infile, 'evayamls', diagdir)
+
+files = os.listdir('evayamls')
+for file in files:
+    infile = os.path.join('evayamls', file)
+    print('running eva on', infile)
+    subprocess.run(['eva', infile], check=True)
+
+#######################################
+# calculate diag statistics
+#######################################
+
+# As of 11/12/2024 not working
+# diag_statistics.get_diag_stats()

utils/soca/fig_gallery/run_marine_analysis_vrfy_manual.job

@@ -0,0 +1,45 @@
+#!/bin/bash
+#SBATCH --job-name=marine_vrfy    # Assign a name to the job (customize as needed)
+#SBATCH --account=da-cpu
+#SBATCH --qos=debug
+#SBATCH -A da-cpu
+#SBATCH --output=run_marine_vrfy_analysis.out
+#SBATCH --nodes=1                 # Request 1 node
+#SBATCH --ntasks=40               # Request 40 total tasks (processors across nodes)
+#SBATCH --partition=hera          # Specify the partition (cluster) named "hera"
+#SBATCH --cpus-per-task=1         # Set 1 CPU per task (equivalent to ppn=40 and tpp=1)
+#SBATCH --mem=24GB                # Request 24GB of memory
+#SBATCH --time=00:30:00           # Set the walltime limit to 30 minutes
+
+# Define HOMEgfs
+export HOMEgfs="/scratch1/NCEPDEV/da/Mindo.Choi/workflow_11122024/global-workflow/"
+
+# Load EVA module 
+module use ${HOMEgfs}sorc/gdas.cd/modulefiles
+module load EVA/hera
+
+# Set PYTHONPATH using HOMEgfs
+export PYTHONPATH="${HOMEgfs}sorc/gdas.cd/ush/:\
+${HOMEgfs}sorc/gdas.cd/ush/eva/:\
+${HOMEgfs}sorc/gdas.cd/ush/soca/:\
+$PYTHONPATH"
+
+# Set flags to control plotConfig in the Python script
+export RUN_ENSENBLE_ANALYSIS=OFF       # Check if ensemble run is ON
+export RUN_BACKGROUND_ERROR_ANALYSIS=ON
+export RUN_BACKGROUND_ANALYSIS=ON
+export RUN_INCREMENT_ANLYSIS=ON
+
+# Define and export the environment variables
+export cyc="00"
+export RUN="gdas"
+export PSLOT="gdas_test"
+export PDY="20210827"
+
+# Define and export environment variables with paths
+export COM_OCEAN_ANALYSIS="/scratch1/NCEPDEV/da/Mindo.Choi/sandbox/marine_vrfy/gdas.20210827/00/analysis/ocean"
+export COM_ICE_HISTORY_PREV="/scratch1/NCEPDEV/da/Mindo.Choi/sandbox/marine_vrfy/gdas.20210826/18/model/ice/history"
+export COM_OCEAN_HISTORY_PREV="/scratch1/NCEPDEV/da/Mindo.Choi/sandbox/marine_vrfy/gdas.20210826/18/model/ocean/history"
+
+# Excute Marine Verify Analysis
+python3 ${HOMEgfs}sorc/gdas.cd/scripts/exgdas_global_marine_analysis_vrfy_manual.py