Skip to content

Commit

Permalink
Refactoring of vrfy - transferred bulk of plotting to sub-script (#538)
Browse files Browse the repository at this point in the history 
* initial changes for refactoring

* python coding norms

* removed paleolithic header

* python coding norms

* further removal of plotting guts from script script

* python style

* python style

* added plotting class

* rest of plot suites added, cleanup

* python code style

* python code style

* correcting camel case

* python coding norms

* python style thing in something completely unrelated
  • Loading branch information
@AndrewEichmann-NOAA
AndrewEichmann-NOAA authored Aug 3, 2023
1 parent 75ed6fe commit 12c3d4f
Show file tree
Hide file tree
Showing 3 changed files with 166 additions and 152 deletions.
242 changes: 102 additions & 140 deletions scripts/exgdas_global_marine_analysis_vrfy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,189 +22,147 @@
import gen_eva_obs_yaml
import marine_eva_post
import diag_statistics
from soca_vrfy import plot_config, plot_horizontal_slice, plot_zonal_slice
from soca_vrfy import statePlotter, plotConfig
import subprocess
from datetime import datetime, timedelta

comout = os.getenv('COM_OCEAN_ANALYSIS')
com_ice_history = os.getenv('COM_ICE_HISTORY_PREV')
com_ocean_history = os.getenv('COM_OCEAN_HISTORY_PREV')
data = os.getenv('DATA')
pdy = os.getenv('PDY')
cyc = os.getenv('cyc')
RUN = os.getenv('RUN')
bcyc = str((int(cyc) - 3) % 24).zfill(2)
gcyc = str((int(cyc) - 6) % 24).zfill(2)
gcdate = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') - timedelta(hours=int(os.getenv('assim_freq')))

bcyc = str((int(cyc) - 3) % 24).zfill(2)
grid_file = os.path.join(comout, f'{RUN}.t'+bcyc+'z.ocngrid.nc')

# for eva
diagdir = os.path.join(comout, 'diags')
HOMEgfs = os.getenv('HOMEgfs')


#######################################
# INCREMENT
# ocean increment
#######################################
incr_cmap = 'RdBu'

data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocninc.nc')
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap=incr_cmap,
comout=os.path.join(comout, 'vrfy', 'incr'))
config = plotConfig(grid_file=grid_file,
data_file=data_file,
lats=np.arange(-60, 60, 10),
variables_zonal=['Temp', 'Salt'],
variables_horiz=['Temp', 'Salt', 'ave_ssh'],
allbounds={'Temp': [-0.5, 0.5],
'Salt': [-0.1, 0.1],
'ave_ssh': [-0.1, 0.1]},
colormap='RdBu',
comout=os.path.join(comout, 'vrfy', 'incr'))
ocnIncPlotter = statePlotter(config)
ocnIncPlotter.plot()

#######################################
# zonal slices

for lat in np.arange(-60, 60, 10):

for max_depth in [700.0, 5000.0]:
config['lats'] = [lat]
config['max depth'] = max_depth

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-.5, .5]})
plot_zonal_slice(config)

# Salinity
config.update({'variable': 'Salt', 'levels': [1], 'bounds': [-.1, .1]})
plot_zonal_slice(config)

# sea ice increment
#######################################
# Horizontal slices

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-1, 1]})
plot_horizontal_slice(config)

# Salinity
config.update({'variable': 'Salt', 'bounds': [-0.1, 0.1]})
plot_horizontal_slice(config)

# Sea surface height
config.update({'variable': 'ave_ssh', 'bounds': [-0.1, 0.1]})
plot_horizontal_slice(config)

#######################################
# Sea ice
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ice.incr.nc')
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap=incr_cmap,
comout=os.path.join(comout, 'vrfy', 'incr'))

for proj in ['North', 'South']:
# concentration
config.update({'variable': 'aicen', 'bounds': [-0.2, 0.2], 'proj': proj})
plot_horizontal_slice(config)

# thickness
config.update({'variable': 'hicen', 'bounds': [-0.5, 0.5], 'proj': proj})
plot_horizontal_slice(config)
config = plotConfig(grid_file=grid_file,
data_file=data_file,
lats=np.arange(-60, 60, 10),
variables_horiz=['aicen', 'hicen', 'hsnon'],
allbounds={'aicen': [-0.2, 0.2],
'hicen': [-0.5, 0.5],
'hsnon': [-0.1, 0.1]},
colormap='RdBu',
projs=['North', 'South'],
comout=os.path.join(comout, 'vrfy', 'incr'))
iceIncPlotter = statePlotter(config)
iceIncPlotter.plot()

# snow depth
config.update({'variable': 'hsnon', 'bounds': [-0.1, 0.1], 'proj': proj})
plot_horizontal_slice(config)

#######################################
# Analysis/Background
#######################################

# sea ice analysis
#######################################
# Sea ice
data_files = [os.path.join(comout, f'{RUN}.t{cyc}z.iceana.nc'),
os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc')]
dirs_out = ['ana', 'bkg']
ice_vars = {'bkg': ['aice_h', 'hs_h', 'hi_h'], 'ana': ['aicen', 'hicen', 'hsnon']}
for data_file, dir_out in zip(data_files, dirs_out):
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap='jet',
comout=os.path.join(comout, 'vrfy', dir_out))

