Source code for pyTMD.predict

#!/usr/bin/env python
u"""
predict.py
Written by Tyler Sutterley (05/2025)
Prediction routines for ocean, load, equilibrium and solid earth tides

REFERENCES:
    G. D. Egbert and S. Erofeeva, "Efficient Inverse Modeling of Barotropic
        Ocean Tides", Journal of Atmospheric and Oceanic Technology, (2002).

PYTHON DEPENDENCIES:
    numpy: Scientific Computing Tools For Python
        https://numpy.org
        https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
    timescale: Python tools for time and astronomical calculations
        https://pypi.org/project/timescale/

PROGRAM DEPENDENCIES:
    arguments.py: loads nodal corrections for tidal constituents
    astro.py: computes the basic astronomical mean longitudes
    crs.py: Coordinate Reference System (CRS) routines
    spatial.py: utilities for working with geospatial data

UPDATE HISTORY:
    Updated 05/2025: pass keyword arguments to nodal corrections functions
    Updated 03/2025: changed argument for method calculating mean longitudes
    Updated 02/2025: verify dimensions of harmonic constants
    Updated 11/2024: use Love numbers for long-period tides when inferring
        move body tide Love/Shida numbers to arguments module
    Updated 10/2024: use PREM as the default Earth model for Love numbers
        more descriptive error message if cannot infer minor constituents
        updated calculation of long-period equilibrium tides
        added option to use Munk-Cartwright admittance interpolation for minor
    Updated 09/2024: verify order of minor constituents to infer
        fix to use case insensitive assertions of string argument values
        split infer minor function into short and long period calculations
        add two new functions to infer semi-diurnal and diurnal tides separately
    Updated 08/2024: minor nodal angle corrections in radians to match arguments
        include inference of eps2 and eta2 when predicting from GOT models
        add keyword argument to allow inferring specific minor constituents
    	use nodal arguments for all non-OTIS model type cases
        add load pole tide function that exports in cartesian coordinates
        add ocean pole tide function that exports in cartesian coordinates
    Updated 07/2024: use normalize_angle from pyTMD astro module
        make number of days to convert tide time to MJD a variable
    Updated 02/2024: changed class name for ellipsoid parameters to datum
    Updated 01/2024: moved minor arguments calculation into new function
        moved constituent parameters function from predict to arguments
    Updated 12/2023: phase_angles function renamed to doodson_arguments
    Updated 09/2023: moved constituent parameters function within this module
    Updated 08/2023: changed ESR netCDF4 format to TMD3 format
    Updated 04/2023: using renamed astro mean_longitudes function
        using renamed arguments function for nodal corrections
        adding prediction routine for solid earth tides
        output solid earth tide corrections as combined XYZ components
    Updated 03/2023: add basic variable typing to function inputs
    Updated 12/2022: merged prediction functions into a single module
    Updated 05/2022: added ESR netCDF4 formats to list of model types
    Updated 04/2022: updated docstrings to numpy documentation format
    Updated 02/2021: replaced numpy bool to prevent deprecation warning
    Updated 09/2020: append output mask over each constituent
    Updated 08/2020: change time variable names to not overwrite functions
    Updated 07/2020: added function docstrings
    Updated 11/2019: output as numpy masked arrays instead of nan-filled arrays
    Updated 09/2019: added netcdf option to CORRECTIONS option
    Updated 08/2018: added correction option ATLAS for localized OTIS solutions
    Updated 07/2018: added option to use GSFC GOT nodal corrections
    Updated 09/2017: Rewritten in Python
"""
from __future__ import annotations

import logging
import numpy as np
import pyTMD.arguments
import pyTMD.astro
import pyTMD.math
import pyTMD.spatial
import timescale.time

__all__ = [
    "map",
    "drift",
    "time_series",
    "infer_minor",
    "_infer_short_period",
    "_infer_semi_diurnal",
    "_infer_diurnal",
    "_infer_long_period",
    "equilibrium_tide",
    "load_pole_tide",
    "ocean_pole_tide",
    "solid_earth_tide",
    "_out_of_phase_diurnal",
    "_out_of_phase_semidiurnal",
    "_latitude_dependence",
    "_frequency_dependence_diurnal",
    "_frequency_dependence_long_period",
    "_free_to_mean"
]

# number of days between the Julian day epoch and MJD
_jd_mjd = 2400000.5
# number of days between MJD and the tide epoch (1992-01-01T00:00:00)
_mjd_tide = 48622.0
# number of days between the Julian day epoch and the tide epoch
_jd_tide = _jd_mjd + _mjd_tide

