Plot Antarctic Tidal Currents¶

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

# available model list
model_list = sorted(pyTMD.io.model.global_current() + pyTMD.io.model.antarctic_current())
# display widgets for setting directory and model
TMDwidgets = pyTMD.tools.widgets()
TMDwidgets.model.options = model_list
TMDwidgets.model.value = 'CATS2008'
TMDwidgets.VBox([
    TMDwidgets.directory,
    TMDwidgets.model,
    TMDwidgets.atlas,
    TMDwidgets.compress,
    TMDwidgets.datepick
])

# get model parameters
model = pyTMD.io.model(TMDwidgets.directory.value,
    format=TMDwidgets.atlas.value,
    compressed=TMDwidgets.compress.value
   ).current(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())

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

# save tide currents
tide = {}
# iterate over u and v currents
for TYPE in model.type:
    # read tidal constants and interpolate to grid points
    if model.format in ('OTIS','ATLAS','TMD3'):
        amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,
            model.model_file['u'], model.projection, type=TYPE, crop=True,
            method='spline', grid=model.format)
        DELTAT = np.zeros_like(ts.tide)
    elif (model.format == 'netcdf'):
        amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file,
            model.model_file[TYPE], type=TYPE, crop=True, method='spline',
            scale=model.scale, compressed=model.compressed)
        DELTAT = np.zeros_like(ts.tide)
    elif (model.format == 'GOT'):
        amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file[TYPE],
            crop=True, method='spline', scale=model.scale,
            compressed=model.compressed)
        # delta time (TT - UT1)
        DELTAT = ts.tt_ut1
    elif (model.format == 'FES'):
        amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file[TYPE],
            type=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 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 = pyTMD.predict.map(ts.tide[hour], hc, c, deltat=DELTAT[hour],
            corrections=model.format)
        MINOR = pyTMD.predict.infer_minor(ts.tide[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_title(f'{model.name} Tidal Velocity', fontsize=13,
    rotation=0, y=-1.65, va='top')
cbar.ax.set_xlabel('cm/s', fontsize=13, rotation=0, va='center')
cbar.ax.xaxis.set_label_coords(1.075, 0.5)
# 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, y=0.97, 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 model.type:
        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())