diff --git a/scripts/exgdas_global_marine_analysis_vrfy.py b/scripts/exgdas_global_marine_analysis_vrfy.py index cd302e1c0..2a61994d5 100755 --- a/scripts/exgdas_global_marine_analysis_vrfy.py +++ b/scripts/exgdas_global_marine_analysis_vrfy.py @@ -22,7 +22,7 @@ import gen_eva_obs_yaml import marine_eva_post import diag_statistics -from soca_vrfy import StatePlotter, plot_config, plot_increment, plot_analysis +from soca_vrfy import StatePlotter, plot_config import subprocess from datetime import datetime, timedelta @@ -40,27 +40,111 @@ diagdir = os.path.join(comout, 'diags') HOMEgfs = os.getenv('HOMEgfs') +####################################### +# ocean increment +####################################### + +data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocninc.nc') +config = plot_config(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() ####################################### -# INCREMENT +# sea ice increment ####################################### -plot_increment(comout, cyc, RUN, grid_file) +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, + 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() ####################################### -# Analysis/Background +# sea ice analysis ####################################### -plot_analysis(comout, - com_ice_history, - com_ocean_history, - cyc, - RUN, - grid_file, - gcyc) +data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.iceana.nc') +config = plot_config(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() +####################################### +# sea ice background +####################################### +data_file = os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc') +config = plot_config(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() +####################################### +# ocean surface analysis +####################################### + +data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc') +config = plot_config(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() + +####################################### +# ocean surface background +####################################### + +data_file = os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc') +config = plot_config(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() + +####################################### +# background error +####################################### data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc') config = plot_config(grid_file = grid_file, @@ -76,10 +160,9 @@ 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') @@ -108,5 +191,6 @@ ####################################### # calculate diag statistics +####################################### diag_statistics.get_diag_stats() diff --git a/ush/soca/soca_vrfy.py b/ush/soca/soca_vrfy.py index cfb4e7336..5e8c9d492 100755 --- a/ush/soca/soca_vrfy.py +++ b/ush/soca/soca_vrfy.py @@ -28,27 +28,30 @@ def plot_config(grid_file=[], variables_horiz = [], variables_zonal = [], lat=np.nan, - lats=np.arange(-60, 60, 10)): + lats=np.arange(-60, 60, 10), + proj='set me', + projs=['Global']): """ Prepares the configuration for the plotting functions below """ config = {} + config['comout'] = comout # output directory config['grid file'] = grid_file config['fields file'] = data_file - config['variable'] = variable config['levels'] = [1] - config['bounds'] = bounds - config['all bounds'] = allbounds config['colormap'] = colormap + config['all bounds'] = allbounds + config['bounds'] = bounds config['lats'] = lats # all the lats to plot config['lat'] = lat # the lat being currently plotted - config['comout'] = comout - config['max depth'] = max_depth - config['max depths'] = max_depths - config['horiz variables'] = variables_horiz - config['zonal variables'] = variables_zonal - config['proj'] = 'Global' + config['max depths'] = max_depths # all the max depths to plot + config['max depth'] = max_depth # the max depth currently plotted + config['horiz variables'] = variables_horiz # all the vars for horiz plots + config['zonal variables'] = variables_zonal # all the vars for zonal plots + config['variable'] = variable # the variable currently plotted + config['projs'] = projs # all the projections etc. + config['proj'] = proj return config @@ -120,194 +123,10 @@ def plot_zonal_slice(config): plt.savefig(figname, bbox_inches='tight', dpi=600) -def plot_increment(comout, cyc, RUN, grid_file): - - 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')) - - ####################################### - # zonal slices - - for lat in np.arange(-60, 60, 10): - - for max_depth in [700.0, 5000.0]: - config['lat'] = 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) - - ####################################### - # 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) - - # snow depth - config.update({'variable': 'hsnon', 'bounds': [-0.1, 0.1], 'proj': proj}) - plot_horizontal_slice(config) - - -def plot_analysis(comout, - com_ice_history, - com_ocean_history, - cyc, - RUN, - grid_file, - gcyc): - - ####################################### - # 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) - - ####################################### - # 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) - - # sss - config.update({'variable': 'Salt', 'bounds': [30, 38], 'proj': proj, 'levels': [1]}) - plot_horizontal_slice(config) - -#def plot_error(comout, -# com_ice_history, -# com_ocean_history, -# cyc, -# RUN, -# grid_file, -# gcyc): - -# ####################################### -# # Std Bkg. Error -# ####################################### -# 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['lat'] = 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) -# -# ####################################### -# # Horizontal slices -# -# # Temperature -# config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 2]}) -# plot_horizontal_slice(config) -# -# # Salinity -# config.update({'variable': 'Salt', 'bounds': [0, 0.2]}) -# plot_horizontal_slice(config) -# -# # Sea surface height -# config.update({'variable': 'ave_ssh', 'bounds': [0, 0.1]}) -# plot_horizontal_slice(config) -# class StatePlotter: def __init__(self, config_dict): self.config = config_dict -# self.grid_file = config['grid file'] -# self.data_file = config['fields file'] -# self.variable = config['variable'] -# self.levels = config['levels'] -# self.bounds = config['bounds'] -# self.colormap = config['colormap'] -# self.lats = config['lats'] -# self.comout = config['comout'] -# self.maxdepth = 5000.0 -# self.proj = 'Global' -# - ... - # Where config contains the bounds, filename and whatever else is needed to plot stuff def plot(self): # Loop over variables, slices (horiz and vertical) and projections ... and whatever else is needed @@ -322,33 +141,14 @@ def plot(self): self.config['max depth'] = max_depth for variable in self.config['zonal variables']: - # Temperature bounds = self.config['all bounds'][variable] self.config.update({'variable': variable, 'bounds': bounds}) plot_zonal_slice(self.config) - # Salinity -# self.config.update({'variable': 'Salt', 'levels': [1], 'bounds': [0, .2]}) -# plot_zonal_slice(config) - ####################################### # Horizontal slices - - for variable in self.config['horiz variables']: - bounds = self.config['all bounds'][variable] - self.config.update({'variable': variable, 'bounds': bounds}) - plot_horizontal_slice(self.config) - -# # Temperature -# self.config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 2]}) -# plot_horizontal_slice(config) -# -# # Salinity -# self.config.update({'variable': 'Salt', 'bounds': [0, 0.2]}) -# plot_horizontal_slice(config) -# -# # Sea surface height -# self.config.update({'variable': 'ave_ssh', 'bounds': [0, 0.1]}) -# plot_horizontal_slice(config) - - + for proj in self.config['projs']: + for variable in self.config['horiz variables']: + bounds = self.config['all bounds'][variable] + self.config.update({'variable': variable, 'bounds': bounds, 'proj': proj}) + plot_horizontal_slice(self.config)