diff --git a/scripts/exgdas_global_marine_analysis_vrfy.py b/scripts/exgdas_global_marine_analysis_vrfy.py index 88a33f49b..49ecf93d5 100755 --- a/scripts/exgdas_global_marine_analysis_vrfy.py +++ b/scripts/exgdas_global_marine_analysis_vrfy.py @@ -48,11 +48,18 @@ 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]}, + ### + lons=np.arange(-280, 80, 10), + ### + 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]}, + ### Added lines for meridional + variables_meridional={'Temp': [-0.5, 0.5], + 'Salt': [-0.1, 0.1]}, + ### colormap='RdBu', comout=os.path.join(comout, 'vrfy', 'incr')) ocnIncPlotter = statePlotter(config) @@ -66,10 +73,9 @@ 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]}, + variables_horiz={'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')) @@ -83,10 +89,9 @@ 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]}, + variables_horiz={'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')) @@ -100,10 +105,9 @@ 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]}, + variables_horiz={'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')) @@ -117,10 +121,9 @@ 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]}, + variables_horiz={'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) @@ -133,10 +136,9 @@ 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]}, + variables_horiz={'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) @@ -150,16 +152,25 @@ 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]}, + ### + lons=np.arange(-280, 80, 10), + ### + variables_zonal={'Temp': [0, 2], + 'Salt': [0, 0.2]}, + variables_horiz={'Temp': [0, 2], + 'Salt': [0, 0.2], + 'ave_ssh': [0, 0.1]}, + ### Added lines for meridional + variables_meridional={'Temp': [0, 2], #### + 'Salt': [0, 0.2]}, #### + ### colormap='jet', comout=os.path.join(comout, 'vrfy', 'bkgerr')) bkgErrPlotter = statePlotter(config) bkgErrPlotter.plot() - +# +quit() # for quick vrty working with ctest +# ####################################### # eva plots ####################################### diff --git a/ush/soca/soca_vrfy.py b/ush/soca/soca_vrfy.py index 248505711..1bba0f7f9 100755 --- a/ush/soca/soca_vrfy.py +++ b/ush/soca/soca_vrfy.py @@ -19,16 +19,20 @@ 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=[], + variables_horiz={}, + variables_zonal={}, + #-------------------------- + variables_meridional={}, + #-------------------------- lat=np.nan, lats=np.arange(-60, 60, 10), + lon=np.nan, + lons=np.arange(-280, 80, 10), proj='set me', projs=['Global']): @@ -41,14 +45,18 @@ def plotConfig(grid_file=[], config['fields file'] = data_file config['levels'] = [1] 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['lons'] = lons + config['lon'] = lon 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['meridional variables'] = variables_meridional # all the vars for meridional plots + #-------------------------- config['variable'] = variable # the variable currently plotted config['projs'] = projs # all the projections etc. config['proj'] = proj @@ -86,7 +94,7 @@ def plotHorizontalSlice(config): cmap=config['colormap']) plt.colorbar(label=label_colorbar, shrink=0.5, orientation='horizontal') - ax.coastlines() # TODO: make this work on hpc + #ax.coastlines() # TODO: make this work on hpc ax.gridlines(draw_labels=True) if config['proj'] == 'South': ax.set_extent([-180, 180, -90, -50], ccrs.PlateCarree()) @@ -94,6 +102,7 @@ def plotHorizontalSlice(config): 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.close(fig) def plotZonalSlice(config): @@ -121,7 +130,42 @@ def plotZonalSlice(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) + plt.close(fig) +########################################################################################## +# +def plotMeridionalSlice(config): + """ + pcolormesh of a Meridional slice of an ocean field + """ + lon = float(config['lon']) # <- changed lat to lon + 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))) + lon_index = np.argmin(np.array(np.abs(np.squeeze(grid.lon)[:, 0]-lon))) # <- changed lat to lon + # + slice_data = np.squeeze(np.array(data[config['variable']]))[:, lon_index, :] # <- changed lat to lon0 + depth = np.squeeze(np.array(grid['h']))[:, lon_index, :] # <- changed lat to lon + depth[np.where(np.abs(depth) > 10000.0)] = 0.0 + depth = np.cumsum(depth, axis=0) + bounds = config['bounds'] + # + y = np.tile(np.squeeze(grid.lon[:, lon_index]), (np.shape(depth)[0], 1)) # <- changed lat to lon NEED to be attention + # + fig, ax = plt.subplots(figsize=(8, 5)) + plt.pcolormesh(y, -depth, slice_data, + vmin=bounds[0], vmax=bounds[1], + cmap=config['colormap']) + plt.colorbar(label=config['variable']+' Lon '+str(lon), shrink=0.5, orientation='horizontal') + ax.set_ylim(-config['max depth'], 0) + dirname = os.path.join(config['comout'], config['variable']) + 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.close(fig) +# +########################################################################################## class statePlotter: @@ -140,15 +184,34 @@ def plot(self): 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] + variableBounds = self.config['zonal variables'] + for variable in variableBounds.keys(): + bounds = variableBounds[variable] self.config.update({'variable': variable, 'bounds': bounds}) plotZonalSlice(self.config) + + ###################################### + # Meridional slices + + for lon in self.config['lons']: + self.config['lon'] = lon + + for max_depth in self.config['max depths']: + self.config['max depth'] = max_depth + + variableBounds = self.config['meridional variables'] + for variable in variableBounds.keys(): + bounds = variableBounds[variable] + self.config.update({'variable': variable, 'bounds': bounds}) + plotMeridionalSlice(self.config) + ####################################### # Horizontal slices for proj in self.config['projs']: - for variable in self.config['horiz variables']: - bounds = self.config['all bounds'][variable] + + variableBounds = self.config['horiz variables'] + for variable in variableBounds.keys(): + bounds = variableBounds[variable] self.config.update({'variable': variable, 'bounds': bounds, 'proj': proj}) plotHorizontalSlice(self.config)