Demonstrates plotting hourly tidal currents around Antarctica
OTIS format tidal solutions provided by Ohio State University and ESR
Finite Element Solution (FES) provided by AVISO
calc_astrol_longitudes.py
: computes the basic astronomical mean longitudescalc_delta_time.py
: calculates difference between universal and dynamic timeconvert_ll_xy.py
: convert lat/lon points to and from projected coordinatesload_constituent.py
: loads parameters for a given tidal constituentload_nodal_corrections.py
: load the nodal corrections for tidal constituentsinfer_minor_corrections.py
: return corrections for minor constituentsread_tide_model.py
: extract tidal harmonic constants from OTIS tide modelsread_netcdf_model.py
: extract tidal harmonic constants from netcdf modelsread_FES_model.py
: extract tidal harmonic constants from FES tide modelspredict_tide.py
: predict tidal elevation at a single time using harmonic constantsThis notebook uses Jupyter widgets to set parameters for calculating the tidal maps.
The widgets can be installed as described below.
pip3 install --user ipywidgets
jupyter nbextension install --user --py widgetsnbextension
jupyter nbextension enable --user --py widgetsnbextension
jupyter-notebook
import os
import pyproj
import datetime
import numpy as np
import matplotlib
matplotlib.rcParams['axes.linewidth'] = 2.0
matplotlib.rcParams["animation.html"] = "jshtml"
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import cartopy.crs as ccrs
import ipywidgets as widgets
from IPython.display import HTML
import pyTMD.time
from pyTMD.calc_delta_time import calc_delta_time
from pyTMD.read_tide_model import extract_tidal_constants
from pyTMD.read_netcdf_model import extract_netcdf_constants
from pyTMD.read_GOT_model import extract_GOT_constants
from pyTMD.read_FES_model import extract_FES_constants
from pyTMD.infer_minor_corrections import infer_minor_corrections
from pyTMD.predict_tide import predict_tide
#-- autoreload
%load_ext autoreload
%autoreload 2
#-- set the directory with tide models
dirText = widgets.Text(
value=os.getcwd(),
description='Directory:',
disabled=False
)
#-- dropdown menu for setting tide model
model_list = ['CATS0201','CATS2008','TPXO9-atlas','TPXO9-atlas-v2','TPXO9.1',
'TPXO8-atlas','TPXO7.2','FES2014']
modelDropdown = widgets.Dropdown(
options=model_list,
value='CATS2008',
description='Model:',
disabled=False,
)
#-- date picker widget for setting time
datepick = widgets.DatePicker(
description='Date:',
value = datetime.date.today(),
disabled=False
)
#-- display widgets for setting directory, model and date
widgets.VBox([dirText,modelDropdown,datepick])
#-- directory with tide models
tide_dir = os.path.expanduser(dirText.value)
MODEL = modelDropdown.value
if (MODEL == 'CATS0201'):
grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS')
model_file = os.path.join(tide_dir,'cats0201_tmd','UV0_CATS02_01')
model_format = 'OTIS'
EPSG = '4326'
TYPES = ['u','v']
elif (MODEL == 'CATS2008'):
grid_file = os.path.join(tide_dir,'CATS2008','grid_CATS2008')
model_file = os.path.join(tide_dir,'CATS2008','uv.CATS2008.out')
model_format = 'OTIS'
EPSG = 'CATS2008'
TYPES = ['u','v']
elif (MODEL == 'TPXO9-atlas'):
model_directory = os.path.join(tide_dir,'TPXO9_atlas')
grid_file = 'grid_tpxo9_atlas.nc.gz'
model_files = {}
model_files['u'] = ['u_q1_tpxo9_atlas_30.nc.gz','u_o1_tpxo9_atlas_30.nc.gz',
'u_p1_tpxo9_atlas_30.nc.gz','u_k1_tpxo9_atlas_30.nc.gz',
'u_n2_tpxo9_atlas_30.nc.gz','u_m2_tpxo9_atlas_30.nc.gz',
'u_s2_tpxo9_atlas_30.nc.gz','u_k2_tpxo9_atlas_30.nc.gz',
'u_m4_tpxo9_atlas_30.nc.gz','u_ms4_tpxo9_atlas_30.nc.gz',
'u_mn4_tpxo9_atlas_30.nc.gz','u_2n2_tpxo9_atlas_30.nc.gz']
model_files['v'] = ['v_q1_tpxo9_atlas_30.nc.gz','v_o1_tpxo9_atlas_30.nc.gz',
'v_p1_tpxo9_atlas_30.nc.gz','v_k1_tpxo9_atlas_30.nc.gz',
'v_n2_tpxo9_atlas_30.nc.gz','v_m2_tpxo9_atlas_30.nc.gz',
'v_s2_tpxo9_atlas_30.nc.gz','v_k2_tpxo9_atlas_30.nc.gz',
'v_m4_tpxo9_atlas_30.nc.gz','v_ms4_tpxo9_atlas_30.nc.gz',
'v_mn4_tpxo9_atlas_30.nc.gz','v_2n2_tpxo9_atlas_30.nc.gz']
model_format = 'netcdf'
TYPES = ['u','v']
SCALE = 1.0/100.0
GZIP = True
elif (MODEL == 'TPXO9-atlas-v2'):
model_directory = os.path.join(tide_dir,'TPXO9_atlas_v2')
grid_file = 'grid_tpxo9_atlas_30_v2.nc.gz'
model_files = {}
model_files['u'] = ['u_q1_tpxo9_atlas_30_v2.nc.gz','u_o1_tpxo9_atlas_30_v2.nc.gz',
'u_p1_tpxo9_atlas_30_v2.nc.gz','u_k1_tpxo9_atlas_30_v2.nc.gz',
'u_n2_tpxo9_atlas_30_v2.nc.gz','u_m2_tpxo9_atlas_30_v2.nc.gz',
'u_s2_tpxo9_atlas_30_v2.nc.gz','u_k2_tpxo9_atlas_30_v2.nc.gz',
'u_m4_tpxo9_atlas_30_v2.nc.gz','u_ms4_tpxo9_atlas_30_v2.nc.gz',
'u_mn4_tpxo9_atlas_30_v2.nc.gz','u_2n2_tpxo9_atlas_30_v2.nc.gz']
model_files['v'] = ['v_q1_tpxo9_atlas_30_v2.nc.gz','v_o1_tpxo9_atlas_30_v2.nc.gz',
'v_p1_tpxo9_atlas_30_v2.nc.gz','v_k1_tpxo9_atlas_30_v2.nc.gz',
'v_n2_tpxo9_atlas_30_v2.nc.gz','v_m2_tpxo9_atlas_30_v2.nc.gz',
'v_s2_tpxo9_atlas_30_v2.nc.gz','v_k2_tpxo9_atlas_30_v2.nc.gz',
'v_m4_tpxo9_atlas_30_v2.nc.gz','v_ms4_tpxo9_atlas_30_v2.nc.gz',
'v_mn4_tpxo9_atlas_30_v2.nc.gz','v_2n2_tpxo9_atlas_30_v2.nc.gz']
model_format = 'netcdf'
TYPES = ['u','v']
SCALE = 1.0/100.0
GZIP = True
elif (MODEL == 'TPXO9.1'):
grid_file = os.path.join(tide_dir,'TPXO9.1','DATA','grid_tpxo9')
model_file = os.path.join(tide_dir,'TPXO9.1','DATA','u_tpxo9.v1')
model_format = 'OTIS'
EPSG = '4326'
TYPES = ['u','v']
elif (MODEL == 'TPXO8-atlas'):
grid_file = os.path.join(tide_dir,'tpxo8_atlas','grid_tpxo8atlas_30_v1')
model_file = os.path.join(tide_dir,'tpxo8_atlas','uv.tpxo8_atlas_30_v1')
model_format = 'ATLAS'
EPSG = '4326'
TYPES = ['u','v']
elif (MODEL == 'TPXO7.2'):
grid_file = os.path.join(tide_dir,'TPXO7.2_tmd','grid_tpxo7.2')
model_file = os.path.join(tide_dir,'TPXO7.2_tmd','u_tpxo7.2')
model_format = 'OTIS'
EPSG = '4326'
TYPES = ['u','v']
elif (MODEL == 'FES2014'):
model_directory = {}
model_directory['u'] = os.path.join(tide_dir,'fes2014','eastward_velocity')
model_directory['v'] = os.path.join(tide_dir,'fes2014','northward_velocity')
model_files = ['2n2.nc.gz','eps2.nc.gz','j1.nc.gz','k1.nc.gz',
'k2.nc.gz','l2.nc.gz','la2.nc.gz','m2.nc.gz','m3.nc.gz','m4.nc.gz',
'm6.nc.gz','m8.nc.gz','mf.nc.gz','mks2.nc.gz','mm.nc.gz',
'mn4.nc.gz','ms4.nc.gz','msf.nc.gz','msqm.nc.gz','mtm.nc.gz',
'mu2.nc.gz','n2.nc.gz','n4.nc.gz','nu2.nc.gz','o1.nc.gz','p1.nc.gz',
'q1.nc.gz','r2.nc.gz','s1.nc.gz','s2.nc.gz','s4.nc.gz','sa.nc.gz',
'ssa.nc.gz','t2.nc.gz']
c = ['2n2','eps2','j1','k1','k2','l2','lambda2','m2','m3','m4','m6',
'm8','mf','mks2','mm','mn4','ms4','msf','msqm','mtm','mu2','n2',
'n4','nu2','o1','p1','q1','r2','s1','s2','s4','sa','ssa','t2']
model_format = 'FES'
TYPES = ['u','v']
#-- 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 = np.int((xlimits[1]-xlimits[0])/spacing[0])+1
ny = np.int((ylimits[0]-ylimits[1])/spacing[1])+1
#-- convert image coordinates from polar stereographic to latitude/longitude
crs1 = pyproj.CRS.from_string("epsg:{0:d}".format(3031))
crs2 = pyproj.CRS.from_string("epsg:{0:d}".format(4326))
transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)
lon,lat = transformer.transform(xgrid.flatten(), ygrid.flatten())
#-- convert from calendar date to days relative to Jan 1, 1992 (48622 MJD)
YMD = datepick.value
tide_time = pyTMD.time.convert_calendar_dates(YMD.year, YMD.month,
YMD.day, hour=np.arange(24))
#-- save tide currents
tide = {}
#-- iterate over u and v currents
for TYPE in TYPES:
#-- read tidal constants and interpolate to grid points
if model_format in ('OTIS','ATLAS'):
amp,ph,D,c = extract_tidal_constants(lon, lat, grid_file, model_file,
EPSG, TYPE=TYPE, METHOD='spline', GRID=model_format)
DELTAT = np.zeros_like(tide_time)
elif (model_format == 'netcdf'):
amp,ph,D,c = extract_netcdf_constants(lon, lat, model_directory,
grid_file, model_files[TYPE], TYPE=TYPE, METHOD='spline',
SCALE=SCALE, GZIP=GZIP)
DELTAT = np.zeros_like(tide_time)
elif (model_format == 'GOT'):
amp,ph = extract_GOT_constants(lon, lat, model_directory, model_files,
METHOD='spline', SCALE=SCALE)
#-- interpolate delta times from calendar dates to tide time
delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
DELTAT = calc_delta_time(delta_file, tide_time)
elif (model_format == 'FES'):
amp,ph = extract_FES_constants(lon, lat, model_directory[TYPE], model_files,
TYPE=TYPE, VERSION=MODEL, METHOD=METHOD, SCALE=SCALE)
#-- interpolate delta times from calendar dates to tide time
delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
DELTAT = calc_delta_time(delta_file, tide_time)
#-- calculate complex phase in radians for Euler's
cph = -1j*ph*np.pi/180.0
#-- calculate constituent oscillation
hc = amp*np.exp(cph)
#-- allocate for tide current map calculated every hour
tide[TYPE] = np.ma.zeros((ny,nx,24))
for hour in range(24):
#-- predict tidal elevations at time and infer minor corrections
TIDE = predict_tide(tide_time[hour], hc, c, DELTAT=DELTAT[hour],
CORRECTIONS=model_format)
MINOR = infer_minor_corrections(tide_time[hour], hc, c,
DELTAT=DELTAT[hour], CORRECTIONS=model_format)
#-- add major and minor components and reform grid
tide[TYPE][:,:,hour] = np.reshape((TIDE+MINOR),(ny,nx))
#-- output Antarctic Tidal Current Animation
projection = ccrs.Stereographic(central_longitude=0.0,
central_latitude=-90.0,true_scale_latitude=-71.0)
#-- figure axis and image objects
ax1,im = ({},{})
fig, (ax1['u'],ax1['v']) = plt.subplots(num=1, ncols=2,
figsize=(11.5,7), subplot_kw=dict(projection=projection))
vmin = np.min([tide['u'].min(),tide['v'].min()])
vmax = np.max([tide['u'].max(),tide['v'].max()])
extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])
for TYPE,ax in ax1.items():
#-- plot tidal currents
im[TYPE] = ax.imshow(np.zeros((ny,nx)),
interpolation='nearest', vmin=vmin, vmax=vmax,
transform=projection, extent=extent, origin='upper',
animated=True)
#-- add high resolution cartopy coastlines
ax.coastlines('10m')
#-- 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')
#-- Add colorbar with a colorbar axis
#-- Add an axes at position rect [left, bottom, width, height]
cbar_ax = fig.add_axes([0.085, 0.075, 0.83, 0.035])
#-- extend = add extension triangles to upper and lower bounds
#-- options: neither, both, min, max
cbar = fig.colorbar(im['u'], cax=cbar_ax, extend='both',
extendfrac=0.0375, drawedges=False, orientation='horizontal')
#-- rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
#-- Add label to the colorbar
cbar.ax.set_xlabel('{0} Tidal Velocity'.format(MODEL), fontsize=13)
cbar.ax.set_ylabel('cm/s', fontsize=13, rotation=0)
cbar.ax.yaxis.set_label_coords(1.04, 0.15)
#-- ticks lines all the way across
cbar.ax.tick_params(which='both', width=1, length=18,
labelsize=13, direction='in')
#-- add title (date and time)
ttl = fig.suptitle(None, fontsize=13)
# adjust subplot within figure
fig.subplots_adjust(left=0.02,right=0.98,bottom=0.1,top=0.98,wspace=0.04)
#-- animate each map
def animate_maps(hour):
#-- set map data iterating over u and v currents
for TYPE in TYPES:
im[TYPE].set_data(tide[TYPE][:,:,hour])
#-- set title
args = (YMD.year,YMD.month,YMD.day,hour)
ttl.set_text('{0:4d}-{1:02d}-{2:02d}T{3:02d}:00:00'.format(*args))
#-- set animation
anim = animation.FuncAnimation(fig, animate_maps, frames=24)
%matplotlib inline
HTML(anim.to_jshtml())