diff --git a/diagnostics/physics/NWA12/plot_gulf_stream.py b/diagnostics/physics/NWA12/plot_gulf_stream.py new file mode 100644 index 000000000..5a869746f --- /dev/null +++ b/diagnostics/physics/NWA12/plot_gulf_stream.py @@ -0,0 +1,157 @@ +""" +Plot of the Gulf Stream position and index, +Uses whatever model data can be found within the directory pp_root, +and does not try to match the model and observed time periods. +How to use: +python plot_gulf_stream.py -p /archive/acr/fre/NWA/2023_04/NWA12_COBALT_2023_04_kpo4-coastatten-physics/gfdl.ncrc5-intel22-prod +""" +import xarray +import xesmf +import pandas as pd +import numpy as np +import cartopy.feature as feature +import cartopy.crs as ccrs +from cartopy.mpl.geoaxes import GeoAxes +import matplotlib.gridspec as gridspec +from matplotlib.lines import Line2D +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import AxesGrid + +# Need to append physics dir to path to access plot common +import sys +sys.path.append("..") +from plot_common import open_var, add_ticks, save_figure + +def compute_gs(ssh, data_grid=None): + lons = np.arange(360-72, 360-51.9, 1) + lats = np.arange(36, 42, 0.1) + target_grid = {'lat': lats, 'lon': lons} + + if data_grid is None: + data_grid = {'lat': ssh.lat, 'lon': ssh.lon} + + ssh_to_grid = xesmf.Regridder( + data_grid, + target_grid, + method='bilinear' + ) + + # Interpolate the SSH data onto the index grid. + regridded = ssh_to_grid(ssh) + + # Find anomalies relative to the calendar month mean SSH over the full model run. + anom = regridded.groupby('time.month') - regridded.groupby('time.month').mean('time') + + # For each longitude point, the Gulf Stream is located at the latitude with the maximum SSH anomaly variance. + stdev = anom.std('time') + amax = stdev.argmax('lat').compute() + gs_points = stdev.lat.isel(lat=amax).compute() + + # The index is the mean latitude of the Gulf Stream, divided by the standard deviation of the mean latitude of the Gulf Stream. + index = ((anom.isel(lat=amax).mean('lon')) / anom.isel(lat=amax).mean('lon').std('time')).compute() + + # Move times to the beginning of the month to match observations. + monthly_index = index.to_pandas().resample('1MS').first() + return monthly_index, gs_points + +def plot_gulf_stream(pp_root, label): + + # Load Natural Earth Shapefiles + _LAND_50M = feature.NaturalEarthFeature( + 'physical', 'land', '50m', + edgecolor='face', + facecolor='#999999' + ) + + # Get model grid + model_grid = xarray.open_dataset( '../../data/geography/ocean_static.nc' ) + + # Get model thetao data TODO: maki this comment better + model_thetao = open_var(pp_root, 'ocean_monthly_z', 'thetao') + + if '01_l' in model_thetao.coords: + model_thetao = model_thetao.rename({'01_l': 'z_l'}) + + model_t200 = model_thetao.interp(z_l=200).mean('time') + + # Ideally would use SSH, but some diag_tables only saved zos + try: + model_ssh = open_var(pp_root, 'ocean_monthly', 'ssh') + except: + print('Using zos') + model_ssh = open_var(pp_root, 'ocean_monthly', 'zos') + + model_ssh_index, model_ssh_points = compute_gs( + model_ssh, + data_grid=model_grid[['geolon', 'geolat']].rename({'geolon': 'lon', 'geolat': 'lat'}) + ) + + # Get Glorys data + glorys_t200 = xarray.open_dataarray('../../data/diagnostics/glorys_T200.nc') + + # Get satellite points + #satellite_ssh_index, satellite_ssh_points = compute_gs(satellite['adt']) + #satellite_ssh_points.to_netcdf('../data/obs/satellite_ssh_points.nc') + #satellite_ssh_index.to_pickle('../data/obs/satellite_ssh_index.pkl') + #read pre-calculate satellite_ssh_index and points + satellite_ssh_points = xarray.open_dataset('../../data/obs/satellite_ssh_points.nc') + satellite_ssh_index = pd.read_pickle('../../data/obs/satellite_ssh_index.pkl') + satellite_rolled = satellite_ssh_index.rolling(25, center=True, min_periods=25).mean().dropna() + + #satellite = xarray.open_mfdataset([f'/net2/acr/altimetry/SEALEVEL_GLO_PHY_L4_MY_008_047/adt_{y}_{m:02d}.nc' for y in range(1993, 2020) for m in range(1, 13)]) + #satellite = satellite.rename({'longitude': 'lon', 'latitude': 'lat'}) + #satellite = satellite.resample(time='1MS').mean('time') + + # Get rolling averages and correlations + model_rolled = model_ssh_index.rolling(25, center=True, min_periods=25).mean().dropna() + corr = pd.concat((model_ssh_index, satellite_ssh_index), axis=1).corr().iloc[0, 1] + corr_rolled = pd.concat((model_rolled, satellite_rolled), axis=1).corr().iloc[0, 1] + + # Plot of Gulf Stream position and index based on SSH, + # plus position based on T200. + fig = plt.figure(figsize=(10, 6), tight_layout=True) + gs = gridspec.GridSpec(2, 2, hspace=.25) + + # Set projection of input data files so that data is correctly tranformed when plotting + proj = ccrs.PlateCarree() + + ax = fig.add_subplot(gs[0, 0], projection = proj) + ax.add_feature(_LAND_50M) + ax.contour(model_grid.geolon, model_grid.geolat, model_t200, levels=[15], colors='r') + ax.contour(glorys_t200.longitude, glorys_t200.latitude, glorys_t200, levels=[15], colors='k') + add_ticks(ax, xlabelinterval=5) + ax.set_extent([-82, -50, 25, 41]) + ax.set_title('(a) Gulf Stream position based on T200') + custom_lines = [Line2D([0], [0], color=c, lw=2) for c in ['r', 'k']] + ax.legend(custom_lines, ['Model', 'GLORYS12'], loc='lower right', frameon=False) + + ax = fig.add_subplot(gs[0, 1], projection = proj) + ax.add_feature(_LAND_50M) + ax.plot(model_ssh_points.lon-360, model_ssh_points, c='r') + ax.plot(satellite_ssh_points.lon-360, satellite_ssh_points['__xarray_dataarray_variable__'], c='k') + add_ticks(ax, xlabelinterval=5) + ax.set_extent([-82, -50, 25, 41]) + ax.set_title('(b) Gulf Stream position based on SSH variance') + ax.legend(custom_lines, ['Model', 'Altimetry'], loc='lower right', frameon=False) + + ax = fig.add_subplot(gs[1, :]) + model_ssh_index.plot(ax=ax, c='#ffbbbb', label='Model') + satellite_ssh_index.plot(ax=ax, c='#bbbbbb', label=f'Altimetry (r={corr:2.2f})') + model_rolled.plot(ax=ax, c='r', label='Model rolling mean') + satellite_rolled.plot(ax=ax, c='k', label=f'Altimetry rolling mean (r={corr_rolled:2.2f})') + ax.set_title('(c) Gulf Stream index based on SSH variance') + ax.set_xlabel('') + ax.set_ylim(-3, 3) + ax.set_ylabel('Index (positive north)') + ax.legend(ncol=4, loc='lower right', frameon=False, fontsize=8) + + # default to saving figures in current dir instead of dedicated figures dir + save_figure('gulfstream_eval', label=label, pdf=True, output_dir='../figures/') + +if __name__ == '__main__': + from argparse import ArgumentParser + parser = ArgumentParser() + parser.add_argument('-p','--pp_root', help='Path to postprocessed data (up to but not including /pp/)', required = True) + parser.add_argument('-l', '--label', help='Label to add to figure file names', type=str, default='') + args = parser.parse_args() + plot_gulf_stream(args.pp_root, args.label) diff --git a/diagnostics/physics/config.yaml b/diagnostics/physics/config.yaml index 1134591e9..fb2c97a58 100644 --- a/diagnostics/physics/config.yaml +++ b/diagnostics/physics/config.yaml @@ -1,5 +1,6 @@ figures_dir: 'figures/' glorys: '/work/acr/mom6/diagnostics/glorys/glorys_sfc.nc' +glorys_zos: '/work/acr/glorys/GLOBAL_MULTIYEAR_PHY_001_030/monthly/glorys_monthly_z_fine_*.nc' model_grid: '../data/geography/ocean_static.nc' # Variables to rename @@ -68,10 +69,13 @@ bias_max: 2.1 bias_min_trends: -1.5 bias_max_trends: 1.51 bias_step: 0.25 - ticks: [-2, -1, 0, 1, 2] - # SST Trends Settings start_year: "2005" end_year: "2019" + +# Colorbar for ssh plots +ssh_levels_min: -1.1 +ssh_levels_max: .8 +ssh_levels_step: .1 diff --git a/diagnostics/physics/plot_common.py b/diagnostics/physics/plot_common.py index 73ba2feef..ea454526f 100644 --- a/diagnostics/physics/plot_common.py +++ b/diagnostics/physics/plot_common.py @@ -55,7 +55,7 @@ def get_map_norm(cmap, levels, no_offset=True): norm = BoundaryNorm(levels, ncolors=nlev, clip=False) return cmap, norm -def annotate_skill(model, obs, ax, dim=['yh', 'xh'], x0=-98.5, y0=54, yint=4, xint=4, weights=None, cols=1, proj = ccrs.PlateCarree(), plot_lat=False,**kwargs): +def annotate_skill(model, obs, ax, dim=['yh', 'xh'], x0=-98.5, y0=54, yint=4, xint=4, weights=None, cols=1, proj = ccrs.PlateCarree(), plot_lat=False, **kwargs): """ Annotate an axis with model vs obs skill metrics """ @@ -65,6 +65,7 @@ def annotate_skill(model, obs, ax, dim=['yh', 'xh'], x0=-98.5, y0=54, yint=4, xi medae = xskillscore.median_absolute_error(model, obs, dim=dim, skipna=True) ax.text(x0, y0, f'Bias: {float(bias):2.2f}', transform=proj, **kwargs) + # Set plot_lat=True in order to plot skill along a line of latitude. Otherwise, plot along longitude if plot_lat: ax.text(x0-xint, y0, f'RMSE: {float(rmse):2.2f}', transform=proj, **kwargs) @@ -113,20 +114,20 @@ def autoextend_colorbar(ax, plot, plot_array=None, **kwargs): extend = 'neither' return ax.colorbar(plot, extend=extend, **kwargs) -def add_ticks(ax, xticks=np.arange(-100, -31, 1), yticks=np.arange(2, 61, 1), xlabelinterval=2, ylabelinterval=2, fontsize=10, **kwargs): +def add_ticks(ax, xticks=np.arange(-100, -31, 1), yticks=np.arange(2, 61, 1), xlabelinterval=2, ylabelinterval=2, fontsize=10, projection = ccrs.PlateCarree(), **kwargs): """ Add lat and lon ticks and labels to a plot axis. By default, tick at 1 degree intervals for x and y, and label every other tick. Additional kwargs are passed to LongitudeFormatter and LatitudeFormatter. """ ax.yaxis.tick_right() - ax.set_xticks(xticks, crs=ccrs.PlateCarree()) + ax.set_xticks(xticks, crs = projection) if xlabelinterval == 0: plt.setp(ax.get_xticklabels(), visible=False) else: plt.setp([l for i, l in enumerate(ax.get_xticklabels()) if i % xlabelinterval != 0], visible=False, fontsize=fontsize) plt.setp([l for i, l in enumerate(ax.get_xticklabels()) if i % xlabelinterval == 0], fontsize=fontsize) - ax.set_yticks(yticks, crs=ccrs.PlateCarree()) + ax.set_yticks(yticks, crs = projection) if ylabelinterval == 0: plt.setp(ax.get_yticklabels(), visible=False) else: diff --git a/diagnostics/physics/ssh_eval.py b/diagnostics/physics/ssh_eval.py index 149d74782..8ee5f6e9b 100644 --- a/diagnostics/physics/ssh_eval.py +++ b/diagnostics/physics/ssh_eval.py @@ -1,163 +1,72 @@ """ Compare model SSH with 1993--2019 data from GLORYS. -Includes a plot of the Gulf Stream position and index, -and a comparison of time mean SSH. Uses whatever model data can be found within the directory pp_root, and does not try to match the model and observed time periods. How to use: -python ssh_eval.py /archive/acr/fre/NWA/2023_04/NWA12_COBALT_2023_04_kpo4-coastatten-physics/gfdl.ncrc5-intel22-prod +python ssh_eval.py -p /archive/acr/fre/NWA/2023_04/NWA12_COBALT_2023_04_kpo4-coastatten-physics/gfdl.ncrc5-intel22-prod -c config.yaml """ -import cartopy.feature as feature import cartopy.crs as ccrs from cartopy.mpl.geoaxes import GeoAxes -import matplotlib.gridspec as gridspec -from matplotlib.lines import Line2D import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import AxesGrid -import numpy as np -import pandas as pd import xarray import xesmf +import numpy as np +import logging -from plot_common import autoextend_colorbar, corners, get_map_norm, open_var, add_ticks, annotate_skill, save_figure - - - -_LAND_50M = feature.NaturalEarthFeature( - 'physical', 'land', '50m', - edgecolor='face', - facecolor='#999999' -) - -PC = ccrs.PlateCarree() - -def compute_gs(ssh, data_grid=None): - lons = np.arange(360-72, 360-51.9, 1) - lats = np.arange(36, 42, 0.1) - target_grid = {'lat': lats, 'lon': lons} - - if data_grid is None: - data_grid = {'lat': ssh.lat, 'lon': ssh.lon} - - ssh_to_grid = xesmf.Regridder( - data_grid, - target_grid, - method='bilinear' - ) - - # Interpolate the SSH data onto the index grid. - regridded = ssh_to_grid(ssh) - - # Find anomalies relative to the calendar month mean SSH over the full model run. - anom = regridded.groupby('time.month') - regridded.groupby('time.month').mean('time') - - # For each longitude point, the Gulf Stream is located at the latitude with the maximum SSH anomaly variance. - stdev = anom.std('time') - amax = stdev.argmax('lat').compute() - gs_points = stdev.lat.isel(lat=amax).compute() - - # The index is the mean latitude of the Gulf Stream, divided by the standard deviation of the mean latitude of the Gulf Stream. - index = ((anom.isel(lat=amax).mean('lon')) / anom.isel(lat=amax).mean('lon').std('time')).compute() +from plot_common import autoextend_colorbar, corners, get_map_norm, open_var, add_ticks, annotate_skill, save_figure, load_config - # Move times to the beginning of the month to match observations. - monthly_index = index.to_pandas().resample('1MS').first() - return monthly_index, gs_points +# Configure logging for ssh_eval +logger = logging.getLogger(__name__) +logging.basicConfig(filename="ssh_eval.log", format='%(asctime)s %(levelname)s:%(name)s: %(message)s',level=logging.INFO) -def plot_ssh_eval(pp_root, label): +def plot_ssh_eval(pp_root, config, label): # Ideally would use SSH, but some diag_tables only saved zos try: model_ssh = open_var(pp_root, 'ocean_monthly', 'ssh') except: print('Using zos') model_ssh = open_var(pp_root, 'ocean_monthly', 'zos') - model_thetao = open_var(pp_root, 'ocean_monthly_z', 'thetao') - - if '01_l' in model_thetao.coords: - model_thetao = model_thetao.rename({'01_l': 'z_l'}) - model_grid = xarray.open_dataset('../data/geography/ocean_static.nc') + model_grid = xarray.open_dataset( config['model_grid'] ) model_ssh_ave = model_ssh.mean('time') - model_t200 = model_thetao.interp(z_l=200).mean('time') + logger.info("MODEL_GRID: %s",model_grid) + logger.info("MODEL_SSH_AVE: %s",model_ssh_ave) + logger.info("Successfully opened model grid and took mean over time") - glorys_zos = xarray.open_mfdataset('/work/acr/glorys/GLOBAL_MULTIYEAR_PHY_001_030/monthly/glorys_monthly_z_fine_*.nc').zos + # Verify that xh/yh are set as coordinates, then make sure model coordinates match grid data + model_grid = model_grid.assign_coords( {'xh':model_grid.xh, 'yh':model_grid.yh } ) + model_ssh_ave = xarray.align(model_grid, model_ssh_ave,join='override')[1] + + glorys_zos = xarray.open_mfdataset( config['glorys_zos'] ).zos glorys_zos = glorys_zos.rename({'latitude': 'lat', 'longitude': 'lon'}) glorys_zos_ave = glorys_zos.mean('time').load() - glorys_t200 = xarray.open_dataarray('../data/diagnostics/glorys_T200.nc') glorys_lonc, glorys_latc = corners(glorys_zos.lon, glorys_zos.lat) + logger.info("GLORYS_ZOS_AVE: %s",glorys_zos_ave) - #satellite = xarray.open_mfdataset([f'/net2/acr/altimetry/SEALEVEL_GLO_PHY_L4_MY_008_047/adt_{y}_{m:02d}.nc' for y in range(1993, 2020) for m in range(1, 13)]) - #satellite = satellite.rename({'longitude': 'lon', 'latitude': 'lat'}) - #satellite = satellite.resample(time='1MS').mean('time') - - model_ssh_index, model_ssh_points = compute_gs( - model_ssh, - data_grid=model_grid[['geolon', 'geolat']].rename({'geolon': 'lon', 'geolat': 'lat'}) - ) - #satellite_ssh_index, satellite_ssh_points = compute_gs(satellite['adt']) - #satellite_ssh_points.to_netcdf('../data/obs/satellite_ssh_points.nc') - #satellite_ssh_index.to_pickle('../data/obs/satellite_ssh_index.pkl') - #read pre-calculate satellite_ssh_index and points - satellite_ssh_points = xarray.open_dataset('../data/obs/satellite_ssh_points.nc') - satellite_ssh_index = pd.read_pickle('../data/obs/satellite_ssh_index.pkl') - - model_rolled = model_ssh_index.rolling(25, center=True, min_periods=25).mean().dropna() - satellite_rolled = satellite_ssh_index.rolling(25, center=True, min_periods=25).mean().dropna() - - corr = pd.concat((model_ssh_index, satellite_ssh_index), axis=1).corr().iloc[0, 1] - corr_rolled = pd.concat((model_rolled, satellite_rolled), axis=1).corr().iloc[0, 1] - - # Plot of Gulf Stream position and index based on SSH, - # plus position based on T200. - fig = plt.figure(figsize=(10, 6), tight_layout=True) - gs = gridspec.GridSpec(2, 2, hspace=.25) - - ax = fig.add_subplot(gs[0, 0], projection=PC) - ax.add_feature(_LAND_50M) - ax.contour(model_grid.geolon, model_grid.geolat, model_t200, levels=[15], colors='r') - ax.contour(glorys_t200.longitude, glorys_t200.latitude, glorys_t200, levels=[15], colors='k') - add_ticks(ax, xlabelinterval=5) - ax.set_extent([-82, -50, 25, 41]) - ax.set_title('(a) Gulf Stream position based on T200') - custom_lines = [Line2D([0], [0], color=c, lw=2) for c in ['r', 'k']] - ax.legend(custom_lines, ['Model', 'GLORYS12'], loc='lower right', frameon=False) - - ax = fig.add_subplot(gs[0, 1], projection=PC) - ax.add_feature(_LAND_50M) - ax.plot(model_ssh_points.lon-360, model_ssh_points, c='r') - ax.plot(satellite_ssh_points.lon-360, satellite_ssh_points['__xarray_dataarray_variable__'], c='k') - add_ticks(ax, xlabelinterval=5) - ax.set_extent([-82, -50, 25, 41]) - ax.set_title('(b) Gulf Stream position based on SSH variance') - ax.legend(custom_lines, ['Model', 'Altimetry'], loc='lower right', frameon=False) - - ax = fig.add_subplot(gs[1, :]) - model_ssh_index.plot(ax=ax, c='#ffbbbb', label='Model') - satellite_ssh_index.plot(ax=ax, c='#bbbbbb', label=f'Altimetry (r={corr:2.2f})') - model_rolled.plot(ax=ax, c='r', label='Model rolling mean') - satellite_rolled.plot(ax=ax, c='k', label=f'Altimetry rolling mean (r={corr_rolled:2.2f})') - ax.set_title('(c) Gulf Stream index based on SSH variance') - ax.set_xlabel('') - ax.set_ylim(-3, 3) - ax.set_ylabel('Index (positive north)') - ax.legend(ncol=4, loc='lower right', frameon=False, fontsize=8) - - save_figure('gulfstream_eval', label=label, pdf=True) + # Set projection of each grid in the plot + # For now, sst_eval.py will only support a projection for the arctic and a projection for all other domains + if config['projection_grid'] == 'NorthPolarStereo': + p = ccrs.NorthPolarStereo() + else: + p = ccrs.PlateCarree() # Plot of time-mean SSH. fig = plt.figure(figsize=(6, 8)) grid = AxesGrid(fig, 111, nrows_ncols=(2, 1), - axes_class = (GeoAxes, dict(projection=PC)), + axes_class = ( GeoAxes, dict( projection = p ) ), axes_pad=0.55, cbar_location='right', cbar_mode='edge', cbar_pad=0.5, cbar_size='5%', - label_mode='' + label_mode='keep' ) + logger.info("Successfully created grid") - levels = np.arange(-1.1, .8, .1) + levels = np.arange( config['ssh_levels_min'], config['ssh_levels_max'], config['ssh_levels_step']) try: import colorcet except ModuleNotFoundError: @@ -166,22 +75,34 @@ def plot_ssh_eval(pp_root, label): cm = 'cet_CET_L8' cmap, norm = get_map_norm(cm, levels=levels) - xl, xr = -90, -35 - yb, yt = 20, 50 + # Set projection of input data files so that data is correctly tranformed when plotting + # For now, sst_eval.py will only support a projection for the arctic and a projection for all other domains + if config['projection_data'] == 'NorthPolarStereo': + proj = ccrs.NorthPolarStereo() + else: + proj = ccrs.PlateCarree() target_grid = model_grid[['geolon', 'geolat']].rename({'geolon': 'lon', 'geolat': 'lat'}) glorys_to_mom = xesmf.Regridder(glorys_zos_ave, target_grid, method='bilinear', unmapped_to_nan=True) glorys_rg = glorys_to_mom(glorys_zos_ave) - mod_mask = (model_grid.geolon >= xl) & (model_grid.geolon <= xr) & (model_grid.geolat >= yb) & (model_grid.geolat <= yt) + # TODO: Find more robust solution for model_grids using [-180,180] coordinates + mod_mask = ( (model_grid.geolon >= (config['lon']['west'] - 360.0)) + & (model_grid.geolon <= (config['lon']['east'] - 360.0)) + & (model_grid.geolat >= config['lat']['south']) + & (model_grid.geolat <= config['lat']['north']) + ) - p = grid[0].pcolormesh(model_grid.geolon_c, model_grid.geolat_c, model_ssh_ave, cmap=cmap, norm=norm) - cbar = autoextend_colorbar(grid.cbar_axes[0], p) + # MODEL + p0 = grid[0].pcolormesh(model_grid.geolon_c, model_grid.geolat_c, model_ssh_ave, cmap=cmap, norm=norm, transform = proj ) + cbar = autoextend_colorbar(grid.cbar_axes[0], p0) cbar.ax.set_title('SSH (m)', fontsize=10) grid[0].set_title('(a) Model mean SSH') + logger.info("Successfully plotted model data") - p = grid[1].pcolormesh(glorys_lonc, glorys_latc, glorys_zos_ave, cmap=cmap, norm=norm) - cbar = autoextend_colorbar(grid.cbar_axes[1], p) + # GLORYS + p1 = grid[1].pcolormesh(glorys_lonc, glorys_latc, glorys_zos_ave, cmap=cmap, norm=norm, transform = proj ) + cbar = autoextend_colorbar(grid.cbar_axes[1], p1) cbar.ax.set_title('SSH (m)', fontsize=10) grid[1].set_title('(b) GLORYS12 mean SSH') annotate_skill( @@ -189,21 +110,33 @@ def plot_ssh_eval(pp_root, label): glorys_rg.where(mod_mask), grid[1], weights=model_grid.areacello.where(mod_mask).fillna(0), - x0=-89, - y0=48, - yint=2, - fontsize=8 + fontsize=8, + x0=config['text_x'], + y0=config['text_y'], + xint=config['text_xint'], + plot_lat=config['plot_lat'] ) + logger.info("Successfully plotted glorys") for ax in grid: - add_ticks(ax, xlabelinterval=2, ylabelinterval=2, xticks=np.arange(-100, -31, 5), yticks=np.arange(0, 61, 5)) - ax.set_xlim(xl, xr) - ax.set_ylim(yb, yt) - ax.set_xlabel('') - ax.set_ylabel('') - ax.set_facecolor('#bbbbbb') - for s in ax.spines.values(): - s.set_visible(False) + # The set_xticks function called by add_ticks does not currently support non-rectangular projections. + # Until support is added, simply skip function if non-rectangular projection is used. Assume + # NorthPolarStereo will be the only non-rectangular projection inputted + if config['projection_grid'] != 'NorthPolarStereo': + add_ticks(ax, xlabelinterval=2, ylabelinterval=2, + xticks=np.arange( config['lon']['west'], config['lon']['east'], 5 ), + yticks=np.arange( config['lat']['south'], config['lat']['north'], 5 ), + projection = proj + ) + else: + logger.warning("WARNING: Cannot set tick marks in non-rectangular projection. Will not call add_ticks") + + ax.set_extent([ config['x']['min'], config['x']['max'], config['y']['min'], config['y']['max'] ], crs=proj ) + ax.set_xlabel('') + ax.set_ylabel('') + ax.set_facecolor('#bbbbbb') + for s in ax.spines.values(): + s.set_visible(False) save_figure('mean_ssh_eval', label=label) @@ -211,7 +144,9 @@ def plot_ssh_eval(pp_root, label): if __name__ == '__main__': from argparse import ArgumentParser parser = ArgumentParser() - parser.add_argument('pp_root', help='Path to postprocessed data (up to but not including /pp/)') + parser.add_argument('-p','--pp_root', help='Path to postprocessed data (up to but not including /pp/)', required = True) + parser.add_argument('-c','--config', help='Path to config.yaml file', required = True) parser.add_argument('-l', '--label', help='Label to add to figure file names', type=str, default='') args = parser.parse_args() - plot_ssh_eval(args.pp_root, args.label) + config = load_config( args.config ) + plot_ssh_eval(args.pp_root, config, args.label)