|
| 1 | +""" |
| 2 | +Plot of the Gulf Stream position and index, |
| 3 | +Uses whatever model data can be found within the directory pp_root, |
| 4 | +and does not try to match the model and observed time periods. |
| 5 | +How to use: |
| 6 | +python plot_gulf_stream.py -p /archive/acr/fre/NWA/2023_04/NWA12_COBALT_2023_04_kpo4-coastatten-physics/gfdl.ncrc5-intel22-prod |
| 7 | +""" |
| 8 | +import xarray |
| 9 | +import xesmf |
| 10 | +import pandas as pd |
| 11 | +import numpy as np |
| 12 | +import cartopy.feature as feature |
| 13 | +import cartopy.crs as ccrs |
| 14 | +from cartopy.mpl.geoaxes import GeoAxes |
| 15 | +import matplotlib.gridspec as gridspec |
| 16 | +from matplotlib.lines import Line2D |
| 17 | +import matplotlib.pyplot as plt |
| 18 | +from mpl_toolkits.axes_grid1 import AxesGrid |
| 19 | + |
| 20 | +# Need to append physics dir to path to access plot common |
| 21 | +import sys |
| 22 | +sys.path.append("..") |
| 23 | +from plot_common import open_var, add_ticks, save_figure |
| 24 | + |
| 25 | +def compute_gs(ssh, data_grid=None): |
| 26 | + lons = np.arange(360-72, 360-51.9, 1) |
| 27 | + lats = np.arange(36, 42, 0.1) |
| 28 | + target_grid = {'lat': lats, 'lon': lons} |
| 29 | + |
| 30 | + if data_grid is None: |
| 31 | + data_grid = {'lat': ssh.lat, 'lon': ssh.lon} |
| 32 | + |
| 33 | + ssh_to_grid = xesmf.Regridder( |
| 34 | + data_grid, |
| 35 | + target_grid, |
| 36 | + method='bilinear' |
| 37 | + ) |
| 38 | + |
| 39 | + # Interpolate the SSH data onto the index grid. |
| 40 | + regridded = ssh_to_grid(ssh) |
| 41 | + |
| 42 | + # Find anomalies relative to the calendar month mean SSH over the full model run. |
| 43 | + anom = regridded.groupby('time.month') - regridded.groupby('time.month').mean('time') |
| 44 | + |
| 45 | + # For each longitude point, the Gulf Stream is located at the latitude with the maximum SSH anomaly variance. |
| 46 | + stdev = anom.std('time') |
| 47 | + amax = stdev.argmax('lat').compute() |
| 48 | + gs_points = stdev.lat.isel(lat=amax).compute() |
| 49 | + |
| 50 | + # The index is the mean latitude of the Gulf Stream, divided by the standard deviation of the mean latitude of the Gulf Stream. |
| 51 | + index = ((anom.isel(lat=amax).mean('lon')) / anom.isel(lat=amax).mean('lon').std('time')).compute() |
| 52 | + |
| 53 | + # Move times to the beginning of the month to match observations. |
| 54 | + monthly_index = index.to_pandas().resample('1MS').first() |
| 55 | + return monthly_index, gs_points |
| 56 | + |
| 57 | +def plot_gulf_stream(pp_root, label): |
| 58 | + |
| 59 | + # Load Natural Earth Shapefiles |
| 60 | + _LAND_50M = feature.NaturalEarthFeature( |
| 61 | + 'physical', 'land', '50m', |
| 62 | + edgecolor='face', |
| 63 | + facecolor='#999999' |
| 64 | + ) |
| 65 | + |
| 66 | + # Get model grid |
| 67 | + model_grid = xarray.open_dataset( '../../data/geography/ocean_static.nc' ) |
| 68 | + |
| 69 | + # Get model thetao data TODO: maki this comment better |
| 70 | + model_thetao = open_var(pp_root, 'ocean_monthly_z', 'thetao') |
| 71 | + |
| 72 | + if '01_l' in model_thetao.coords: |
| 73 | + model_thetao = model_thetao.rename({'01_l': 'z_l'}) |
| 74 | + |
| 75 | + model_t200 = model_thetao.interp(z_l=200).mean('time') |
| 76 | + |
| 77 | + # Ideally would use SSH, but some diag_tables only saved zos |
| 78 | + try: |
| 79 | + model_ssh = open_var(pp_root, 'ocean_monthly', 'ssh') |
| 80 | + except: |
| 81 | + print('Using zos') |
| 82 | + model_ssh = open_var(pp_root, 'ocean_monthly', 'zos') |
| 83 | + |
| 84 | + model_ssh_index, model_ssh_points = compute_gs( |
| 85 | + model_ssh, |
| 86 | + data_grid=model_grid[['geolon', 'geolat']].rename({'geolon': 'lon', 'geolat': 'lat'}) |
| 87 | + ) |
| 88 | + |
| 89 | + # Get Glorys data |
| 90 | + glorys_t200 = xarray.open_dataarray('../../data/diagnostics/glorys_T200.nc') |
| 91 | + |
| 92 | + # Get satellite points |
| 93 | + #satellite_ssh_index, satellite_ssh_points = compute_gs(satellite['adt']) |
| 94 | + #satellite_ssh_points.to_netcdf('../data/obs/satellite_ssh_points.nc') |
| 95 | + #satellite_ssh_index.to_pickle('../data/obs/satellite_ssh_index.pkl') |
| 96 | + #read pre-calculate satellite_ssh_index and points |
| 97 | + satellite_ssh_points = xarray.open_dataset('../../data/obs/satellite_ssh_points.nc') |
| 98 | + satellite_ssh_index = pd.read_pickle('../../data/obs/satellite_ssh_index.pkl') |
| 99 | + satellite_rolled = satellite_ssh_index.rolling(25, center=True, min_periods=25).mean().dropna() |
| 100 | + |
| 101 | + #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)]) |
| 102 | + #satellite = satellite.rename({'longitude': 'lon', 'latitude': 'lat'}) |
| 103 | + #satellite = satellite.resample(time='1MS').mean('time') |
| 104 | + |
| 105 | + # Get rolling averages and correlations |
| 106 | + model_rolled = model_ssh_index.rolling(25, center=True, min_periods=25).mean().dropna() |
| 107 | + corr = pd.concat((model_ssh_index, satellite_ssh_index), axis=1).corr().iloc[0, 1] |
| 108 | + corr_rolled = pd.concat((model_rolled, satellite_rolled), axis=1).corr().iloc[0, 1] |
| 109 | + |
| 110 | + # Plot of Gulf Stream position and index based on SSH, |
| 111 | + # plus position based on T200. |
| 112 | + fig = plt.figure(figsize=(10, 6), tight_layout=True) |
| 113 | + gs = gridspec.GridSpec(2, 2, hspace=.25) |
| 114 | + |
| 115 | + # Set projection of input data files so that data is correctly tranformed when plotting |
| 116 | + proj = ccrs.PlateCarree() |
| 117 | + |
| 118 | + ax = fig.add_subplot(gs[0, 0], projection = proj) |
| 119 | + ax.add_feature(_LAND_50M) |
| 120 | + ax.contour(model_grid.geolon, model_grid.geolat, model_t200, levels=[15], colors='r') |
| 121 | + ax.contour(glorys_t200.longitude, glorys_t200.latitude, glorys_t200, levels=[15], colors='k') |
| 122 | + add_ticks(ax, xlabelinterval=5) |
| 123 | + ax.set_extent([-82, -50, 25, 41]) |
| 124 | + ax.set_title('(a) Gulf Stream position based on T200') |
| 125 | + custom_lines = [Line2D([0], [0], color=c, lw=2) for c in ['r', 'k']] |
| 126 | + ax.legend(custom_lines, ['Model', 'GLORYS12'], loc='lower right', frameon=False) |
| 127 | + |
| 128 | + ax = fig.add_subplot(gs[0, 1], projection = proj) |
| 129 | + ax.add_feature(_LAND_50M) |
| 130 | + ax.plot(model_ssh_points.lon-360, model_ssh_points, c='r') |
| 131 | + ax.plot(satellite_ssh_points.lon-360, satellite_ssh_points['__xarray_dataarray_variable__'], c='k') |
| 132 | + add_ticks(ax, xlabelinterval=5) |
| 133 | + ax.set_extent([-82, -50, 25, 41]) |
| 134 | + ax.set_title('(b) Gulf Stream position based on SSH variance') |
| 135 | + ax.legend(custom_lines, ['Model', 'Altimetry'], loc='lower right', frameon=False) |
| 136 | + |
| 137 | + ax = fig.add_subplot(gs[1, :]) |
| 138 | + model_ssh_index.plot(ax=ax, c='#ffbbbb', label='Model') |
| 139 | + satellite_ssh_index.plot(ax=ax, c='#bbbbbb', label=f'Altimetry (r={corr:2.2f})') |
| 140 | + model_rolled.plot(ax=ax, c='r', label='Model rolling mean') |
| 141 | + satellite_rolled.plot(ax=ax, c='k', label=f'Altimetry rolling mean (r={corr_rolled:2.2f})') |
| 142 | + ax.set_title('(c) Gulf Stream index based on SSH variance') |
| 143 | + ax.set_xlabel('') |
| 144 | + ax.set_ylim(-3, 3) |
| 145 | + ax.set_ylabel('Index (positive north)') |
| 146 | + ax.legend(ncol=4, loc='lower right', frameon=False, fontsize=8) |
| 147 | + |
| 148 | + # default to saving figures in current dir instead of dedicated figures dir |
| 149 | + save_figure('gulfstream_eval', label=label, pdf=True, output_dir='../figures/') |
| 150 | + |
| 151 | +if __name__ == '__main__': |
| 152 | + from argparse import ArgumentParser |
| 153 | + parser = ArgumentParser() |
| 154 | + parser.add_argument('-p','--pp_root', help='Path to postprocessed data (up to but not including /pp/)', required = True) |
| 155 | + parser.add_argument('-l', '--label', help='Label to add to figure file names', type=str, default='') |
| 156 | + args = parser.parse_args() |
| 157 | + plot_gulf_stream(args.pp_root, args.label) |
0 commit comments