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

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

@@ -62,6 +62,7 @@
 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
@@ -225,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
@@ -291,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
 
@@ -357,7 +361,7 @@ def tide_elevations(
     amp, ph, c = model.extract_constants(lon, lat, type=model.type,
         crop=CROP, bounds=BOUNDS, method=METHOD,
         extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
-        apply_flexure=APPLY_FLEXURE)
+        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

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/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:
         """

pyTMD/io/model.py

@@ -12,6 +12,7 @@
 
 UPDATE HISTORY:
     Updated 10/2024: add wrapper functions to read and interpolate constants
+        add functions to append node tide equilibrium values to amplitudes
     Updated 09/2024: use JSON database for known model parameters
         drop support for the ascii definition file format
         add file_format and nodal correction attributes
@@ -977,6 +978,8 @@ def extract_constants(self,
             Longitude point in degrees
         lat: np.ndarray
             Latitude point in degrees
+        append_node: bool, default False
+            Append equilibrium amplitudes for node tides
         kwargs: dict
             Keyword arguments for extracting constants
 
@@ -993,6 +996,7 @@ def extract_constants(self,
         from pyTMD.io import OTIS, ATLAS, GOT, FES
         # set default keyword arguments
         kwargs.setdefault('type', self.type)
+        kwargs.setdefault('append_node', False)
         kwargs.setdefault('crop', False)
         kwargs.setdefault('bounds', None)
         kwargs.setdefault('method', 'spline')
@@ -1041,6 +1045,15 @@ def extract_constants(self,
                 **kwargs)
             # available model constituents
             c = self.constituents
+        # append node equilibrium tide if not in constituents list
+        if kwargs['append_node'] and ('node' not in c):
+            # calculate node equilibrium tide
+            anode, pnode = self.node_equilibrium(lat)
+            # concatenate to output
+            amp = np.ma.concatenate((amp, anode[:,None]), axis=1)
+            ph = np.ma.concatenate((ph, pnode[:,None]), axis=1)
+            # convert to complex amplitude and append
+            c.append('node')
         # return the amplitude, phase, and constituents
         return (amp, ph, c)
 
@@ -1050,13 +1063,16 @@ def read_constants(self, **kwargs):
 
         Parameters
         ----------
+        append_node: bool, default False
+            Append equilibrium amplitudes for node tides
         kwargs: dict
             Keyword arguments for reading constants
         """
         # import tide model functions
         from pyTMD.io import OTIS, ATLAS, GOT, FES
         # set default keyword arguments
         kwargs.setdefault('type', self.type)
+        kwargs.setdefault('append_node', False)
         kwargs.setdefault('crop', False)
         kwargs.setdefault('bounds', None)
         kwargs.setdefault('apply_flexure', False)
@@ -1096,6 +1112,21 @@ def read_constants(self, **kwargs):
             c = FES.read_constants(model_file,
                 version=self.version, compressed=self.compressed,
                 **kwargs)
+        # append node equilibrium tide if not in constituents list
+        if kwargs['append_node'] and ('node' not in c.fields):
+            # calculate node equilibrium tide
+            amp, ph = self.node_equilibrium(c.latitude)
+            # broadcast to shape of constituents
+            if (np.shape(c.latitude) != c.shape):
+                amp = np.broadcast_to(amp[:,None], c.shape, subok=True)
+                ph = np.broadcast_to(ph[:,None], c.shape, subok=True)
+                amp.mask = np.zeros_like(amp.data, dtype=bool)
+                ph.mask = np.zeros_like(ph.data, dtype=bool)
+            # calculate complex phase in radians for Euler's
+            cph = -1j*ph*np.pi/180.0
+            hc = amp*np.exp(cph)
+            # append constituent
+            c.append('node', hc)
         # return the tidal constituents
         self._constituents = c
         return c
@@ -1153,6 +1184,53 @@ def interpolate_constants(self,
         # return the amplitude and phase
         return (amp, ph)
 
+    @staticmethod
+    def node_equilibrium(lat: np.ndarray):
+        """
+        Compute the complex amplitude of the 18.6 year equilibrium node tide
+        for inferrence when not supplied as a constituent [1]_ [2]_
+
+        Parameters
+        ----------
+        lat: np.ndarray
+            latitude (degrees north)
+
+        Returns
+        -------
+        amp: np.ndarray
+            Tidal amplitude in meters
+        phase: np.ndarray
+            Tidal phase in degrees
+
+        References
+        ----------
+        .. [1] D. E. Cartwright and R. J. Tayler,
+            "New Computations of the Tide-generating Potential,"
+            *Geophysical Journal of the Royal Astronomical Society*,
+            23(1), 45--73. (1971). `doi: 10.1111/j.1365-246X.1971.tb01803.x
+            <https://doi.org/10.1111/j.1365-246X.1971.tb01803.x>`_
+        .. [2] D. E. Cartwright and A. C. Edden,
+            "Corrected Tables of Tidal Harmonics,"
+            *Geophysical Journal of the Royal Astronomical Society*,
+            33(3), 253--264, (1973). `doi: 10.1111/j.1365-246X.1973.tb03420.x
+            <https://doi.org/10.1111/j.1365-246X.1973.tb03420.x>`_
+        """
+        # Cartwright and Edden potential amplitude
+        amajor = 0.027929# node
+        # love numbers
+        k2 = 0.302
+        h2 = 0.609
+        gamma_2 = (1.0 + k2 - h2)
+        # 2nd degree Legendre polynomials
+        P20 = 0.5*(3.0*np.sin(lat*np.pi/180.0)**2 - 1.0)
+        # calculate equilibrium node constants
+        amp = np.ma.zeros_like(lat, dtype=np.float64)
+        amp.data[:] = amajor*gamma_2*P20*np.sqrt((4.0 + 1.0)/(4.0*np.pi))
+        amp.mask = np.zeros_like(lat, dtype=bool)
+        phase = np.ma.zeros_like(amp)
+        phase.data[:] = 180.0
+        return (amp, phase)
+
     def __str__(self):
         """String representation of the ``io.model`` object
         """