Plot Arctic Ocean Map
This (notebook
) demonstrates creating an animation of hourly tidal elevations for the Arctic Ocean
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
arguments.py
: load the nodal corrections for tidal constituentsastro.py
: computes the basic astronomical mean longitudescrs.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 modelsio.GOT.py
: extract tidal harmonic constants from GOT tide modelspredict.py
: predict tidal values using harmonic constantstime.py
: utilities for calculating time operationsutilities.py
: download and management utilities for files
This notebook uses Jupyter widgets to set parameters for calculating the tidal maps.
Load modules
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
from IPython.display import HTML
# import tide programs
import pyTMD.io
import pyTMD.predict
import pyTMD.tools
import pyTMD.utilities
import timescale.time
# autoreload
%load_ext autoreload
%autoreload 2
Set parameters for program
Model directory
Tide model
Date to run
# 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 = 'TPXO8-atlas'
TMDwidgets.VBox([
TMDwidgets.directory,
TMDwidgets.model,
TMDwidgets.compress,
TMDwidgets.datepick
])
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 the Arctic Ocean
# use NSIDC Polar Stereographic definitions
# https://nsidc.org/data/polar-stereo/ps_grids.html
xlimits = [-3850000,3750000]
ylimits = [-5350000,5850000]
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(3413)
crs2 = pyproj.CRS.from_epsg(4326)
transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)
lon,lat = transformer.transform(xgrid.flatten(), ygrid.flatten())
Calculate tide map
# convert from calendar date to days relative to Jan 1, 1992 (48622 MJD)
YMD = TMDwidgets.datepick.value
ts = timescale.time.Timescale().from_calendar(YMD.year, YMD.month,
YMD.day, hour=np.arange(24))
# 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)
DELTAT = np.zeros_like(ts.tide)
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)
DELTAT = np.zeros_like(ts.tide)
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)
# delta time (TT - UT1)
DELTAT = ts.tt_ut1
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
# delta time (TT - UT1)
DELTAT = ts.tt_ut1
# 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 map calculated every hour
tide_cm = np.ma.zeros((ny,nx,24))
for hour in range(24):
# predict tidal elevations at time and infer minor corrections
TIDE = pyTMD.predict.map(ts.tide[hour], hc, c, deltat=DELTAT[hour],
corrections=model.corrections)
MINOR = pyTMD.predict.infer_minor(ts.tide[hour], hc, c,
deltat=DELTAT[hour], corrections=model.corrections)
# add major and minor components and reform grid
# convert from meters to centimeters
tide_cm[:,:,hour] = 100.0*np.reshape((TIDE+MINOR),(ny,nx))
Create animation of hourly tidal oscillation
# output Arctic Ocean Tide Animation
projection = ccrs.Stereographic(central_longitude=-45.0,
central_latitude=+90.0,true_scale_latitude=+70.0)
fig, ax = plt.subplots(num=1, figsize=(8,9),
subplot_kw=dict(projection=projection))
# plot tide height
vmin,vmax = (np.min(tide_cm), np.max(tide_cm))
extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])
im = ax.imshow(np.zeros((ny,nx)), interpolation='nearest',
vmin=vmin, vmax=vmax, transform=projection,
extent=extent, origin='upper', animated=True)
# add 50m resolution cartopy coastlines
ax.coastlines('50m')
# 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='both',
extendfrac=0.0375, shrink=0.90, aspect=25.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 Height', 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=21,
labelsize=13, direction='in')
# add title (date and time)
ttl = ax.set_title(None, fontsize=13)
# 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.95)
# animate each map
def animate_maps(hour):
# set map data
im.set_data(tide_cm[:,:,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())