Demonstrates plotting the tidal amplitudes around Antarctica
OTIS format tidal solutions provided by Oregon State University and ESR
Global Tide Model (GOT) solutions provided by Richard Ray at GSFC
Finite Element Solution (FES) provided by AVISO
crs.py
: Coordinate Reference System (CRS) routinesio.model.py
: retrieves tide model parameters for named tide modelsio.OTIS.py
: extract tidal harmonic constants from OTIS tide modelsio.ATLAS.py
: extract tidal harmonic constants from ATLAS netcdf modelsio.FES.py
: extract tidal harmonic constants from FES tide modelsThis notebook uses Jupyter widgets to set parameters for calculating the tidal maps.
import os
import pyproj
import numpy as np
import matplotlib
matplotlib.rcParams['axes.linewidth'] = 2.0
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# import tide programs
import pyTMD.io
import pyTMD.tools
# autoreload
%load_ext autoreload
%autoreload 2
# available model list
model_list = sorted(pyTMD.io.model.ocean_elevation())
# display widgets for setting directory and model
TMDwidgets = pyTMD.tools.widgets()
TMDwidgets.model.options = model_list
TMDwidgets.model.value = 'GOT4.10'
TMDwidgets.VBox([
TMDwidgets.directory,
TMDwidgets.model,
TMDwidgets.atlas,
TMDwidgets.compress
])
# get model parameters
model = pyTMD.io.model(TMDwidgets.directory.value,
format=TMDwidgets.atlas.value,
compressed=TMDwidgets.compress.value
).elevation(TMDwidgets.model.value)
# create an image around Antarctica
xlimits = [-560.*5e3,560.*5e3]
ylimits = [-560.*5e3,560.*5e3]
spacing = [20e3,-20e3]
# x and y coordinates
x = np.arange(xlimits[0],xlimits[1]+spacing[0],spacing[0])
y = np.arange(ylimits[1],ylimits[0]+spacing[1],spacing[1])
xgrid,ygrid = np.meshgrid(x,y)
# x and y dimensions
nx = int((xlimits[1]-xlimits[0])/spacing[0])+1
ny = int((ylimits[0]-ylimits[1])/spacing[1])+1
# convert image coordinates from polar stereographic to latitude/longitude
crs1 = pyproj.CRS.from_epsg(3031)
crs2 = pyproj.CRS.from_epsg(4326)
transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)
lon,lat = transformer.transform(xgrid.flatten(), ygrid.flatten())
def infer_minor_amplitudes(zmajor,constituents):
# number of constituents
npts,nc = np.shape(zmajor)
cindex = ['q1','o1','p1','k1','n2','m2','s2','k2']
# re-order zmajor to correspond to cindex
z8 = np.ma.zeros((npts,8))
ni = 0
for i,c in enumerate(cindex):
j = [j for j,val in enumerate(constituents) if val == c]
if j:
j1, = j
z8[:,i] = zmajor[:,j1]
ni += 1
# list of minor constituents
minor = ['2q1','sigma1','rho1','m1','m1','chi1','pi1','phi1','theta1',
'j1','oo1','2n2','mu2','nu2','lambda2','l2','l2','t2']
# only add minor constituents that are not on the list of major values
minor_flag = [m not in constituents for m in minor]
# estimate minor constituents
zmin = np.zeros((npts,18))
zmin[:,0] = 0.263*z8[:,0] - 0.0252*z8[:,1]# 2Q1
zmin[:,1] = 0.297*z8[:,0] - 0.0264*z8[:,1]# sigma1
zmin[:,2] = 0.164*z8[:,0] + 0.0048*z8[:,1]# rho1
zmin[:,3] = 0.0140*z8[:,1] + 0.0101*z8[:,3]# M1
zmin[:,4] = 0.0389*z8[:,1] + 0.0282*z8[:,3]# M1
zmin[:,5] = 0.0064*z8[:,1] + 0.0060*z8[:,3]# chi1
zmin[:,6] = 0.0030*z8[:,1] + 0.0171*z8[:,3]# pi1
zmin[:,7] = -0.0015*z8[:,1] + 0.0152*z8[:,3]# phi1
zmin[:,8] = -0.0065*z8[:,1] + 0.0155*z8[:,3]# theta1
zmin[:,9] = -0.0389*z8[:,1] + 0.0836*z8[:,3]# J1
zmin[:,10] = -0.0431*z8[:,1] + 0.0613*z8[:,3]# OO1
zmin[:,11] = 0.264*z8[:,4] - 0.0253*z8[:,5]# 2N2
zmin[:,12] = 0.298*z8[:,4] - 0.0264*z8[:,5]# mu2
zmin[:,13] = 0.165*z8[:,4] + 0.00487*z8[:,5]# nu2
zmin[:,14] = 0.0040*z8[:,5] + 0.0074*z8[:,6]# lambda2
zmin[:,15] = 0.0131*z8[:,5] + 0.0326*z8[:,6]# L2
zmin[:,16] = 0.0033*z8[:,5] + 0.0082*z8[:,6]# L2
zmin[:,17] = 0.0585*z8[:,6]# t2
# only add minor constituents that are not major
return np.where(minor_flag, np.abs(zmin), 0.0)
# read tidal constants and interpolate to grid points
if model.format in ('OTIS','ATLAS-compact','TMD3'):
amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,
model.model_file, model.projection, type=model.type, crop=True,
method='spline', grid=model.file_format)
elif (model.format == 'ATLAS-netcdf'):
amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file,
model.model_file, type=model.type, crop=True, method='spline',
scale=model.scale, compressed=model.compressed)
elif model.format in ('GOT-ascii', 'GOT-netcdf'):
amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file,
grid=model.file_format, crop=True, method='spline',
scale=model.scale, compressed=model.compressed)
elif (model.format == 'FES-netcdf'):
amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file,
type=model.type, version=model.version, crop=True,
method='spline', scale=model.scale, compressed=model.compressed)
c = model.constituents
# calculate minor constituent amplitudes
minor_amp = infer_minor_amplitudes(amp,c)
# calculate sum of major and minor amplitudes
tide_range = np.sum(amp,axis=1) + np.sum(minor_amp,axis=1)
# convert from meters to centimeters
tide_cm = 100.0*np.reshape(tide_range,(ny,nx))
# plot Antarctic tide range
projection = ccrs.Stereographic(central_longitude=0.0,
central_latitude=-90,true_scale_latitude=-71.0)
fig, ax = plt.subplots(num=1, figsize=(9,8),
subplot_kw=dict(projection=projection))
# plot tide height
vmin,vmax = (0, np.max(tide_cm))
extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])
im = ax.imshow(tide_cm, interpolation='nearest',
vmin=vmin, vmax=vmax, transform=projection,
extent=extent, origin='upper')
# add high resolution cartopy coastlines
ax.coastlines('10m')
# Add colorbar and adjust size
# pad = distance from main plot axis
# extend = add extension triangles to upper and lower bounds
# options: neither, both, min, max
# shrink = percent size of colorbar
# aspect = lengthXwidth aspect of colorbar
cbar = plt.colorbar(im, ax=ax, pad=0.025, extend='max',
extendfrac=0.0375, shrink=0.85, aspect=22.5, drawedges=False)
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
cbar.ax.set_ylabel(f'{model.name} Tide Range', labelpad=10, fontsize=13)
cbar.ax.set_title('cm', fontsize=13, va='bottom')
# ticks lines all the way across
cbar.ax.tick_params(which='both', width=1, length=20,
labelsize=13, direction='in')
# set x and y limits
ax.set_xlim(xlimits)
ax.set_ylim(ylimits)
# stronger linewidth on frame
ax.spines['geo'].set_linewidth(2.0)
ax.spines['geo'].set_capstyle('projecting')
# adjust subplot within figure
fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)
# show the plot
plt.show()