Plot Antarctic Tide Range

This (notebook) demonstrates plotting the total combined tidal amplitude for all constituents 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

Python Dependencies

Program Dependencies

  • crs.py: Coordinate Reference System (CRS) routines

  • io.model.py: retrieves tide model parameters for named tide models

  • io.OTIS.py: extract tidal harmonic constants from OTIS tide models

  • io.ATLAS.py: extract tidal harmonic constants from ATLAS netcdf models

  • io.GOT.py: extract tidal harmonic constants from GOT tide models

  • io.FES.py: extract tidal harmonic constants from FES tide models

This notebook uses Jupyter widgets to set parameters for calculating the tidal maps.

Load modules

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

Set parameters for program

  • Model directory

  • Tide model

# 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.compress
])

Setup tide model parameters

# get model parameters
model = pyTMD.io.model(TMDwidgets.directory.value,
    compressed=TMDwidgets.compress.value
   ).elevation(TMDwidgets.model.value)

Setup coordinates for calculating tides

# 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())

Infer minor amplitudes from the major constituents

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)

Calculate tide map

# 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))

Create plot of tidal range

# 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()