From 939c7a01ee5e9bcf55f347fa3aee6b2f6e842552 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Mon, 15 Jan 2024 11:55:29 +1300 Subject: [PATCH 1/7] refactor: made the inferrence of minor constituents an option --- pyTMD/compute_tide_corrections.py | 36 +++++++++++++++------- pyTMD/io/OTIS.py | 4 +-- scripts/compute_tidal_currents.py | 40 ++++++++++++++++++------- scripts/compute_tidal_elevations.py | 46 +++++++++++++++++++++-------- 4 files changed, 91 insertions(+), 35 deletions(-) diff --git a/pyTMD/compute_tide_corrections.py b/pyTMD/compute_tide_corrections.py index af4420b2..a834abe8 100644 --- a/pyTMD/compute_tide_corrections.py +++ b/pyTMD/compute_tide_corrections.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tide_corrections.py -Written by Tyler Sutterley (12/2023) +Written by Tyler Sutterley (01/2024) Calculates tidal elevations for correcting elevation or imagery data Ocean and Load Tides @@ -59,6 +59,7 @@ interpolate.py: interpolation routines for spatial data UPDATE HISTORY: + Updated 01/2024: made the inferrence of minor constituents an option Updated 12/2023: use new crs class for coordinate reprojection Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 05/2023: use timescale class for time conversion operations @@ -184,6 +185,7 @@ def compute_tide_corrections( METHOD: str = 'spline', EXTRAPOLATE: bool = False, CUTOFF: int | float=10.0, + INFER_MINOR: bool = True, APPLY_FLEXURE: bool = False, FILL_VALUE: float = np.nan, **kwargs @@ -245,8 +247,10 @@ def compute_tide_corrections( Extrapolation cutoff in kilometers Set to ``np.inf`` to extrapolate for all points + INFER_MINOR: bool, default True + Infer the height values for minor tidal constituents APPLY_FLEXURE: bool, default False - Apply ice flexure scaling factor to height constituents + Apply ice flexure scaling factor to height values Only valid for models containing flexure fields FILL_VALUE: float, default np.nan @@ -352,8 +356,12 @@ def compute_tide_corrections( for i in range(nt): TIDE = pyTMD.predict.map(timescale.tide[i], hc, c, deltat=deltat[i], corrections=model.format) - MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, - deltat=deltat[i], corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid tide[:,:,i] = np.reshape((TIDE+MINOR), (ny,nx)) tide.mask[:,:,i] = np.reshape((TIDE.mask | MINOR.mask), (ny,nx)) @@ -362,18 +370,26 @@ def compute_tide_corrections( tide.mask = np.any(hc.mask,axis=1) tide.data[:] = pyTMD.predict.drift(timescale.tide, hc, c, deltat=deltat, corrections=model.format) - minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, - deltat=deltat, corrections=model.format) - tide.data[:] += minor.data[:] + # calculate values for minor constituents by inferrence + if INFER_MINOR: + minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + tide.data[:] += minor.data[:] elif (TYPE.lower() == 'time series'): nstation = len(x) tide = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) tide.mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): - TIDE = pyTMD.predict.time_series(timescale.tide, hc[s,None,:], c, - deltat=deltat, corrections=model.format) - MINOR = pyTMD.predict.infer_minor(timescale.tide, hc[s,None,:], c, + HC = hc[s,None,:] + TIDE = pyTMD.predict.time_series(timescale.tide, HC, c, deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components tide.data[s,:] = TIDE.data[:] + MINOR.data[:] tide.mask[s,:] = (TIDE.mask | MINOR.mask) # replace invalid values with fill value diff --git a/pyTMD/io/OTIS.py b/pyTMD/io/OTIS.py index b7f5e122..4f141179 100644 --- a/pyTMD/io/OTIS.py +++ b/pyTMD/io/OTIS.py @@ -186,7 +186,7 @@ def extract_constants( - ``'OTIS'``: combined global or local solution - ``'TMD3'``: combined global or local netCDF4 solution apply_flexure: bool, default False - Apply ice flexure scaling factor to height constituents + Apply ice flexure scaling factor to height values Returns ------- @@ -491,7 +491,7 @@ def read_constants( - ``'OTIS'``: combined global or local solution - ``'TMD3'``: combined global or local netCDF4 solution apply_flexure: bool, default False - Apply ice flexure scaling factor to height constituents + Apply ice flexure scaling factor to height values Returns ------- diff --git a/scripts/compute_tidal_currents.py b/scripts/compute_tidal_currents.py index 254dbe4a..a9c742d8 100755 --- a/scripts/compute_tidal_currents.py +++ b/scripts/compute_tidal_currents.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tidal_currents.py -Written by Tyler Sutterley (12/2023) +Written by Tyler Sutterley (01/2024) Calculates zonal and meridional tidal currents for an input file Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -56,6 +56,7 @@ -E X, --extrapolate X: Extrapolate with nearest-neighbors -c X, --cutoff X: Extrapolation cutoff in kilometers set to inf to extrapolate for all points + --infer-minor: Infer the current values for minor constituents -f X, --fill-value X: Invalid value for spatial fields -V, --verbose: Verbose output of processing run -M X, --mode X: Permission mode of output file @@ -94,6 +95,7 @@ predict.py: predict tidal values using harmonic constants UPDATE HISTORY: + Updated 01/2024: made the inferrence of minor constituents an option Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 08/2023: changed ESR netCDF4 format to TMD3 format @@ -195,6 +197,7 @@ def compute_tidal_currents(tide_dir, input_file, output_file, METHOD='spline', EXTRAPOLATE=False, CUTOFF=None, + INFER_MINOR=False, FILL_VALUE=-9999.0, VERBOSE=False, MODE=0o775): @@ -305,15 +308,19 @@ def compute_tidal_currents(tide_dir, input_file, output_file, # calculate constituent oscillation hc = amp*np.exp(cph) - # predict tidal currents at time and infer minor corrections + # predict tidal currents at time if (TYPE == 'grid'): tide[t] = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) tide[t].mask = np.zeros((ny,nx,nt),dtype=bool) for i in range(nt): TIDE = pyTMD.predict.map(timescale.tide[i], hc, c, deltat=deltat[i], corrections=model.format) - MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, - deltat=deltat[i], corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid tide[t][:,:,i] = np.reshape((TIDE+MINOR), (ny,nx)) tide[t].mask[:,:,i] = np.reshape((TIDE.mask | MINOR.mask), @@ -323,18 +330,26 @@ def compute_tidal_currents(tide_dir, input_file, output_file, tide[t].mask = np.any(hc.mask,axis=1) tide[t].data[:] = pyTMD.predict.drift(timescale.tide, hc, c, deltat=deltat, corrections=model.format) - minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, - deltat=deltat, corrections=model.format) - tide[t].data[:] += minor.data[:] + # calculate values for minor constituents by inferrence + if INFER_MINOR: + minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + tide[t].data[:] += minor.data[:] elif (TYPE == 'time series'): tide[t] = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) tide[t].mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): # calculate constituent oscillation for station - TIDE = pyTMD.predict.time_series(timescale.tide, hc[s,None,:], c, - deltat=deltat, corrections=model.format) - MINOR = pyTMD.predict.infer_minor(timescale.tide, hc[s,None,:], c, + HC = hc[s,None,:] + TIDE = pyTMD.predict.time_series(timescale.tide, HC, c, deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components tide[t].data[s,:] = TIDE.data[:] + MINOR.data[:] tide[t].mask[s,:] = (TIDE.mask | MINOR.mask) # replace invalid values with fill value @@ -498,6 +513,10 @@ def arguments(): parser.add_argument('--cutoff','-c', type=np.float64, default=10.0, help='Extrapolation cutoff in kilometers') + # infer minor constituents from major + parser.add_argument('--infer-minor', + default=False, action='store_true', + help='Infer the current values for minor constituents') # fill value for output spatial fields parser.add_argument('--fill-value','-f', type=float, default=-9999.0, @@ -543,6 +562,7 @@ def main(): METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, CUTOFF=args.cutoff, + INFER_MINOR=args.infer_minor, FILL_VALUE=args.fill_value, VERBOSE=args.verbose, MODE=args.mode) diff --git a/scripts/compute_tidal_elevations.py b/scripts/compute_tidal_elevations.py index 4e68bc29..c869997d 100755 --- a/scripts/compute_tidal_elevations.py +++ b/scripts/compute_tidal_elevations.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tidal_elevations.py -Written by Tyler Sutterley (12/2023) +Written by Tyler Sutterley (01/2024) Calculates tidal elevations for an input file Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -57,7 +57,8 @@ -E X, --extrapolate X: Extrapolate with nearest-neighbors -c X, --cutoff X: Extrapolation cutoff in kilometers set to inf to extrapolate for all points - --apply-flexure: Apply ice flexure scaling factor to height constituents + --infer-minor: Infer the height values for minor constituents + --apply-flexure: Apply ice flexure scaling factor to height values Only valid for models containing flexure fields -f X, --fill-value X: Invalid value for spatial fields -V, --verbose: Verbose output of processing run @@ -97,6 +98,7 @@ predict.py: predict tidal values using harmonic constants UPDATE HISTORY: + Updated 01/2024: made the inferrence of minor constituents an option Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 08/2023: changed ESR netCDF4 format to TMD3 format @@ -199,6 +201,7 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, METHOD='spline', EXTRAPOLATE=False, CUTOFF=None, + INFER_MINOR=False, APPLY_FLEXURE=False, FILL_VALUE=-9999.0, VERBOSE=False, @@ -312,15 +315,19 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, # calculate constituent oscillation hc = amp*np.exp(cph) - # predict tidal elevations at time and infer minor corrections + # predict tidal elevations at time if (TYPE == 'grid'): tide = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) tide.mask = np.zeros((ny,nx,nt),dtype=bool) for i in range(nt): TIDE = pyTMD.predict.map(timescale.tide[i], hc, c, deltat=deltat[i], corrections=model.format) - MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, - deltat=deltat[i], corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) # add major and minor components and reform grid tide[:,:,i] = np.reshape((TIDE+MINOR), (ny,nx)) tide.mask[:,:,i] = np.reshape((TIDE.mask | MINOR.mask), (ny,nx)) @@ -329,18 +336,26 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, tide.mask = np.any(hc.mask,axis=1) tide.data[:] = pyTMD.predict.drift(timescale.tide, hc, c, deltat=deltat, corrections=model.format) - minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, - deltat=deltat, corrections=model.format) - tide.data[:] += minor.data[:] + # calculate values for minor constituents by inferrence + if INFER_MINOR: + minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + tide.data[:] += minor.data[:] elif (TYPE == 'time series'): tide = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) tide.mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): # calculate constituent oscillation for station - TIDE = pyTMD.predict.time_series(timescale.tide, hc[s,None,:], c, - deltat=deltat, corrections=model.format) - MINOR = pyTMD.predict.infer_minor(timescale.tide, hc[s,None,:], c, + HC = hc[s,None,:] + TIDE = pyTMD.predict.time_series(timescale.tide, HC, c, deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components tide.data[s,:] = TIDE.data[:] + MINOR.data[:] tide.mask[s,:] = (TIDE.mask | MINOR.mask) # replace invalid values with fill value @@ -494,10 +509,14 @@ def arguments(): parser.add_argument('--cutoff','-c', type=np.float64, default=10.0, help='Extrapolation cutoff in kilometers') - # apply flexure scaling factors to height constituents + # infer minor constituents from major + parser.add_argument('--infer-minor', + default=False, action='store_true', + help='Infer the height values for minor constituents') + # apply flexure scaling factors to height values parser.add_argument('--apply-flexure', default=False, action='store_true', - help='Apply ice flexure scaling factor to height constituents') + help='Apply ice flexure scaling factor to height values') # fill value for output spatial fields parser.add_argument('--fill-value','-f', type=float, default=-9999.0, @@ -545,6 +564,7 @@ def main(): EXTRAPOLATE=args.extrapolate, CUTOFF=args.cutoff, APPLY_FLEXURE=args.apply_flexure, + INFER_MINOR=args.infer_minor, FILL_VALUE=args.fill_value, VERBOSE=args.verbose, MODE=args.mode) From 02ef63d8f6b12ba1a45e2f21b216a1994a45899e Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 18 Jan 2024 11:01:44 +1300 Subject: [PATCH 2/7] refactor: 1-liners in compute.py refactor: lunisolar ephemerides functions feat: added more constituent parameters for OTIS/ATLAS predictions fix: add option to return None and not raise error for Doodson numbers\ docs: add more definitions to the glossary --- doc/source/api_reference/astro.rst | 4 + doc/source/api_reference/compute.rst | 56 + .../compute_tide_corrections.rst | 41 - .../getting_started/Getting-Started.rst | 19 +- doc/source/getting_started/Glossary.rst | 48 +- doc/source/index.rst | 2 +- pyTMD/__init__.py | 1 + pyTMD/arguments.py | 40 +- pyTMD/astro.py | 87 +- pyTMD/compute.py | 1274 +++++++++++++++++ pyTMD/compute_tide_corrections.py | 924 +----------- pyTMD/predict.py | 25 +- scripts/compute_SET_displacements.py | 13 +- test/test_atlas_read.py | 3 +- test/test_download_and_read.py | 7 +- test/test_perth3_read.py | 9 +- test/test_solid_earth.py | 23 +- 17 files changed, 1577 insertions(+), 999 deletions(-) create mode 100644 doc/source/api_reference/compute.rst delete mode 100644 doc/source/api_reference/compute_tide_corrections.rst create mode 100644 pyTMD/compute.py diff --git a/doc/source/api_reference/astro.rst b/doc/source/api_reference/astro.rst index ad9001bc..e299e62e 100644 --- a/doc/source/api_reference/astro.rst +++ b/doc/source/api_reference/astro.rst @@ -31,10 +31,14 @@ Calling Sequence .. autofunction:: pyTMD.astro.solar_ecef +.. autofunction:: pyTMD.astro.solar_approximate + .. autofunction:: pyTMD.astro.solar_ephemerides .. autofunction:: pyTMD.astro.lunar_ecef +.. autofunction:: pyTMD.astro.lunar_approximate + .. autofunction:: pyTMD.astro.lunar_ephemerides .. autofunction:: pyTMD.astro.gast diff --git a/doc/source/api_reference/compute.rst b/doc/source/api_reference/compute.rst new file mode 100644 index 00000000..f5df85a2 --- /dev/null +++ b/doc/source/api_reference/compute.rst @@ -0,0 +1,56 @@ +======= +compute +======= + +- Calculates tidal elevations at points and times + + * Can use OTIS format tidal solutions provided by Ohio State University and ESR + * Can use Global Tide Model (GOT) solutions provided by Richard Ray at GSFC + * Can use Finite Element Solution (FES) models provided by AVISO +- Calculates tidal currents at points and times + + * Can use OTIS format tidal solutions provided by Ohio State University and ESR + * Can use Finite Element Solution (FES) models provided by AVISO +- Calculates long-period equilibrium tides (LPET) at points and times + + * Uses the summation of fifteen tidal spectral lines from `Cartwright and Edden, (1973) `_ +- Calculates radial pole load tides (LPT) at points and times + + * Following `IERS Convention (2010) guidelines `_ +- Calculates radial ocean pole load tides (OPT) at points and times + + * Following `IERS Convention (2010) guidelines `_ +- Calculates radial solid earth tides (SET) at points and times + + * Following `IERS Convention (2010) guidelines `_ + +Calling Sequence +---------------- + +.. code-block:: python + + import pyTMD + tide_h = pyTMD.compute.tide_elevations(x, y, delta_time, + DIRECTORY=DIRECTORY, MODEL=MODEL, EPOCH=(2000,1,1,0,0,0), + EPSG=3031, TYPE='drift') + tide_uv = pyTMD.compute.tide_currents(x, y, delta_time, + DIRECTORY=DIRECTORY, MODEL=MODEL, EPOCH=(2000,1,1,0,0,0), + EPSG=3031, TYPE='drift') + +`Source code`__ + +.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/compute.py + +.. autofunction:: pyTMD.compute.corrections + +.. autofunction:: pyTMD.compute.tide_elevations + +.. autofunction:: pyTMD.compute.tide_currents + +.. autofunction:: pyTMD.compute.LPET_elevations + +.. autofunction:: pyTMD.compute.LPT_displacements + +.. autofunction:: pyTMD.compute.OPT_displacements + +.. autofunction:: pyTMD.compute.SET_displacements diff --git a/doc/source/api_reference/compute_tide_corrections.rst b/doc/source/api_reference/compute_tide_corrections.rst deleted file mode 100644 index 438f3250..00000000 --- a/doc/source/api_reference/compute_tide_corrections.rst +++ /dev/null @@ -1,41 +0,0 @@ -======================== -compute_tide_corrections -======================== - -- Calculates tidal elevations for correcting elevation or imagery data - - * Can use OTIS format tidal solutions provided by Ohio State University and ESR - * Can use Global Tide Model (GOT) solutions provided by Richard Ray at GSFC - * Can use Finite Element Solution (FES) models provided by AVISO -- Calculates long-period equilibrium tides (LPET) for correcting elevation or imagery data - - * Uses the summation of fifteen tidal spectral lines from `Cartwright and Edden, (1973) `_ -- Calculates radial pole load tides (LPT) for correcting elevation or imagery data following `IERS Convention (2010) guidelines `_ -- Calculates radial ocean pole load tides (OPT) for correcting elevation or imagery data following `IERS Convention (2010) guidelines `_ -- Calculates radial solid earth tides (SET) for correcting elevation or imagery data following `IERS Convention (2010) guidelines `_ - -Calling Sequence ----------------- - -.. code-block:: python - - import pyTMD - tide = pyTMD.compute_tide_corrections(x, y, delta_time, - DIRECTORY=DIRECTORY, MODEL=MODEL, EPOCH=(2000,1,1,0,0,0), - EPSG=3031, TYPE='drift') - -`Source code`__ - -.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/compute_tide_corrections.py - -.. autofunction:: pyTMD.compute_corrections - -.. autofunction:: pyTMD.compute_tide_corrections - -.. autofunction:: pyTMD.compute_LPET_corrections - -.. autofunction:: pyTMD.compute_LPT_corrections - -.. autofunction:: pyTMD.compute_OPT_corrections - -.. autofunction:: pyTMD.compute_SET_corrections diff --git a/doc/source/getting_started/Getting-Started.rst b/doc/source/getting_started/Getting-Started.rst index 3324ae31..fee8aecb 100644 --- a/doc/source/getting_started/Getting-Started.rst +++ b/doc/source/getting_started/Getting-Started.rst @@ -83,14 +83,23 @@ Programs ######## For users wanting to compute tide corrections for use with numpy arrays or pandas dataframes -`compute_tide_corrections() `_ -is the place to start. It is a function that takes ``x``, ``y``, and ``time`` coordinates and -computes the corresponding tidal elevation. +`pyTMD.compute `_ +is the place to start. +These are a series of functions that take ``x``, ``y``, and ``time`` coordinates and +compute the corresponding tidal elevation or currents. .. code-block:: python - tide = compute_tide_corrections(x, y, delta_time, DIRECTORY=path_to_tide_models, - MODEL='CATS2008', EPSG=3031, EPOCH=(2000,1,1,0,0,0), TYPE='drift', TIME='GPS', + import pyTMD + tide_h = pyTMD.compute.tide_elevations(x, y, delta_time, + DIRECTORY=path_to_tide_models, + MODEL='CATS2008', EPSG=3031, EPOCH=(2000,1,1,0,0,0), + TYPE='drift', TIME='GPS', + METHOD='spline', FILL_VALUE=np.nan) + tide_uv = pyTMD.compute.tide_currents(x, y, delta_time, + DIRECTORY=path_to_tide_models, + MODEL='CATS2008', EPSG=3031, EPOCH=(2000,1,1,0,0,0), + TYPE='drift', TIME='GPS', METHOD='spline', FILL_VALUE=np.nan) diff --git a/doc/source/getting_started/Glossary.rst b/doc/source/getting_started/Glossary.rst index adfa4bc1..5b3d0a3a 100644 --- a/doc/source/getting_started/Glossary.rst +++ b/doc/source/getting_started/Glossary.rst @@ -4,12 +4,25 @@ Glossary .. glossary:: + Aliasing + spurious tidal frequency appearing in an analysis due to inadequate temporal sampling + + Amplitude + one half of the range of a tidal constituent + Anomaly angular distance between the perihelion and the position of a celestial body + Apogee + point of an orbit where a celestial body is furthest to Earth + + see :term:`Perigee` + Aphelion point of an orbit where a celestial body is furthest from the sun + see :term:`Perihelion` + Ascending Node point of an orbit where a celestial body intersects the ecliptic, and the latitude coordinate is increasing @@ -23,7 +36,10 @@ Glossary small, semi-periodic deviations in the motion of the pole of rotation Constituent - tidal oscillation corresponding with a specific periodic motion of the Earth, Sun and Moon + see :term:`Harmonic Constituents` + + Declination + angular distance of an astronomical body north or south of the celestial equator Descending Node point of an orbit where a celestial body intersects the ecliptic, and the latitude coordinate is decreasing @@ -43,6 +59,27 @@ Glossary Epoch fixed point in time used as a reference value + Equilibrium Tide + hypothetical tide produced solely by lunisolar gravitational forces under equilibrium theory + + Frequency + number of cycles in a unit time + + Geoid + equipotential surface coinciding with the ocean surface in the absence of astronomical or dynamical effects + + Harmonic Analysis + mathematical process by which the tides are separated into :term:`Harmonic Constituents` + + Harmonic Constants + amplitude and phase of the :term:`Harmonic Constituents` + + Harmonic Constituents + harmonic elements of the tide-producing force corresponding with a periodic change of relative position of the Earth, Sun and Moon + + Harmonic Prediction + Method of estimating tidal elevations and currents through a combination of the :term:`Harmonic Constituents` + Load Tide elastic deformation of the solid Earth due to ocean and atmospheric tides @@ -55,6 +92,9 @@ Glossary Mean Tide model with both direct and indirect permanent tidal effects retained + Nodal Corrections + adjustments to the amplitudes and phases of harmonic constituents to allow for periodic modulations + Nutation short-period oscillations in the motion of the pole of rotation @@ -67,9 +107,13 @@ Glossary Perigee point of an orbit where a celestial body is closest to Earth + see :term:`Apogee` + Perihelion point of an orbit where a celestial body is closest to the sun + see :term:`Aphelion` + Period time it takes to make one complete revolution @@ -92,7 +136,7 @@ Glossary Tidal Species classification of tidal constituents based on period - + see :term:`Semi-diurnal Tide`, :term:`Diurnal Tide`, and :term:`Long Period Tide` Tidal Stream diff --git a/doc/source/index.rst b/doc/source/index.rst index d39ae9f1..0badaea4 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -33,7 +33,7 @@ ocean, load, solid Earth and pole tides api_reference/arguments.rst api_reference/astro.rst api_reference/check_points.rst - api_reference/compute_tide_corrections.rst + api_reference/compute.rst api_reference/constants.rst api_reference/crs.rst api_reference/ellipse.rst diff --git a/pyTMD/__init__.py b/pyTMD/__init__.py index 55141db0..3ffa1dcc 100644 --- a/pyTMD/__init__.py +++ b/pyTMD/__init__.py @@ -13,6 +13,7 @@ """ import pyTMD.arguments import pyTMD.astro +import pyTMD.compute import pyTMD.eop import pyTMD.interpolate import pyTMD.predict diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index e4b7f2d2..49fd173e 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -42,6 +42,7 @@ add function to calculate the arguments for minor constituents include multiples of 90 degrees as variable following Ray 2017 add function to calculate the Doodson numbers for constituents + add option to return None and not raise error for Doodson numbers Updated 12/2023: made keyword argument for selecting M1 coefficients Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 04/2023: using renamed astro mean_longitudes function @@ -691,6 +692,8 @@ def doodson_number( - ``'Cartwright'`` - ``'Doodson'`` + raise_error: bool, default True + Raise exception if constituent is unsupported Returns ------- @@ -708,6 +711,7 @@ def doodson_number( # set default keyword arguments kwargs.setdefault('corrections', 'OTIS') kwargs.setdefault('formalism', 'Doodson') + kwargs.setdefault('raise_error', True) # validate inputs assert kwargs['formalism'].title() in ('Cartwright', 'Doodson'), \ f'Unknown formalism {kwargs["formalism"]}' @@ -722,9 +726,12 @@ def doodson_number( # get the table of coefficients coefficients = _arguments_table(**kwargs) if isinstance(constituents, str): + # check that given constituents is supported in tidal program + if (constituents.lower() not in cindex) and kwargs['raise_error']: + raise ValueError(f'Unsupported constituent {constituents}') + elif (constituents.lower() not in cindex): + return None # map between given constituents and supported in tidal program - assert constituents.lower() in cindex, \ - f'Unsupported constituent {constituents}' j, = [j for j,val in enumerate(cindex) if (val == constituents.lower())] # extract identifier in formalism if (kwargs['formalism'] == 'Cartwright'): @@ -732,14 +739,19 @@ def doodson_number( numbers = np.array(coefficients[:6,j]) elif (kwargs['formalism'] == 'Doodson'): # convert from coefficients to Doodson number - numbers = _to_doodson_number(coefficients[:,j]) + numbers = _to_doodson_number(coefficients[:,j], **kwargs) else: # output dictionary with Doodson numbers numbers = {} # for each input constituent for i,c in enumerate(constituents): + # check that given constituents is supported in tidal program + if (c.lower() not in cindex) and kwargs['raise_error']: + raise ValueError(f'Unsupported constituent {c}') + elif (c.lower() not in cindex): + numbers[c] = None + continue # map between given constituents and supported in tidal program - assert c.lower() in cindex, f'Unsupported constituent {c}' j, = [j for j,val in enumerate(cindex) if (val == c.lower())] # convert from coefficients to Doodson number if (kwargs['formalism'] == 'Cartwright'): @@ -747,7 +759,7 @@ def doodson_number( numbers[c] = np.array(coefficients[:6,j]) elif (kwargs['formalism'] == 'Doodson'): # convert from coefficients to Doodson number - numbers[c] = _to_doodson_number(coefficients[:,j]) + numbers[c] = _to_doodson_number(coefficients[:,j], **kwargs) # return the Doodson or Cartwright number return numbers @@ -908,7 +920,7 @@ def _minor_table(**kwargs): # return the coefficient table return coef -def _to_doodson_number(coef: list | np.ndarray): +def _to_doodson_number(coef: list | np.ndarray, **kwargs): """ Converts Cartwright numbers into a Doodson number @@ -916,19 +928,29 @@ def _to_doodson_number(coef: list | np.ndarray): ---------- coef: list or np.ndarray Doodson coefficients (Cartwright numbers) for constituent + raise_error: bool, default True + Raise exception if constituent is unsupported Returns ------- DO: float Doodson number for constituent """ + # default keyword arguments + kwargs.setdefault('raise_error', True) # assert length and verify array coef = np.array(coef[:6]) # add 5 to values following Doodson convention (prevent negatives) coef[1:] += 5 - # convert to single number and round off floating point errors - DO = np.sum([v*10**(2-o) for o,v in enumerate(coef)]) - return np.round(DO, decimals=3) + # check for unsupported constituents + if (np.any(coef < 0) or np.any(coef > 10)) and kwargs['raise_error']: + raise ValueError('Unsupported constituent') + elif (np.any(coef < 0) or np.any(coef > 10)): + return None + else: + # convert to single number and round off floating point errors + DO = np.sum([v*10**(2-o) for o,v in enumerate(coef)]) + return np.round(DO, decimals=3) def _from_doodson_number(DO: float | np.ndarray): """ diff --git a/pyTMD/astro.py b/pyTMD/astro.py index 10c3a17c..9a5e9833 100644 --- a/pyTMD/astro.py +++ b/pyTMD/astro.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" astro.py -Written by Tyler Sutterley (12/2023) +Written by Tyler Sutterley (01/2024) Astronomical and nutation routines PYTHON DEPENDENCIES: @@ -16,6 +16,7 @@ Oliver Montenbruck, Practical Ephemeris Calculations, 1989. UPDATE HISTORY: + Updated 01/2024: refactored lunisolar ephemerides functions Updated 12/2023: refactored phase_angles function to doodson_arguments added option to compute mean lunar time using equinox method Updated 05/2023: add wrapper function for nutation angles @@ -430,7 +431,47 @@ def mean_obliquity(MJD: np.ndarray): return atr*polynomial_sum(epsilon0, T) # PURPOSE: compute coordinates of the sun in an ECEF frame -def solar_ecef(MJD: np.ndarray): +def solar_ecef(MJD: np.ndarray, **kwargs): + """ + Wrapper function for calculating the positional coordinates + of the sun in an Earth-centric, Earth-Fixed (ECEF) frame + [1]_ [2]_ [3]_ + + Parameters + ---------- + MJD: np.ndarray + Modified Julian Day (MJD) of input date + ephemerides: str, default 'approximate' + Method for calculating solar ephemerides + + - ``'approximate'``: low-resolution ephemerides + - ``'JPL'``: computed ephemerides from JPL kernels + kwargs: dict + Keyword options for ephemeris calculation + + Returns + ------- + X, Y, Z: np.ndarray + ECEF coordinates of the sun (meters) + + References + ---------- + .. [1] J. Meeus, *Astronomical Algorithms*, 2nd edition, 477 pp., (1998). + .. [2] O. Montenbruck, *Practical Ephemeris Calculations*, + 146 pp., (1989). + .. [3] R. S. Park, W. M. Folkner, and J. G. Williams, and D. H. Boggs, + "The JPL Planetary and Lunar Ephemerides DE440 and DE441", + *The Astronomical Journal*, 161(3), 105, (2021). + `doi: 10.3847/1538-3881/abd414 + `_ + """ + kwargs.setdefault('ephemerides', 'approximate') + if (kwargs['ephemerides'].lower() == 'approximate'): + return solar_approximate(MJD, **kwargs) + elif (kwargs['ephemerides'].upper() == 'JPL'): + return solar_ephemerides(MJD, **kwargs) + +def solar_approximate(MJD, **kwargs): """ Computes approximate positional coordinates of the sun in an Earth-centric, Earth-Fixed (ECEF) frame [1]_ [2]_ @@ -537,7 +578,47 @@ def solar_ephemerides(MJD: np.ndarray, **kwargs): return (X, Y, Z) # PURPOSE: compute coordinates of the moon in an ECEF frame -def lunar_ecef(MJD: np.ndarray): +def lunar_ecef(MJD: np.ndarray, **kwargs): + """ + Wrapper function for calculating the positional coordinates + of the moon in an Earth-centric, Earth-Fixed (ECEF) frame + [1]_ [2]_ [3]_ + + Parameters + ---------- + MJD: np.ndarray + Modified Julian Day (MJD) of input date + ephemerides: str, default 'approximate' + Method for calculating lunar ephemerides + + - ``'approximate'``: low-resolution ephemerides + - ``'JPL'``: computed ephemerides from JPL kernels + kwargs: dict + Keyword options for ephemeris calculation + + Returns + ------- + X, Y, Z: np.ndarray + ECEF coordinates of the moon (meters) + + References + ---------- + .. [1] J. Meeus, *Astronomical Algorithms*, 2nd edition, 477 pp., (1998). + .. [2] O. Montenbruck, *Practical Ephemeris Calculations*, + 146 pp., (1989). + .. [3] R. S. Park, W. M. Folkner, and J. G. Williams, and D. H. Boggs, + "The JPL Planetary and Lunar Ephemerides DE440 and DE441", + *The Astronomical Journal*, 161(3), 105, (2021). + `doi: 10.3847/1538-3881/abd414 + `_ + """ + kwargs.setdefault('ephemerides', 'approximate') + if (kwargs['ephemerides'].lower() == 'approximate'): + return lunar_approximate(MJD, **kwargs) + elif (kwargs['ephemerides'].upper() == 'JPL'): + return lunar_ephemerides(MJD, **kwargs) + +def lunar_approximate(MJD, **kwargs): """ Computes approximate positional coordinates of the moon in an Earth-centric, Earth-Fixed (ECEF) frame [1]_ [2]_ diff --git a/pyTMD/compute.py b/pyTMD/compute.py new file mode 100644 index 00000000..097deed0 --- /dev/null +++ b/pyTMD/compute.py @@ -0,0 +1,1274 @@ +#!/usr/bin/env python +u""" +compute.py +Written by Tyler Sutterley (01/2024) +Calculates tidal elevations for correcting elevation or imagery data +Calculates tidal currents at locations and times + +Ocean and Load Tides +Uses OTIS format tidal solutions provided by Ohio State University and ESR + http://volkov.oce.orst.edu/tides/region.html + https://www.esr.org/research/polar-tide-models/list-of-polar-tide-models/ + ftp://ftp.esr.org/pub/datasets/tmd/ +Global Tide Model (GOT) solutions provided by Richard Ray at GSFC +or Finite Element Solution (FES) models provided by AVISO + +Long-Period Equilibrium Tides (LPET) +Calculates long-period equilibrium tidal elevations for correcting +elevation or imagery data from the summation of fifteen spectral lines + https://doi.org/10.1111/j.1365-246X.1973.tb03420.x + +Load Pole Tides (LPT) +Calculates radial load pole tide displacements following IERS Convention +(2010) guidelines for correcting elevation or imagery data + https://iers-conventions.obspm.fr/chapter7.php + +Ocean Pole Tides (OPT) +Calculates radial ocean pole load tide displacements following IERS Convention +(2010) guidelines for correcting elevation or imagery data + https://iers-conventions.obspm.fr/chapter7.php + +Ocean Pole Tides (SET) +Calculates radial Solid Earth tide displacements following IERS Convention +(2010) guidelines for correcting elevation or imagery data + https://iers-conventions.obspm.fr/chapter7.php + +PYTHON DEPENDENCIES: + numpy: Scientific Computing Tools For Python + https://numpy.org + https://numpy.org/doc/stable/user/numpy-for-matlab-users.html + scipy: Scientific Tools for Python + https://docs.scipy.org/doc/ + netCDF4: Python interface to the netCDF C library + https://unidata.github.io/netcdf4-python/netCDF4/index.html + pyproj: Python interface to PROJ library + https://pypi.org/project/pyproj/ + +PROGRAM DEPENDENCIES: + time.py: utilities for calculating time operations + spatial: utilities for reading, writing and operating on spatial data + utilities.py: download and management utilities for syncing files + arguments.py: load the nodal corrections for tidal constituents + astro.py: computes the basic astronomical mean longitudes + crs.py: Coordinate Reference System (CRS) routines + predict.py: predict tide values using harmonic constants + 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 netcdf models + io/GOT.py: extract tidal harmonic constants from GSFC GOT models + io/FES.py: extract tidal harmonic constants from FES tide models + interpolate.py: interpolation routines for spatial data + +UPDATE HISTORY: + Updated 01/2024: made the inferrence of minor constituents an option + refactored lunisolar ephemerides functions + renamed module to compute and added tidal currents function + Updated 12/2023: use new crs class for coordinate reprojection + Updated 08/2023: changed ESR netCDF4 format to TMD3 format + Updated 05/2023: use timescale class for time conversion operations + use defaults from eop module for pole tide and EOP files + add option for using higher resolution ephemerides from JPL + Updated 04/2023: added function for radial solid earth tides + using pathlib to define and expand paths + Updated 03/2023: add basic variable typing to function inputs + added function for long-period equilibrium tides + added function for radial load pole tides + added function for radial ocean pole tides + Updated 12/2022: refactored tide read and prediction programs + Updated 11/2022: place some imports within try/except statements + use f-strings for formatting verbose or ascii output + Updated 05/2022: added ESR netCDF4 formats to list of model types + updated keyword arguments to read tide model programs + added option to apply flexure to heights for applicable models + Updated 04/2022: updated docstrings to numpy documentation format + Updated 12/2021: added function to calculate a tidal time series + verify coordinate dimensions for each input data type + added option for converting from LORAN times to UTC + Updated 09/2021: refactor to use model class for files and attributes + Updated 07/2021: can use numpy datetime arrays as input time variable + added function for determining the input spatial variable type + added check that tide model directory is accessible + Updated 06/2021: added new Gr1km-v2 1km Greenland model from ESR + add try/except for input projection strings + Updated 05/2021: added option for extrapolation cutoff in kilometers + Updated 03/2021: added TPXO9-atlas-v4 in binary OTIS format + simplified netcdf inputs to be similar to binary OTIS read program + Updated 02/2021: replaced numpy bool to prevent deprecation warning + Updated 12/2020: added valid data extrapolation with nearest_extrap + Updated 11/2020: added model constituents from TPXO9-atlas-v3 + Updated 08/2020: using builtin time operations. + calculate difference in leap seconds from start of epoch + using conversion protocols following pyproj-2 updates + Updated 07/2020: added function docstrings, FES2014 and TPXO9-atlas-v2 + use merged delta time files combining biannual, monthly and daily files + Updated 03/2020: added TYPE, TIME, FILL_VALUE and METHOD options + Written 03/2020 +""" +from __future__ import print_function, annotations + +import logging +import pathlib +import numpy as np +import scipy.interpolate +import pyTMD.constants +import pyTMD.crs +import pyTMD.eop +import pyTMD.io +import pyTMD.time +import pyTMD.io.model +import pyTMD.predict +import pyTMD.spatial +import pyTMD.utilities + +# attempt imports +try: + import pyproj +except (AttributeError, ImportError, ModuleNotFoundError) as exc: + logging.critical("pyproj not available") + +# PURPOSE: wrapper function for computing values +def corrections( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + CORRECTION: str = 'ocean', + **kwargs + ): + """ + Wrapper function to compute tide corrections at points and times + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + CORRECTION: str, default 'ocean' + Correction type to compute + + - ``'ocean'``: ocean tide from model constituents + - ``'load'``: load tide from model constituents + - ``'LPET'``: long-period equilibrium tide + - ``'LPT'``: solid earth load pole tide + - ``'OPT'``: ocean pole tide + - ``'SET'``: solid earth tide + **kwargs: dict + keyword arguments for correction functions + + Returns + ------- + values: np.ndarray + tidal correction at coordinates and time in meters + """ + if CORRECTION.lower() in ('ocean', 'load'): + return tide_elevations(x, y, delta_time, **kwargs) + elif (CORRECTION.upper() == 'LPET'): + return LPET_elevations(x, y, delta_time, **kwargs) + elif (CORRECTION.upper() == 'LPT'): + return LPT_displacements(x, y, delta_time, **kwargs) + elif (CORRECTION.upper() == 'OPT'): + return OPT_displacements(x, y, delta_time, **kwargs) + elif (CORRECTION.upper() == 'SET'): + return SET_displacements(x, y, delta_time, **kwargs) + else: + raise ValueError(f'Unrecognized correction type: {CORRECTION}') + +# PURPOSE: compute tides at points and times using tide model algorithms +def tide_elevations( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + DIRECTORY: str | pathlib.Path | None = None, + MODEL: str | None = None, + ATLAS_FORMAT: str = 'netcdf', + GZIP: bool = False, + DEFINITION_FILE: str | pathlib.Path | None = None, + EPSG: str | int = 3031, + EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), + TYPE: str or None = 'drift', + TIME: str = 'UTC', + METHOD: str = 'spline', + EXTRAPOLATE: bool = False, + CUTOFF: int | float=10.0, + INFER_MINOR: bool = True, + APPLY_FLEXURE: bool = False, + FILL_VALUE: float = np.nan, + **kwargs + ): + """ + Compute ocean or load tides at points and times from + model constituents + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + DIRECTORY: str or NoneType, default None + working data directory for tide models + MODEL: str or NoneType, default None + Tide model to use in correction + ATLAS_FORMAT: str, default 'netcdf' + ATLAS tide model format + + - ``'OTIS'`` + - ``'netcdf'`` + GZIP: bool, default False + Tide model files are gzip compressed + DEFINITION_FILE: str or NoneType, default None + Tide model definition file for use + EPSG: int, default: 3031 (Polar Stereographic South, WGS84) + Input coordinate system + EPOCH: tuple, default (2000,1,1,0,0,0) + Time period for calculating delta times + TYPE: str or NoneType, default 'drift' + Input data type + + - ``None``: determined from input variable dimensions + - ``'drift'``: drift buoys or satellite/airborne altimetry + - ``'grid'``: spatial grids or images + - ``'time series'``: time series at a single point + TIME: str, default 'UTC' + Time type if need to compute leap seconds to convert to UTC + + - ``'GPS'``: leap seconds needed + - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) + - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) + - ``'UTC'``: no leap seconds needed + - ``'datetime'``: numpy datatime array in UTC + METHOD: str + Interpolation method + + - ```bilinear```: quick bilinear interpolation + - ```spline```: scipy bivariate spline interpolation + - ```linear```, ```nearest```: scipy regular grid interpolations + + EXTRAPOLATE: bool, default False + Extrapolate with nearest-neighbors + CUTOFF: int or float, default 10.0 + Extrapolation cutoff in kilometers + + Set to ``np.inf`` to extrapolate for all points + INFER_MINOR: bool, default True + Infer the height values for minor tidal constituents + APPLY_FLEXURE: bool, default False + Apply ice flexure scaling factor to height values + + Only valid for models containing flexure fields + FILL_VALUE: float, default np.nan + Output invalid value + + Returns + ------- + tide: np.ndarray + tidal elevation at coordinates and time in meters + """ + + # check that tide directory is accessible + if DIRECTORY is not None: + DIRECTORY = pathlib.Path(DIRECTORY).expanduser() + if not DIRECTORY.exists(): + raise FileNotFoundError("Invalid tide directory") + + # validate input arguments + assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') + assert METHOD in ('bilinear', 'spline', 'linear', 'nearest') + + # get parameters for tide model + if DEFINITION_FILE is not None: + model = pyTMD.io.model(DIRECTORY).from_file( + pathlib.Path(DEFINITION_FILE).expanduser()) + else: + model = pyTMD.io.model(DIRECTORY, format=ATLAS_FORMAT, + compressed=GZIP).elevation(MODEL) + + # determine input data type based on variable dimensions + if not TYPE: + TYPE = pyTMD.spatial.data_type(x, y, delta_time) + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon, lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + timescale = pyTMD.time.timescale().from_datetime( + delta_time.flatten()) + else: + timescale = pyTMD.time.timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + # number of time points + nt = len(timescale) + + # read tidal constants and interpolate to grid points + if model.format in ('OTIS', 'ATLAS', 'TMD3'): + amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file, + model.model_file, model.projection, type=model.type, + method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + grid=model.format, apply_flexure=APPLY_FLEXURE) + # use delta time at 2000.0 to match TMD outputs + deltat = np.zeros((nt), dtype=np.float64) + elif (model.format == 'netcdf'): + amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file, + model.model_file, type=model.type, method=METHOD, + extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, + compressed=model.compressed) + # use delta time at 2000.0 to match TMD outputs + deltat = np.zeros((nt), dtype=np.float64) + elif (model.format == 'GOT'): + amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file, + method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + scale=model.scale, compressed=model.compressed) + # delta time (TT - UT1) + deltat = timescale.tt_ut1 + elif (model.format == 'FES'): + amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file, + type=model.type, version=model.version, method=METHOD, + extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, + compressed=model.compressed) + # available model constituents + c = model.constituents + # delta time (TT - UT1) + deltat = timescale.tt_ut1 + + # calculate complex phase in radians for Euler's + cph = -1j*ph*np.pi/180.0 + # calculate constituent oscillation + hc = amp*np.exp(cph) + + # predict tidal elevations at time + if (TYPE.lower() == 'grid'): + ny,nx = np.shape(x) + tide = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) + tide.mask = np.zeros((ny,nx,nt),dtype=bool) + for i in range(nt): + TIDE = pyTMD.predict.map(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components and reform grid + tide[:,:,i] = np.reshape((TIDE+MINOR), (ny,nx)) + tide.mask[:,:,i] = np.reshape((TIDE.mask | MINOR.mask), (ny,nx)) + elif (TYPE.lower() == 'drift'): + tide = np.ma.zeros((nt), fill_value=FILL_VALUE) + tide.mask = np.any(hc.mask,axis=1) + tide.data[:] = pyTMD.predict.drift(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + tide.data[:] += minor.data[:] + elif (TYPE.lower() == 'time series'): + nstation = len(x) + tide = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) + tide.mask = np.zeros((nstation,nt),dtype=bool) + for s in range(nstation): + HC = hc[s,None,:] + TIDE = pyTMD.predict.time_series(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components + tide.data[s,:] = TIDE.data[:] + MINOR.data[:] + tide.mask[s,:] = (TIDE.mask | MINOR.mask) + # replace invalid values with fill value + tide.data[tide.mask] = tide.fill_value + + # return the ocean or load tide correction + return tide + +# PURPOSE: compute tides at points and times using tide model algorithms +def tide_currents( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + DIRECTORY: str | pathlib.Path | None = None, + MODEL: str | None = None, + ATLAS_FORMAT: str = 'netcdf', + GZIP: bool = False, + DEFINITION_FILE: str | pathlib.Path | None = None, + EPSG: str | int = 3031, + EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), + TYPE: str or None = 'drift', + TIME: str = 'UTC', + METHOD: str = 'spline', + EXTRAPOLATE: bool = False, + CUTOFF: int | float=10.0, + INFER_MINOR: bool = True, + FILL_VALUE: float = np.nan, + **kwargs + ): + """ + Compute ocean tide currents at points and times from + model constituents + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + DIRECTORY: str or NoneType, default None + working data directory for tide models + MODEL: str or NoneType, default None + Tide model to use in correction + ATLAS_FORMAT: str, default 'netcdf' + ATLAS tide model format + + - ``'OTIS'`` + - ``'netcdf'`` + GZIP: bool, default False + Tide model files are gzip compressed + DEFINITION_FILE: str or NoneType, default None + Tide model definition file for use + EPSG: int, default: 3031 (Polar Stereographic South, WGS84) + Input coordinate system + EPOCH: tuple, default (2000,1,1,0,0,0) + Time period for calculating delta times + TYPE: str or NoneType, default 'drift' + Input data type + + - ``None``: determined from input variable dimensions + - ``'drift'``: drift buoys or satellite/airborne altimetry + - ``'grid'``: spatial grids or images + - ``'time series'``: time series at a single point + TIME: str, default 'UTC' + Time type if need to compute leap seconds to convert to UTC + + - ``'GPS'``: leap seconds needed + - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) + - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) + - ``'UTC'``: no leap seconds needed + - ``'datetime'``: numpy datatime array in UTC + METHOD: str + Interpolation method + + - ```bilinear```: quick bilinear interpolation + - ```spline```: scipy bivariate spline interpolation + - ```linear```, ```nearest```: scipy regular grid interpolations + + EXTRAPOLATE: bool, default False + Extrapolate with nearest-neighbors + CUTOFF: int or float, default 10.0 + Extrapolation cutoff in kilometers + + Set to ``np.inf`` to extrapolate for all points + INFER_MINOR: bool, default True + Infer the height values for minor tidal constituents + FILL_VALUE: float, default np.nan + Output invalid value + + Returns + ------- + tide: dict + tidal currents at coordinates and times + """ + + # check that tide directory is accessible + if DIRECTORY is not None: + DIRECTORY = pathlib.Path(DIRECTORY).expanduser() + if not DIRECTORY.exists(): + raise FileNotFoundError("Invalid tide directory") + + # validate input arguments + assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') + assert METHOD in ('bilinear', 'spline', 'linear', 'nearest') + + # get parameters for tide model + if DEFINITION_FILE is not None: + model = pyTMD.io.model(DIRECTORY).from_file( + pathlib.Path(DEFINITION_FILE).expanduser()) + else: + model = pyTMD.io.model(DIRECTORY, format=ATLAS_FORMAT, + compressed=GZIP).current(MODEL) + + # determine input data type based on variable dimensions + if not TYPE: + TYPE = pyTMD.spatial.data_type(x, y, delta_time) + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon, lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + timescale = pyTMD.time.timescale().from_datetime( + delta_time.flatten()) + else: + timescale = pyTMD.time.timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + # number of time points + nt = len(timescale) + + # python dictionary with tide model data + tide = {} + # iterate over u and v currents + for t in model.type: + # read tidal constants and interpolate to grid points + if model.format in ('OTIS', 'ATLAS', 'TMD3'): + amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file, + model.model_file['u'], model.projection, type=t, + method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, + grid=model.format) + # use delta time at 2000.0 to match TMD outputs + deltat = np.zeros((nt), dtype=np.float64) + elif (model.format == 'netcdf'): + amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file, + model.model_file[t], type=t, method=METHOD, + extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, + compressed=model.compressed) + # use delta time at 2000.0 to match TMD outputs + deltat = np.zeros((nt), dtype=np.float64) + elif (model.format == 'FES'): + amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file[t], + type=t, version=model.version, method=METHOD, + extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, + compressed=model.compressed) + # available model constituents + c = model.constituents + # delta time (TT - UT1) + deltat = timescale.tt_ut1 + + # calculate complex phase in radians for Euler's + cph = -1j*ph*np.pi/180.0 + # calculate constituent oscillation + hc = amp*np.exp(cph) + + # predict tidal currents at time + if (TYPE.lower() == 'grid'): + ny,nx = np.shape(x) + tide[t] = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) + tide[t].mask = np.zeros((ny,nx,nt),dtype=bool) + for i in range(nt): + TIDE = pyTMD.predict.map(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, + deltat=deltat[i], corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components and reform grid + tide[t][:,:,i] = np.reshape((TIDE+MINOR), (ny,nx)) + tide[t].mask[:,:,i] = np.reshape((TIDE.mask | MINOR.mask), (ny,nx)) + elif (TYPE.lower() == 'drift'): + tide[t] = np.ma.zeros((nt), fill_value=FILL_VALUE) + tide[t].mask = np.any(hc.mask,axis=1) + tide[t].data[:] = pyTMD.predict.drift(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, + deltat=deltat, corrections=model.format) + tide[t].data[:] += minor.data[:] + elif (TYPE.lower() == 'time series'): + nstation = len(x) + tide[t] = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) + tide[t].mask = np.zeros((nstation,nt),dtype=bool) + for s in range(nstation): + HC = hc[s,None,:] + TIDE = pyTMD.predict.time_series(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + # calculate values for minor constituents by inferrence + if INFER_MINOR: + MINOR = pyTMD.predict.infer_minor(timescale.tide, HC, c, + deltat=deltat, corrections=model.format) + else: + MINOR = np.ma.zeros_like(TIDE) + # add major and minor components + tide[t].data[s,:] = TIDE.data[:] + MINOR.data[:] + tide[t].mask[s,:] = (TIDE.mask | MINOR.mask) + # replace invalid values with fill value + tide[t].data[tide[t].mask] = tide[t].fill_value + + # return the ocean tide currents + return tide + +# PURPOSE: compute long-period equilibrium tidal elevations +def LPET_elevations( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + EPSG: str | int = 3031, + EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), + TYPE: str or None = 'drift', + TIME: str = 'UTC', + **kwargs + ): + """ + Compute long-period equilibrium tidal elevations at points and times + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + EPSG: int, default: 3031 (Polar Stereographic South, WGS84) + Input coordinate system + EPOCH: tuple, default (2000,1,1,0,0,0) + Time period for calculating delta times + TYPE: str or NoneType, default 'drift' + Input data type + + - ``None``: determined from input variable dimensions + - ``'drift'``: drift buoys or satellite/airborne altimetry + - ``'grid'``: spatial grids or images + - ``'time series'``: time series at a single point + TIME: str, default 'UTC' + Time type if need to compute leap seconds to convert to UTC + + - ``'GPS'``: leap seconds needed + - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) + - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) + - ``'UTC'``: no leap seconds needed + - ``'datetime'``: numpy datatime array in UTC + FILL_VALUE: float, default np.nan + Output invalid value + + Returns + ------- + tide_lpe: np.ndarray + long-period equilibrium tide at coordinates and time in meters + """ + + # validate input arguments + assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') + # determine input data type based on variable dimensions + if not TYPE: + TYPE = pyTMD.spatial.data_type(x, y, delta_time) + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon, lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + timescale = pyTMD.time.timescale().from_datetime( + delta_time.flatten()) + else: + timescale = pyTMD.time.timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + # number of time points + nt = len(timescale) + # convert tide times to dynamic time + tide_time = timescale.tide + timescale.tt_ut1 + + # predict long-period equilibrium tides at time + if (TYPE == 'grid'): + ny,nx = np.shape(x) + tide_lpe = np.zeros((ny,nx,nt)) + for i in range(nt): + lpet = pyTMD.predict.equilibrium_tide(tide_time[i], lat) + tide_lpe[:,:,i] = np.reshape(lpet,(ny,nx)) + elif (TYPE == 'drift'): + tide_lpe = pyTMD.predict.equilibrium_tide(tide_time, lat) + elif (TYPE == 'time series'): + nstation = len(x) + tide_lpe = np.zeros((nstation,nt)) + for s in range(nstation): + lpet = pyTMD.predict.equilibrium_tide(tide_time, lat[s]) + tide_lpe[s,:] = np.copy(lpet) + + # return the long-period equilibrium tide elevations + return tide_lpe + +# PURPOSE: compute radial load pole tide displacements +# following IERS Convention (2010) guidelines +def LPT_displacements( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + h: float | np.ndarray = 0.0, + EPSG: str | int = 3031, + EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), + TYPE: str or None = 'drift', + TIME: str = 'UTC', + ELLIPSOID: str = 'WGS84', + CONVENTION: str = '2018', + FILL_VALUE: float = np.nan, + **kwargs + ): + """ + Compute radial load pole tide displacements at points and times + following IERS Convention (2010) guidelines + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + h: float or np.ndarray, default 0.0 + height coordinates in meters + EPSG: int, default: 3031 (Polar Stereographic South, WGS84) + Input coordinate system + EPOCH: tuple, default (2000,1,1,0,0,0) + Time period for calculating delta times + TYPE: str or NoneType, default 'drift' + Input data type + + - ``None``: determined from input variable dimensions + - ``'drift'``: drift buoys or satellite/airborne altimetry + - ``'grid'``: spatial grids or images + - ``'time series'``: time series at a single point + TIME: str, default 'UTC' + Time type if need to compute leap seconds to convert to UTC + + - ``'GPS'``: leap seconds needed + - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) + - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) + - ``'UTC'``: no leap seconds needed + - ``'datetime'``: numpy datatime array in UTC + ELLIPSOID: str, default 'WGS84' + Ellipsoid for calculating Earth parameters + CONVENTION: str, default '2018' + IERS Mean or Secular Pole Convention + + - ``'2003'`` + - ``'2010'`` + - ``'2015'`` + - ``'2018'`` + FILL_VALUE: float, default np.nan + Output invalid value + + Returns + ------- + Srad: np.ndarray + solid earth pole tide at coordinates and time in meters + """ + + # validate input arguments + assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') + assert ELLIPSOID in pyTMD._ellipsoids + assert CONVENTION in pyTMD.eop._conventions + # determine input data type based on variable dimensions + if not TYPE: + TYPE = pyTMD.spatial.data_type(x, y, delta_time) + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon,lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + timescale = pyTMD.time.timescale().from_datetime( + delta_time.flatten()) + else: + timescale = pyTMD.time.timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + + # convert dynamic time to Modified Julian Days (MJD) + MJD = timescale.tt - 2400000.5 + # convert Julian days to calendar dates + Y,M,D,h,m,s = pyTMD.time.convert_julian(timescale.tt, format='tuple') + # calculate time in year-decimal format + time_decimal = pyTMD.time.convert_calendar_decimal(Y, M, day=D, + hour=h, minute=m, second=s) + # number of time points + nt = len(time_decimal) + + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # earth and physical parameters for ellipsoid + units = pyTMD.constants(ELLIPSOID) + # tidal love number appropriate for the load tide + hb2 = 0.6207 + + # flatten heights + h = np.array(h).flatten() if np.ndim(h) else h + # convert from geodetic latitude to geocentric latitude + # calculate X, Y and Z from geodetic latitude and longitude + X,Y,Z = pyTMD.spatial.to_cartesian(lon.flatten(), lat.flatten(), h=h, + a_axis=units.a_axis, flat=units.flat) + rr = np.sqrt(X**2.0 + Y**2.0 + Z**2.0) + # calculate geocentric latitude and convert to degrees + latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))/dtr + # geocentric colatitude and longitude in radians + theta = dtr*(90.0 - latitude_geocentric) + phi = dtr*lon.flatten() + + # compute normal gravity at spatial location and elevation of points. + # Normal gravity at height h. p. 82, Eqn.(2-215) + gamma_h = units.gamma_h(theta, h) + dfactor = -hb2*atr*(units.omega**2*rr**2)/(2.0*gamma_h) + + # calculate angular coordinates of mean/secular pole at time + mpx, mpy, fl = pyTMD.eop.iers_mean_pole(time_decimal, convention=CONVENTION) + # read and interpolate IERS daily polar motion values + px, py = pyTMD.eop.iers_polar_motion(MJD, k=3, s=0) + # calculate differentials from mean/secular pole positions + mx = px - mpx + my = -(py - mpy) + + # calculate radial displacement at time + if (TYPE == 'grid'): + ny,nx = np.shape(x) + Srad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) + Srad.mask = np.zeros((ny,nx,nt),dtype=bool) + for i in range(nt): + SRAD = dfactor*np.sin(2.0*theta)*(mx[i]*np.cos(phi)+my[i]*np.sin(phi)) + # reform grid + Srad.data[:,:,i] = np.reshape(SRAD, (ny,nx)) + Srad.mask[:,:,i] = np.isnan(Srad.data[:,:,i]) + elif (TYPE == 'drift'): + Srad = np.ma.zeros((nt), fill_value=FILL_VALUE) + Srad.data[:] = dfactor*np.sin(2.0*theta)*(mx*np.cos(phi)+my*np.sin(phi)) + Srad.mask = np.isnan(Srad.data) + elif (TYPE == 'time series'): + nstation = len(x) + Srad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) + Srad.mask = np.zeros((nstation,nt),dtype=bool) + for s in range(nstation): + SRAD = dfactor[s]*np.sin(2.0*theta[s])*(mx*np.cos(phi[s])+my*np.sin(phi[s])) + Srad.data[s,:] = np.copy(SRAD) + Srad.mask[s,:] = np.isnan(Srad.data[s,:]) + # replace invalid data with fill values + Srad.data[Srad.mask] = Srad.fill_value + + # return the solid earth pole tide displacements + return Srad + +# PURPOSE: compute radial load pole tide displacements +# following IERS Convention (2010) guidelines +def OPT_displacements( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + h: float | np.ndarray = 0.0, + EPSG: str | int = 3031, + EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), + TYPE: str or None = 'drift', + TIME: str = 'UTC', + ELLIPSOID: str = 'WGS84', + CONVENTION: str = '2018', + METHOD: str = 'spline', + FILL_VALUE: float = np.nan, + **kwargs + ): + """ + Compute radial ocean pole tide displacements at points and times + following IERS Convention (2010) guidelines + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + h: float or np.ndarray, default 0.0 + height coordinates in meters + EPSG: int, default: 3031 (Polar Stereographic South, WGS84) + Input coordinate system + EPOCH: tuple, default (2000,1,1,0,0,0) + Time period for calculating delta times + TYPE: str or NoneType, default 'drift' + Input data type + + - ``None``: determined from input variable dimensions + - ``'drift'``: drift buoys or satellite/airborne altimetry + - ``'grid'``: spatial grids or images + - ``'time series'``: time series at a single point + TIME: str, default 'UTC' + Time type if need to compute leap seconds to convert to UTC + + - ``'GPS'``: leap seconds needed + - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) + - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) + - ``'UTC'``: no leap seconds needed + - ``'datetime'``: numpy datatime array in UTC + ELLIPSOID: str, default 'WGS84' + Ellipsoid for calculating Earth parameters + CONVENTION: str, default '2018' + IERS Mean or Secular Pole Convention + + - ``'2003'`` + - ``'2010'`` + - ``'2015'`` + - ``'2018'`` + METHOD: str + Interpolation method + + - ```bilinear```: quick bilinear interpolation + - ```spline```: scipy bivariate spline interpolation + - ```linear```, ```nearest```: scipy regular grid interpolations + FILL_VALUE: float, default np.nan + Output invalid value + + Returns + ------- + Urad: np.ndarray + ocean pole tide at coordinates and time in meters + """ + + # validate input arguments + assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') + assert ELLIPSOID in pyTMD._ellipsoids + assert CONVENTION in pyTMD.eop._conventions + assert METHOD in ('bilinear', 'spline', 'linear', 'nearest') + # determine input data type based on variable dimensions + if not TYPE: + TYPE = pyTMD.spatial.data_type(x, y, delta_time) + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon,lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + timescale = pyTMD.time.timescale().from_datetime( + delta_time.flatten()) + else: + timescale = pyTMD.time.timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + + # convert dynamic time to Modified Julian Days (MJD) + MJD = timescale.tt - 2400000.5 + # convert Julian days to calendar dates + Y,M,D,h,m,s = pyTMD.time.convert_julian(timescale.tt, format='tuple') + # calculate time in year-decimal format + time_decimal = pyTMD.time.convert_calendar_decimal(Y, M, day=D, + hour=h, minute=m, second=s) + # number of time points + nt = len(time_decimal) + + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # earth and physical parameters for ellipsoid + units = pyTMD.constants(ELLIPSOID) + # mean equatorial gravitational acceleration [m/s^2] + ge = 9.7803278 + # density of sea water [kg/m^3] + rho_w = 1025.0 + # tidal love number differential (1 + kl - hl) for pole tide frequencies + gamma = 0.6870 + 0.0036j + + # flatten heights + h = np.array(h).flatten() if np.ndim(h) else h + # convert from geodetic latitude to geocentric latitude + # calculate X, Y and Z from geodetic latitude and longitude + X,Y,Z = pyTMD.spatial.to_cartesian(lon.flatten(), lat.flatten(), h=h, + a_axis=units.a_axis, flat=units.flat) + rr = np.sqrt(X**2.0 + Y**2.0 + Z**2.0) + # calculate geocentric latitude and convert to degrees + latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))/dtr + + # pole tide displacement scale factor + Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM + K = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis/(3.0*ge) + K1 = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis**3/(3.0*units.GM) + + # calculate angular coordinates of mean/secular pole at time + mpx, mpy, fl = pyTMD.eop.iers_mean_pole(time_decimal, convention=CONVENTION) + # read and interpolate IERS daily polar motion values + px, py = pyTMD.eop.iers_polar_motion(MJD, k=3, s=0) + # calculate differentials from mean/secular pole positions + mx = px - mpx + my = -(py - mpy) + + # read ocean pole tide map from Desai (2002) + iur, iun, iue, ilon, ilat = pyTMD.io.ocean_pole_tide() + # interpolate ocean pole tide map from Desai (2002) + if (METHOD == 'spline'): + # use scipy bivariate splines to interpolate to output points + f1 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], + iur[:,::-1].real, kx=1, ky=1) + f2 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], + iur[:,::-1].imag, kx=1, ky=1) + UR = np.zeros((len(latitude_geocentric)), dtype=np.longcomplex) + UR.real = f1.ev(lon.flatten(), latitude_geocentric) + UR.imag = f2.ev(lon.flatten(), latitude_geocentric) + else: + # use scipy regular grid to interpolate values for a given method + r1 = scipy.interpolate.RegularGridInterpolator((ilon,ilat[::-1]), + iur[:,::-1], method=METHOD) + UR = r1.__call__(np.c_[lon.flatten(), latitude_geocentric]) + + # calculate radial displacement at time + if (TYPE == 'grid'): + ny,nx = np.shape(x) + Urad = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) + Urad.mask = np.zeros((ny,nx,nt),dtype=bool) + for i in range(nt): + URAD = K*atr*np.real((mx[i]*gamma.real + my[i]*gamma.imag)*UR.real + + (my[i]*gamma.real - mx[i]*gamma.imag)*UR.imag) + # reform grid + Urad.data[:,:,i] = np.reshape(URAD, (ny,nx)) + Urad.mask[:,:,i] = np.isnan(Urad.data[:,:,i]) + elif (TYPE == 'drift'): + Urad = np.ma.zeros((nt),fill_value=FILL_VALUE) + Urad.data[:] = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real + + (my*gamma.real - mx*gamma.imag)*UR.imag) + Urad.mask = np.isnan(Urad.data) + elif (TYPE == 'time series'): + nstation = len(x) + Urad = np.ma.zeros((nstation,nt),fill_value=FILL_VALUE) + Urad.mask = np.zeros((nstation,nt),dtype=bool) + for s in range(nstation): + URAD = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real[s] + + (my*gamma.real - mx*gamma.imag)*UR.imag[s]) + Urad.data[s,:] = np.copy(URAD) + Urad.mask[s,:] = np.isnan(Urad.data[s,:]) + # replace invalid data with fill values + Urad.data[Urad.mask] = Urad.fill_value + + # return the ocean pole tide displacements + return Urad + +# PURPOSE: compute solid earth tidal elevations +def SET_displacements( + x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, + h: float | np.ndarray = 0.0, + EPSG: str | int = 3031, + EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), + TYPE: str or None = 'drift', + TIME: str = 'UTC', + ELLIPSOID: str = 'WGS84', + TIDE_SYSTEM='tide_free', + EPHEMERIDES='approximate', + **kwargs + ): + """ + Compute solid earth tidal elevations at points and times + following IERS Convention (2010) guidelines + + Parameters + ---------- + x: np.ndarray + x-coordinates in projection EPSG + y: np.ndarray + y-coordinates in projection EPSG + delta_time: np.ndarray + seconds since EPOCH or datetime array + h: float or np.ndarray, default 0.0 + height coordinates in meters + EPSG: int, default: 3031 (Polar Stereographic South, WGS84) + Input coordinate system + EPOCH: tuple, default (2000,1,1,0,0,0) + Time period for calculating delta times + TYPE: str or NoneType, default 'drift' + Input data type + + - ``None``: determined from input variable dimensions + - ``'drift'``: drift buoys or satellite/airborne altimetry + - ``'grid'``: spatial grids or images + - ``'time series'``: time series at a single point + TIME: str, default 'UTC' + Time type if need to compute leap seconds to convert to UTC + + - ``'GPS'``: leap seconds needed + - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) + - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) + - ``'UTC'``: no leap seconds needed + - ``'datetime'``: numpy datatime array in UTC + ELLIPSOID: str, default 'WGS84' + Ellipsoid for calculating Earth parameters + 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) + EPHEMERIDES: str, default 'approximate' + Ephemerides for calculating Earth parameters + + - ``'approximate'``: approximate lunisolar parameters + - ``'JPL'``: computed from JPL ephmerides kernel + + Returns + ------- + tide_se: np.ndarray + solid earth tide at coordinates and time in meters + """ + + # validate input arguments + assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') + assert TIDE_SYSTEM in ('mean_tide', 'tide_free') + assert EPHEMERIDES in ('approximate', 'JPL') + # determine input data type based on variable dimensions + if not TYPE: + TYPE = pyTMD.spatial.data_type(x, y, delta_time) + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon, lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + timescale = pyTMD.time.timescale().from_datetime( + delta_time.flatten()) + else: + timescale = pyTMD.time.timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + # convert tide times to dynamical time + tide_time = timescale.tide + timescale.tt_ut1 + # number of time points + nt = len(timescale) + + # earth and physical parameters for ellipsoid + units = pyTMD.constants(ELLIPSOID) + + # convert input coordinates to cartesian + X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, h=h, + a_axis=units.a_axis, flat=units.flat) + # compute ephemerides for lunisolar coordinates + SX, SY, SZ = pyTMD.astro.solar_ecef(timescale.MJD, ephemerides=EPHEMERIDES) + LX, LY, LZ = pyTMD.astro.lunar_ecef(timescale.MJD, ephemerides=EPHEMERIDES) + + # calculate radial displacement at time + if (TYPE == 'grid'): + ny,nx = np.shape(x) + tide_se = np.zeros((ny,nx,nt)) + # convert coordinates to column arrays + XYZ = np.c_[X, Y, Z] + for i in range(nt): + # reshape time to match spatial + t = tide_time[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=TIDE_SYSTEM) + # calculate radial component of solid earth tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # remove effects of original topography + # (if added when computing cartesian coordinates) + tide_se[:,:,i] = np.reshape(drad - h, (ny,nx)) + elif (TYPE == 'drift'): + # convert coordinates to column arrays + XYZ = np.c_[X, Y, Z] + SXYZ = np.c_[SX, SY, SZ] + LXYZ = np.c_[LX, LY, LZ] + # predict solid earth tides (cartesian) + dxi = pyTMD.predict.solid_earth_tide(tide_time, + XYZ, SXYZ, LXYZ, a_axis=units.a_axis, + tide_system=TIDE_SYSTEM) + # calculate radial component of solid earth tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # remove effects of original topography + # (if added when computing cartesian coordinates) + tide_se = drad - h + elif (TYPE == 'time series'): + nstation = len(x) + tide_se = np.zeros((nstation,nt)) + # convert coordinates to column arrays + SXYZ = np.c_[SX, SY, SZ] + LXYZ = np.c_[LX, LY, LZ] + for s in range(nstation): + # convert coordinates to column arrays + XYZ = np.repeat(np.c_[X[s], Y[s], Z[s]], nt, axis=0) + # predict solid earth tides (cartesian) + dxi = pyTMD.predict.solid_earth_tide(tide_time, + XYZ, SXYZ, LXYZ, a_axis=units.a_axis, + tide_system=TIDE_SYSTEM) + # calculate radial component of solid earth tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # remove effects of original topography + # (if added when computing cartesian coordinates) + tide_se[s,:] = drad - h + + # return the solid earth tide displacements + return tide_se diff --git a/pyTMD/compute_tide_corrections.py b/pyTMD/compute_tide_corrections.py index a834abe8..8170cb06 100644 --- a/pyTMD/compute_tide_corrections.py +++ b/pyTMD/compute_tide_corrections.py @@ -60,6 +60,8 @@ UPDATE HISTORY: Updated 01/2024: made the inferrence of minor constituents an option + refactored lunisolar ephemerides functions + deprecated in favor of refactored compute.py module Updated 12/2023: use new crs class for coordinate reprojection Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 05/2023: use timescale class for time conversion operations @@ -103,25 +105,9 @@ """ from __future__ import print_function, annotations -import logging -import pathlib +import warnings import numpy as np -import scipy.interpolate -import pyTMD.constants -import pyTMD.crs -import pyTMD.eop -import pyTMD.io -import pyTMD.time -import pyTMD.io.model -import pyTMD.predict -import pyTMD.spatial -import pyTMD.utilities - -# attempt imports -try: - import pyproj -except (AttributeError, ImportError, ModuleNotFoundError) as exc: - logging.critical("pyproj not available") +import pyTMD.compute # PURPOSE: wrapper function for computing corrections def compute_corrections( @@ -157,900 +143,32 @@ def compute_corrections( correction: np.ndarray tidal correction at coordinates and time in meters """ + warnings.warn("Deprecated. Please use pyTMD.compute instead", + DeprecationWarning) if CORRECTION.lower() in ('ocean', 'load'): - return compute_tide_corrections(x, y, delta_time, **kwargs) + return pyTMD.compute.tide_elevations(x, y, delta_time, **kwargs) elif (CORRECTION.upper() == 'LPET'): - return compute_LPET_corrections(x, y, delta_time, **kwargs) + return pyTMD.compute.LPET_elevations(x, y, delta_time, **kwargs) elif (CORRECTION.upper() == 'LPT'): - return compute_LPT_corrections(x, y, delta_time, **kwargs) + return pyTMD.compute.LPT_displacements(x, y, delta_time, **kwargs) elif (CORRECTION.upper() == 'OPT'): - return compute_OPT_corrections(x, y, delta_time, **kwargs) + return pyTMD.compute.OPT_displacements(x, y, delta_time, **kwargs) elif (CORRECTION.upper() == 'SET'): - return compute_SET_corrections(x, y, delta_time, **kwargs) + return pyTMD.compute.SET_displacements(x, y, delta_time, **kwargs) else: raise ValueError(f'Unrecognized correction type: {CORRECTION}') -# PURPOSE: compute tides at points and times using tide model algorithms -def compute_tide_corrections( - x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - DIRECTORY: str | pathlib.Path | None = None, - MODEL: str | None = None, - ATLAS_FORMAT: str = 'netcdf', - GZIP: bool = False, - DEFINITION_FILE: str | pathlib.Path | None = None, - EPSG: str | int = 3031, - EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), - TYPE: str or None = 'drift', - TIME: str = 'UTC', - METHOD: str = 'spline', - EXTRAPOLATE: bool = False, - CUTOFF: int | float=10.0, - INFER_MINOR: bool = True, - APPLY_FLEXURE: bool = False, - FILL_VALUE: float = np.nan, - **kwargs - ): - """ - Compute ocean or load tides at points and times from - model constituents - - Parameters - ---------- - x: np.ndarray - x-coordinates in projection EPSG - y: np.ndarray - y-coordinates in projection EPSG - delta_time: np.ndarray - seconds since EPOCH or datetime array - DIRECTORY: str or NoneType, default None - working data directory for tide models - MODEL: str or NoneType, default None - Tide model to use in correction - ATLAS_FORMAT: str, default 'netcdf' - ATLAS tide model format - - - ``'OTIS'`` - - ``'netcdf'`` - GZIP: bool, default False - Tide model files are gzip compressed - DEFINITION_FILE: str or NoneType, default None - Tide model definition file for use - EPSG: int, default: 3031 (Polar Stereographic South, WGS84) - Input coordinate system - EPOCH: tuple, default (2000,1,1,0,0,0) - Time period for calculating delta times - TYPE: str or NoneType, default 'drift' - Input data type - - - ``None``: determined from input variable dimensions - - ``'drift'``: drift buoys or satellite/airborne altimetry - - ``'grid'``: spatial grids or images - - ``'time series'``: time series at a single point - TIME: str, default 'UTC' - Time type if need to compute leap seconds to convert to UTC - - - ``'GPS'``: leap seconds needed - - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) - - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) - - ``'UTC'``: no leap seconds needed - - ``'datetime'``: numpy datatime array in UTC - METHOD: str - Interpolation method - - - ```bilinear```: quick bilinear interpolation - - ```spline```: scipy bivariate spline interpolation - - ```linear```, ```nearest```: scipy regular grid interpolations - - EXTRAPOLATE: bool, default False - Extrapolate with nearest-neighbors - CUTOFF: int or float, default 10.0 - Extrapolation cutoff in kilometers - - Set to ``np.inf`` to extrapolate for all points - INFER_MINOR: bool, default True - Infer the height values for minor tidal constituents - APPLY_FLEXURE: bool, default False - Apply ice flexure scaling factor to height values - - Only valid for models containing flexure fields - FILL_VALUE: float, default np.nan - Output invalid value - - Returns - ------- - tide: np.ndarray - tidal elevation at coordinates and time in meters - """ - - # check that tide directory is accessible - if DIRECTORY is not None: - DIRECTORY = pathlib.Path(DIRECTORY).expanduser() - if not DIRECTORY.exists(): - raise FileNotFoundError("Invalid tide directory") - - # validate input arguments - assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') - assert METHOD in ('bilinear', 'spline', 'linear', 'nearest') - - # get parameters for tide model - if DEFINITION_FILE is not None: - model = pyTMD.io.model(DIRECTORY).from_file( - pathlib.Path(DEFINITION_FILE).expanduser()) - else: - model = pyTMD.io.model(DIRECTORY, format=ATLAS_FORMAT, - compressed=GZIP).elevation(MODEL) - - # determine input data type based on variable dimensions - if not TYPE: - TYPE = pyTMD.spatial.data_type(x, y, delta_time) - # reform coordinate dimensions for input grids - # or verify coordinate dimension shapes - if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): - x,y = np.meshgrid(np.copy(x),np.copy(y)) - elif (TYPE.lower() == 'grid'): - x = np.atleast_2d(x) - y = np.atleast_2d(y) - elif TYPE.lower() in ('time series', 'drift'): - x = np.atleast_1d(x) - y = np.atleast_1d(y) - - # converting x,y from EPSG to latitude/longitude - crs1 = pyTMD.crs().from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - lon, lat = transformer.transform(x.flatten(), y.flatten()) - - # assert delta time is an array - delta_time = np.atleast_1d(delta_time) - # convert delta times or datetimes objects to timescale - if (TIME.lower() == 'datetime'): - timescale = pyTMD.time.timescale().from_datetime( - delta_time.flatten()) - else: - timescale = pyTMD.time.timescale().from_deltatime(delta_time, - epoch=EPOCH, standard=TIME) - # number of time points - nt = len(timescale) - - # read tidal constants and interpolate to grid points - if model.format in ('OTIS', 'ATLAS', 'TMD3'): - amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file, - model.model_file, model.projection, type=model.type, - method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, - grid=model.format, apply_flexure=APPLY_FLEXURE) - # use delta time at 2000.0 to match TMD outputs - deltat = np.zeros((nt), dtype=np.float64) - elif (model.format == 'netcdf'): - amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file, - model.model_file, type=model.type, method=METHOD, - extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, - compressed=model.compressed) - # use delta time at 2000.0 to match TMD outputs - deltat = np.zeros((nt), dtype=np.float64) - elif (model.format == 'GOT'): - amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file, - method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF, - scale=model.scale, compressed=model.compressed) - # delta time (TT - UT1) - deltat = timescale.tt_ut1 - elif (model.format == 'FES'): - amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file, - type=model.type, version=model.version, method=METHOD, - extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale, - compressed=model.compressed) - # available model constituents - c = model.constituents - # delta time (TT - UT1) - deltat = timescale.tt_ut1 - - # calculate complex phase in radians for Euler's - cph = -1j*ph*np.pi/180.0 - # calculate constituent oscillation - hc = amp*np.exp(cph) - - # predict tidal elevations at time and infer minor corrections - if (TYPE.lower() == 'grid'): - ny,nx = np.shape(x) - tide = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) - tide.mask = np.zeros((ny,nx,nt),dtype=bool) - for i in range(nt): - TIDE = pyTMD.predict.map(timescale.tide[i], hc, c, - deltat=deltat[i], corrections=model.format) - # calculate values for minor constituents by inferrence - if INFER_MINOR: - MINOR = pyTMD.predict.infer_minor(timescale.tide[i], hc, c, - deltat=deltat[i], corrections=model.format) - else: - MINOR = np.ma.zeros_like(TIDE) - # add major and minor components and reform grid - tide[:,:,i] = np.reshape((TIDE+MINOR), (ny,nx)) - tide.mask[:,:,i] = np.reshape((TIDE.mask | MINOR.mask), (ny,nx)) - elif (TYPE.lower() == 'drift'): - tide = np.ma.zeros((nt), fill_value=FILL_VALUE) - tide.mask = np.any(hc.mask,axis=1) - tide.data[:] = pyTMD.predict.drift(timescale.tide, hc, c, - deltat=deltat, corrections=model.format) - # calculate values for minor constituents by inferrence - if INFER_MINOR: - minor = pyTMD.predict.infer_minor(timescale.tide, hc, c, - deltat=deltat, corrections=model.format) - tide.data[:] += minor.data[:] - elif (TYPE.lower() == 'time series'): - nstation = len(x) - tide = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) - tide.mask = np.zeros((nstation,nt),dtype=bool) - for s in range(nstation): - HC = hc[s,None,:] - TIDE = pyTMD.predict.time_series(timescale.tide, HC, c, - deltat=deltat, corrections=model.format) - # calculate values for minor constituents by inferrence - if INFER_MINOR: - MINOR = pyTMD.predict.infer_minor(timescale.tide, HC, c, - deltat=deltat, corrections=model.format) - else: - MINOR = np.ma.zeros_like(TIDE) - # add major and minor components - tide.data[s,:] = TIDE.data[:] + MINOR.data[:] - tide.mask[s,:] = (TIDE.mask | MINOR.mask) - # replace invalid values with fill value - tide.data[tide.mask] = tide.fill_value - - # return the ocean or load tide correction - return tide - -# PURPOSE: compute long-period equilibrium tidal elevations -def compute_LPET_corrections( - x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - EPSG: str | int = 3031, - EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), - TYPE: str or None = 'drift', - TIME: str = 'UTC', - **kwargs - ): - """ - Compute long-period equilibrium tidal elevations at points and times - - Parameters - ---------- - x: np.ndarray - x-coordinates in projection EPSG - y: np.ndarray - y-coordinates in projection EPSG - delta_time: np.ndarray - seconds since EPOCH or datetime array - EPSG: int, default: 3031 (Polar Stereographic South, WGS84) - Input coordinate system - EPOCH: tuple, default (2000,1,1,0,0,0) - Time period for calculating delta times - TYPE: str or NoneType, default 'drift' - Input data type - - - ``None``: determined from input variable dimensions - - ``'drift'``: drift buoys or satellite/airborne altimetry - - ``'grid'``: spatial grids or images - - ``'time series'``: time series at a single point - TIME: str, default 'UTC' - Time type if need to compute leap seconds to convert to UTC - - - ``'GPS'``: leap seconds needed - - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) - - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) - - ``'UTC'``: no leap seconds needed - - ``'datetime'``: numpy datatime array in UTC - FILL_VALUE: float, default np.nan - Output invalid value - - Returns - ------- - tide_lpe: np.ndarray - long-period equilibrium tide at coordinates and time in meters - """ - - # validate input arguments - assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') - # determine input data type based on variable dimensions - if not TYPE: - TYPE = pyTMD.spatial.data_type(x, y, delta_time) - # reform coordinate dimensions for input grids - # or verify coordinate dimension shapes - if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): - x,y = np.meshgrid(np.copy(x),np.copy(y)) - elif (TYPE.lower() == 'grid'): - x = np.atleast_2d(x) - y = np.atleast_2d(y) - elif TYPE.lower() in ('time series', 'drift'): - x = np.atleast_1d(x) - y = np.atleast_1d(y) - - # converting x,y from EPSG to latitude/longitude - crs1 = pyTMD.crs().from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - lon, lat = transformer.transform(x.flatten(), y.flatten()) - - # assert delta time is an array - delta_time = np.atleast_1d(delta_time) - # convert delta times or datetimes objects to timescale - if (TIME.lower() == 'datetime'): - timescale = pyTMD.time.timescale().from_datetime( - delta_time.flatten()) - else: - timescale = pyTMD.time.timescale().from_deltatime(delta_time, - epoch=EPOCH, standard=TIME) - # number of time points - nt = len(timescale) - # convert tide times to dynamic time - tide_time = timescale.tide + timescale.tt_ut1 - - # predict long-period equilibrium tides at time - if (TYPE == 'grid'): - ny,nx = np.shape(x) - tide_lpe = np.zeros((ny,nx,nt)) - for i in range(nt): - lpet = pyTMD.predict.equilibrium_tide(tide_time[i], lat) - tide_lpe[:,:,i] = np.reshape(lpet,(ny,nx)) - elif (TYPE == 'drift'): - tide_lpe = pyTMD.predict.equilibrium_tide(tide_time, lat) - elif (TYPE == 'time series'): - nstation = len(x) - tide_lpe = np.zeros((nstation,nt)) - for s in range(nstation): - lpet = pyTMD.predict.equilibrium_tide(tide_time, lat[s]) - tide_lpe[s,:] = np.copy(lpet) - - # return the long-period equilibrium tide corrections - return tide_lpe - -# PURPOSE: compute radial load pole tide displacements -# following IERS Convention (2010) guidelines -def compute_LPT_corrections( - x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - h: float | np.ndarray = 0.0, - EPSG: str | int = 3031, - EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), - TYPE: str or None = 'drift', - TIME: str = 'UTC', - ELLIPSOID: str = 'WGS84', - CONVENTION: str = '2018', - FILL_VALUE: float = np.nan, - **kwargs - ): - """ - Compute radial load pole tide displacements at points and times - following IERS Convention (2010) guidelines - - Parameters - ---------- - x: np.ndarray - x-coordinates in projection EPSG - y: np.ndarray - y-coordinates in projection EPSG - delta_time: np.ndarray - seconds since EPOCH or datetime array - h: float or np.ndarray, default 0.0 - height coordinates in meters - EPSG: int, default: 3031 (Polar Stereographic South, WGS84) - Input coordinate system - EPOCH: tuple, default (2000,1,1,0,0,0) - Time period for calculating delta times - TYPE: str or NoneType, default 'drift' - Input data type - - - ``None``: determined from input variable dimensions - - ``'drift'``: drift buoys or satellite/airborne altimetry - - ``'grid'``: spatial grids or images - - ``'time series'``: time series at a single point - TIME: str, default 'UTC' - Time type if need to compute leap seconds to convert to UTC - - - ``'GPS'``: leap seconds needed - - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) - - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) - - ``'UTC'``: no leap seconds needed - - ``'datetime'``: numpy datatime array in UTC - ELLIPSOID: str, default 'WGS84' - Ellipsoid for calculating Earth parameters - CONVENTION: str, default '2018' - IERS Mean or Secular Pole Convention - - - ``'2003'`` - - ``'2010'`` - - ``'2015'`` - - ``'2018'`` - FILL_VALUE: float, default np.nan - Output invalid value - - Returns - ------- - Srad: np.ndarray - solid earth pole tide at coordinates and time in meters - """ - - # validate input arguments - assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') - assert ELLIPSOID in pyTMD._ellipsoids - assert CONVENTION in pyTMD.eop._conventions - # determine input data type based on variable dimensions - if not TYPE: - TYPE = pyTMD.spatial.data_type(x, y, delta_time) - # reform coordinate dimensions for input grids - # or verify coordinate dimension shapes - if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): - x,y = np.meshgrid(np.copy(x),np.copy(y)) - elif (TYPE.lower() == 'grid'): - x = np.atleast_2d(x) - y = np.atleast_2d(y) - elif TYPE.lower() in ('time series', 'drift'): - x = np.atleast_1d(x) - y = np.atleast_1d(y) - - # converting x,y from EPSG to latitude/longitude - crs1 = pyTMD.crs().from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - lon,lat = transformer.transform(x.flatten(), y.flatten()) - - # assert delta time is an array - delta_time = np.atleast_1d(delta_time) - # convert delta times or datetimes objects to timescale - if (TIME.lower() == 'datetime'): - timescale = pyTMD.time.timescale().from_datetime( - delta_time.flatten()) - else: - timescale = pyTMD.time.timescale().from_deltatime(delta_time, - epoch=EPOCH, standard=TIME) - - # convert dynamic time to Modified Julian Days (MJD) - MJD = timescale.tt - 2400000.5 - # convert Julian days to calendar dates - Y,M,D,h,m,s = pyTMD.time.convert_julian(timescale.tt, format='tuple') - # calculate time in year-decimal format - time_decimal = pyTMD.time.convert_calendar_decimal(Y, M, day=D, - hour=h, minute=m, second=s) - # number of time points - nt = len(time_decimal) - - # degrees and arcseconds to radians - dtr = np.pi/180.0 - atr = np.pi/648000.0 - # earth and physical parameters for ellipsoid - units = pyTMD.constants(ELLIPSOID) - # tidal love number appropriate for the load tide - hb2 = 0.6207 - - # flatten heights - h = np.array(h).flatten() if np.ndim(h) else h - # convert from geodetic latitude to geocentric latitude - # calculate X, Y and Z from geodetic latitude and longitude - X,Y,Z = pyTMD.spatial.to_cartesian(lon.flatten(), lat.flatten(), h=h, - a_axis=units.a_axis, flat=units.flat) - rr = np.sqrt(X**2.0 + Y**2.0 + Z**2.0) - # calculate geocentric latitude and convert to degrees - latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))/dtr - # geocentric colatitude and longitude in radians - theta = dtr*(90.0 - latitude_geocentric) - phi = dtr*lon.flatten() - - # compute normal gravity at spatial location and elevation of points. - # Normal gravity at height h. p. 82, Eqn.(2-215) - gamma_h = units.gamma_h(theta, h) - dfactor = -hb2*atr*(units.omega**2*rr**2)/(2.0*gamma_h) - - # calculate angular coordinates of mean/secular pole at time - mpx, mpy, fl = pyTMD.eop.iers_mean_pole(time_decimal, convention=CONVENTION) - # read and interpolate IERS daily polar motion values - px, py = pyTMD.eop.iers_polar_motion(MJD, k=3, s=0) - # calculate differentials from mean/secular pole positions - mx = px - mpx - my = -(py - mpy) - - # calculate radial displacement at time - if (TYPE == 'grid'): - ny,nx = np.shape(x) - Srad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) - Srad.mask = np.zeros((ny,nx,nt),dtype=bool) - for i in range(nt): - SRAD = dfactor*np.sin(2.0*theta)*(mx[i]*np.cos(phi)+my[i]*np.sin(phi)) - # reform grid - Srad.data[:,:,i] = np.reshape(SRAD, (ny,nx)) - Srad.mask[:,:,i] = np.isnan(Srad.data[:,:,i]) - elif (TYPE == 'drift'): - Srad = np.ma.zeros((nt), fill_value=FILL_VALUE) - Srad.data[:] = dfactor*np.sin(2.0*theta)*(mx*np.cos(phi)+my*np.sin(phi)) - Srad.mask = np.isnan(Srad.data) - elif (TYPE == 'time series'): - nstation = len(x) - Srad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) - Srad.mask = np.zeros((nstation,nt),dtype=bool) - for s in range(nstation): - SRAD = dfactor[s]*np.sin(2.0*theta[s])*(mx*np.cos(phi[s])+my*np.sin(phi[s])) - Srad.data[s,:] = np.copy(SRAD) - Srad.mask[s,:] = np.isnan(Srad.data[s,:]) - # replace invalid data with fill values - Srad.data[Srad.mask] = Srad.fill_value - - # return the solid earth pole tide corrections - return Srad - -# PURPOSE: compute radial load pole tide displacements -# following IERS Convention (2010) guidelines -def compute_OPT_corrections( - x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - h: float | np.ndarray = 0.0, - EPSG: str | int = 3031, - EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), - TYPE: str or None = 'drift', - TIME: str = 'UTC', - ELLIPSOID: str = 'WGS84', - CONVENTION: str = '2018', - METHOD: str = 'spline', - FILL_VALUE: float = np.nan, - **kwargs - ): - """ - Compute radial ocean pole tide displacements at points and times - following IERS Convention (2010) guidelines - - Parameters - ---------- - x: np.ndarray - x-coordinates in projection EPSG - y: np.ndarray - y-coordinates in projection EPSG - delta_time: np.ndarray - seconds since EPOCH or datetime array - h: float or np.ndarray, default 0.0 - height coordinates in meters - EPSG: int, default: 3031 (Polar Stereographic South, WGS84) - Input coordinate system - EPOCH: tuple, default (2000,1,1,0,0,0) - Time period for calculating delta times - TYPE: str or NoneType, default 'drift' - Input data type - - - ``None``: determined from input variable dimensions - - ``'drift'``: drift buoys or satellite/airborne altimetry - - ``'grid'``: spatial grids or images - - ``'time series'``: time series at a single point - TIME: str, default 'UTC' - Time type if need to compute leap seconds to convert to UTC - - - ``'GPS'``: leap seconds needed - - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) - - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) - - ``'UTC'``: no leap seconds needed - - ``'datetime'``: numpy datatime array in UTC - ELLIPSOID: str, default 'WGS84' - Ellipsoid for calculating Earth parameters - CONVENTION: str, default '2018' - IERS Mean or Secular Pole Convention - - - ``'2003'`` - - ``'2010'`` - - ``'2015'`` - - ``'2018'`` - METHOD: str - Interpolation method - - - ```bilinear```: quick bilinear interpolation - - ```spline```: scipy bivariate spline interpolation - - ```linear```, ```nearest```: scipy regular grid interpolations - FILL_VALUE: float, default np.nan - Output invalid value - - Returns - ------- - Urad: np.ndarray - ocean pole tide at coordinates and time in meters - """ - - # validate input arguments - assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') - assert ELLIPSOID in pyTMD._ellipsoids - assert CONVENTION in pyTMD.eop._conventions - assert METHOD in ('bilinear', 'spline', 'linear', 'nearest') - # determine input data type based on variable dimensions - if not TYPE: - TYPE = pyTMD.spatial.data_type(x, y, delta_time) - # reform coordinate dimensions for input grids - # or verify coordinate dimension shapes - if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): - x,y = np.meshgrid(np.copy(x),np.copy(y)) - elif (TYPE.lower() == 'grid'): - x = np.atleast_2d(x) - y = np.atleast_2d(y) - elif TYPE.lower() in ('time series', 'drift'): - x = np.atleast_1d(x) - y = np.atleast_1d(y) - - # converting x,y from EPSG to latitude/longitude - crs1 = pyTMD.crs().from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - lon,lat = transformer.transform(x.flatten(), y.flatten()) - - # assert delta time is an array - delta_time = np.atleast_1d(delta_time) - # convert delta times or datetimes objects to timescale - if (TIME.lower() == 'datetime'): - timescale = pyTMD.time.timescale().from_datetime( - delta_time.flatten()) - else: - timescale = pyTMD.time.timescale().from_deltatime(delta_time, - epoch=EPOCH, standard=TIME) - - # convert dynamic time to Modified Julian Days (MJD) - MJD = timescale.tt - 2400000.5 - # convert Julian days to calendar dates - Y,M,D,h,m,s = pyTMD.time.convert_julian(timescale.tt, format='tuple') - # calculate time in year-decimal format - time_decimal = pyTMD.time.convert_calendar_decimal(Y, M, day=D, - hour=h, minute=m, second=s) - # number of time points - nt = len(time_decimal) - - # degrees and arcseconds to radians - dtr = np.pi/180.0 - atr = np.pi/648000.0 - # earth and physical parameters for ellipsoid - units = pyTMD.constants(ELLIPSOID) - # mean equatorial gravitational acceleration [m/s^2] - ge = 9.7803278 - # density of sea water [kg/m^3] - rho_w = 1025.0 - # tidal love number differential (1 + kl - hl) for pole tide frequencies - gamma = 0.6870 + 0.0036j - - # flatten heights - h = np.array(h).flatten() if np.ndim(h) else h - # convert from geodetic latitude to geocentric latitude - # calculate X, Y and Z from geodetic latitude and longitude - X,Y,Z = pyTMD.spatial.to_cartesian(lon.flatten(), lat.flatten(), h=h, - a_axis=units.a_axis, flat=units.flat) - rr = np.sqrt(X**2.0 + Y**2.0 + Z**2.0) - # calculate geocentric latitude and convert to degrees - latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))/dtr - - # pole tide displacement scale factor - Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM - K = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis/(3.0*ge) - K1 = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis**3/(3.0*units.GM) - - # calculate angular coordinates of mean/secular pole at time - mpx, mpy, fl = pyTMD.eop.iers_mean_pole(time_decimal, convention=CONVENTION) - # read and interpolate IERS daily polar motion values - px, py = pyTMD.eop.iers_polar_motion(MJD, k=3, s=0) - # calculate differentials from mean/secular pole positions - mx = px - mpx - my = -(py - mpy) - - # read ocean pole tide map from Desai (2002) - iur, iun, iue, ilon, ilat = pyTMD.io.ocean_pole_tide() - # interpolate ocean pole tide map from Desai (2002) - if (METHOD == 'spline'): - # use scipy bivariate splines to interpolate to output points - f1 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], - iur[:,::-1].real, kx=1, ky=1) - f2 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], - iur[:,::-1].imag, kx=1, ky=1) - UR = np.zeros((len(latitude_geocentric)), dtype=np.longcomplex) - UR.real = f1.ev(lon.flatten(), latitude_geocentric) - UR.imag = f2.ev(lon.flatten(), latitude_geocentric) - else: - # use scipy regular grid to interpolate values for a given method - r1 = scipy.interpolate.RegularGridInterpolator((ilon,ilat[::-1]), - iur[:,::-1], method=METHOD) - UR = r1.__call__(np.c_[lon.flatten(), latitude_geocentric]) - - # calculate radial displacement at time - if (TYPE == 'grid'): - ny,nx = np.shape(x) - Urad = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) - Urad.mask = np.zeros((ny,nx,nt),dtype=bool) - for i in range(nt): - URAD = K*atr*np.real((mx[i]*gamma.real + my[i]*gamma.imag)*UR.real + - (my[i]*gamma.real - mx[i]*gamma.imag)*UR.imag) - # reform grid - Urad.data[:,:,i] = np.reshape(URAD, (ny,nx)) - Urad.mask[:,:,i] = np.isnan(Urad.data[:,:,i]) - elif (TYPE == 'drift'): - Urad = np.ma.zeros((nt),fill_value=FILL_VALUE) - Urad.data[:] = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real + - (my*gamma.real - mx*gamma.imag)*UR.imag) - Urad.mask = np.isnan(Urad.data) - elif (TYPE == 'time series'): - nstation = len(x) - Urad = np.ma.zeros((nstation,nt),fill_value=FILL_VALUE) - Urad.mask = np.zeros((nstation,nt),dtype=bool) - for s in range(nstation): - URAD = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real[s] + - (my*gamma.real - mx*gamma.imag)*UR.imag[s]) - Urad.data[s,:] = np.copy(URAD) - Urad.mask[s,:] = np.isnan(Urad.data[s,:]) - # replace invalid data with fill values - Urad.data[Urad.mask] = Urad.fill_value - - # return the ocean pole tide corrections - return Urad - -# PURPOSE: compute solid earth tidal elevations -def compute_SET_corrections( - x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - h: float | np.ndarray = 0.0, - EPSG: str | int = 3031, - EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), - TYPE: str or None = 'drift', - TIME: str = 'UTC', - ELLIPSOID: str = 'WGS84', - TIDE_SYSTEM='tide_free', - EPHEMERIDES='approximate', - **kwargs - ): - """ - Compute solid earth tidal elevations at points and times - following IERS Convention (2010) guidelines - - Parameters - ---------- - x: np.ndarray - x-coordinates in projection EPSG - y: np.ndarray - y-coordinates in projection EPSG - delta_time: np.ndarray - seconds since EPOCH or datetime array - h: float or np.ndarray, default 0.0 - height coordinates in meters - EPSG: int, default: 3031 (Polar Stereographic South, WGS84) - Input coordinate system - EPOCH: tuple, default (2000,1,1,0,0,0) - Time period for calculating delta times - TYPE: str or NoneType, default 'drift' - Input data type - - - ``None``: determined from input variable dimensions - - ``'drift'``: drift buoys or satellite/airborne altimetry - - ``'grid'``: spatial grids or images - - ``'time series'``: time series at a single point - TIME: str, default 'UTC' - Time type if need to compute leap seconds to convert to UTC - - - ``'GPS'``: leap seconds needed - - ``'LORAN'``: leap seconds needed (LORAN = GPS + 9 seconds) - - ``'TAI'``: leap seconds needed (TAI = GPS + 19 seconds) - - ``'UTC'``: no leap seconds needed - - ``'datetime'``: numpy datatime array in UTC - ELLIPSOID: str, default 'WGS84' - Ellipsoid for calculating Earth parameters - 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) - EPHEMERIDES: str, default 'approximate' - Ephemerides for calculating Earth parameters - - - ``'approximate'``: approximate lunisolar parameters - - ``'JPL'``: computed from JPL ephmerides kernel - - Returns - ------- - tide_se: np.ndarray - solid earth tide at coordinates and time in meters - """ - - # validate input arguments - assert TIME in ('GPS', 'LORAN', 'TAI', 'UTC', 'datetime') - assert TIDE_SYSTEM in ('mean_tide', 'tide_free') - assert EPHEMERIDES in ('approximate', 'JPL') - # determine input data type based on variable dimensions - if not TYPE: - TYPE = pyTMD.spatial.data_type(x, y, delta_time) - # reform coordinate dimensions for input grids - # or verify coordinate dimension shapes - if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): - x,y = np.meshgrid(np.copy(x),np.copy(y)) - elif (TYPE.lower() == 'grid'): - x = np.atleast_2d(x) - y = np.atleast_2d(y) - elif TYPE.lower() in ('time series', 'drift'): - x = np.atleast_1d(x) - y = np.atleast_1d(y) - - # converting x,y from EPSG to latitude/longitude - crs1 = pyTMD.crs().from_input(EPSG) - crs2 = pyproj.CRS.from_epsg(4326) - transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - lon, lat = transformer.transform(x.flatten(), y.flatten()) - - # assert delta time is an array - delta_time = np.atleast_1d(delta_time) - # convert delta times or datetimes objects to timescale - if (TIME.lower() == 'datetime'): - timescale = pyTMD.time.timescale().from_datetime( - delta_time.flatten()) - else: - timescale = pyTMD.time.timescale().from_deltatime(delta_time, - epoch=EPOCH, standard=TIME) - # convert tide times to dynamical time - tide_time = timescale.tide + timescale.tt_ut1 - # number of time points - nt = len(timescale) +def compute_tide_corrections(*args, **kwargs): + return pyTMD.compute.tide_elevations(*args, **kwargs) - # earth and physical parameters for ellipsoid - units = pyTMD.constants(ELLIPSOID) +def compute_LPET_corrections(*args, **kwargs): + return pyTMD.compute.LPET_elevations(*args, **kwargs) - # convert input coordinates to cartesian - X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, h=h, - a_axis=units.a_axis, flat=units.flat) - # compute ephemerides for lunisolar coordinates - if (EPHEMERIDES.lower() == 'approximate'): - # get low-resolution solar and lunar ephemerides - SX, SY, SZ = pyTMD.astro.solar_ecef(timescale.MJD) - LX, LY, LZ = pyTMD.astro.lunar_ecef(timescale.MJD) - elif (EPHEMERIDES.upper() == 'JPL'): - # compute solar and lunar ephemerides from JPL kernel - SX, SY, SZ = pyTMD.astro.solar_ephemerides(timescale.MJD) - LX, LY, LZ = pyTMD.astro.lunar_ephemerides(timescale.MJD) +def compute_LPT_corrections(*args, **kwargs): + return pyTMD.compute.LPT_displacements(*args, **kwargs) - # calculate radial displacement at time - if (TYPE == 'grid'): - ny,nx = np.shape(x) - tide_se = np.zeros((ny,nx,nt)) - # convert coordinates to column arrays - XYZ = np.c_[X, Y, Z] - for i in range(nt): - # reshape time to match spatial - t = tide_time[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=TIDE_SYSTEM) - # calculate radial component of solid earth tides - dln,dlt,drad = pyTMD.spatial.to_geodetic( - X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], - a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se[:,:,i] = np.reshape(drad - h, (ny,nx)) - elif (TYPE == 'drift'): - # convert coordinates to column arrays - XYZ = np.c_[X, Y, Z] - SXYZ = np.c_[SX, SY, SZ] - LXYZ = np.c_[LX, LY, LZ] - # predict solid earth tides (cartesian) - dxi = pyTMD.predict.solid_earth_tide(tide_time, - XYZ, SXYZ, LXYZ, a_axis=units.a_axis, - tide_system=TIDE_SYSTEM) - # calculate radial component of solid earth tides - dln,dlt,drad = pyTMD.spatial.to_geodetic( - X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], - a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se = drad - h - elif (TYPE == 'time series'): - nstation = len(x) - tide_se = np.zeros((nstation,nt)) - # convert coordinates to column arrays - SXYZ = np.c_[SX, SY, SZ] - LXYZ = np.c_[LX, LY, LZ] - for s in range(nstation): - # convert coordinates to column arrays - XYZ = np.repeat(np.c_[X[s], Y[s], Z[s]], nt, axis=0) - # predict solid earth tides (cartesian) - dxi = pyTMD.predict.solid_earth_tide(tide_time, - XYZ, SXYZ, LXYZ, a_axis=units.a_axis, - tide_system=TIDE_SYSTEM) - # calculate radial component of solid earth tides - dln,dlt,drad = pyTMD.spatial.to_geodetic( - X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], - a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se[s,:] = drad - h +def compute_OPT_corrections(*args, **kwargs): + return pyTMD.compute.OPT_displacements(*args, **kwargs) - # return the solid earth tide corrections - return tide_se +def compute_SET_corrections(*args, **kwargs): + return pyTMD.compute.SET_displacements(*args, **kwargs) diff --git a/pyTMD/predict.py b/pyTMD/predict.py index 14d8b9cc..46482cf8 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -21,6 +21,7 @@ UPDATE HISTORY: Updated 01/2024: moved minor arguments calculation into new function + added more constituent parameters for OTIS/ATLAS predictions 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 @@ -582,7 +583,7 @@ def solid_earth_tide( # return the solid earth tide return dxt -def _constituent_parameters(c): +def _constituent_parameters(c, **kwargs): """ Loads parameters for a given tidal constituent @@ -590,6 +591,8 @@ def _constituent_parameters(c): ---------- c: list tidal constituent ID + raise_error: bool, default False + Raise exception if constituent is unsupported Returns ------- @@ -613,26 +616,30 @@ def _constituent_parameters(c): .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 """ + # default keyword arguments + kwargs.setdefault('raise_error', False) # constituents array that are included in tidal program cindex = ['m2', 's2', 'k1', 'o1', 'n2', 'p1', 'k2', 'q1', '2n2', 'mu2', 'nu2', 'l2', 't2', 'j1', 'm1', 'oo1', 'rho1', 'mf', 'mm', 'ssa', - 'm4', 'ms4', 'mn4', 'm6', 'm8', 'mk3', 's6', '2sm2', '2mk3'] + 'm4', 'ms4', 'mn4', 'm6', 'm8', 'mk3', 's6', '2sm2', '2mk3', + 'msf', 'sa', 'mt', '2q1'] # species type (spherical harmonic dependence of quadrupole potential) _species = np.array([2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, - 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]) # Load Love numbers # alpha = correction factor for first order load tides _alpha = np.array([0.693, 0.693, 0.736, 0.695, 0.693, 0.706, 0.693, 0.695, 0.693, 0.693, 0.693, 0.693, 0.693, 0.695, 0.695, 0.695, 0.695, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, - 0.693, 0.693]) + 0.693, 0.693, 0.693, 0.693, 0.693, 0.693]) # omega: angular frequency of constituent, in radians _omega = np.array([1.405189e-04, 1.454441e-04, 7.292117e-05, 6.759774e-05, 1.378797e-04, 7.252295e-05, 1.458423e-04, 6.495854e-05, 1.352405e-04, 1.355937e-04, 1.382329e-04, 1.431581e-04, 1.452450e-04, 7.556036e-05, 7.028195e-05, 7.824458e-05, 6.531174e-05, 0.053234e-04, 0.026392e-04, 0.003982e-04, 2.810377e-04, 2.859630e-04, 2.783984e-04, 4.215566e-04, - 5.620755e-04, 2.134402e-04, 4.363323e-04, 1.503693e-04, 2.081166e-04]) + 5.620755e-04, 2.134402e-04, 4.363323e-04, 1.503693e-04, 2.081166e-04, + 4.925200e-06, 1.990970e-07, 7.962619e-06, 6.231934e-05]) # Astronomical arguments (relative to t0 = 1 Jan 0:00 1992) # phases for each constituent are referred to the time when the phase of # the forcing for that constituent is zero on the Greenwich meridian @@ -641,13 +648,15 @@ def _constituent_parameters(c): 3.463115091, 5.427136701, 0.553986502, 0.052841931, 2.137025284, 2.436575100, 1.929046130, 5.254133027, 1.756042456, 1.964021610, 3.487600001, 3.463115091, 1.731557546, 1.499093481, 5.194672637, - 6.926230184, 1.904561220, 0.000000000, 4.551627762, 3.809122439]) + 6.926230184, 1.904561220, 0.000000000, 4.551627762, 3.809122439, + 4.551627762, 6.232786837, 3.720064066, 3.91369596]) # amplitudes of equilibrium tide in meters # _amplitude = np.array([0.242334,0.112743,0.141565,0.100661,0.046397, _amplitude = np.array([0.2441, 0.112743, 0.141565, 0.100661, 0.046397, 0.046848, 0.030684, 0.019273, 0.006141, 0.007408, 0.008811, 0.006931, 0.006608, 0.007915, 0.007915, 0.004338, 0.003661, 0.042041, 0.022191, - 0.019567, 0., 0., 0., 0., 0., 0., 0., 0., 0.]) + 0.019567, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.003681, 0.003104, + 0.008044, 0.002565]) # map between input constituent and cindex j = [j for j,val in enumerate(cindex) if (val == c.lower())] @@ -658,6 +667,8 @@ def _constituent_parameters(c): omega, = _omega[j] alpha, = _alpha[j] species, = _species[j] + elif kwargs['raise_error']: + raise ValueError(f'Unsupported constituent {c}') else: amplitude = 0.0 phase = 0.0 diff --git a/scripts/compute_SET_displacements.py b/scripts/compute_SET_displacements.py index c70fdfdf..7adef650 100644 --- a/scripts/compute_SET_displacements.py +++ b/scripts/compute_SET_displacements.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_SET_displacements.py -Written by Tyler Sutterley (12/2023) +Written by Tyler Sutterley (01/2024) Calculates radial solid earth tide displacements for an input file following IERS Convention (2010) guidelines https://iers-conventions.obspm.fr/chapter7.php @@ -97,6 +97,7 @@ doi: 10.1111/j.1365-246X.1981.tb02690.x UPDATE HISTORY: + Updated 01/2024: refactored lunisolar ephemerides functions Updated 12/2023: use new crs class to get projection information Updated 10/2023: can write datetime as time column for csv files Updated 05/2023: use timescale class for time conversion operations @@ -239,14 +240,8 @@ def compute_SET_displacements(input_file, output_file, X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, h=h, a_axis=units.a_axis, flat=units.flat) # compute ephemerides for lunisolar coordinates - if (EPHEMERIDES.lower() == 'approximate'): - # get low-resolution solar and lunar ephemerides - SX, SY, SZ = pyTMD.astro.solar_ecef(timescale.MJD) - LX, LY, LZ = pyTMD.astro.lunar_ecef(timescale.MJD) - elif (EPHEMERIDES.upper() == 'JPL'): - # compute solar and lunar ephemerides from JPL kernel - SX, SY, SZ = pyTMD.astro.solar_ephemerides(timescale.MJD) - LX, LY, LZ = pyTMD.astro.lunar_ephemerides(timescale.MJD) + SX, SY, SZ = pyTMD.astro.solar_ecef(timescale.MJD, ephemerides=EPHEMERIDES) + LX, LY, LZ = pyTMD.astro.lunar_ecef(timescale.MJD, ephemerides=EPHEMERIDES) # calculate radial displacement at time if (TYPE == 'grid'): diff --git a/test/test_atlas_read.py b/test/test_atlas_read.py index 291ded8e..2098e0c3 100644 --- a/test/test_atlas_read.py +++ b/test/test_atlas_read.py @@ -17,6 +17,7 @@ UPDATE HISTORY: Updated 01/2024: test doodson and cartwright numbers of each constituent + refactored compute functions into compute.py Updated 04/2023: using pathlib to define and expand paths Updated 12/2022: add check for read and interpolate constants Updated 11/2022: use f-strings for formatting verbose or ascii output @@ -389,7 +390,7 @@ def test_Ross_Ice_Shelf(MODEL, METHOD, EXTRAPOLATE): # time dimension delta_time = 0.0 # calculate tide map - tide = pyTMD.compute_tide_corrections(xgrid, ygrid, delta_time, + tide = pyTMD.compute.tide_elevations(xgrid, ygrid, delta_time, DIRECTORY=filepath, MODEL=MODEL, GZIP=True, EPOCH=(2000,1,1,0,0,0), TYPE='grid', TIME='TAI', EPSG=3031, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE) diff --git a/test/test_download_and_read.py b/test/test_download_and_read.py index 041fa6f3..777618b6 100644 --- a/test/test_download_and_read.py +++ b/test/test_download_and_read.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -test_download_and_read.py (04/2023) +test_download_and_read.py (01/2024) Tests that CATS2008 data can be downloaded from the US Antarctic Program (USAP) Tests that AOTIM-5-2018 data can be downloaded from the NSF ArcticData server Tests the read program to verify that constituents are being extracted @@ -19,6 +19,7 @@ https://boto3.amazonaws.com/v1/documentation/api/latest/index.html UPDATE HISTORY: + Updated 01/2024: refactored compute functions into compute.py Updated 04/2023: using pathlib to define and expand paths Updated 12/2022: add check for read and interpolate constants Updated 11/2022: added encoding for writing ascii files @@ -690,7 +691,7 @@ def test_Ross_Ice_Shelf(self, METHOD, EXTRAPOLATE): # time dimension delta_time = np.zeros((24))*3600 # calculate tide drift corrections - tide = pyTMD.compute_tide_corrections(x, y, delta_time, + tide = pyTMD.compute.tide_elevations(x, y, delta_time, DIRECTORY=filepath, MODEL='CATS2008', GZIP=False, EPOCH=pyTMD.time._j2000_epoch, TYPE='drift', TIME='UTC', EPSG=3031, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE) @@ -967,7 +968,7 @@ def test_Arctic_Ocean(self, METHOD, EXTRAPOLATE): # time dimension delta_time = 0.0 # calculate tide map - tide = pyTMD.compute_tide_corrections(xgrid, ygrid, delta_time, + tide = pyTMD.compute.tide_elevations(xgrid, ygrid, delta_time, DIRECTORY=filepath, MODEL='AOTIM-5-2018', GZIP=False, EPOCH=pyTMD.time._j2000_epoch, TYPE='grid', TIME='UTC', EPSG=3413, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE) diff --git a/test/test_perth3_read.py b/test/test_perth3_read.py index 9ee0aca2..3f43d13e 100644 --- a/test/test_perth3_read.py +++ b/test/test_perth3_read.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -test_perth3_read.py (04/2023) +test_perth3_read.py (01/2024) Tests that GOT4.7 data can be downloaded from AWS S3 bucket Tests the read program to verify that constituents are being extracted Tests that interpolated results are comparable to NASA PERTH3 program @@ -15,6 +15,7 @@ https://boto3.amazonaws.com/v1/documentation/api/latest/index.html UPDATE HISTORY: + Updated 01/2024: refactored compute functions into compute.py Updated 04/2023: using pathlib to define and expand paths Updated 12/2022: add check for read and interpolate constants Updated 09/2021: added test for model definition files @@ -38,8 +39,8 @@ import pyTMD.time import pyTMD.io.model import pyTMD.utilities +import pyTMD.compute import pyTMD.predict -import pyTMD.compute_tide_corrections import pyTMD.check_points # current file path @@ -238,7 +239,7 @@ def test_Ross_Ice_Shelf(METHOD, EXTRAPOLATE): # time dimension delta_time = 0.0 # calculate tide map - tide = pyTMD.compute_tide_corrections(xgrid, ygrid, delta_time, + tide = pyTMD.compute.tide_elevations(xgrid, ygrid, delta_time, DIRECTORY=filepath, MODEL='GOT4.7', GZIP=True, EPOCH=pyTMD.time._atlas_sdp_epoch, TYPE='grid', TIME='GPS', EPSG=3031, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE) @@ -277,5 +278,5 @@ def test_extend_array(): def test_unlisted_model(): ermsg = "Unlisted tide model" with pytest.raises(Exception, match=ermsg): - pyTMD.compute_tide_corrections(None, None, None, + pyTMD.compute.tide_elevations(None, None, None, DIRECTORY=filepath, MODEL='invalid') diff --git a/test/test_solid_earth.py b/test/test_solid_earth.py index 5c5d77f9..d4ab4c2c 100644 --- a/test/test_solid_earth.py +++ b/test/test_solid_earth.py @@ -1,5 +1,5 @@ """ -test_solid_earth.py (12/2023) +test_solid_earth.py (01/2024) Tests the steps for calculating the solid earth tides PYTHON DEPENDENCIES: @@ -8,6 +8,7 @@ https://numpy.org/doc/stable/user/numpy-for-matlab-users.html UPDATE HISTORY: + Updated 01/2024: refactored lunisolar ephemerides functions Updated 12/2023: phase_angles function renamed to doodson_arguments Updated 04/2023: added test for using JPL ephemerides for positions Written 04/2023 @@ -15,10 +16,10 @@ import pytest import numpy as np import pyTMD.astro +import pyTMD.compute import pyTMD.predict import pyTMD.time import pyTMD.utilities -from pyTMD.compute_tide_corrections import compute_SET_corrections def test_out_of_phase_diurnal(): """Test out-of-phase diurnal corrections with IERS outputs @@ -262,10 +263,10 @@ def test_solid_earth_radial(EPHEMERIDES): -0.09726749,-0.09726755,-0.11400376,-0.11400391, -0.11400412,-0.11400434]) # predict radial solid earth tides - tide_free = compute_SET_corrections(longitudes, latitudes, times, + tide_free = pyTMD.compute.SET_displacements(longitudes, latitudes, times, EPSG=4326, TYPE='drift', TIME='datetime', ELLIPSOID='WGS84', EPHEMERIDES=EPHEMERIDES) - tide_mean = compute_SET_corrections(longitudes, latitudes, times, + tide_mean = pyTMD.compute.SET_displacements(longitudes, latitudes, times, EPSG=4326, TYPE='drift', TIME='datetime', ELLIPSOID='WGS84', TIDE_SYSTEM='mean_tide', EPHEMERIDES=EPHEMERIDES) # as using estimated ephemerides, assert within 1/2 mm @@ -298,12 +299,12 @@ def download_jpl_ephemerides(): def test_solar_ecef(): """Test solar ECEF coordinates with ephemeride predictions """ - # calculate solar ephemerides MJD = 55414.0 - x1, y1, z1 = pyTMD.astro.solar_ecef(MJD) + # calculate approximate solar ephemerides + x1, y1, z1 = pyTMD.astro.solar_ecef(MJD, ephemerides='approximate') r1 = np.sqrt(x1**2 + y1**2 + z1**2) # predict solar ephemerides - x2, y2, z2 = pyTMD.astro.solar_ephemerides(MJD) + x2, y2, z2 = pyTMD.astro.solar_ecef(MJD, ephemerides='JPL') r2 = np.sqrt(x2**2 + y2**2 + z2**2) # test distances assert np.isclose(np.c_[x1,y1,z1], np.c_[x2,y2,z2], atol=1e9).all() @@ -313,12 +314,12 @@ def test_solar_ecef(): def test_lunar_ecef(): """Test lunar ECEF coordinates with ephemeride predictions """ - # calculate solar ephemerides MJD = 55414.0 - x1, y1, z1 = pyTMD.astro.lunar_ecef(MJD) + # calculate approximate lunar ephemerides + x1, y1, z1 = pyTMD.astro.lunar_ecef(MJD, ephemerides='approximate') r1 = np.sqrt(x1**2 + y1**2 + z1**2) - # predict solar ephemerides - x2, y2, z2 = pyTMD.astro.lunar_ephemerides(MJD) + # predict lunar ephemerides + x2, y2, z2 = pyTMD.astro.lunar_ecef(MJD, ephemerides='JPL') r2 = np.sqrt(x2**2 + y2**2 + z2**2) # test distances assert np.isclose(np.c_[x1,y1,z1], np.c_[x2,y2,z2], atol=5e6).all() From be5907b68d419ff218b51a3d237d6754262e7183 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 18 Jan 2024 11:29:23 +1300 Subject: [PATCH 3/7] refactor: moved constituent parameters function from predict to arguments feat: add functions for tide generating forces and potentials --- doc/source/api_reference/arguments.rst | 2 + doc/source/api_reference/predict.rst | 2 - pyTMD/arguments.py | 97 ++++++++++++++++ pyTMD/load_constituent.py | 7 +- pyTMD/predict.py | 106 ++---------------- pyTMD/solve.py | 147 ++++++++++++++++++++++++- 6 files changed, 252 insertions(+), 109 deletions(-) diff --git a/doc/source/api_reference/arguments.rst b/doc/source/api_reference/arguments.rst index 6e6de179..07794f3b 100644 --- a/doc/source/api_reference/arguments.rst +++ b/doc/source/api_reference/arguments.rst @@ -28,6 +28,8 @@ Calling Sequence .. autofunction:: pyTMD.arguments._minor_table +.. autofunction:: pyTMD.arguments._constituent_parameters + .. autofunction:: pyTMD.arguments._to_doodson_number .. autofunction:: pyTMD.arguments._from_doodson_number diff --git a/doc/source/api_reference/predict.rst b/doc/source/api_reference/predict.rst index 0a0638fb..aabb32e9 100644 --- a/doc/source/api_reference/predict.rst +++ b/doc/source/api_reference/predict.rst @@ -36,8 +36,6 @@ Calling Sequence .. autofunction:: pyTMD.predict.solid_earth_tide -.. autofunction:: pyTMD.predict._constituent_parameters - .. autofunction:: pyTMD.predict._out_of_phase_diurnal .. autofunction:: pyTMD.predict._out_of_phase_semidiurnal diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index 49fd173e..33e0f8ab 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -43,6 +43,8 @@ include multiples of 90 degrees as variable following Ray 2017 add function to calculate the Doodson numbers for constituents add option to return None and not raise error for Doodson numbers + moved constituent parameters function from predict to arguments + added more constituent parameters for OTIS/ATLAS predictions Updated 12/2023: made keyword argument for selecting M1 coefficients Updated 08/2023: changed ESR netCDF4 format to TMD3 format Updated 04/2023: using renamed astro mean_longitudes function @@ -920,6 +922,101 @@ def _minor_table(**kwargs): # return the coefficient table return coef +def _constituent_parameters(c, **kwargs): + """ + Loads parameters for a given tidal constituent + + Parameters + ---------- + c: list + tidal constituent ID + raise_error: bool, default False + Raise exception if constituent is unsupported + + Returns + ------- + amplitude: float + amplitude of equilibrium tide for tidal constituent (meters) + phase: float + phase of tidal constituent (radians) + omega: float + angular frequency of constituent (radians) + alpha: float + load love number of tidal constituent + species: float + spherical harmonic dependence of quadrupole potential + + 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`__ + + .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + """ + # default keyword arguments + kwargs.setdefault('raise_error', False) + # constituents array that are included in tidal program + cindex = ['m2', 's2', 'k1', 'o1', 'n2', 'p1', 'k2', 'q1', '2n2', 'mu2', + 'nu2', 'l2', 't2', 'j1', 'm1', 'oo1', 'rho1', 'mf', 'mm', 'ssa', + 'm4', 'ms4', 'mn4', 'm6', 'm8', 'mk3', 's6', '2sm2', '2mk3', + 'msf', 'sa', 'mt', '2q1'] + # species type (spherical harmonic dependence of quadrupole potential) + _species = np.array([2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]) + # Load Love numbers + # alpha = correction factor for first order load tides + _alpha = np.array([0.693, 0.693, 0.736, 0.695, 0.693, 0.706, 0.693, + 0.695, 0.693, 0.693, 0.693, 0.693, 0.693, 0.695, 0.695, 0.695, 0.695, + 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, + 0.693, 0.693, 0.693, 0.693, 0.693, 0.693]) + # omega: angular frequency of constituent, in radians + _omega = np.array([1.405189e-04, 1.454441e-04, 7.292117e-05, 6.759774e-05, + 1.378797e-04, 7.252295e-05, 1.458423e-04, 6.495854e-05, 1.352405e-04, + 1.355937e-04, 1.382329e-04, 1.431581e-04, 1.452450e-04, 7.556036e-05, + 7.028195e-05, 7.824458e-05, 6.531174e-05, 0.053234e-04, 0.026392e-04, + 0.003982e-04, 2.810377e-04, 2.859630e-04, 2.783984e-04, 4.215566e-04, + 5.620755e-04, 2.134402e-04, 4.363323e-04, 1.503693e-04, 2.081166e-04, + 4.925200e-06, 1.990970e-07, 7.962619e-06, 6.231934e-05]) + # Astronomical arguments (relative to t0 = 1 Jan 0:00 1992) + # phases for each constituent are referred to the time when the phase of + # the forcing for that constituent is zero on the Greenwich meridian + _phase = np.array([1.731557546, 0.000000000, 0.173003674, 1.558553872, + 6.050721243, 6.110181633, 3.487600001, 5.877717569, 4.086699633, + 3.463115091, 5.427136701, 0.553986502, 0.052841931, 2.137025284, + 2.436575100, 1.929046130, 5.254133027, 1.756042456, 1.964021610, + 3.487600001, 3.463115091, 1.731557546, 1.499093481, 5.194672637, + 6.926230184, 1.904561220, 0.000000000, 4.551627762, 3.809122439, + 4.551627762, 6.232786837, 3.720064066, 3.91369596]) + # amplitudes of equilibrium tide in meters + # _amplitude = np.array([0.242334,0.112743,0.141565,0.100661,0.046397, + _amplitude = np.array([0.2441, 0.112743, 0.141565, 0.100661, 0.046397, + 0.046848, 0.030684, 0.019273, 0.006141, 0.007408, 0.008811, 0.006931, + 0.006608, 0.007915, 0.007915, 0.004338, 0.003661, 0.042041, 0.022191, + 0.019567, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.003681, 0.003104, + 0.008044, 0.002565]) + + # map between input constituent and cindex + j = [j for j,val in enumerate(cindex) if (val == c.lower())] + # set the values for the constituent + if j: + amplitude, = _amplitude[j] + phase, = _phase[j] + omega, = _omega[j] + alpha, = _alpha[j] + species, = _species[j] + elif kwargs['raise_error']: + raise ValueError(f'Unsupported constituent {c}') + else: + amplitude = 0.0 + phase = 0.0 + omega = 0.0 + alpha = 0.0 + species = 0 + # return the values for the constituent + return (amplitude, phase, omega, alpha, species) + def _to_doodson_number(coef: list | np.ndarray, **kwargs): """ Converts Cartwright numbers into a Doodson number diff --git a/pyTMD/load_constituent.py b/pyTMD/load_constituent.py index 69f331c2..7d6805b5 100644 --- a/pyTMD/load_constituent.py +++ b/pyTMD/load_constituent.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -load_constituent.py (09/2023) +load_constituent.py (01/2024) Loads parameters for a given tidal constituent CALLING SEQUENCE: @@ -26,13 +26,14 @@ https://numpy.org/doc/stable/user/numpy-for-matlab-users.html UPDATE HISTORY: + Updated 01/2024: moved constituent parameters function to arguments Updated 09/2023: deprecated in favor of pyTMD.predict function Updated 04/2022: updated docstrings to numpy documentation format Updated 07/2020: add more constituents from OTPSnc and function docstrings Updated 09/2017: Rewritten in Python """ import warnings -import pyTMD.predict +import pyTMD.arguments def load_constituent(c): """ @@ -69,4 +70,4 @@ def load_constituent(c): warnings.warn("Deprecated. Please use pyTMD.predict instead", DeprecationWarning) # call updated function to not break current workflows - return pyTMD.predict._constituent_parameters(c) + return pyTMD.arguments._constituent_parameters(c) diff --git a/pyTMD/predict.py b/pyTMD/predict.py index 46482cf8..66ca4af8 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -21,7 +21,7 @@ UPDATE HISTORY: Updated 01/2024: moved minor arguments calculation into new function - added more constituent parameters for OTIS/ATLAS predictions + 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 @@ -103,7 +103,8 @@ def map(t: float | np.ndarray, for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent - amp, ph, omega, alpha, species = _constituent_parameters(c) + 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] elif corrections in ('GOT', 'FES'): @@ -168,7 +169,8 @@ def drift(t: float | np.ndarray, for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent - amp, ph, omega, alpha, species = _constituent_parameters(c) + 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] elif corrections in ('GOT', 'FES'): @@ -233,7 +235,8 @@ def time_series(t: float | np.ndarray, for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): # load parameters for each constituent - amp, ph, omega, alpha, species = _constituent_parameters(c) + 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] elif corrections in ('GOT', 'FES'): @@ -583,101 +586,6 @@ def solid_earth_tide( # return the solid earth tide return dxt -def _constituent_parameters(c, **kwargs): - """ - Loads parameters for a given tidal constituent - - Parameters - ---------- - c: list - tidal constituent ID - raise_error: bool, default False - Raise exception if constituent is unsupported - - Returns - ------- - amplitude: float - amplitude of equilibrium tide for tidal constituent (meters) - phase: float - phase of tidal constituent (radians) - omega: float - angular frequency of constituent (radians) - alpha: float - load love number of tidal constituent - species: float - spherical harmonic dependence of quadrupole potential - - 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`__ - - .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 - """ - # default keyword arguments - kwargs.setdefault('raise_error', False) - # constituents array that are included in tidal program - cindex = ['m2', 's2', 'k1', 'o1', 'n2', 'p1', 'k2', 'q1', '2n2', 'mu2', - 'nu2', 'l2', 't2', 'j1', 'm1', 'oo1', 'rho1', 'mf', 'mm', 'ssa', - 'm4', 'ms4', 'mn4', 'm6', 'm8', 'mk3', 's6', '2sm2', '2mk3', - 'msf', 'sa', 'mt', '2q1'] - # species type (spherical harmonic dependence of quadrupole potential) - _species = np.array([2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, - 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]) - # Load Love numbers - # alpha = correction factor for first order load tides - _alpha = np.array([0.693, 0.693, 0.736, 0.695, 0.693, 0.706, 0.693, - 0.695, 0.693, 0.693, 0.693, 0.693, 0.693, 0.695, 0.695, 0.695, 0.695, - 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, 0.693, - 0.693, 0.693, 0.693, 0.693, 0.693, 0.693]) - # omega: angular frequency of constituent, in radians - _omega = np.array([1.405189e-04, 1.454441e-04, 7.292117e-05, 6.759774e-05, - 1.378797e-04, 7.252295e-05, 1.458423e-04, 6.495854e-05, 1.352405e-04, - 1.355937e-04, 1.382329e-04, 1.431581e-04, 1.452450e-04, 7.556036e-05, - 7.028195e-05, 7.824458e-05, 6.531174e-05, 0.053234e-04, 0.026392e-04, - 0.003982e-04, 2.810377e-04, 2.859630e-04, 2.783984e-04, 4.215566e-04, - 5.620755e-04, 2.134402e-04, 4.363323e-04, 1.503693e-04, 2.081166e-04, - 4.925200e-06, 1.990970e-07, 7.962619e-06, 6.231934e-05]) - # Astronomical arguments (relative to t0 = 1 Jan 0:00 1992) - # phases for each constituent are referred to the time when the phase of - # the forcing for that constituent is zero on the Greenwich meridian - _phase = np.array([1.731557546, 0.000000000, 0.173003674, 1.558553872, - 6.050721243, 6.110181633, 3.487600001, 5.877717569, 4.086699633, - 3.463115091, 5.427136701, 0.553986502, 0.052841931, 2.137025284, - 2.436575100, 1.929046130, 5.254133027, 1.756042456, 1.964021610, - 3.487600001, 3.463115091, 1.731557546, 1.499093481, 5.194672637, - 6.926230184, 1.904561220, 0.000000000, 4.551627762, 3.809122439, - 4.551627762, 6.232786837, 3.720064066, 3.91369596]) - # amplitudes of equilibrium tide in meters - # _amplitude = np.array([0.242334,0.112743,0.141565,0.100661,0.046397, - _amplitude = np.array([0.2441, 0.112743, 0.141565, 0.100661, 0.046397, - 0.046848, 0.030684, 0.019273, 0.006141, 0.007408, 0.008811, 0.006931, - 0.006608, 0.007915, 0.007915, 0.004338, 0.003661, 0.042041, 0.022191, - 0.019567, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.003681, 0.003104, - 0.008044, 0.002565]) - - # map between input constituent and cindex - j = [j for j,val in enumerate(cindex) if (val == c.lower())] - # set the values for the constituent - if j: - amplitude, = _amplitude[j] - phase, = _phase[j] - omega, = _omega[j] - alpha, = _alpha[j] - species, = _species[j] - elif kwargs['raise_error']: - raise ValueError(f'Unsupported constituent {c}') - else: - amplitude = 0.0 - phase = 0.0 - omega = 0.0 - alpha = 0.0 - species = 0 - # return the values for the constituent - return (amplitude, phase, omega, alpha, species) - def _out_of_phase_diurnal( XYZ: np.ndarray, SXYZ: np.ndarray, diff --git a/pyTMD/solve.py b/pyTMD/solve.py index 437748c9..69e7ed6b 100644 --- a/pyTMD/solve.py +++ b/pyTMD/solve.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" solve.py -Written by Tyler Sutterley (12/2023) +Written by Tyler Sutterley (01/2024) Routines for estimating the harmonic constants for ocean tides REFERENCES: @@ -21,6 +21,7 @@ predict.py: predict tidal values using harmonic constants UPDATE HISTORY: + Updated 01/2024: add functions for tide generating forces and potentials Written 12/2023 """ @@ -28,8 +29,11 @@ import numpy as np import scipy.linalg -import pyTMD.predict -from pyTMD.arguments import arguments +import pyTMD.arguments +from pyTMD.constants import constants + +# get WGS84 ellipsoid parameters in meters-kg-seconds +_wgs84 = constants(ellipsoid='WGS84', units='MKS') def constants(t: float | np.ndarray, ht: np.ndarray, @@ -88,7 +92,7 @@ def constants(t: float | np.ndarray, # load the nodal corrections # convert time to Modified Julian Days (MJD) - pu, pf, G = arguments(t + 48622.0, constituents, + pu, pf, G = pyTMD.arguments.arguments(t + 48622.0, constituents, deltat=deltat, corrections=corrections) # create design matrix @@ -99,7 +103,7 @@ def constants(t: float | np.ndarray, for k,c in enumerate(constituents): if corrections in ('OTIS', 'ATLAS', 'TMD3', 'netcdf'): amp, ph, omega, alpha, species = \ - pyTMD.predict._constituent_parameters(c) + pyTMD.arguments._constituent_parameters(c) th = omega*t*86400.0 + ph + pu[:,k] elif corrections in ('GOT', 'FES'): th = G[:,k]*np.pi/180.0 + pu[:,k] @@ -129,3 +133,136 @@ def constants(t: float | np.ndarray, # return the amplitude and phase return (amp, phase) + +# PURPOSE: calculate the astronomical tide generating force +def _generating_force( + lon: np.ndarray, + lat: np.ndarray, + c: str, + ): + """ + Computes the astronomical tide generating force for a tidal + constituent [1]_ + + Parameters + ---------- + lon: np.ndarray + longitudes of grid centers (degrees east) + lat: np.ndarray + latitudes of grid centers (degrees north) + c: str + tidal constituent ID + + 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`__ + + .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + """ + + # grid spacing in radians + dph = (lon[1] - lon[0])*np.pi/180.0 + dth = (lat[1] - lat[0])*np.pi/180.0 + # calculate meshgrid from latitude and longitude + gridlon, gridlat = np.meshgrid(lon, lat) + + # longitude (in radians) of u and v nodes + phi_v = gridlon*np.pi/180.0 + phi_u = phi_u - dph/2.0 + + # colatitudes (in radians) of u and v nodes + th_u = (90.0 - gridlat)*np.pi/180.0 + th_v = th_v + dth/2.0 + + # load parameters for each constituent + amp, ph, omega, alpha, species = \ + pyTMD.arguments._constituent_parameters(c) + # uniform zonal dependence of astronomical forcing + ph = 0.0 + + # loading effects + c = (alpha*_wgs84.gamma*amp/_wgs84.rad_e) + # calculate forcing for constituent + Fu = c*np.exp(1j*phi_u*species + 1j*ph) + Fv = c*np.exp(1j*phi_v*species + 1j*ph) + # calculate latitudinal dependence of forcing + # for a given spherical harmonic dependence + if (species == 1): + # diurnal species + Fu *= 2j*np.cos(th_u) + Fv *= 2.0*(2.0*np.sin(th_v)**2 - 1.0) + elif (species == 2): + # semidiurnal species + Fu *= 2j*np.sin(th_u) + Fv *= -2.0*(np.sin(th_v)*np.cos(th_v)) + else: + # long-period species + Fu *= 0.0 + 0j + Fv *= -3.0*(np.sin(th_v)*np.cos(th_v)) + # return the generating forces + return (Fu, Fv) + +# PURPOSE: calculate the astronomical tide generating potential +# converted to an equilibrium tide elevation +def _generating_potential( + lon: np.ndarray, + lat: np.ndarray, + c: str, + ): + """ + Computes the astronomical tide generating potential for a tidal + constituent converted to an equilibrium tide elevation [1]_ + + Parameters + ---------- + lon: np.ndarray + longitudes of grid centers (degrees east) + lat: np.ndarray + latitudes of grid centers (degrees north) + c: str + tidal constituent ID + + 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`__ + + .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 + """ + + # calculate meshgrid from latitude and longitude + gridlon, gridlat = np.meshgrid(lon, lat) + + # longitude (in radians) + phi_h = gridlon*np.pi/180.0 + # colatitudes (in radians) + th_h = (90.0 - gridlat)*np.pi/180.0 + + # load parameters for each constituent + amp, ph, omega, alpha, species = \ + pyTMD.arguments._constituent_parameters(c) + # uniform zonal dependence of astronomical potential + ph = 0.0 + + # calculate potential for constituent + c = (alpha*amp) + zeta = c*np.exp(1j*phi_h*species + 1j*ph) + + # calculate latitudinal dependence of potential + # for a given spherical harmonic dependence + if (species == 1): + # diurnal species + zeta *= 2*0*np.sin(th_h)*np.cos(th_h) + elif (species == 2): + # semidiurnal species + zeta *= np.sin(th_h)**2 + else: + # long-period species + zeta *= -1.0*(3.0/2.0*np.cos(th_h)**2 - 1.0/2.0) + # return the generating potential + return zeta From bb68c24abbbca097cdb9b9851d8bc1d07d1bbbf4 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 18 Jan 2024 11:33:48 +1300 Subject: [PATCH 4/7] fix: variable typing for `c` in `_constituent_parameters` --- pyTMD/arguments.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index 33e0f8ab..71e3255c 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -922,13 +922,13 @@ def _minor_table(**kwargs): # return the coefficient table return coef -def _constituent_parameters(c, **kwargs): +def _constituent_parameters(c: str, **kwargs): """ Loads parameters for a given tidal constituent Parameters ---------- - c: list + c: str tidal constituent ID raise_error: bool, default False Raise exception if constituent is unsupported From be1c2685e1c18f21db827b293597ff31219d7b3a Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 18 Jan 2024 11:39:59 +1300 Subject: [PATCH 5/7] ci: omit deprecated functions in coverage report --- setup.cfg | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/setup.cfg b/setup.cfg index 5623e39c..bc2cde0c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -14,6 +14,12 @@ omit = setup.py conf.py scripts/* + pyTMD/calc_astrol_longitudes.py + pyTMD/compute_tide_corrections.py + pyTMD/convert_crs.py + pyTMD/convert_ll_xy.py + pyTMD/load_constituent.py + pyTMD/load_nodal_corrections.py [coverage:report] show_missing = true From 0b74d425f37d3931ffd0b01c4f38777c2398d050 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 18 Jan 2024 11:46:15 +1300 Subject: [PATCH 6/7] fix: flake8 check --- pyTMD/io/OTIS.py | 8 ++++---- pyTMD/solve.py | 6 ++++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/pyTMD/io/OTIS.py b/pyTMD/io/OTIS.py index 4f141179..a4232734 100644 --- a/pyTMD/io/OTIS.py +++ b/pyTMD/io/OTIS.py @@ -292,7 +292,7 @@ def extract_constants( # masks zero values mask = (hu == 0) | mu.astype(bool) bathymetry = np.ma.array(hu, mask=mask) - # x coordinates for u transports + # x-coordinates for u transports xi -= dx/2.0 elif kwargs['type'] in ('v','V'): # create current masks and bathymetry estimates @@ -307,7 +307,7 @@ def extract_constants( # masks zero values mask = (hv == 0) | mv.astype(bool) bathymetry = np.ma.array(hv, mask=mask) - # y coordinates for v transports + # y-coordinates for v transports yi -= dy/2.0 # interpolate bathymetry and mask to output points @@ -557,7 +557,7 @@ def read_constants( # masks zero values mask = (hu == 0) | mu.astype(bool) bathymetry = np.ma.array(hu, mask=mask) - # x coordinates for u transports + # x-coordinates for u transports xi -= dx/2.0 elif kwargs['type'] in ('v','V'): # create current masks and bathymetry estimates @@ -572,7 +572,7 @@ def read_constants( # masks zero values mask = (hv == 0) | mv.astype(bool) bathymetry = np.ma.array(hv, mask=mask) - # y coordinates for v transports + # y-coordinates for v transports yi -= dy/2.0 # read each constituent diff --git a/pyTMD/solve.py b/pyTMD/solve.py index 69e7ed6b..1ebd57cf 100644 --- a/pyTMD/solve.py +++ b/pyTMD/solve.py @@ -171,11 +171,13 @@ def _generating_force( # longitude (in radians) of u and v nodes phi_v = gridlon*np.pi/180.0 - phi_u = phi_u - dph/2.0 + # x-coordinates for u transports + phi_u = phi_v - dph/2.0 # colatitudes (in radians) of u and v nodes th_u = (90.0 - gridlat)*np.pi/180.0 - th_v = th_v + dth/2.0 + # y-coordinates for v transports + th_v = th_u + dth/2.0 # load parameters for each constituent amp, ph, omega, alpha, species = \ From f4bee55a93ae8077588b6b4842c3074d205c97cf Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 18 Jan 2024 13:24:11 +1300 Subject: [PATCH 7/7] docs: add toctree for `io` subdirectory test: add quick test for currents wrapper function --- doc/source/api_reference/io/ATLAS.rst | 9 +++---- doc/source/api_reference/io/FES.rst | 6 ++--- doc/source/api_reference/io/GOT.rst | 8 +++---- doc/source/api_reference/io/OTIS.rst | 12 ++++++---- doc/source/api_reference/io/constituents.rst | 6 ++--- doc/source/api_reference/io/io.rst | 14 +++++++++++ doc/source/api_reference/io/model.rst | 6 ++--- .../api_reference/io/ocean_pole_tide.rst | 6 ++--- doc/source/index.rst | 8 +------ pyTMD/load_constituent.py | 2 +- pyTMD/solve.py | 4 ++-- test/test_download_and_read.py | 24 +++++++++++++++++++ 12 files changed, 70 insertions(+), 35 deletions(-) create mode 100644 doc/source/api_reference/io/io.rst diff --git a/doc/source/api_reference/io/ATLAS.rst b/doc/source/api_reference/io/ATLAS.rst index 29c3fd85..57bd46fe 100644 --- a/doc/source/api_reference/io/ATLAS.rst +++ b/doc/source/api_reference/io/ATLAS.rst @@ -1,6 +1,6 @@ -======== -io.ATLAS -======== +===== +ATLAS +===== - Reads netCDF format tidal solutions provided by Ohio State University and ESR - Spatially interpolates tidal constituents to input coordinates @@ -12,7 +12,8 @@ Calling Sequence import pyTMD.io amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(ilon, ilat, grid_file, model_files, - type='z', method='spline') + type='z', + method='spline') `Source code`__ diff --git a/doc/source/api_reference/io/FES.rst b/doc/source/api_reference/io/FES.rst index 9c34b637..82ba61ca 100644 --- a/doc/source/api_reference/io/FES.rst +++ b/doc/source/api_reference/io/FES.rst @@ -1,6 +1,6 @@ -====== -io.FES -====== +=== +FES +=== - Reads files for Finite Element Solution (FES), Empirical Ocean Tide (EOT), and Hamburg direct data Assimilation Methods for Tides (HAMTIDE) models diff --git a/doc/source/api_reference/io/GOT.rst b/doc/source/api_reference/io/GOT.rst index f18b3263..4a578e09 100644 --- a/doc/source/api_reference/io/GOT.rst +++ b/doc/source/api_reference/io/GOT.rst @@ -1,6 +1,6 @@ -====== -io.GOT -====== +=== +GOT +=== - Reads files for Richard Ray's Global Ocean Tide (GOT) models - Spatially interpolates tidal constituents to input coordinates @@ -11,7 +11,7 @@ Calling Sequence .. code-block:: python import pyTMD.io - amp,ph,c = pyTMD.io.GOT.extract_constants(ilon,i lat, model_files, + amp,ph,c = pyTMD.io.GOT.extract_constants(ilon, ilat, model_files, method='spline') `Source code`__ diff --git a/doc/source/api_reference/io/OTIS.rst b/doc/source/api_reference/io/OTIS.rst index 7d1a7a09..6085350a 100644 --- a/doc/source/api_reference/io/OTIS.rst +++ b/doc/source/api_reference/io/OTIS.rst @@ -1,6 +1,6 @@ -======= -io.OTIS -======= +==== +OTIS +==== - Reads OTIS format tidal solutions provided by Ohio State University and ESR @@ -16,8 +16,10 @@ Calling Sequence .. code-block:: python import pyTMD.io - amp,ph,D,c = pyTMD.io.OTIS.extract_constants(ilon, ilat, grid_file, model_file, - EPSG, type='z', method='spline', grid='OTIS') + amp,ph,D,c = pyTMD.io.OTIS.extract_constants(ilon, ilat, grid_file, model_file, EPSG, + type='z', + method='spline', + grid='OTIS') `Source code`__ diff --git a/doc/source/api_reference/io/constituents.rst b/doc/source/api_reference/io/constituents.rst index 2fba5cad..f25d744e 100644 --- a/doc/source/api_reference/io/constituents.rst +++ b/doc/source/api_reference/io/constituents.rst @@ -1,6 +1,6 @@ -=============== -io.constituents -=============== +============ +constituents +============ Basic tide model constituent class diff --git a/doc/source/api_reference/io/io.rst b/doc/source/api_reference/io/io.rst new file mode 100644 index 00000000..17523d74 --- /dev/null +++ b/doc/source/api_reference/io/io.rst @@ -0,0 +1,14 @@ +== +io +== + +.. toctree:: + :maxdepth: 1 + + ATLAS.rst + FES.rst + GOT.rst + OTIS.rst + constituents.rst + model.rst + ocean_pole_tide.rst diff --git a/doc/source/api_reference/io/model.rst b/doc/source/api_reference/io/model.rst index ed426a74..8706dc8b 100644 --- a/doc/source/api_reference/io/model.rst +++ b/doc/source/api_reference/io/model.rst @@ -1,6 +1,6 @@ -======== -io.model -======== +===== +model +===== Retrieves tide model parameters for named tide models and from model definition files diff --git a/doc/source/api_reference/io/ocean_pole_tide.rst b/doc/source/api_reference/io/ocean_pole_tide.rst index 370928f4..1e6991af 100644 --- a/doc/source/api_reference/io/ocean_pole_tide.rst +++ b/doc/source/api_reference/io/ocean_pole_tide.rst @@ -1,6 +1,6 @@ -================== -io.ocean_pole_tide -================== +=============== +ocean_pole_tide +=============== - Reads ocean pole load tide coefficients provided by IERS as computed by `Desai et al. (2002) `_ and `Desai et al. (2015) `_ - See `materials from Chapter 7 of the IERS Conventions `_ diff --git a/doc/source/index.rst b/doc/source/index.rst index 0badaea4..01ee5540 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -39,13 +39,7 @@ ocean, load, solid Earth and pole tides api_reference/ellipse.rst api_reference/eop.rst api_reference/interpolate.rst - api_reference/io/ATLAS.rst - api_reference/io/FES.rst - api_reference/io/GOT.rst - api_reference/io/OTIS.rst - api_reference/io/constituents.rst - api_reference/io/model.rst - api_reference/io/ocean_pole_tide.rst + api_reference/io/io.rst api_reference/predict.rst api_reference/solve.rst api_reference/spatial.rst diff --git a/pyTMD/load_constituent.py b/pyTMD/load_constituent.py index 7d6805b5..9e251ce2 100644 --- a/pyTMD/load_constituent.py +++ b/pyTMD/load_constituent.py @@ -67,7 +67,7 @@ def load_constituent(c): .. __: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2 """ # raise warning for deprecated function call - warnings.warn("Deprecated. Please use pyTMD.predict instead", + warnings.warn("Deprecated. Please use pyTMD.arguments instead", DeprecationWarning) # call updated function to not break current workflows return pyTMD.arguments._constituent_parameters(c) diff --git a/pyTMD/solve.py b/pyTMD/solve.py index 1ebd57cf..db3bb35d 100644 --- a/pyTMD/solve.py +++ b/pyTMD/solve.py @@ -141,8 +141,8 @@ def _generating_force( c: str, ): """ - Computes the astronomical tide generating force for a tidal - constituent [1]_ + Computes the astronomical tide generating force (the horizontal gradient + of the tide generating potential) for a tidal constituent [1]_ Parameters ---------- diff --git a/test/test_download_and_read.py b/test/test_download_and_read.py index 777618b6..2ab7be04 100644 --- a/test/test_download_and_read.py +++ b/test/test_download_and_read.py @@ -52,6 +52,7 @@ import numpy as np import pyTMD.io import pyTMD.io.model +import pyTMD.compute import pyTMD.predict import pyTMD.time import pyTMD.utilities @@ -697,6 +698,29 @@ def test_Ross_Ice_Shelf(self, METHOD, EXTRAPOLATE): EPSG=3031, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE) assert np.any(tide) + # parameterize interpolation method + # only use fast interpolation routines + @pytest.mark.parametrize("METHOD", ['spline','nearest']) + @pytest.mark.parametrize("EXTRAPOLATE", [True]) + # PURPOSE: test the tide currents wrapper function + def test_Ross_Ice_Shelf_currents(self, METHOD, EXTRAPOLATE): + # create a drift track along the Ross Ice Shelf + xlimits = np.array([-740000,520000]) + ylimits = np.array([-1430000,-300000]) + # x and y coordinates + x = np.linspace(xlimits[0],xlimits[1],24) + y = np.linspace(ylimits[0],ylimits[1],24) + # time dimension + delta_time = np.zeros((24))*3600 + # calculate tide drift corrections + tide = pyTMD.compute.tide_currents(x, y, delta_time, + DIRECTORY=filepath, MODEL='CATS2008', GZIP=False, + EPOCH=pyTMD.time._j2000_epoch, TYPE='drift', TIME='UTC', + EPSG=3031, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE) + # iterate over zonal and meridional currents + for key,val in tide.items(): + assert np.any(val) + # PURPOSE: test definition file functionality @pytest.mark.parametrize("MODEL", ['CATS2008']) def test_definition_file(self, MODEL):