diff --git a/scripts/exgdas_global_marine_analysis_vrfy.py b/scripts/exgdas_global_marine_analysis_vrfy.py index ed68fb100..88a33f49b 100755 --- a/scripts/exgdas_global_marine_analysis_vrfy.py +++ b/scripts/exgdas_global_marine_analysis_vrfy.py @@ -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') @@ -231,4 +189,8 @@ print('running eva on', infile) subprocess.run(['eva', infile], check=True) +####################################### +# calculate diag statistics +####################################### + diag_statistics.get_diag_stats() diff --git a/ush/diag_statistics.py b/ush/soca/diag_statistics.py similarity index 100% rename from ush/diag_statistics.py rename to ush/soca/diag_statistics.py diff --git a/ush/soca_vrfy.py b/ush/soca/soca_vrfy.py similarity index 59% rename from ush/soca_vrfy.py rename to ush/soca/soca_vrfy.py index f06e4cd4f..248505711 100755 --- a/ush/soca_vrfy.py +++ b/ush/soca/soca_vrfy.py @@ -15,26 +15,47 @@ 'Global': ccrs.Mollweide(central_longitude=-150)} -def plot_config(grid_file=[], data_file=[], - variable=[], levels=[], bounds=[], colormap=[], comout=[], lats=[]): +def plotConfig(grid_file=[], + data_file=[], + variable=[], + levels=[], + allbounds=[], + bounds=[], + colormap=[], + max_depth=np.nan, + max_depths=[700.0, 5000.0], + comout=[], + variables_horiz=[], + variables_zonal=[], + lat=np.nan, + 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'] = levels - config['bounds'] = bounds + config['levels'] = [1] config['colormap'] = colormap - config['lats'] = lats - config['comout'] = comout - config['max depth'] = 5000.0 - config['proj'] = 'Global' + config['all bounds'] = allbounds + config['bounds'] = bounds + config['lats'] = lats # all the lats to plot + config['lat'] = lat # the lat being currently plotted + 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 -def plot_horizontal_slice(config): +def plotHorizontalSlice(config): """ pcolormesh of a horizontal slice of an ocean field """ @@ -75,11 +96,11 @@ def plot_horizontal_slice(config): plt.savefig(figname, bbox_inches='tight', dpi=600) -def plot_zonal_slice(config): +def plotZonalSlice(config): """ pcolormesh of a zonal slice of an ocean field """ - lat = float(config['lats'][0]) + lat = float(config['lat']) grid = xr.open_dataset(config['grid file']) data = xr.open_dataset(config['fields file']) lat_index = np.argmin(np.array(np.abs(np.squeeze(grid.lat)[:, 0]-lat))) @@ -100,3 +121,34 @@ def plot_zonal_slice(config): 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) + + +class statePlotter: + + def __init__(self, config_dict): + self.config = config_dict + + def plot(self): + # Loop over variables, slices (horiz and vertical) and projections ... and whatever else is needed + + ####################################### + # zonal slices + + for lat in self.config['lats']: + self.config['lat'] = lat + + for max_depth in self.config['max depths']: + self.config['max depth'] = max_depth + + for variable in self.config['zonal variables']: + bounds = self.config['all bounds'][variable] + self.config.update({'variable': variable, 'bounds': bounds}) + plotZonalSlice(self.config) + + ####################################### + # Horizontal slices + 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}) + plotHorizontalSlice(self.config)