Compare with Tide Gauges

This (notebook) compares model predictions with tide gauge data provided by the University of Hawaii Sea Level Center, and computes a local solution from the tide gauge data

OTIS format tidal solutions provided by Oregon State University and ESR

Global Tide Model (GOT) solutions provided by Richard Ray at GSFC

https://earth.gsfc.nasa.gov/geo/data/ocean-tide-models

Finite Element Solution (FES) provided by AVISO

https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html

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
solve.py: estimates the harmonic constants for ocean tides
time.py: utilities for calculating time operations

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

Load modules

from __future__ import print_function

import netCDF4
import ipywidgets
import numpy as np
import scipy.signal
import IPython.display
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

# import tide programs
import pyTMD.tools
import pyTMD.compute
import pyTMD.solve
import pyTMD.utilities
import timescale
# 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,
])

Find Tide Gauge Data

# remote directory with tide gauge data
HOST = 'https://uhslc.soest.hawaii.edu/data/netcdf/rqds/global/hourly/'
f = pyTMD.utilities.uhslc_list(HOST, pattern=r'(.*?).nc', sort=True)
# create dropdown with all tide gauge data
TMDwidgets.gauges = ipywidgets.Dropdown(
    options=f,
    value=f[0],
    description='Tide Gauges:',
    disabled=False,
    style=TMDwidgets.style,
)
display(TMDwidgets.gauges)

Download and Read Tide Gauge Data

# open tide gauge data
fid = pyTMD.utilities.from_http([HOST,TMDwidgets.gauges.value])
with netCDF4.Dataset(fid.filename, memory=fid.read()) as fileID:
    # read time and station name
    delta_time = fileID.variables['time'][:].squeeze()
    date_string = fileID.variables['time'].units
    station_name, = netCDF4.chartostring(fileID.variables['station_name'][:])
    print(f'{fid.filename}: {station_name}')
    # get station latitude and longitude
    lat, = fileID.variables['lat'][:]
    lon, = fileID.variables['lon'][:]
    # get sea level heights
    sea_level = fileID.variables['sea_level'][:].squeeze()
# reduce to valid points
valid = np.logical_not(sea_level.mask | np.isnan(sea_level))
# convert time
epoch, to_sec = timescale.time.parse_date_string(date_string)
ts = timescale.from_deltatime(delta_time[valid]*to_sec, epoch=epoch)
print(f'{ts.min().to_datetime()[0]}')
print(f'{ts.max().to_datetime()[0]}')
deltat = ts.tt_ut1
# remove the mean and convert to meters
h = (sea_level.compressed() - sea_level[valid].mean())/1000.0

%matplotlib widget
# check coordinates on tide grid
f1 = plt.figure(num=1, clear=True, figsize=(8.25,5.25))
ax1 = f1.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
point, = ax1.plot(lon, lat, 'r*', transform=ccrs.PlateCarree())
# add title
ax1.set_title(f'{station_name} ({lat:0.2f}\u00B0N, {lon:0.2f}\u00B0E)')
# add coastlines
ax1.coastlines()
ax1.add_feature(cfeature.LAND, facecolor='0.85')
# axis = equal
ax1.set_aspect('equal', adjustable='box')
# no ticks on the x and y axes
ax1.get_xaxis().set_ticks([])
ax1.get_yaxis().set_ticks([])
ax1.set_extent([-180, 180, -90, 90])
# stronger linewidth on frame
ax1.spines['geo'].set_linewidth(2.0)
ax1.spines['geo'].set_capstyle('projecting')
# adjust subplot within figure
f1.tight_layout()
plt.show()

Predict tides at measurement times

# get model parameters
model = pyTMD.io.model(verify=False).elevation(TMDwidgets.model.value)
model.parse_constituents()
# calculate tide elevations
tide = pyTMD.compute.tide_elevations(lon, lat, ts.to_datetime(),
    DIRECTORY=TMDwidgets.directory.value, TYPE='time series', 
    MODEL=TMDwidgets.model.value, GZIP=TMDwidgets.compress.value,
    EPSG=4326, TIME='datetime', EXTRAPOLATE=True, CUTOFF=20).squeeze()

Adjust global solution for regional data

