Skip to content

Commit

Permalink
implemented ctest for meridional plotting capability
Browse files Browse the repository at this point in the history
  • Loading branch information
Mindo Choi authored and Mindo Choi committed Aug 22, 2023
1 parent 09757ce commit 1698c15
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 40 deletions.
73 changes: 42 additions & 31 deletions scripts/exgdas_global_marine_analysis_vrfy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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'))
Expand All @@ -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'))
Expand All @@ -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'))
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
#######################################
Expand Down
81 changes: 72 additions & 9 deletions ush/soca/soca_vrfy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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']):

Expand All @@ -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
Expand Down Expand Up @@ -86,14 +94,15 @@ 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())
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.close(fig)


def plotZonalSlice(config):
Expand Down Expand Up @@ -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:

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

0 comments on commit 1698c15

Please sign in to comment.