# PURPOSE: Predict tides at single times
[docs] def map(t: float | np.ndarray, hc: np.ndarray, constituents: list | np.ndarray, deltat: float | np.ndarray = 0.0, corrections: str = 'OTIS', **kwargs ): """ Predict tides at a single time using harmonic constants :cite:p:`Egbert:2002ge` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 hc: np.ndarray harmonic constant vector constituents: list or np.ndarray tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models **kwargs: dict keyword arguments for nodal corrections functions Returns ------- ht: np.ndarray tide values reconstructed using the nodal corrections """ # number of points and number of constituents npts, nc = np.shape(hc) # verify dimensions of harmonic constants hc = np.ma.atleast_2d(hc) # load the nodal corrections # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide, constituents, deltat=deltat, corrections=corrections, **kwargs ) # allocate for output tidal elevation ht = np.ma.zeros((npts)) ht.mask = np.zeros((npts), dtype=bool) # for each constituent for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent amp, ph, omega, alpha, species = \ pyTMD.arguments._constituent_parameters(c) # add component for constituent to output tidal elevation th = omega*t*86400.0 + ph + pu[0,k] else: th = G[0,k]*np.pi/180.0 + pu[0,k] # sum over all tides ht.data[:] += pf[0,k]*hc.real[:,k]*np.cos(th) - \ pf[0,k]*hc.imag[:,k]*np.sin(th) ht.mask[:] |= (hc.real.mask[:,k] | hc.imag.mask[:,k]) # return the tidal elevation after removing singleton dimensions return np.squeeze(ht)
# PURPOSE: Predict tides at drift bouys or altimetry points
[docs] def drift(t: float | np.ndarray, hc: np.ndarray, constituents: list | np.ndarray, deltat: float | np.ndarray = 0.0, corrections: str = 'OTIS', **kwargs ): """ Predict tides at multiple times and locations using harmonic constants :cite:p:`Egbert:2002ge` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 hc: np.ndarray harmonic constant vector constituents: list or np.ndarray tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models **kwargs: dict keyword arguments for nodal corrections functions Returns ------- ht: np.ndarray tidal time series reconstructed using the nodal corrections """ # number of points nt = len(t) # verify dimensions of harmonic constants hc = np.ma.atleast_2d(hc) # load the nodal corrections # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide, constituents, deltat=deltat, corrections=corrections, **kwargs ) # allocate for output time series ht = np.ma.zeros((nt)) ht.mask = np.zeros((nt), dtype=bool) # for each constituent for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent amp, ph, omega, alpha, species = \ pyTMD.arguments._constituent_parameters(c) # add component for constituent to output tidal elevation th = omega*t*86400.0 + ph + pu[:,k] else: th = G[:,k]*np.pi/180.0 + pu[:,k] # sum over all tides ht.data[:] += pf[:,k]*hc.real[:,k]*np.cos(th) - \ pf[:,k]*hc.imag[:,k]*np.sin(th) ht.mask[:] |= (hc.real.mask[:,k] | hc.imag.mask[:,k]) # return tides return ht
# PURPOSE: Predict a tidal time series at a location
[docs] def time_series(t: float | np.ndarray, hc: np.ndarray, constituents: list | np.ndarray, deltat: float | np.ndarray = 0.0, corrections: str = 'OTIS', **kwargs ): """ Predict tidal time series at a single location using harmonic constants :cite:p:`Egbert:2002ge` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 hc: np.ndarray harmonic constant vector constituents: list or np.ndarray tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models **kwargs: dict keyword arguments for nodal corrections functions Returns ------- ht: np.ndarray tidal time series reconstructed using the nodal corrections """ # number of time points nt = len(t) # verify dimensions of harmonic constants hc = np.ma.atleast_2d(hc) # load the nodal corrections # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide, constituents, deltat=deltat, corrections=corrections, **kwargs ) # allocate for output time series ht = np.ma.zeros((nt)) ht.mask = np.zeros((nt), dtype=bool) # for each constituent for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent amp, ph, omega, alpha, species = \ pyTMD.arguments._constituent_parameters(c) # add component for constituent to output tidal time series th = omega*t*86400.0 + ph + pu[:,k] else: th = G[:,k]*np.pi/180.0 + pu[:,k] # sum over all tides at location ht.data[:] += pf[:,k]*hc.real[0,k]*np.cos(th) - \ pf[:,k]*hc.imag[0,k]*np.sin(th) ht.mask[:] |= np.any(hc.real.mask[0,k] | hc.imag.mask[0,k]) # return the tidal time series return ht
# PURPOSE: infer the minor corrections from the major constituents
[docs] def infer_minor( t: float | np.ndarray, zmajor: np.ndarray, constituents: list | np.ndarray, **kwargs ): """ Infer the tidal values for minor constituents using their relation with major constituents :cite:p:`Doodson:1941td,Schureman:1958ty,Foreman:1989dt,Egbert:2002ge` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 zmajor: np.ndarray Complex HC for given constituents/points constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models raise_exception: bool, default False Raise a ``ValueError`` if major constituents are not found minor: list or None, default None tidal constituent IDs of minor constituents for inference raise_exception: bool, default False Raise a ``ValueError`` if major constituents are not found Returns ------- dh: np.ndarray tidal time series for minor constituents """ # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'OTIS') kwargs.setdefault('raise_exception', False) # list of minor constituents kwargs.setdefault('minor', None) # infer the minor tidal constituents dh = 0.0 # infer short-period tides for minor constituents if kwargs['corrections'] in ('GOT',): dh += _infer_semi_diurnal(t, zmajor, constituents, **kwargs) dh += _infer_diurnal(t, zmajor, constituents, **kwargs) else: dh += _infer_short_period(t, zmajor, constituents, **kwargs) # infer long-period tides for minor constituents dh += _infer_long_period(t, zmajor, constituents, **kwargs) # return the inferred values return dh
# PURPOSE: infer short-period tides for minor constituents
[docs] def _infer_short_period( t: float | np.ndarray, zmajor: np.ndarray, constituents: list | np.ndarray, **kwargs ): """ Infer the tidal values for short-period minor constituents using their relation with major constituents :cite:p:`Egbert:2002ge,Ray:1999vm` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 zmajor: np.ndarray Complex HC for given constituents/points constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models minor: list or None, default None tidal constituent IDs of minor constituents for inference raise_exception: bool, default False Raise a ``ValueError`` if major constituents are not found Returns ------- dh: np.ndarray tidal time series for minor constituents """ # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'OTIS') kwargs.setdefault('raise_exception', False) # list of minor constituents kwargs.setdefault('minor', None) # number of constituents zmajor = np.ma.atleast_2d(zmajor) npts, nc = np.shape(zmajor) nt = len(np.atleast_1d(t)) # number of data points to calculate if running time series/drift/map n = nt if ((npts == 1) & (nt > 1)) else npts # allocate for output elevation correction dh = np.ma.zeros((n)) # major constituents used for inferring minor tides cindex = ['q1', 'o1', 'p1', 'k1', 'n2', 'm2', 's2', 'k2', '2n2'] # re-order major tides to correspond to order of cindex z = np.ma.zeros((n,len(cindex)), dtype=np.complex64) nz = 0 for i,c in enumerate(cindex): j = [j for j,val in enumerate(constituents) if (val.lower() == c)] if j: j1, = j z[:,i] = zmajor[:,j1] nz += 1 # raise exception or log error msg = 'Not enough constituents for inference of short-period tides' if (nz < 6) and kwargs['raise_exception']: raise Exception(msg) elif (nz < 6): logging.debug(msg) return 0.0 # complete list of minor constituents minor_constituents = ['2q1', 'sigma1', 'rho1', 'm1b', 'm1', 'chi1', 'pi1', 'phi1', 'theta1', 'j1', 'oo1', '2n2', 'mu2', 'nu2', 'lambda2', 'l2', 'l2b', 't2', 'eps2', 'eta2'] # possibly reduced list of minor constituents minor = kwargs['minor'] or minor_constituents # only add minor constituents that are not on the list of major values minor_indices = [i for i,m in enumerate(minor_constituents) if (m not in constituents) and (m in minor)] # if there are no constituents to infer msg = 'No short-period tidal constituents to infer' if not np.any(minor_indices): logging.debug(msg) return 0.0 # relationship between major and minor constituent amplitude and phase zmin = np.zeros((n, 20), dtype=np.complex64) zmin[:,0] = 0.263*z[:,0] - 0.0252*z[:,1]# 2Q1 zmin[:,1] = 0.297*z[:,0] - 0.0264*z[:,1]# sigma1 zmin[:,2] = 0.164*z[:,0] + 0.0048*z[:,1]# rho1 zmin[:,3] = 0.0140*z[:,1] + 0.0101*z[:,3]# M12 zmin[:,4] = 0.0389*z[:,1] + 0.0282*z[:,3]# M11 zmin[:,5] = 0.0064*z[:,1] + 0.0060*z[:,3]# chi1 zmin[:,6] = 0.0030*z[:,1] + 0.0171*z[:,3]# pi1 zmin[:,7] = -0.0015*z[:,1] + 0.0152*z[:,3]# phi1 zmin[:,8] = -0.0065*z[:,1] + 0.0155*z[:,3]# theta1 zmin[:,9] = -0.0389*z[:,1] + 0.0836*z[:,3]# J1 zmin[:,10] = -0.0431*z[:,1] + 0.0613*z[:,3]# OO1 zmin[:,11] = 0.264*z[:,4] - 0.0253*z[:,5]# 2N2 zmin[:,12] = 0.298*z[:,4] - 0.0264*z[:,5]# mu2 zmin[:,13] = 0.165*z[:,4] + 0.00487*z[:,5]# nu2 zmin[:,14] = 0.0040*z[:,5] + 0.0074*z[:,6]# lambda2 zmin[:,15] = 0.0131*z[:,5] + 0.0326*z[:,6]# L2 zmin[:,16] = 0.0033*z[:,5] + 0.0082*z[:,6]# L2 zmin[:,17] = 0.0585*z[:,6]# t2 # additional coefficients for FES and GOT models if kwargs['corrections'] in ('FES',): # spline coefficients for admittances mu2 = [0.069439968323, 0.351535557706, -0.046278307672] nu2 = [-0.006104695053, 0.156878802427, 0.006755704028] l2 = [0.077137765667, -0.051653455134, 0.027869916824] t2 = [0.180480173707, -0.020101177502, 0.008331518844] lda2 = [0.016503557465, -0.013307812292, 0.007753383202] zmin[:,12] = mu2[0]*z[:,7] + mu2[1]*z[:,4] + mu2[2]*z[:,5]# mu2 zmin[:,13] = nu2[0]*z[:,7] + nu2[1]*z[:,4] + nu2[2]*z[:,5]# nu2 zmin[:,14] = lda2[0]*z[:,7] + lda2[1]*z[:,4] + lda2[2]*z[:,5]# lambda2 zmin[:,16] = l2[0]*z[:,7] + l2[1]*z[:,4] + l2[2]*z[:,5]# L2 zmin[:,17] = t2[0]*z[:,7] + t2[1]*z[:,4] + t2[2]*z[:,5]# t2 zmin[:,18] = 0.53285*z[:,8] - 0.03304*z[:,4]# eps2 zmin[:,19] = -0.0034925*z[:,5] + 0.0831707*z[:,7]# eta2 # load the nodal corrections for minor constituents # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.minor_arguments(t + _mjd_tide, deltat=kwargs['deltat'], corrections=kwargs['corrections'] ) # sum over the minor tidal constituents of interest for k in minor_indices: th = G[:,k]*np.pi/180.0 + pu[:,k] dh += zmin.real[:,k]*pf[:,k]*np.cos(th) - \ zmin.imag[:,k]*pf[:,k]*np.sin(th) # return the inferred values return dh
[docs] def _infer_semi_diurnal( t: float | np.ndarray, zmajor: np.ndarray, constituents: list | np.ndarray, **kwargs ): """ Infer the tidal values for semi-diurnal minor constituents using their relation with major constituents :cite:p:`Munk:1966go,Ray:1999vm,Cartwright:1971iz` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 zmajor: np.ndarray Complex HC for given constituents/points constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) minor: list or None, default None tidal constituent IDs of minor constituents for inference method: str, default 'linear' method for interpolating between major constituents * 'linear': linear interpolation * 'admittance': Munk-Cartwright admittance interpolation raise_exception: bool, default False Raise a ``ValueError`` if major constituents are not found Returns ------- dh: np.ndarray tidal time series for minor constituents """ # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'GOT') kwargs.setdefault('method', 'linear') kwargs.setdefault('raise_exception', False) # list of minor constituents kwargs.setdefault('minor', None) # validate interpolation method assert kwargs['method'].lower() in ('linear', 'admittance') # number of constituents zmajor = np.ma.atleast_2d(zmajor) npts, nc = np.shape(zmajor) nt = len(np.atleast_1d(t)) # number of data points to calculate if running time series/drift/map n = nt if ((npts == 1) & (nt > 1)) else npts # allocate for output elevation correction dh = np.ma.zeros((n)) # major constituents used for inferring semi-diurnal minor tides # pivot waves listed in Table 6.7 of the 2010 IERS Conventions cindex = ['n2', 'm2', 's2'] # angular frequencies for major constituents omajor = pyTMD.arguments.frequency(cindex, **kwargs) # Cartwright and Edden potential amplitudes for major constituents amajor = np.zeros((3)) amajor[0] = 0.121006# n2 amajor[1] = 0.631931# m2 amajor[2] = 0.294019# s2 # re-order major tides to correspond to order of cindex z = np.ma.zeros((n,len(cindex)), dtype=np.complex64) nz = 0 for i,c in enumerate(cindex): j = [j for j,val in enumerate(constituents) if (val.lower() == c)] if j: j1, = j # "normalize" tide values z[:,i] = zmajor[:,j1]/amajor[i] nz += 1 # raise exception or log error msg = 'Not enough constituents for inference of semi-diurnal tides' if (nz < 3) and kwargs['raise_exception']: raise Exception(msg) elif (nz < 3): logging.debug(msg) return 0.0 # complete list of minor constituents minor_constituents = ['eps2', '2n2', 'mu2', 'nu2', 'gamma2', 'alpha2', 'beta2', 'delta2', 'lambda2', 'l2', 't2', 'r2', 'k2', 'eta2'] # possibly reduced list of minor constituents minor = kwargs['minor'] or minor_constituents # only add minor constituents that are not on the list of major values minor_indices = [i for i,m in enumerate(minor_constituents) if (m not in constituents) and (m in minor)] # if there are no constituents to infer msg = 'No semi-diurnal tidal constituents to infer' if not np.any(minor_indices): logging.debug(msg) return 0.0 # angular frequencies for inferred constituents omega = pyTMD.arguments.frequency(minor_constituents, **kwargs) # Cartwright and Edden potential amplitudes for inferred constituents amin = np.zeros((14)) amin[0] = 0.004669# eps2 amin[1] = 0.016011# 2n2 amin[2] = 0.019316# mu2 amin[3] = 0.022983# nu2 amin[4] = 0.001902# gamma2 amin[5] = 0.002178# alpha2 amin[6] = 0.001921# beta2 amin[7] = 0.000714# delta2 amin[8] = 0.004662# lambda2 amin[9] = 0.017862# l2 amin[10] = 0.017180# t2 amin[11] = 0.002463# r2 amin[12] = 0.079924# k2 amin[13] = 0.004467# eta # load the nodal corrections for minor constituents # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide, minor_constituents, deltat=kwargs['deltat'], corrections=kwargs['corrections'] ) # sum over the minor tidal constituents of interest for k in minor_indices: # interpolate from major constituents if (kwargs['method'].lower() == 'linear'): # linearly interpolate between major constituents if (omajor[0] < omajor[1]) and (omega[k] < omajor[1]): slope = (z[:,1] - z[:,0])/(omajor[1] - omajor[0]) zmin = amin[k]*(z[:,0] + slope*(omega[k] - omajor[0])) else: slope = (z[:,2] - z[:,1])/(omajor[2] - omajor[1]) zmin = amin[k]*(z[:,1] + slope*(omega[k] - omajor[1])) elif (kwargs['method'].lower() == 'admittance'): # admittance interpolation using Munk-Cartwright approach Ainv = np.array([[3.3133, -4.2538, 1.9405], [-3.3133, 4.2538, -0.9405], [1.5018, -3.2579, 1.7561]]) coef = np.dot(Ainv, z.T) # convert frequency to radians per 48 hours # following Munk and Cartwright (1966) f = 2.0*omega[k]*86400.0 # calculate interpolated values for constituent interp = coef[0,:] + coef[1,:]*np.cos(f) + coef[2,:]*np.sin(f) # rescale tide values zmin = amin[k]*interp # sum over all tides th = G[:,k]*np.pi/180.0 + pu[:,k] dh += zmin.real*pf[:,k]*np.cos(th) - \ zmin.imag*pf[:,k]*np.sin(th) # return the inferred values return dh
[docs] def _infer_diurnal( t: float | np.ndarray, zmajor: np.ndarray, constituents: list | np.ndarray, **kwargs ): """ Infer the tidal values for diurnal minor constituents using their relation with major constituents taking into account resonance due to free core nutation :cite:p:`Munk:1966go,Ray:2017jx,Wahr:1981if,Cartwright:1973em` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 zmajor: np.ndarray Complex HC for given constituents/points constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) minor: list or None, default None tidal constituent IDs of minor constituents for inference method: str, default 'linear' method for interpolating between major constituents * 'linear': linear interpolation * 'admittance': Munk-Cartwright admittance interpolation raise_exception: bool, default False Raise a ``ValueError`` if major constituents are not found Returns ------- dh: np.ndarray tidal time series for minor constituents """ # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'GOT') kwargs.setdefault('method', 'linear') kwargs.setdefault('raise_exception', False) # list of minor constituents kwargs.setdefault('minor', None) # validate interpolation method assert kwargs['method'].lower() in ('linear', 'admittance') # number of constituents zmajor = np.ma.atleast_2d(zmajor) npts, nc = np.shape(zmajor) nt = len(np.atleast_1d(t)) # number of data points to calculate if running time series/drift/map n = nt if ((npts == 1) & (nt > 1)) else npts # allocate for output elevation correction dh = np.ma.zeros((n)) # major constituents used for inferring diurnal minor tides # pivot waves listed in Table 6.7 of the 2010 IERS Conventions cindex = ['q1', 'o1', 'k1'] # angular frequencies for major constituents omajor = pyTMD.arguments.frequency(cindex, **kwargs) # Cartwright and Edden potential amplitudes for major constituents amajor = np.zeros((3)) amajor[0] = 0.050184# q1 amajor[1] = 0.262163# o1 amajor[2] = 0.368731# k1 # re-order major tides to correspond to order of cindex z = np.ma.zeros((n,len(cindex)), dtype=np.complex64) nz = 0 for i,c in enumerate(cindex): j = [j for j,val in enumerate(constituents) if (val.lower() == c)] if j: j1, = j # Love numbers of degree 2 for constituent h2, k2, l2 = pyTMD.arguments._love_numbers(omajor[i]) # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # "normalize" tide values z[:,i] = zmajor[:,j1]/(amajor[i]*gamma_2) nz += 1 # raise exception or log error msg = 'Not enough constituents for inference of diurnal tides' if (nz < 3) and kwargs['raise_exception']: raise Exception(msg) elif (nz < 3): logging.debug(msg) return 0.0 # complete list of minor constituents minor_constituents = ['2q1', 'sigma1', 'rho1', 'tau1', 'beta1', 'm1a', 'm1b', 'chi1', 'pi1', 'p1', 'psi1', 'phi1', 'theta1', 'j1', 'so1', 'oo1', 'ups1'] # possibly reduced list of minor constituents minor = kwargs['minor'] or minor_constituents # only add minor constituents that are not on the list of major values minor_indices = [i for i,m in enumerate(minor_constituents) if (m not in constituents) and (m in minor)] # if there are no constituents to infer msg = 'No diurnal tidal constituents to infer' if not np.any(minor_indices): logging.debug(msg) return 0.0 # angular frequencies for inferred constituents omega = pyTMD.arguments.frequency(minor_constituents, **kwargs) # Cartwright and Edden potential amplitudes for inferred constituents amin = np.zeros((17)) amin[0] = 0.006638# 2q1 amin[1] = 0.008023# sigma1 amin[2] = 0.009540# rho1 amin[3] = 0.003430# tau1 amin[4] = 0.001941# beta1 amin[5] = 0.020604# m1a amin[6] = 0.007420# m1b amin[7] = 0.003925# chi1 amin[8] = 0.007125# pi1 amin[9] = 0.122008# p1 amin[10] = 0.002929# psi1 amin[11] = 0.005247# phi1 amin[12] = 0.003966# theta1 amin[13] = 0.020618# j1 amin[14] = 0.003417# so1 amin[15] = 0.011293# oo1 amin[16] = 0.002157# ups1 # load the nodal corrections for minor constituents # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide, minor_constituents, deltat=kwargs['deltat'], corrections=kwargs['corrections'] ) # sum over the minor tidal constituents of interest for k in minor_indices: # Love numbers of degree 2 for constituent h2, k2, l2 = pyTMD.arguments._love_numbers(omega[k]) # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # interpolate from major constituents if (kwargs['method'].lower() == 'linear'): # linearly interpolate between major constituents if (omajor[0] < omajor[1]) and (omega[k] < omajor[1]): slope = (z[:,1] - z[:,0])/(omajor[1] - omajor[0]) zmin = amin[k]*gamma_2*(z[:,0] + slope*(omega[k] - omajor[0])) else: slope = (z[:,2] - z[:,1])/(omajor[2] - omajor[1]) zmin = amin[k]*gamma_2*(z[:,1] + slope*(omega[k] - omajor[1])) elif (kwargs['method'].lower() == 'admittance'): # admittance interpolation using Munk-Cartwright approach Ainv = np.array([[3.1214, -3.8494, 1.728], [-3.1727, 3.9559, -0.7832], [1.438, -3.0297, 1.5917]]) coef = np.dot(Ainv, z.T) # convert frequency to radians per 48 hours # following Munk and Cartwright (1966) f = 2.0*omega[k]*86400.0 # calculate interpolated values for constituent interp = coef[0,:] + coef[1,:]*np.cos(f) + coef[2,:]*np.sin(f) # rescale tide values zmin = amin[k]*gamma_2*interp # sum over all tides th = G[:,k]*np.pi/180.0 + pu[:,k] dh += zmin.real*pf[:,k]*np.cos(th) - \ zmin.imag*pf[:,k]*np.sin(th) # return the inferred values return dh
# PURPOSE: infer long-period tides for minor constituents
[docs] def _infer_long_period( t: float | np.ndarray, zmajor: np.ndarray, constituents: list | np.ndarray, **kwargs ): """ Infer the tidal values for long-period minor constituents using their relation with major constituents :cite:p:`Ray:1999vm,Ray:2014fu,Cartwright:1973em` Parameters ---------- t: float or np.ndarray days relative to 1992-01-01T00:00:00 zmajor: np.ndarray Complex HC for given constituents/points constituents: list tidal constituent IDs deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) minor: list or None, default None tidal constituent IDs of minor constituents for inference raise_exception: bool, default False Raise a ``ValueError`` if major constituents are not found Returns ------- dh: np.ndarray tidal time series for minor constituents """ # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'OTIS') kwargs.setdefault('raise_exception', False) # list of minor constituents kwargs.setdefault('minor', None) # number of constituents zmajor = np.ma.atleast_2d(zmajor) npts, nc = np.shape(zmajor) nt = len(np.atleast_1d(t)) # number of data points to calculate if running time series/drift/map n = nt if ((npts == 1) & (nt > 1)) else npts # allocate for output elevation correction dh = np.ma.zeros((n)) # major constituents used for inferring long period minor tides # pivot waves listed in Table 6.7 of the 2010 IERS Conventions cindex = ['node', 'mm', 'mf'] # angular frequencies for major constituents omajor = pyTMD.arguments.frequency(cindex, **kwargs) # Cartwright and Edden potential amplitudes for major constituents amajor = np.zeros((3)) amajor[0] = 0.027929# node amajor[1] = 0.035184# mm amajor[2] = 0.066607# mf # re-order major tides to correspond to order of cindex z = np.ma.zeros((n,len(cindex)), dtype=np.complex64) nz = 0 for i,c in enumerate(cindex): j = [j for j,val in enumerate(constituents) if (val.lower() == c)] if j: j1, = j z[:,i] = zmajor[:,j1]/amajor[i] nz += 1 # raise exception or log error msg = 'Not enough constituents for inference of long-period tides' if (nz < 3) and kwargs['raise_exception']: raise Exception(msg) elif (nz < 3): logging.debug(msg) return 0.0 # complete list of minor constituents minor_constituents = ['sa', 'ssa', 'sta', 'msm', 'msf', 'mst', 'mt', 'msqm', 'mq'] # possibly reduced list of minor constituents minor = kwargs['minor'] or minor_constituents # only add minor constituents that are not on the list of major values minor_indices = [i for i,m in enumerate(minor_constituents) if (m not in constituents) and (m in minor)] # if there are no constituents to infer msg = 'No long-period tidal constituents to infer' if not np.any(minor_indices): logging.debug(msg) return 0.0 # angular frequencies for inferred constituents omega = pyTMD.arguments.frequency(minor_constituents, **kwargs) # Cartwright and Edden potential amplitudes for inferred constituents amin = np.zeros((9)) amin[0] = 0.004922# sa amin[1] = 0.030988# ssa amin[2] = 0.001809# sta amin[3] = 0.006728# msm amin[4] = 0.005837# msf amin[5] = 0.002422# mst amin[6] = 0.012753# mt amin[7] = 0.002037# msqm amin[8] = 0.001687# mq # load the nodal corrections for minor constituents # convert time to Modified Julian Days (MJD) pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide, minor_constituents, deltat=kwargs['deltat'], corrections=kwargs['corrections'] ) # sum over the minor tidal constituents of interest for k in minor_indices: # linearly interpolate between major constituents if (omajor[0] < omajor[1]) and (omega[k] < omajor[1]): slope = (z[:,1] - z[:,0])/(omajor[1] - omajor[0]) zmin = amin[k]*(z[:,0] + slope*(omega[k] - omajor[0])) else: slope = (z[:,2] - z[:,1])/(omajor[2] - omajor[1]) zmin = amin[k]*(z[:,1] + slope*(omega[k] - omajor[1])) # sum over all tides th = G[:,k]*np.pi/180.0 + pu[:,k] dh += zmin.real*pf[:,k]*np.cos(th) - \ zmin.imag*pf[:,k]*np.sin(th) # return the inferred values return dh
# PURPOSE: estimate long-period equilibrium tides
[docs] def equilibrium_tide( t: np.ndarray, lat: np.ndarray, **kwargs ): """ Compute the long-period equilibrium tides the summation of fifteen tidal spectral lines from Cartwright-Tayler-Edden tables :cite:p:`Cartwright:1971iz,Cartwright:1973em` Parameters ---------- t: np.ndarray time (days relative to January 1, 1992) lat: np.ndarray latitude (degrees north) deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models constituents: list long-period tidal constituent IDs Returns ------- lpet: np.ndarray long-period equilibrium tide in meters """ # set default keyword arguments cindex = ['node', 'sa', 'ssa', 'msm', '065.445', 'mm', '065.465', 'msf', '075.355', 'mf', 'mf+', '075.575', 'mst', 'mt', '085.465'] kwargs.setdefault('constituents', cindex) kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'OTIS') # number of input points nt = len(np.atleast_1d(t)) nlat = len(np.atleast_1d(lat)) # set function for astronomical longitudes # use ASTRO5 routines if not using an OTIS type model if kwargs['corrections'] in ('OTIS','ATLAS','TMD3','netcdf'): method = 'Cartwright' else: method = 'ASTRO5' # convert from Modified Julian Dates into Ephemeris Time MJD = t + _mjd_tide # compute principal mean longitudes s, h, p, N, pp = pyTMD.astro.mean_longitudes(MJD + kwargs['deltat'], method=method) # convert to negative mean longitude of the ascending node (N') n = pyTMD.math.normalize_angle(360.0 - N) # determine equilibrium arguments fargs = np.c_[s, h, p, n, pp] # Cartwright and Edden potential amplitudes (centimeters) # assemble long-period tide potential from 15 CTE terms greater than 1 mm amajor = np.zeros((15)) # group 0,0 # nodal term is included but not the constant term. amajor[0] = 2.7929# node amajor[1] = -0.4922# sa amajor[2] = -3.0988# ssa # group 0,1 amajor[3] = -0.6728# msm amajor[4] = 0.231 amajor[5] = -3.5184# mm amajor[6] = 0.228 # group 0,2 amajor[7] = -0.5837# msf amajor[8] = -0.288 amajor[9] = -6.6607# mf amajor[10] = -2.763# mf+ amajor[11] = -0.258 # group 0,3 amajor[12] = -0.2422# mst amajor[13] = -1.2753# mt amajor[14] = -0.528 # set constituents to be iterable and lower case if isinstance(kwargs['constituents'], str): constituents = [kwargs['constituents'].lower()] else: constituents = [c.lower() for c in kwargs['constituents']] # reduce potential amplitudes to constituents CTE = np.zeros((15)) for i,c in enumerate(cindex): if c in constituents: CTE[i] = amajor[i] # Doodson coefficients for 15 long-period terms coef = np.zeros((5, 15)) # group 0,0 coef[:,0] = [0.0, 0.0, 0.0, 1.0, 0.0]# node coef[:,1] = [0.0, 1.0, 0.0, 0.0, -1.0]# sa coef[:,2] = [0.0, 2.0, 0.0, 0.0, 0.0]# ssa # group 0,1 coef[:,3] = [1.0, -2.0, 1.0, 0.0, 0.0]# msm coef[:,4] = [1.0, 0.0, -1.0, -1.0, 0.0] coef[:,5] = [1.0, 0.0, -1.0, 0.0, 0.0]# mm coef[:,6] = [1.0, 0.0, -1.0, 1.0, 0.0] # group 0,2 coef[:,7] = [2.0, -2.0, 0.0, 0.0, 0.0]# msf coef[:,8] = [2.0, 0.0, -2.0, 0.0, 0.0] coef[:,9] = [2.0, 0.0, 0.0, 0.0, 0.0]# mf coef[:,10] = [2.0, 0.0, 0.0, 1.0, 0.0]# mf+ coef[:,11] = [2.0, 0.0, 0.0, 2.0, 0.0] # group 0,3 coef[:,12] = [3.0, -2.0, 1.0, 0.0, 0.0]# mst coef[:,13] = [3.0, 0.0, -1.0, 0.0, 0.0]# mt coef[:,14] = [3.0, 0.0, -1.0, 1.0, 0.0] # determine equilibrium arguments G = np.dot(fargs, coef) Z = np.inner(np.cos(G*np.pi/180.0), CTE) # Love numbers for long-period tides (Wahr, 1981) k2 = 0.299 h2 = 0.606 # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # 2nd degree Legendre polynomials P20 = 0.5*(3.0*np.sin(lat*np.pi/180.0)**2 - 1.0) # calculate long-period equilibrium tide and convert to meters # Multiply by gamma_2 * normalization * P20(lat) if (nlat != nt): lpet = gamma_2*np.sqrt((4.0+1.0)/(4.0*np.pi))*np.outer(P20,Z/100.0) else: lpet = gamma_2*np.sqrt((4.0+1.0)/(4.0*np.pi))*P20*(Z/100.0) # return the long-period equilibrium tides return lpet
# PURPOSE: estimate load pole tides in Cartesian coordinates
[docs] def load_pole_tide( t: np.ndarray, XYZ: np.ndarray, deltat: float = 0.0, gamma_0: float = 9.80665, omega: float = 7.2921151467e-5, h2: float = 0.6207, l2: float = 0.0836, convention: str = '2018' ): """ Estimate load pole tide displacements in Cartesian coordinates :cite:p:`Petit:2010tp` Parameters ---------- t: np.ndarray Time (days relative to January 1, 1992) XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) gamma_0: float, default 9.80665 Normal gravity (m/s^2) omega: float, default 7.2921151467e-5 Earth's rotation rate (radians/second) h2: float, default 0.6207 Degree-2 Love number of vertical displacement l2: float, default 0.0836 Degree-2 Love (Shida) number of horizontal displacement convention: str, default '2018' IERS Mean or Secular Pole Convention - ``'2003'`` - ``'2010'`` - ``'2015'`` - ``'2018'`` Returns ------- dxt: np.ndarray Load pole tide displacements in meters in Cartesian coordinates """ # degrees and arcseconds to radians dtr = np.pi/180.0 atr = np.pi/648000.0 # convert time to Terrestial Time (TT) tt = t + _jd_tide + deltat # convert time to Modified Julian Days (MJD) MJD = tt - _jd_mjd # convert Julian days to calendar dates Y,M,D,h,m,s = timescale.time.convert_julian(tt, format='tuple') # calculate time in year-decimal format time_decimal = timescale.time.convert_calendar_decimal(Y, M, day=D, hour=h, minute=m, second=s) # radius of the Earth radius = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2 + XYZ[:,2]**2) # geocentric latitude (radians) latitude = np.arctan(XYZ[:,2] / np.sqrt(XYZ[:,0]**2.0 + XYZ[:,1]**2.0)) # geocentric colatitude (radians) theta = (np.pi/2.0 - latitude) # calculate longitude (radians) phi = np.arctan2(XYZ[:,1], XYZ[:,0]) # calculate angular coordinates of mean/secular pole at time mpx, mpy, fl = timescale.eop.iers_mean_pole(time_decimal, convention=convention) # read and interpolate IERS daily polar motion values px, py = timescale.eop.iers_polar_motion(MJD, k=3, s=0) # calculate differentials from mean/secular pole positions # using the latest definition from IERS Conventions (2010) mx = px - mpx my = -(py - mpy) # number of points n = np.maximum(len(time_decimal), len(theta)) # conversion factors in latitude, longitude, and radial directions dfactor = np.zeros((n, 3)) dfactor[:,0] = -l2*atr*(omega**2 * radius**2)/(gamma_0) dfactor[:,1] = l2*atr*(omega**2 * radius**2)/(gamma_0) dfactor[:,2] = -h2*atr*(omega**2 * radius**2)/(2.0*gamma_0) # calculate pole tide displacements (meters) S = np.zeros((n, 3)) # pole tide displacements in latitude, longitude, and radial directions S[:,0] = dfactor[:,0]*np.cos(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi)) S[:,1] = dfactor[:,1]*np.cos(theta)*(mx*np.sin(phi) - my*np.cos(phi)) S[:,2] = dfactor[:,2]*np.sin(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi)) # rotation matrix R = np.zeros((3, 3, n)) R[0,0,:] = np.cos(phi)*np.cos(theta) R[0,1,:] = -np.sin(phi) R[0,2,:] = np.cos(phi)*np.sin(theta) R[1,0,:] = np.sin(phi)*np.cos(theta) R[1,1,:] = np.cos(phi) R[1,2,:] = np.sin(phi)*np.sin(theta) R[2,0,:] = -np.sin(theta) R[2,2,:] = np.cos(theta) # rotate displacements to ECEF coordinates dxt = np.einsum('ti...,jit...->tj...', S, R) # return the pole tide displacements # in Cartesian coordinates return dxt
# PURPOSE: estimate ocean pole tides in Cartesian coordinates
[docs] def ocean_pole_tide( t: np.ndarray, XYZ: np.ndarray, UXYZ: np.ndarray, deltat: float = 0.0, gamma_0: float = 9.780325, a_axis: float = 6378136.3, GM: float = 3.986004418e14, omega: float = 7.2921151467e-5, rho_w: float = 1025.0, g2: complex = 0.6870 + 0.0036j, convention: str = '2018' ): """ Estimate ocean pole tide displacements in Cartesian coordinates :cite:p:`Desai:2002ev,Desai:2015jr,Petit:2010tp` Parameters ---------- t: np.ndarray Time (days relative to January 1, 1992) XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) UXYZ: np.ndarray Ocean pole tide values from Desai (2002) deltat: float or np.ndarray, default 0.0 time correction for converting to Ephemeris Time (days) a_axis: float, default 6378136.3 Semi-major axis of the Earth (meters) gamma_0: float, default 9.780325 Normal gravity (m/s^2) GM: float, default 3.986004418e14 Earth's gravitational constant [m^3/s^2] omega: float, default 7.2921151467e-5 Earth's rotation rate (radians/second) rho_w: float, default 1025.0 Density of sea water [kg/m^3] g2: complex, default 0.6870 + 0.0036j Degree-2 Love number differential (1 + k2 - h2) convention: str, default '2018' IERS Mean or Secular Pole Convention - ``'2003'`` - ``'2010'`` - ``'2015'`` - ``'2018'`` Returns ------- dxt: np.ndarray Load pole tide displacements in meters in Cartesian coordinates """ # degrees and arcseconds to radians dtr = np.pi/180.0 atr = np.pi/648000.0 # convert time to Terrestial Time (TT) tt = t + _jd_tide + deltat # convert time to Modified Julian Days (MJD) MJD = tt - _jd_mjd # convert Julian days to calendar dates Y,M,D,h,m,s = timescale.time.convert_julian(tt, format='tuple') # calculate time in year-decimal format time_decimal = timescale.time.convert_calendar_decimal(Y, M, day=D, hour=h, minute=m, second=s) # radius of the Earth radius = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2 + XYZ[:,2]**2) # geocentric latitude (radians) latitude = np.arctan(XYZ[:,2] / np.sqrt(XYZ[:,0]**2.0 + XYZ[:,1]**2.0)) # geocentric colatitude (radians) theta = (np.pi/2.0 - latitude) # calculate longitude (radians) phi = np.arctan2(XYZ[:,1], XYZ[:,0]) # universal gravitational constant [N*m^2/kg^2] G = 6.67430e-11 # calculate angular coordinates of mean/secular pole at time mpx, mpy, fl = timescale.eop.iers_mean_pole(time_decimal, convention=convention) # read and interpolate IERS daily polar motion values px, py = timescale.eop.iers_polar_motion(MJD, k=3, s=0) # calculate differentials from mean/secular pole positions # using the latest definition from IERS Conventions (2010) mx = px - mpx my = -(py - mpy) # pole tide displacement factors Hp = np.sqrt(8.0*np.pi/15.0)*(omega**2 * a_axis**4)/GM K = 4.0*np.pi*G*rho_w*Hp*a_axis/(3.0*gamma_0) # number of points n = np.maximum(len(time_decimal), len(theta)) # calculate ocean pole tide displacements (meters) dxt = np.zeros((n, 3)) for i in range(3): dxt[:,i] = K*atr*np.real( (mx*g2.real + my*g2.imag)*UXYZ.real[:,i] + (my*g2.real - mx*g2.imag)*UXYZ.imag[:,i]) # return the ocean pole tide displacements # in Cartesian coordinates return dxt
# get IERS parameters _iers = pyTMD.spatial.datum(ellipsoid='IERS', units='MKS') # PURPOSE: estimate solid Earth tides due to gravitational attraction
[docs] def solid_earth_tide( t: np.ndarray, XYZ: np.ndarray, SXYZ: np.ndarray, LXYZ: np.ndarray, a_axis: float = _iers.a_axis, tide_system: str = 'tide_free', **kwargs ): """ Compute the solid Earth tides due to the gravitational attraction of the moon and sun :cite:p:`Mathews:1991kv,Mathews:1997js,Ries:1992ip,Wahr:1981ea` Parameters ---------- t: np.ndarray Time (days relative to January 1, 1992) XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the moon (meters) a_axis: float, default 6378136.3 Semi-major axis of the Earth (meters) tide_system: str, default 'tide_free' Permanent tide system for the output solid Earth tide - ``'tide_free'``: no permanent direct and indirect tidal potentials - ``'mean_tide'``: permanent tidal potentials (direct and indirect) Returns ------- dxt: np.ndarray Solid Earth tide in meters in Cartesian coordinates """ # set default keyword arguments # nominal Love and Shida numbers kwargs.setdefault('h2', 0.6078) kwargs.setdefault('l2', 0.0847) kwargs.setdefault('h3', 0.292) kwargs.setdefault('l3', 0.015) # mass ratios between earth and sun/moon kwargs.setdefault('mass_ratio_solar', 332946.0482) kwargs.setdefault('mass_ratio_lunar', 0.0123000371) # validate output tide system assert tide_system.lower() in ('tide_free', 'mean_tide') # number of input coordinates nt = len(np.atleast_1d(t)) # convert time to Modified Julian Days (MJD) MJD = t + _mjd_tide # scalar product of input coordinates with sun/moon vectors radius = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2 + XYZ[:,2]**2) solar_radius = np.sqrt(SXYZ[:,0]**2 + SXYZ[:,1]**2 + SXYZ[:,2]**2) lunar_radius = np.sqrt(LXYZ[:,0]**2 + LXYZ[:,1]**2 + LXYZ[:,2]**2) solar_scalar = (XYZ[:,0]*SXYZ[:,0] + XYZ[:,1]*SXYZ[:,1] + XYZ[:,2]*SXYZ[:,2])/(radius*solar_radius) lunar_scalar = (XYZ[:,0]*LXYZ[:,0] + XYZ[:,1]*LXYZ[:,1] + XYZ[:,2]*LXYZ[:,2])/(radius*lunar_radius) # compute new h2 and l2 (Mathews et al., 1997) cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius h2 = kwargs['h2'] - 0.0006*(1.0 - 3.0/2.0*cosphi**2) l2 = kwargs['l2'] + 0.0002*(1.0 - 3.0/2.0*cosphi**2) # compute P2 terms P2_solar = 3.0*(h2/2.0 - l2)*solar_scalar**2 - h2/2.0 P2_lunar = 3.0*(h2/2.0 - l2)*lunar_scalar**2 - h2/2.0 # compute P3 terms P3_solar = 5.0/2.0*(kwargs['h3'] - 3.0*kwargs['l3'])*solar_scalar**3 + \ 3.0/2.0*(kwargs['l3'] - kwargs['h3'])*solar_scalar P3_lunar = 5.0/2.0*(kwargs['h3'] - 3.0*kwargs['l3'])*lunar_scalar**3 + \ 3.0/2.0*(kwargs['l3'] - kwargs['h3'])*lunar_scalar # compute terms in direction of sun/moon vectors X2_solar = 3.0*l2*solar_scalar X2_lunar = 3.0*l2*lunar_scalar X3_solar = 3.0*kwargs['l3']/2.0*(5.0*solar_scalar**2 - 1.0) X3_lunar = 3.0*kwargs['l3']/2.0*(5.0*lunar_scalar**2 - 1.0) # factors for sun and moon using IAU estimates of mass ratios F2_solar = kwargs['mass_ratio_solar']*a_axis*(a_axis/solar_radius)**3 F2_lunar = kwargs['mass_ratio_lunar']*a_axis*(a_axis/lunar_radius)**3 F3_solar = kwargs['mass_ratio_solar']*a_axis*(a_axis/solar_radius)**4 F3_lunar = kwargs['mass_ratio_lunar']*a_axis*(a_axis/lunar_radius)**4 # compute total displacement (Mathews et al. 1997) dxt = np.zeros((nt, 3)) for i in range(3): S2 = F2_solar*(X2_solar*SXYZ[:,i]/solar_radius+P2_solar*XYZ[:,i]/radius) L2 = F2_lunar*(X2_lunar*LXYZ[:,i]/lunar_radius+P2_lunar*XYZ[:,i]/radius) S3 = F3_solar*(X3_solar*SXYZ[:,i]/solar_radius+P3_solar*XYZ[:,i]/radius) L3 = F3_lunar*(X3_lunar*LXYZ[:,i]/lunar_radius+P3_lunar*XYZ[:,i]/radius) dxt[:,i] = S2 + L2 + S3 + L3 # corrections for out-of-phase portions of the Love and Shida numbers dxt += _out_of_phase_diurnal(XYZ, SXYZ, LXYZ, F2_solar, F2_lunar) dxt += _out_of_phase_semidiurnal(XYZ, SXYZ, LXYZ, F2_solar, F2_lunar) # corrections for the latitudinal dependence dxt += _latitude_dependence(XYZ, SXYZ, LXYZ, F2_solar, F2_lunar) # corrections for the frequency dependence dxt += _frequency_dependence_diurnal(XYZ, MJD) dxt += _frequency_dependence_long_period(XYZ, MJD) # convert the permanent tide system if specified if (tide_system.lower() == 'mean_tide'): dxt += _free_to_mean(XYZ, h2, l2) # return the solid earth tide return dxt
[docs] def _out_of_phase_diurnal( XYZ: np.ndarray, SXYZ: np.ndarray, LXYZ: np.ndarray, F2_solar: np.ndarray, F2_lunar: np.ndarray ): """ Computes the out-of-phase corrections induced by mantle anelasticity in the diurnal band Parameters ---------- XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the moon (meters) F2_solar: np.ndarray Factors for the sun F2_lunar: np.ndarray Factors for the moon """ # Love and Shida number corrections dhi = -0.0025 dli = -0.0007 # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius cos2phi = cosphi**2 - sinphi**2 sinla = XYZ[:,1]/cosphi/radius cosla = XYZ[:,0]/cosphi/radius # Compute the normalized position vector of the Sun/Moon solar_radius = np.sqrt(np.sum(SXYZ**2, axis=1)) lunar_radius = np.sqrt(np.sum(LXYZ**2, axis=1)) # calculate offsets dr_solar = -3.0*dhi*sinphi*cosphi*F2_solar*SXYZ[:,2]* \ (SXYZ[:,0]*sinla-SXYZ[:,1]*cosla)/solar_radius**2 dr_lunar = -3.0*dhi*sinphi*cosphi*F2_lunar*LXYZ[:,2]* \ (LXYZ[:,0]*sinla-LXYZ[:,1]*cosla)/lunar_radius**2 dn_solar = -3.0*dli*cos2phi*F2_solar*SXYZ[:,2]* \ (SXYZ[:,0]*sinla-SXYZ[:,1]*cosla)/solar_radius**2 dn_lunar = -3.0*dli*cos2phi*F2_lunar*LXYZ[:,2]* \ (LXYZ[:,0]*sinla-LXYZ[:,1]*cosla)/lunar_radius**2 de_solar = -3.0*dli*sinphi*F2_solar*SXYZ[:,2]* \ (SXYZ[:,0]*cosla+SXYZ[:,1]*sinla)/solar_radius**2 de_lunar = -3.0*dli*sinphi*F2_lunar*LXYZ[:,2]* \ (LXYZ[:,0]*cosla+LXYZ[:,1]*sinla)/lunar_radius**2 # add solar and lunar offsets DR = dr_solar + dr_lunar DN = dn_solar + dn_lunar DE = de_solar + de_lunar # compute corrections DX = DR*cosla*cosphi - DE*sinla - DN*cosla*sinphi DY = DR*sinla*cosphi + DE*cosla - DN*sinla*sinphi DZ = DR*sinphi + DN*cosphi # return the corrections return np.c_[DX, DY, DZ]
[docs] def _out_of_phase_semidiurnal( XYZ: np.ndarray, SXYZ: np.ndarray, LXYZ: np.ndarray, F2_solar: np.ndarray, F2_lunar: np.ndarray ): """ Computes the out-of-phase corrections induced by mantle anelasticity in the semi-diurnal band Parameters ---------- XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the moon (meters) F2_solar: np.ndarray Factors for the sun F2_lunar: np.ndarray Factors for the moon """ # Love and Shida number corrections dhi = -0.0022 dli = -0.0007 # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius sinla = XYZ[:,1]/cosphi/radius cosla = XYZ[:,0]/cosphi/radius cos2la = cosla**2 - sinla**2 sin2la = 2.0*cosla*sinla # Compute the normalized position vector of the Sun/Moon solar_radius = np.sqrt(np.sum(SXYZ**2, axis=1)) lunar_radius = np.sqrt(np.sum(LXYZ**2, axis=1)) # calculate offsets dr_solar = -3.0/4.0*dhi*cosphi**2*F2_solar * \ ((SXYZ[:,0]**2-SXYZ[:,1]**2)*sin2la-2.0*SXYZ[:,0]*SXYZ[:,1]*cos2la) / \ solar_radius**2 dr_lunar = -3.0/4.0*dhi*cosphi**2*F2_lunar * \ ((LXYZ[:,0]**2-LXYZ[:,1]**2)*sin2la-2.0*LXYZ[:,0]*LXYZ[:,1]*cos2la) / \ lunar_radius**2 dn_solar = 3.0/2.0*dli*sinphi*cosphi*F2_solar * \ ((SXYZ[:,0]**2-SXYZ[:,1]**2)*sin2la-2.0*SXYZ[:,0]*SXYZ[:,1]*cos2la) / \ solar_radius**2 dn_lunar = 3.0/2.0*dli*sinphi*cosphi*F2_lunar * \ ((LXYZ[:,0]**2-LXYZ[:,1]**2)*sin2la-2.0*LXYZ[:,0]*LXYZ[:,1]*cos2la) / \ lunar_radius**2 de_solar = -3.0/2.0*dli*cosphi*F2_solar * \ ((SXYZ[:,0]**2-SXYZ[:,1]**2)*cos2la+2.0*SXYZ[:,0]*SXYZ[:,1]*sin2la) / \ solar_radius**2 de_lunar = -3.0/2.0*dli*cosphi*F2_lunar * \ ((LXYZ[:,0]**2-LXYZ[:,1]**2)*cos2la+2.0*LXYZ[:,0]*LXYZ[:,1]*sin2la) / \ lunar_radius**2 # add solar and lunar offsets DR = dr_solar + dr_lunar DN = dn_solar + dn_lunar DE = de_solar + de_lunar # compute corrections DX = DR*cosla*cosphi - DE*sinla - DN*cosla*sinphi DY = DR*sinla*cosphi + DE*cosla - DN*sinla*sinphi DZ = DR*sinphi + DN*cosphi # return the corrections return np.c_[DX, DY, DZ]
[docs] def _latitude_dependence( XYZ: np.ndarray, SXYZ: np.ndarray, LXYZ: np.ndarray, F2_solar: np.ndarray, F2_lunar: np.ndarray ): r""" Computes the corrections induced by the latitude of the dependence given by L\ :sup:`1` Parameters ---------- XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the moon (meters) F2_solar: np.ndarray Factors for the sun F2_lunar: np.ndarray Factors for the moon """ # Love/Shida number corrections (diurnal and semi-diurnal) l1d = 0.0012 l1sd = 0.0024 # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius sinla = XYZ[:,1]/cosphi/radius cosla = XYZ[:,0]/cosphi/radius cos2la = cosla**2 - sinla**2 sin2la = 2.0*cosla*sinla # Compute the normalized position vector of the Sun/Moon solar_radius = np.sqrt(np.sum(SXYZ**2, axis=1)) lunar_radius = np.sqrt(np.sum(LXYZ**2, axis=1)) # calculate offsets for the diurnal band dn_d_solar = -l1d*sinphi**2*F2_solar*SXYZ[:,2] * \ (SXYZ[:,0]*cosla+SXYZ[:,1]*sinla)/solar_radius**2 dn_d_lunar = -l1d*sinphi**2*F2_lunar*LXYZ[:,2] * \ (LXYZ[:,0]*cosla+LXYZ[:,1]*sinla)/lunar_radius**2 de_d_solar = l1d*sinphi*(cosphi**2-sinphi**2)*F2_solar*SXYZ[:,2] * \ (SXYZ[:,0]*sinla-SXYZ[:,1]*cosla)/solar_radius**2 de_d_lunar = l1d*sinphi*(cosphi**2-sinphi**2)*F2_lunar*LXYZ[:,2] * \ (LXYZ[:,0]*sinla-LXYZ[:,1]*cosla)/lunar_radius**2 # calculate offsets for the semi-diurnal band dn_s_solar = -l1sd/2.0*sinphi*cosphi*F2_solar * \ ((SXYZ[:,0]**2-SXYZ[:,1]**2)*cos2la+2.0*SXYZ[:,0]*SXYZ[:,1]*sin2la) / \ solar_radius**2 dn_s_lunar =-l1sd/2.0*sinphi*cosphi*F2_lunar * \ ((LXYZ[:,0]**2-LXYZ[:,1]**2)*cos2la+2.0*LXYZ[:,0]*LXYZ[:,1]*sin2la) / \ lunar_radius**2 de_s_solar =-l1sd/2.0*sinphi**2*cosphi*F2_solar * \ ((SXYZ[:,0]**2-SXYZ[:,1]**2)*sin2la-2.0*SXYZ[:,0]*SXYZ[:,1]*cos2la) / \ solar_radius**2 de_s_lunar =-l1sd/2.0*sinphi**2*cosphi*F2_lunar * \ ((LXYZ[:,0]**2-LXYZ[:,1]**2)*sin2la-2.0*LXYZ[:,0]*LXYZ[:,1]*cos2la) / \ lunar_radius**2 # add solar and lunar offsets (diurnal and semi-diurnal) DN = 3.0*(dn_d_solar + dn_d_lunar + dn_s_solar + dn_s_lunar) DE = 3.0*(de_d_solar + de_d_lunar + de_s_solar + de_s_lunar) # compute combined diurnal and semi-diurnal corrections DX = -DE*sinla - DN*cosla*sinphi DY = DE*cosla - DN*sinla*sinphi DZ = DN*cosphi # return the corrections return np.c_[DX, DY, DZ]
[docs] def _frequency_dependence_diurnal( XYZ: np.ndarray, MJD: np.ndarray ): """ Computes the in-phase and out-of-phase corrections induced by mantle anelasticity in the diurnal band Parameters ---------- XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) MJD: np.ndarray Modified Julian Day (MJD) """ # number of time steps nt = len(np.atleast_1d(MJD)) # Corrections to Diurnal Tides for Frequency Dependence # of Love and Shida Number Parameters # reduced version of table 7.3a from IERS conventions table = np.array([ [-3.0, 0.0, 2.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0], [-3.0, 2.0, 0.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0], [-2.0, 0.0, 1.0, -1.0, 0.0, -0.02, 0.0, 0.0, 0.0], [-2.0, 0.0, 1.0, 0.0, 0.0, -0.08, 0.0, -0.01, 0.01], [-2.0, 2.0, -1.0, 0.0, 0.0, -0.02, 0.0, 0.0, 0.0], [-1.0, 0.0, 0.0,-1.0, 0.0, -0.10, 0.0, 0.0, 0.0], [-1.0, 0.0, 0.0, 0.0, 0.0, -0.51, 0.0, -0.02, 0.03], [-1.0, 2.0, 0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0], [0.0, -2.0, 1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0], [0.0, 0.0, -1.0, 0.0, 0.0, 0.02, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0, 0.06, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 1.0, 0.0, 0.01, 0.0, 0.0, 0.0], [0.0, 2.0, -1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0], [1.0, -3.0, 0.0, 0.0, 1.0, -0.06, 0.0, 0.0, 0.0], [1.0, -2.0, 0.0, -1.0, 0.0, 0.01, 0.0, 0.0, 0.0], [1.0, -2.0, 0.0, 0.0, 0.0, -1.23, -0.07, 0.06, 0.01], [1.0, -1.0, 0.0, 0.0,-1.0, 0.02, 0.0, 0.0, 0.0], [1.0, -1.0, 0.0, 0.0, 1.0, 0.04, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, -1.0, 0.0, -0.22, 0.01, 0.01, 0.0], [1.0, 0.0, 0.0, 0.0, 0.0, 12.00, -0.80, -0.67, -0.03], [1.0, 0.0, 0.0, 1.0, 0.0, 1.73, -0.12, -0.10, 0.0], [1.0, 0.0, 0.0, 2.0, 0.0, -0.04, 0.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0, -1.0, -0.50, -0.01, 0.03, 0.0], [1.0, 1.0, 0.0, 0.0, 1.0, 0.01, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 1.0, -1.0, -0.01, 0.0, 0.0, 0.0], [1.0, 2.0, -2.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0], [1.0, 2.0, 0.0, 0.0, 0.0, -0.11, 0.01, 0.01, 0.0], [2.0, -2.0, 1.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0], [2.0, 0.0,-1.0, 0.0, 0.0, -0.02, 0.0, 0.0, 0.0], [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [3.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0] ]) # get phase angles (Doodson arguments) TAU, S, H, P, ZNS, PS = pyTMD.astro.doodson_arguments(MJD) # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius sinla = XYZ[:,1]/cosphi/radius cosla = XYZ[:,0]/cosphi/radius zla = np.arctan2(XYZ[:,1], XYZ[:,0]) # compute corrections (Mathews et al. 1997) DX = np.zeros((nt)) DY = np.zeros((nt)) DZ = np.zeros((nt)) # iterate over rows in the table for i, row in enumerate(table): thetaf = TAU + S*row[0] + H*row[1] + P*row[2] + \ ZNS*row[3] + PS*row[4] dr = 2.0*row[5]*sinphi*cosphi*np.sin(thetaf + zla) + \ 2.0*row[6]*sinphi*cosphi*np.cos(thetaf + zla) dn = row[7]*(cosphi**2 - sinphi**2)*np.sin(thetaf + zla) + \ row[8]*(cosphi**2 - sinphi**2)*np.cos(thetaf + zla) de = row[7]*sinphi*np.cos(thetaf + zla) - \ row[8]*sinphi*np.sin(thetaf + zla) DX += 1e-3*(dr*cosla*cosphi - de*sinla - dn*cosla*sinphi) DY += 1e-3*(dr*sinla*cosphi + de*cosla - dn*sinla*sinphi) DZ += 1e-3*(dr*sinphi + dn*cosphi) # return the corrections return np.c_[DX, DY, DZ]
[docs] def _frequency_dependence_long_period( XYZ: np.ndarray, MJD: np.ndarray ): """ Computes the in-phase and out-of-phase corrections induced by mantle anelasticity in the long-period band Parameters ---------- XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) MJD: np.ndarray Modified Julian Day (MJD) """ # number of time steps nt = len(np.atleast_1d(MJD)) # Corrections to Long-Period Tides for Frequency Dependence # of Love and Shida Number Parameters # reduced version of table 7.3b from IERS conventions table = np.array([ [0.0, 0.0, 0.0, 1.0, 0.0, 0.47, 0.23, 0.16, 0.07], [0.0, 2.0, 0.0, 0.0, 0.0, -0.20, -0.12, -0.11, -0.05], [1.0, 0.0, -1.0, 0.0, 0.0, -0.11, -0.08, -0.09, -0.04], [2.0, 0.0, 0.0, 0.0, 0.0, -0.13, -0.11, -0.15, -0.07], [2.0, 0.0, 0.0, 1.0, 0.0, -0.05, -0.05, -0.06, -0.03] ]) # get phase angles (Doodson arguments) TAU, S, H, P, ZNS, PS = pyTMD.astro.doodson_arguments(MJD) # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius sinla = XYZ[:,1]/cosphi/radius cosla = XYZ[:,0]/cosphi/radius # compute corrections (Mathews et al. 1997) DX = np.zeros((nt)) DY = np.zeros((nt)) DZ = np.zeros((nt)) # iterate over rows in the table for i, row in enumerate(table): thetaf = S*row[0] + H*row[1] + P*row[2] + ZNS*row[3] + PS*row[4] dr = row[5]*(3.0*sinphi**2 - 1.0)*np.cos(thetaf)/2.0 + \ row[7]*(3.0*sinphi**2 - 1.0)*np.sin(thetaf)/2.0 dn = row[6]*(2.0*cosphi*sinphi)*np.cos(thetaf) + \ row[8]*(2.0*cosphi*sinphi)*np.sin(thetaf) de = 0.0 DX += 1e-3*(dr*cosla*cosphi - de*sinla - dn*cosla*sinphi) DY += 1e-3*(dr*sinla*cosphi + de*cosla - dn*sinla*sinphi) DZ += 1e-3*(dr*sinphi + dn*cosphi) # return the corrections return np.c_[DX, DY, DZ]
[docs] def _free_to_mean( XYZ: np.ndarray, h2: float | np.ndarray, l2: float | np.ndarray ): """ Calculate offsets for converting the permanent tide from a tide-free to a mean-tide state Parameters ---------- XYZ: np.ndarray Cartesian coordinates of the prediction points (meters) h2: float or np.ndarray Degree-2 Love number of vertical displacement l2: float or np.ndarray Degree-2 Love (Shida) number of horizontal displacement """ # Compute the normalized position vector of coordinates radius = np.sqrt(np.sum(XYZ**2, axis=1)) sinphi = XYZ[:,2]/radius cosphi = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2)/radius sinla = XYZ[:,1]/cosphi/radius cosla = XYZ[:,0]/cosphi/radius # time-independent constituent of amplitude (Mathews et al. 1997) H0 = -0.31460 # in Mathews et al. (1997): dR0=-0.1196 m with h2=0.6026 dR0 = np.sqrt(5.0/(4.0*np.pi))*h2*H0 # in Mathews et al. (1997): dN0=-0.0247 m with l2=0.0831 dN0 = np.sqrt(45.0/(16.0*np.pi))*l2*H0 # use double angle formula for sin(2*phi) dr = dR0*(3.0/2.0*sinphi**2 - 1.0/2.0) dn = 2.0*dN0*cosphi*sinphi # compute as an additive correction (Mathews et al. 1997) DX = dr*cosla*cosphi - dn*cosla*sinphi DY = dr*sinla*cosphi - dn*sinla*sinphi DZ = dr*sinphi + dn*cosphi # return the corrections return np.c_[DX, DY, DZ]