diff --git a/doc/source/api_reference/arguments.rst b/doc/source/api_reference/arguments.rst index 55a607a..23bec0e 100644 --- a/doc/source/api_reference/arguments.rst +++ b/doc/source/api_reference/arguments.rst @@ -28,6 +28,8 @@ Calling Sequence .. autofunction:: pyTMD.arguments.nodal +.. autofunction:: pyTMD.arguments.frequency + .. autofunction:: pyTMD.arguments._arguments_table .. autofunction:: pyTMD.arguments._minor_table diff --git a/doc/source/api_reference/predict.rst b/doc/source/api_reference/predict.rst index 737395c..75d389b 100644 --- a/doc/source/api_reference/predict.rst +++ b/doc/source/api_reference/predict.rst @@ -32,6 +32,10 @@ Calling Sequence .. autofunction:: pyTMD.predict.infer_minor +.. autofunction:: pyTMD.predict._infer_short_period + +.. autofunction:: pyTMD.predict._infer_long_period + .. autofunction:: pyTMD.predict.equilibrium_tide .. autofunction:: pyTMD.predict.load_pole_tide diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index 66858b7..741ff31 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" arguments.py -Written by Tyler Sutterley (08/2024) +Written by Tyler Sutterley (09/2024) Calculates the nodal corrections for tidal constituents Modification of ARGUMENTS fortran subroutine by Richard Ray 03/1999 @@ -38,6 +38,7 @@ Ocean Tides", Journal of Atmospheric and Oceanic Technology, (2002). UPDATE HISTORY: + Updated 09/2024: add function to calculate tidal angular frequencies Updated 08/2024: add support for constituents in PERTH5 tables add back nodal arguments from PERTH3 for backwards compatibility Updated 01/2024: add function to create arguments coefficients table @@ -78,6 +79,7 @@ "coefficients_table", "doodson_number", "nodal", + "frequency", "_arguments_table", "_minor_table", "_constituent_parameters", @@ -1060,12 +1062,9 @@ def nodal( model time series," *Advances in Water Resources*, 12(3), 109--120, (1989). `doi: 10.1016/0309-1708(89)90017-1 `_ - .. [4] G. D. Egbert and S. Y. Erofeeva, "Efficient Inverse Modeling of - Barotropic Ocean Tides," *Journal of Atmospheric and Oceanic - Technology*, 19(2), 183--204, (2002). - `doi: 10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2`__ - - .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + .. [4] R. D. Ray, "A global ocean tide model from + Topex/Poseidon altimetry: GOT99.2", + NASA Goddard Space Flight Center, TM-1999-209478, (1999). """ # set default keyword arguments kwargs.setdefault('corrections', 'OTIS') @@ -1724,6 +1723,67 @@ def nodal( # return corrections for constituents return (u, f) +def frequency( + constituents: list | np.ndarray, + **kwargs + ): + """ + Calculates the angular frequency for tidal constituents [1]_ + + Parameters + ---------- + constituents: list + tidal constituent IDs + corrections: str, default 'OTIS' + use nodal corrections from OTIS, FES or GOT models + M1: str, default 'perth5' + coefficients to use for M1 tides + + - ``'Doodson'`` + - ``'Ray'`` + - ``'perth5'`` + + Returns + ------- + omega: np.ndarray + angular frequency in radians per second + + References + ---------- + .. [1] R. D. Ray, "A global ocean tide model from + Topex/Poseidon altimetry: GOT99.2", + NASA Goddard Space Flight Center, TM-1999-209478, (1999). + """ + # set default keyword arguments + kwargs.setdefault('corrections', 'OTIS') + kwargs.setdefault('M1', 'perth5') + # set function for astronomical longitudes + # use ASTRO5 routines if not using an OTIS type model + ASTRO5 = kwargs['corrections'] not in ('OTIS','ATLAS','TMD3','netcdf') + # Modified Julian Dates at J2000 + MJD = np.array([51544.5, 51544.55]) + # time interval in seconds + deltat = 86400.0*(MJD[1] - MJD[0]) + # calculate the mean longitudes of the sun and moon + s, h, p, n, pp = pyTMD.astro.mean_longitudes(MJD, ASTRO5=ASTRO5) + + # number of temporal values + nt = len(np.atleast_1d(MJD)) + # initial time conversions + hour = 24.0*np.mod(MJD, 1) + # convert from hours solar time into mean lunar time in degrees + tau = 15.0*hour - s + h + # variable for multiples of 90 degrees (Ray technical note 2017) + k = 90.0 + np.zeros((nt)) + + # determine equilibrium arguments + fargs = np.c_[tau, s, h, p, n, pp, k] + rates = (fargs[1,:] - fargs[0,:])/deltat + fd = np.dot(rates, coefficients_table(constituents, **kwargs)) + # convert to radians per second + omega = 2.0*np.pi*fd/360.0 + return omega + def _arguments_table(**kwargs): """ Arguments table for tidal constituents [1]_ [2]_ diff --git a/pyTMD/compute.py b/pyTMD/compute.py index 144904e..fe98685 100644 --- a/pyTMD/compute.py +++ b/pyTMD/compute.py @@ -65,6 +65,7 @@ use model class attributes for file format and corrections add keyword argument to select nodal corrections type fix to use case insensitive assertions of string argument values + add model attribute for tide model bulk frequencies Updated 08/2024: allow inferring only specific minor constituents use prediction functions for pole tides in cartesian coordinates use rotation matrix to convert from cartesian to spherical @@ -402,7 +403,7 @@ def tide_elevations( if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide[i], hc, c, deltat=deltat[i], corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid @@ -417,7 +418,7 @@ def tide_elevations( if INFER_MINOR: minor = pyTMD.predict.infer_minor(ts.tide, hc, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) tide.data[:] += minor.data[:] elif (TYPE.lower() == 'time series'): nstation = len(x) @@ -431,7 +432,7 @@ def tide_elevations( if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide, HC, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components @@ -635,7 +636,7 @@ def tide_currents( if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide[i], hc, c, deltat=deltat[i], corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid @@ -650,7 +651,7 @@ def tide_currents( if INFER_MINOR: minor = pyTMD.predict.infer_minor(ts.tide, hc, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) tide[t].data[:] += minor.data[:] elif (TYPE.lower() == 'time series'): nstation = len(x) @@ -664,7 +665,7 @@ def tide_currents( if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide, HC, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components diff --git a/pyTMD/io/model.py b/pyTMD/io/model.py index d9c9340..a1a168e 100644 --- a/pyTMD/io/model.py +++ b/pyTMD/io/model.py @@ -11,6 +11,7 @@ add file_format and nodal correction attributes export database as a dataclass for easier access added variable name and descriptions for long period tides + added attribute for model bulk frequencies for long period tides Updated 08/2024: added attribute for minor constituents to infer allow searching over iterable glob strings in definition files added option to try automatic detection of definition file format @@ -301,6 +302,16 @@ def corrections(self) -> str: else: return part1 + @property + def frequency(self) -> str: + """ + Returns the frequency type for the model + """ + if self.variable in ('tide_lpe',): + return 'long' + else: + return 'short' + @property def file_format(self) -> str: """ diff --git a/pyTMD/predict.py b/pyTMD/predict.py index 0707642..7b85428 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -22,6 +22,7 @@ UPDATE HISTORY: 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 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 @@ -67,6 +68,8 @@ "drift", "time_series", "infer_minor", + "_infer_short_period", + "_infer_long_period", "equilibrium_tide", "load_pole_tide", "ocean_pole_tide", @@ -307,6 +310,11 @@ def infer_minor( time correction for converting to Ephemeris Time (days) corrections: str, default 'OTIS' use nodal corrections from OTIS/ATLAS or GOT/FES models + frequency: str, default 'short' + frequency of tidal constituents to infer + + - 'short': diurnal and semidiurnal constituents + - 'long': long-period constituents minor: list or None, default None tidal constituent IDs @@ -335,9 +343,65 @@ def infer_minor( # set default keyword arguments kwargs.setdefault('deltat', 0.0) kwargs.setdefault('corrections', 'OTIS') + kwargs.setdefault('frequency', 'short') # list of minor constituents kwargs.setdefault('minor', None) + # infer the minor tidal constituents + if kwargs['frequency'] in ('short',): + dh = _infer_short_period(t, zmajor, constituents, **kwargs) + elif kwargs['frequency'] in ('long',): + dh = _infer_long_period(t, zmajor, constituents, **kwargs) + # return the inferred values + return dh + +# PURPOSE: infer short-period tides for minor constituents +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 [1]_ [2]_ + + 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 + + Returns + ------- + dh: np.ndarray + tidal time series for minor constituents + References + ---------- + .. [1] G. D. Egbert and S. Y. Erofeeva, "Efficient Inverse Modeling of + Barotropic Ocean Tides," *Journal of Atmospheric and Oceanic + Technology*, 19(2), 183--204, (2002). + `doi: 10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2`__ + .. [2] R. D. Ray, "A global ocean tide model from + Topex/Poseidon altimetry: GOT99.2", + NASA Goddard Space Flight Center, TM-1999-209478, (1999). + + .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + """ + # set default keyword arguments + kwargs.setdefault('deltat', 0.0) + kwargs.setdefault('corrections', 'OTIS') + # list of minor constituents + kwargs.setdefault('minor', None) # number of constituents npts, nc = np.shape(zmajor) nt = len(np.atleast_1d(t)) @@ -348,7 +412,7 @@ def infer_minor( # 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,9),dtype=np.complex64) + 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)] @@ -424,6 +488,126 @@ def infer_minor( # return the inferred values return dh +# PURPOSE: infer long-period tides for minor constituents +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 [1]_ [2]_ + + 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 + + Returns + ------- + dh: np.ndarray + tidal time series for minor constituents + + References + ---------- + .. [1] R. D. Ray, "A global ocean tide model from + Topex/Poseidon altimetry: GOT99.2", + NASA Goddard Space Flight Center, TM-1999-209478, (1999). + .. [2] R. D. Ray and S. Y. Erofeeva, "Long-period tidal + variations in the length of day", *Journal of Geophysical + Research: Solid Earth*, 119, 1498--1509, (2013). + `doi: 10.1002/2013JB010830 `_ + """ + # set default keyword arguments + kwargs.setdefault('deltat', 0.0) + kwargs.setdefault('corrections', 'OTIS') + # list of minor constituents + kwargs.setdefault('minor', None) + # number of constituents + 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 + 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 + + if (nz < 3): + raise Exception('Not enough constituents for inference') + + # 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)] + + # 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 def equilibrium_tide(t: np.ndarray, lat: np.ndarray): """ diff --git a/scripts/compute_tidal_currents.py b/scripts/compute_tidal_currents.py index 750dcf2..9695359 100755 --- a/scripts/compute_tidal_currents.py +++ b/scripts/compute_tidal_currents.py @@ -103,6 +103,7 @@ drop support for the ascii definition file format use model class attributes for file format and corrections add command line option to select nodal corrections type + use attribute for inferring minor long period constituents Updated 08/2024: allow inferring only specific minor constituents added option to try automatic detection of definition file format changed from 'geotiff' to 'GTiff' and 'cog' formats @@ -345,7 +346,7 @@ def compute_tidal_currents(tide_dir, input_file, output_file, if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide[i], hc, c, deltat=deltat[i], corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid @@ -361,7 +362,7 @@ def compute_tidal_currents(tide_dir, input_file, output_file, if INFER_MINOR: minor = pyTMD.predict.infer_minor(ts.tide, hc, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) tide[t].data[:] += minor.data[:] elif (TYPE == 'time series'): tide[t] = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) @@ -375,7 +376,7 @@ def compute_tidal_currents(tide_dir, input_file, output_file, if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide, HC, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components diff --git a/scripts/compute_tidal_elevations.py b/scripts/compute_tidal_elevations.py index e6e010b..5eae729 100755 --- a/scripts/compute_tidal_elevations.py +++ b/scripts/compute_tidal_elevations.py @@ -106,6 +106,7 @@ drop support for the ascii definition file format use model class attributes for file format and corrections add command line option to select nodal corrections type + use attribute for inferring minor long period constituents Updated 08/2024: allow inferring only specific minor constituents added option to try automatic detection of definition file format changed from 'geotiff' to 'GTiff' and 'cog' formats @@ -353,7 +354,7 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide[i], hc, c, deltat=deltat[i], corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid @@ -368,7 +369,7 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, if INFER_MINOR: minor = pyTMD.predict.infer_minor(ts.tide, hc, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) tide.data[:] += minor.data[:] elif (TYPE == 'time series'): tide = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) @@ -382,7 +383,7 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, if INFER_MINOR: MINOR = pyTMD.predict.infer_minor(ts.tide, HC, c, deltat=deltat, corrections=nodal_corrections, - minor=minor_constituents) + minor=minor_constituents, frequency=model.frequency) else: MINOR = np.ma.zeros_like(TIDE) # add major and minor components