feat: add wrapper functions to read and interpolate constants (#349)

* tests: use simplified wrapper functions
* feat: updated computation of long-period equilibrium tides
* feat: add functions to append node tide equilibrium values to amplitudes
* fix: add messaging if there are no minor constituents to infer
* feat: can convert Doodson numbers formatted as strings
* docs: expand the glossary

doc/source/api_reference/io/model.rst

@@ -15,3 +15,9 @@ General Attributes and Methods
 
 .. autoclass:: pyTMD.io.model
    :members:
+
+.. autofunction:: pyTMD.io.extract_constants
+
+.. autofunction:: pyTMD.io.read_constants
+
+.. autofunction:: pyTMD.io.interpolate_constants

doc/source/getting_started/Glossary.rst

@@ -39,14 +39,17 @@ Glossary
         see :term:`Harmonic Constituents`
 
     Declination
-        angular distance of an astronomical body north or south of the celestial equator
+        angular distance of an astronomical body north or south of the celestial equator, with north being positive and south negative
 
     Descending Node
         point of an orbit where a celestial body intersects the ecliptic, and the latitude coordinate is decreasing
 
     Diurnal Tide
         tidal oscillations with a near-daily period
 
+    Earth Tide
+        tidal variations of the surface of the solid Earth, frequently split into the :term:`Body Tide` and the :term:`Load Tide`
+
     Ecliptic
         plane of the orbit of the Earth around the sun
 
@@ -60,11 +63,17 @@ Glossary
         fixed point in time used as a reference value
 
     Equilibrium Tide
-        hypothetical tide produced solely by lunisolar gravitational forces under equilibrium theory
+        hypothetical tide produced solely by lunisolar gravitational forces under equilibrium theory, in the absence of inertia and currents
+
+    Equinox
+        the ascending node of the Earth's orbit about the sun, and the point on the equatorial plane where it meets the ecliptic plane
 
     Frequency
         number of cycles in a unit time
 
+    Geocentric Tide
+        Tidal variation relative to the Earth's center of figure, which is equal to the sum of the :term:`Ocean Tide` and :term:`Earth Tide`
+
     Geoid
         equipotential surface coinciding with the ocean surface in the absence of astronomical or dynamical effects
 
@@ -75,23 +84,41 @@ Glossary
         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 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`
+        method of estimating tidal elevations and currents through a combination of the :term:`Harmonic Constituents`
+
+    High Water Height
+        height of the maximum vertical elevation of the tide relative to a fixed datum
+
+        see :term:`Low Water Height`
+
+    King Tide
+        see :term:`Spring Tide`
 
     Load Tide
         elastic deformation of the solid Earth due to ocean and atmospheric tides
 
     Long Period Tide
-        tidal oscillations with periods much greater than one day
+        tidal oscillations with periods much greater than one day (typically 9.1 days to 18.6 years)
 
     Love and Shida Numbers
         dimensionless parameters relating the vertical (`h`), horizontal (`l`) and gravitational (`k`) elastic responses to tidal loading
 
+    Low Water Height
+        height of the minimum vertical elevation of the tide relative to a fixed datum
+
+        see :term:`High Water Height`
+
     Mean Tide
         model with both direct and indirect permanent tidal effects retained
 
+    Neap Tide
+        tides of the least amplitude within a 15-day cycle
+
+        see :term:`Spring Tide`
+
     Nodal Corrections
         adjustments to the amplitudes and phases of harmonic constituents to allow for periodic modulations
 
@@ -125,12 +152,20 @@ Glossary
     Pole Tide
         apparent tide due to variations in the Earth's figure axis about its mean
 
+    Range
+        Height difference between the :term:`High Water Height` and the :term:`Low Water Height`
+
     Semi-diurnal Tide
         tidal oscillations with an approximate half-day period
 
     Solid Earth Tide
         deformation of the solid Earth due to gravitational forces
 
+    Spring Tide
+        tides of the greatest amplitude within a 15-day cycle, sometimes called a :term:`King Tide`
+
+        see :term:`Neap Tide`
+
     Tidal Current
         horizontal movement of water due to periodic forces
 

pyTMD/arguments.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 arguments.py
-Written by Tyler Sutterley (09/2024)
+Written by Tyler Sutterley (10/2024)
 Calculates the nodal corrections for tidal constituents
 Modification of ARGUMENTS fortran subroutine by Richard Ray 03/1999
 
@@ -38,6 +38,7 @@
         Ocean Tides", Journal of Atmospheric and Oceanic Technology, (2002).
 
 UPDATE HISTORY:
+    Updated 10/2024: can convert Doodson numbers formatted as strings
     Updated 09/2024: add function to calculate tidal angular frequencies
     Updated 08/2024: add support for constituents in PERTH5 tables
         add back nodal arguments from PERTH3 for backwards compatibility
@@ -1993,13 +1994,13 @@ def _to_doodson_number(coef: list | np.ndarray, **kwargs):
         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):
+def _from_doodson_number(DO: str | float | np.ndarray):
     """
     Converts Doodson numbers into Cartwright numbers
 
     Parameters
     ----------
-    DO: float or np.ndarray
+    DO: str, float or np.ndarray
         Doodson number for constituent
 
     Returns
@@ -2008,6 +2009,8 @@ def _from_doodson_number(DO: float | np.ndarray):
         Doodson coefficients (Cartwright numbers) for constituent
     """
     # convert from Doodson number to Cartwright numbers
+    # verify Doodson numbers are floating point variables
+    DO = np.array(DO).astype(float)
     # multiply by 1000 to prevent floating point errors
     coef = np.array([np.mod(1e3*DO, 10**(6-o))//10**(5-o) for o in range(6)])
     # remove 5 from values following Doodson convention

pyTMD/compute.py

@@ -61,6 +61,8 @@
 
 UPDATE HISTORY:
     Updated 10/2024: compute delta times based on corrections type
+        simplify by using wrapper functions to read and interpolate constants
+        added option to append equilibrium amplitudes for node tides
     Updated 09/2024: use JSON database for known model parameters
         drop support for the ascii definition file format
         use model class attributes for file format and corrections
@@ -224,6 +226,7 @@ def tide_elevations(
         CORRECTIONS: str | None = None,
         INFER_MINOR: bool = True,
         MINOR_CONSTITUENTS: list | None = None,
+        APPEND_NODE: bool = False,
         APPLY_FLEXURE: bool = False,
         FILL_VALUE: float = np.nan,
         **kwargs
@@ -290,6 +293,8 @@ def tide_elevations(
         Infer the height values for minor tidal constituents
     MINOR_CONSTITUENTS: list or None, default None
         Specify constituents to infer
+    APPEND_NODE: bool, default False
+        Append equilibrium amplitudes for node tides
     APPLY_FLEXURE: bool, default False
         Apply ice flexure scaling factor to height values
 
@@ -353,29 +358,10 @@ def tide_elevations(
     nt = len(ts)
 
     # read tidal constants and interpolate to grid points
-    if model.format in ('OTIS', 'ATLAS-compact', 'TMD3'):
-        amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,
-            model.model_file, model.projection, type=model.type,
-            grid=model.file_format, crop=CROP, bounds=BOUNDS, method=METHOD,
-            extrapolate=EXTRAPOLATE, cutoff=CUTOFF, apply_flexure=APPLY_FLEXURE)
-    elif model.format in ('ATLAS-netcdf',):
-        amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file,
-            model.model_file, type=model.type, crop=CROP, bounds=BOUNDS,
-            method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
-            scale=model.scale, compressed=model.compressed)
-    elif model.format in ('GOT-ascii', 'GOT-netcdf'):
-        amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file,
-            grid=model.file_format, crop=CROP, bounds=BOUNDS, method=METHOD,
-            extrapolate=EXTRAPOLATE, cutoff=CUTOFF, scale=model.scale,
-            compressed=model.compressed)
-    elif model.format in ('FES-ascii', 'FES-netcdf'):
-        amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file,
-            type=model.type, version=model.version, crop=CROP, bounds=BOUNDS,
-            method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
-            scale=model.scale, compressed=model.compressed)
-        # available model constituents
-        c = model.constituents
-
+    amp, ph, c = model.extract_constants(lon, lat, type=model.type,
+        crop=CROP, bounds=BOUNDS, method=METHOD,
+        extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
+        append_node=APPEND_NODE, apply_flexure=APPLY_FLEXURE)
     # calculate complex phase in radians for Euler's
     cph = -1j*ph*np.pi/180.0
     # calculate constituent oscillation
@@ -593,24 +579,9 @@ def tide_currents(
     # 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-compact', 'TMD3'):
-            amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,
-                model.model_file['u'], model.projection, type=t,
-                grid=model.file_format, crop=CROP, bounds=BOUNDS,
-                method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF)
-        elif model.format in ('ATLAS-netcdf',):
-            amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file,
-                model.model_file[t], type=t, crop=CROP, bounds=BOUNDS,
-                method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
-                scale=model.scale, compressed=model.compressed)
-        elif model.format in ('FES-ascii', 'FES-netcdf'):
-            amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file[t],
-                type=t, version=model.version, crop=CROP, bounds=BOUNDS,
-                method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
-                scale=model.scale, compressed=model.compressed)
-            # available model constituents
-            c = model.constituents
-
+        amp, ph, c = model.extract_constants(lon, lat, type=t,
+            crop=CROP, bounds=BOUNDS, method=METHOD,
+            extrapolate=EXTRAPOLATE, cutoff=CUTOFF)
         # calculate complex phase in radians for Euler's
         cph = -1j*ph*np.pi/180.0
         # calculate constituent oscillation

pyTMD/io/OTIS.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 OTIS.py
-Written by Tyler Sutterley (09/2024)
+Written by Tyler Sutterley (10/2024)
 
 Reads files for a tidal model and makes initial calculations to run tide program
 Includes functions to extract tidal harmonic constants from OTIS tide models for
@@ -58,6 +58,7 @@
     interpolate.py: interpolation routines for spatial data
 
 UPDATE HISTORY:
+    Updated 10/2024: save latitude and longitude to output constituent object
     Updated 09/2024: using new JSON dictionary format for model projections
     Updated 08/2024: revert change and assume crop bounds are projected
     Updated 07/2024: added crop and bounds keywords for trimming model data
@@ -647,14 +648,19 @@ def read_constants(
         # y-coordinates for v transports
         yi -= dy/2.0
 
+    # calculate geographic coordinates of model grid
+    gridx, gridy = np.meshgrid(xi, yi)
+    lon, lat = crs.transform(gridx, gridy, direction='INVERSE')
+
     # read each constituent
     if isinstance(model_file, list):
         cons = [read_constituents(m)[0].pop() for m in model_file]
     else:
         cons,_ = read_constituents(model_file, grid=kwargs['grid'])
     # save output constituents and coordinate reference system
     constituents = pyTMD.io.constituents(x=xi, y=yi,
-        bathymetry=bathymetry.data, mask=mask, crs=crs)
+        bathymetry=bathymetry.data, mask=mask, crs=crs,
+        longitude=lon, latitude=lat)
 
     # read each model constituent
     for i,c in enumerate(cons):

pyTMD/io/__init__.py

@@ -7,4 +7,10 @@
 from .OTIS import *
 from .IERS import *
 from .constituents import constituents
-from .model import model, load_database
+from .model import (
+    model,
+    load_database,
+    extract_constants,
+    read_constants,
+    interpolate_constants
+)

pyTMD/io/constituents.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 constituents.py
-Written by Tyler Sutterley (09/2024)
+Written by Tyler Sutterley (10/2024)
 Basic tide model constituent class
 
 PYTHON DEPENDENCIES:
@@ -10,6 +10,7 @@
         https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
 
 UPDATE HISTORY:
+    Updated 10/2024: added property for the shape of constituent fields
     Updated 09/2024: add more known constituents to string parser function
     Updated 08/2024: add GOT prime nomenclature for 3rd degree constituents
     Updated 07/2024: add function to parse tidal constituents from strings
@@ -197,6 +198,16 @@ def cartwright_number(self):
         # return the list of Cartwright numbers
         return cartwright_numbers
 
+    @property
+    def shape(self):
+        """Shape of constituent fields
+        """
+        try:
+            field = self.fields[0]
+            return getattr(self, field).shape
+        except:
+            return None
+
     @staticmethod
     def parse(constituent: str) -> str:
         """