From 2bcedf26cea991fc61f640ccb5be39409b1bd931 Mon Sep 17 00:00:00 2001 From: Mindo Choi <141867620+apchoiCMD@users.noreply.github.com> Date: Fri, 15 Nov 2024 15:45:14 -0500 Subject: [PATCH] Enable the manual operation of the marine verification tool (#1373) #### This PR enables the Marine Verification Tool to run outside of the g-w CI workflow by submitting an `sbatch` job manually on Hera Includes, - Vrfy task run by a simple driver in the offline #1345 - Improve cosmetic issues we found #1349 - Bug fixes and more #1314 - ~~Move `exgdas_global_marine_analysis_vrfy.py` to `scripts/old` directory~~ Most up-to-date plots can be found at ``` /scratch1/NCEPDEV/da/Mindo.Choi/sandbox/marine_vrfy/gdas.20210827/00/analysis/ocean/vrfy_final_PR ``` The wall time is as follows: ``` [Mindo.Choi@hfe02 vrfy]$ sacct -j 2477688 --format=JobID,JobName,State,ExitCode,Elapsed JobID JobName State ExitCode Elapsed ------------ ---------- ---------- -------- ---------- 2477688 marine_vr+ COMPLETED 0:0 00:11:54 2477688.bat+ batch COMPLETED 0:0 00:11:54 2477688.ext+ extern COMPLETED 0:0 00:11:54 ``` Additional plotting work will be added by consolidating vrfy task as follows: - SST/SSH time series - Omb time series - Spatial SSH/SST/OHC - HTML (?) Close #1314 , Close #1345 , Close #1349 --------- Co-authored-by: Guillaume Vernieres --- scripts/exgdas_global_marine_analysis_vrfy.py | 0 ush/eva/marine_eva_post.py | 4 +- ush/eva/marine_gdas_plots.yaml | 6 +- ush/soca/soca_vrfy.py | 50 +++-- ...gdas_global_marine_analysis_vrfy_manual.py | 210 ++++++++++++++++++ .../run_marine_analysis_vrfy_manual.job | 45 ++++ 6 files changed, 297 insertions(+), 18 deletions(-) mode change 100755 => 100644 scripts/exgdas_global_marine_analysis_vrfy.py create mode 100644 utils/soca/fig_gallery/exgdas_global_marine_analysis_vrfy_manual.py create mode 100644 utils/soca/fig_gallery/run_marine_analysis_vrfy_manual.job diff --git a/scripts/exgdas_global_marine_analysis_vrfy.py b/scripts/exgdas_global_marine_analysis_vrfy.py old mode 100755 new mode 100644 diff --git a/ush/eva/marine_eva_post.py b/ush/eva/marine_eva_post.py index a355621a1..b537ddb3a 100755 --- a/ush/eva/marine_eva_post.py +++ b/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): diff --git a/ush/eva/marine_gdas_plots.yaml b/ush/eva/marine_gdas_plots.yaml index 5bedd1f69..0a903d0c4 100644 --- a/ush/eva/marine_gdas_plots.yaml +++ b/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}@CHANNELVAR@OmBQC.png 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} diff --git a/ush/soca/soca_vrfy.py b/ush/soca/soca_vrfy.py index 854d7ab69..a4060fecd 100755 --- a/ush/soca/soca_vrfy.py +++ b/ush/soca/soca_vrfy.py @@ -38,6 +38,18 @@ def plotConfig(grid_file=[], proj='set me', projs=['Global']): + # Map variable names to their units + variable_units = { + 'ave_ssh': 'meter', + 'Temp': 'deg C', + 'Salt': 'psu', + 'aice_h': 'meter', + 'hi_h': 'meter', + 'hs_h': 'meter', + 'u': 'm/s', + 'v': 'm/s' + } + """ Prepares the configuration for the plotting functions below """ @@ -64,6 +76,9 @@ def plotConfig(grid_file=[], config['variable'] = variable # the variable currently plotted config['projs'] = projs # all the projections etc. config['proj'] = proj + + # Add units to the config for each variable + config['variable_units'] = variable_units return config @@ -78,6 +93,7 @@ def plotHorizontalSlice(config): os.makedirs(dirname, exist_ok=True) variable = config['variable'] + unit = config['variable_units'].get(config['variable'], 'unknown') exp = config['exp'] PDY = config['PDY'] cyc = config['cyc'] @@ -85,12 +101,12 @@ def plotHorizontalSlice(config): if variable in ['Temp', 'Salt', 'u', 'v']: level = config['levels'][0] slice_data = np.squeeze(data[variable])[level, :, :] - label_colorbar = variable + ' Level ' + str(level) + label_colorbar = f"{variable} ({unit}) Level {level}" figname = os.path.join(dirname, variable + '_Level_' + str(level)) title = f"{exp} {PDY} {cyc} {variable} Level {level}" else: slice_data = np.squeeze(data[variable]) - label_colorbar = variable + label_colorbar = f"{variable} ({unit})" figname = os.path.join(dirname, variable + '_' + config['proj']) title = f"{exp} {PDY} {cyc} {variable}" @@ -99,17 +115,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 +136,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) @@ -138,6 +158,7 @@ def plotZonalSlice(config): Contourf of a zonal slice of an ocean field """ variable = config['variable'] + unit = config['variable_units'].get(config['variable'], 'unknown') exp = config['exp'] PDY = config['PDY'] cyc = config['cyc'] @@ -171,7 +192,7 @@ def plotZonalSlice(config): # Add colorbar for filled contours cbar = fig.colorbar(contourf_plot, ax=ax, shrink=0.5, orientation='horizontal') - cbar.set_label(variable + ' Lat ' + str(lat)) + cbar.set_label(f"{config['variable']} ({unit}) Lat {lat}") # Set the colorbar ticks cbar.set_ticks(contour_levels) @@ -184,7 +205,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) @@ -193,6 +214,7 @@ def plotMeridionalSlice(config): Contourf of a Meridional slice of an ocean field """ variable = config['variable'] + unit = config['variable_units'].get(config['variable'], 'unknown') exp = config['exp'] PDY = config['PDY'] cyc = config['cyc'] @@ -226,7 +248,7 @@ def plotMeridionalSlice(config): # Add colorbar for filled contours cbar = fig.colorbar(contourf_plot, ax=ax, shrink=0.5, orientation='horizontal') - cbar.set_label(variable + ' Lon ' + str(lon)) + cbar.set_label(f"{config['variable']} ({unit}) Lon {lon}") # Set the colorbar ticks cbar.set_ticks(contour_levels) @@ -239,7 +261,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) diff --git a/utils/soca/fig_gallery/exgdas_global_marine_analysis_vrfy_manual.py b/utils/soca/fig_gallery/exgdas_global_marine_analysis_vrfy_manual.py new file mode 100644 index 000000000..7c8efd0a6 --- /dev/null +++ b/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() diff --git a/utils/soca/fig_gallery/run_marine_analysis_vrfy_manual.job b/utils/soca/fig_gallery/run_marine_analysis_vrfy_manual.job new file mode 100644 index 000000000..38ce48ffc --- /dev/null +++ b/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/utils/soca/fig_gallery/exgdas_global_marine_analysis_vrfy_manual.py