Plot Solid Earth Tides
This (notebook) demonstrates creating a global animation of solid Earth tide displacements
Python Dependencies
Program Dependencies
astro.py: computes the basic astronomical mean longitudespredict.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 solid Earth tide maps.
Load modules
from __future__ import print_function
import os
import numpy as np
import ipywidgets
import matplotlib
matplotlib.rcParams['axes.linewidth'] = 2.0
matplotlib.rcParams["animation.html"] = "jshtml"
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.animation as animation
import cartopy.crs as ccrs
from IPython.display import HTML
# import tide programs
import pyTMD.predict
import pyTMD.spatial
import pyTMD.tools
import pyTMD.utilities
import timescale.time
# autoreload
%load_ext autoreload
%autoreload 2
Set parameters for program
# display widgets
TMDwidgets = pyTMD.tools.widgets()
# create dropdown with ephemeride type
TMDwidgets.ephemerides = ipywidgets.Dropdown(
options=['Approximate', 'JPL'],
value='JPL',
description='Ephemerides:',
disabled=False,
style=TMDwidgets.style,
)
# create dropdown with permanent tide type
TMDwidgets.tide_system = ipywidgets.Dropdown(
options=['mean_tide', 'tide_free'],
value='tide_free',
description='Tide System:',
disabled=False,
style=TMDwidgets.style,
)
TMDwidgets.VBox([
TMDwidgets.ephemerides,
TMDwidgets.tide_system,
TMDwidgets.datepick
])
Setup coordinates for calculating tides
# create a global image
xlimits = [-180,180]
ylimits = [-90, 90]
spacing = [1.0, 1.0]
# x and y coordinates
x = np.arange(xlimits[0],xlimits[1]+spacing[0],spacing[0])
y = np.arange(ylimits[0],ylimits[1]+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[1]-ylimits[0])/spacing[1])+1
# flatten latitude and longitude to arrays
lon,lat = xgrid.flatten(), ygrid.flatten()
Calculate and plot solid Earth tides
# convert from calendar date to days relative to Jan 1, 1992 (48622 MJD)
YMD = TMDwidgets.datepick.value
# calculate a forecast every hour
hours = np.arange(24)
# convert time from MJD to days relative to Jan 1, 1992 (48622 MJD)
ts = timescale.time.Timescale().from_calendar(YMD.year, YMD.month,
YMD.day, hour=hours)
nt = len(ts)
# earth and physical parameters for ellipsoid
units = pyTMD.spatial.datum(ellipsoid='IERS', units='MKS')
# compute ephemerides for lunisolar coordinates
EPHEMERIDES = TMDwidgets.ephemerides.value
SX, SY, SZ = pyTMD.astro.solar_ecef(ts.MJD, ephemerides=EPHEMERIDES)
LX, LY, LZ = pyTMD.astro.lunar_ecef(ts.MJD, ephemerides=EPHEMERIDES)
# convert input coordinates to cartesian
X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat,
a_axis=units.a_axis, flat=units.flat)
# convert coordinates to column arrays
XYZ = np.c_[X, Y, Z]
# geocentric latitude (radians)
latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))
npts = len(latitude_geocentric)
# geocentric colatitude (radians)
theta = (np.pi/2.0 - latitude_geocentric)
# calculate longitude (radians)
phi = np.arctan2(Y, X)
# rotation matrix for converting from cartesian coordinates
R = np.zeros((npts, 3, 3))
R[:,0,0] = np.cos(phi)*np.cos(theta)
R[:,1,0] = -np.sin(phi)
R[:,2,0] = np.cos(phi)*np.sin(theta)
R[:,0,1] = np.sin(phi)*np.cos(theta)
R[:,1,1] = np.cos(phi)
R[:,2,1] = np.sin(phi)*np.sin(theta)
R[:,0,2] = -np.sin(theta)
R[:,2,2] = np.cos(theta)
# create a time-varying map
tide_se = np.zeros((ny,nx,nt))
for i in range(nt):
# convert tide times to dynamical time
# reshape time to match spatial
t = ts.tide[i] + ts.tt_ut1[i] + np.ones((ny*nx))
# convert coordinates to column arrays
SXYZ = np.repeat(np.c_[SX[i], SY[i], SZ[i]], ny*nx, axis=0)
LXYZ = np.repeat(np.c_[LX[i], LY[i], LZ[i]], ny*nx, axis=0)
# predict solid earth tides (cartesian)
dxi = pyTMD.predict.solid_earth_tide(t,
XYZ, SXYZ, LXYZ, a_axis=units.a_axis,
tide_system=TMDwidgets.tide_system.value)
# calculate components of solid earth tides
SE = np.einsum('ti...,tji...->tj...', dxi, R)
# reshape to output dimensions and convert to cm
tide_se[:,:,i] = 100.0*np.reshape(SE[:,2], (ny,nx))
# cartopy transform for Equirectangular Projection
projection = ccrs.PlateCarree()
# create figure axis
fig, ax = plt.subplots(num=1, figsize=(5.5,3.5),
subplot_kw=dict(projection=projection))
# plot solid earth tide height
norm = colors.CenteredNorm(halfrange=np.max(np.abs(tide_se)))
extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])
im = ax.imshow(np.zeros((ny,nx)), interpolation='nearest',
norm=norm, cmap='Spectral', transform=projection,
extent=extent, origin='upper', animated=True)
# 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, extend='both',
extendfrac=0.0375, orientation='horizontal', pad=0.025,
shrink=0.90, aspect=22, drawedges=False)
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
cbar.ax.set_title('Radial Tide Displacement', fontsize=13,
rotation=0, y=-2.0, va='top')
cbar.ax.set_xlabel('cm', fontsize=13, rotation=0, va='center')
cbar.ax.xaxis.set_label_coords(1.085, 0.5)
# ticks lines all the way across
cbar.ax.tick_params(which='both', width=1, length=15, labelsize=13,
direction='in')
# add title (date and time)
ttl = ax.set_title(None, fontsize=13)
# axis = equal
ax.set_aspect('equal', adjustable='box')
# set x and y limits
ax.set_xlim(xlimits)
ax.set_ylim(ylimits)
# no ticks on the x and y axes
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
# 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.04,right=0.96,bottom=0.05,top=0.96)
# animate each map
def animate_maps(hour):
# set map data
im.set_data(tide_se[:,:,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())