Skip to content

Commit

Permalink
Enable the manual operation of the marine verification tool (#1373)
Browse files Browse the repository at this point in the history
#### 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 <[email protected]>
  • Loading branch information
apchoiCMD and guillaumevernieres authored Nov 15, 2024
1 parent e514b92 commit 2bcedf2
Show file tree
Hide file tree
Showing 6 changed files with 297 additions and 18 deletions.
Empty file modified scripts/exgdas_global_marine_analysis_vrfy.py
100755 → 100644
Empty file.
4 changes: 3 additions & 1 deletion ush/eva/marine_eva_post.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
6 changes: 3 additions & 3 deletions ush/eva/marine_gdas_plots.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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}

Expand Down
50 changes: 36 additions & 14 deletions ush/soca/soca_vrfy.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""
Expand All @@ -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


Expand All @@ -78,19 +93,20 @@ 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']

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}"

Expand All @@ -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
Expand All @@ -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)


Expand All @@ -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']
Expand Down Expand Up @@ -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)
Expand All @@ -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)


Expand All @@ -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']
Expand Down Expand Up @@ -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)
Expand All @@ -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)


Expand Down
210 changes: 210 additions & 0 deletions utils/soca/fig_gallery/exgdas_global_marine_analysis_vrfy_manual.py
Original file line number Diff line number Diff line change
@@ -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()
Loading

0 comments on commit 2bcedf2

Please sign in to comment.