for proj in ['North', 'South', 'Global']:
# concentration
var = ice_vars[dir_out]
config.update({'variable': var[0], 'bounds': [0.0, 1.0], 'proj': proj})
plot_horizontal_slice(config)

# thickness
config.update({'variable': var[1], 'bounds': [0.0, 4.0], 'proj': proj})
plot_horizontal_slice(config)

# snow depth
config.update({'variable': var[2], 'bounds': [0.0, 0.5], 'proj': proj})
plot_horizontal_slice(config)
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.iceana.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['aicen', 'hicen', 'hsnon'],
allbounds={'aicen': [0.0, 1.0],
'hicen': [0.0, 4.0],
'hsnon': [0.0, 0.5]},
colormap='jet',
projs=['North', 'South', 'Global'],
comout=os.path.join(comout, 'vrfy', 'ana'))
iceAnaPlotter = statePlotter(config)
iceAnaPlotter.plot()

#######################################
# Ocean surface
data_files = [os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc'),
os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc')]
dirs_out = ['ana', 'bkg']
ocn_vars = ['ave_ssh', 'Temp', 'Salt']
for data_file, dir_out in zip(data_files, dirs_out):
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap='jet',
comout=os.path.join(comout, 'vrfy', dir_out))

# ssh
config.update({'variable': 'ave_ssh', 'bounds': [-1.8, 1.3], 'proj': proj, 'levels': [1]})
plot_horizontal_slice(config)

# sst
config.update({'variable': 'Temp', 'bounds': [-1.8, 34.0], 'proj': proj, 'levels': [1]})
plot_horizontal_slice(config)
# sea ice background
#######################################

# sss
config.update({'variable': 'Salt', 'bounds': [30, 38], 'proj': proj, 'levels': [1]})
plot_horizontal_slice(config)
data_file = os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['aice_h', 'hs_h', 'hi_h'],
allbounds={'aice_h': [0.0, 1.0],
'hs_h': [0.0, 4.0],
'hi_h': [0.0, 0.5]},
colormap='jet',
projs=['North', 'South', 'Global'],
comout=os.path.join(comout, 'vrfy', 'bkg'))
iceBkgPlotter = statePlotter(config)
iceBkgPlotter.plot()

#######################################
# Std Bkg. Error
# ocean surface analysis
#######################################
bmat_cmap = 'jet'
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc')
config = plot_config(grid_file=grid_file,
data_file=data_file,
colormap=bmat_cmap,
comout=os.path.join(comout, 'vrfy', 'bkgerr'))

#######################################
# zonal slices

for lat in np.arange(-60, 60, 10):

for max_depth in [700.0, 5000.0]:
config['lats'] = [lat]
config['max depth'] = max_depth

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 1.5]})
plot_zonal_slice(config)

# Salinity
config.update({'variable': 'Salt', 'levels': [1], 'bounds': [0, .2]})
plot_zonal_slice(config)
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['ave_ssh', 'Temp', 'Salt'],
allbounds={'ave_ssh': [-1.8, 1.3],
'Temp': [-1.8, 34.0],
'Salt': [30, 38]},
colormap='jet',
comout=os.path.join(comout, 'vrfy', 'ana'))
ocnAnaPlotter = statePlotter(config)
ocnAnaPlotter.plot()

#######################################
# Horizontal slices
# ocean surface background
#######################################

# Temperature
config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 2]})
plot_horizontal_slice(config)
data_file = os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
variables_horiz=['ave_ssh', 'Temp', 'Salt'],
allbounds={'ave_ssh': [-1.8, 1.3],
'Temp': [-1.8, 34.0],
'Salt': [30, 38]},
colormap='jet',
comout=os.path.join(comout, 'vrfy', 'bkg'))
ocnBkgPlotter = statePlotter(config)
ocnBkgPlotter.plot()

# Salinity
config.update({'variable': 'Salt', 'bounds': [0, 0.2]})
plot_horizontal_slice(config)
#######################################
# background error
#######################################

# Sea surface height
config.update({'variable': 'ave_ssh', 'bounds': [0, 0.1]})
plot_horizontal_slice(config)
data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc')
config = plotConfig(grid_file=grid_file,
data_file=data_file,
lats=np.arange(-60, 60, 10),
variables_zonal=['Temp', 'Salt'],
variables_horiz=['Temp', 'Salt', 'ave_ssh'],
allbounds={'Temp': [0, 2],
'Salt': [0, 0.2],
'ave_ssh': [0, 0.1]},
colormap='jet',
comout=os.path.join(comout, 'vrfy', 'bkgerr'))
bkgErrPlotter = statePlotter(config)
bkgErrPlotter.plot()

#######################################
# eva plots
#######################################

evadir = os.path.join(HOMEgfs, 'sorc', f'{RUN}.cd', 'ush', 'eva')
marinetemplate = os.path.join(evadir, 'marine_gdas_plots.yaml')
Expand All @@ -231,4 +189,8 @@
print('running eva on', infile)
subprocess.run(['eva', infile], check=True)

#######################################
# calculate diag statistics
#######################################

diag_statistics.get_diag_stats()
0 ush/diag_statistics.py → ush/soca/diag_statistics.py
View file
File renamed without changes.
Loading

0 comments on commit 12c3d4f

Please sign in to comment.