# use constituents from original model
c = model.constituents
amp, ph = pyTMD.solve.constants(ts.tide, h-tide.data, c,
    deltat=deltat, corrections=model.corrections)
# calculate complex phase in radians for Euler's
cph = np.atleast_2d(-1j*ph*np.pi/180.0)
# calculate constituent oscillation
hc = np.ma.array(amp*np.exp(cph))
hc.mask = hc.data == hc.fill_value
sol = pyTMD.predict.time_series(ts.tide, hc, c,
    deltat=deltat, corrections=model.corrections)
sol += tide

Calculate Periodograms to Compare Frequencies

# create array of angular frequencies
N = 1000
# frequency range (use Sa to M4)
sa, = pyTMD.arguments.frequency('sa')
m4, = pyTMD.arguments.frequency('m4')
omega = np.linspace(sa, m4, N)
# calculate Lomb-Scargle periodograms
lssa_h = scipy.signal.lombscargle(ts.J2000, h, omega, normalize=True)
lssa_tide = scipy.signal.lombscargle(ts.J2000, tide, omega, normalize=True)
lssa_sol = scipy.signal.lombscargle(ts.J2000, sol, omega, normalize=True)

Compare predictions and measured values

%matplotlib widget
# create figure
f2 = plt.figure(num=2, clear=True, figsize=(9, 6))
ax2 = f2.subplots(ncols=2)
# plot time series results
ax2[0].plot(ts.year, h, label='Tide Gauge', lw=2, color='darkorchid')
ax2[0].plot(ts.year, tide, label='Tide Prediction', color='mediumseagreen')
ax2[0].plot(ts.year, sol, label='Tide Solution', color='darkorange')
# plot Lomb-Scargle power
# convert to cycles per solar day
f = omega*86400.0/(2.0*np.pi)
ax2[1].semilogy(f, lssa_h, label='Tide Gauge', lw=2, color='darkorchid')
ax2[1].semilogy(f, lssa_tide, label='Tide Prediction', color='mediumseagreen')
ax2[1].semilogy(f, lssa_sol, label='Tide Solution', color='darkorange')
# add figure and axes titles
f2.suptitle(f'{station_name} ({lat:0.2f}\u00B0N, {lon:0.2f}\u00B0E)')
ax2[0].set_title('Time Series')
ax2[1].set_title('Lomb-Scargle Periodogram')
# add labels
ax2[0].set_xlabel('Time [yr]')
ax2[0].set_ylabel('Water Height [m]')
ax2[1].set_xlabel('Frequency [cpd]')
ax2[1].set_ylabel('Power')
# add legend
lgd = ax2[1].legend(frameon=False)
lgd.get_frame().set_alpha(1.0)
for line in lgd.get_lines():
    line.set_linewidth(6)
f2.tight_layout()
plt.show()

%matplotlib widget
# create figure
f3 = plt.figure(num=3, clear=True, figsize=(9, 5.5))
ax3 = f3.subplots(ncols=2, sharex=True, sharey=True)
xmin = np.minimum(tide.min(), sol.min())
xmax = np.maximum(tide.max(), sol.max())
ymin = np.minimum(h.min(), tide.min())
ymax = np.maximum(h.max(), tide.max())
ax3[0].hist2d(h, tide, bins=60, range=((xmin, xmax), (ymin, ymax)), cmap='Reds')
ax3[1].hist2d(h, sol, bins=60, range=((xmin, xmax), (ymin, ymax)), cmap='Reds')
ax3[0].axline((0, 0), slope=1, color='0.4', lw=1.0)
ax3[1].axline((0, 0), slope=1, color='0.4', lw=1.0)
ax3[0].set_xlabel('Tide Gauge [m]')
ax3[1].set_xlabel('Tide Gauge [m]')
ax3[0].set_ylabel('Tide Prediction [m]')
ax3[0].set_title(model.name)
ax3[1].set_title('Solution')
ax3[0].set_aspect('equal', adjustable='box')
ax3[1].set_aspect('equal', adjustable='box')
f3.suptitle(f'{station_name} ({lat:0.2f}\u00B0N, {lon:0.2f}\u00B0E)')
# adjust subplot within figure
f3.subplots_adjust(left=0.075,right=0.975,bottom=0.1,top=0.90, wspace=0.075)
f3.show()