Python implementation

.circleci/config.yml

@@ -62,3 +62,5 @@ jobs:
             # run tests
             python ${PYSOLID_HOME}/tests/point.py
             python ${PYSOLID_HOME}/tests/grid.py
+            python ${PYSOLID_HOME}/tests/point_pyimpl.py
+            python ${PYSOLID_HOME}/tests/grid_pyimpl.py

src/pysolid/__init__.py

@@ -11,16 +11,17 @@
 
 
 # top-level functions
-from pysolid.grid import (
+from .grid import (
     calc_solid_earth_tides_grid,
     plot_solid_earth_tides_grid,
 )
-from pysolid.point import (
+from .point import (
     TIDES,
     calc_solid_earth_tides_point,
     plot_solid_earth_tides_point,
     plot_power_spectral_density4tides,
 )
+from . import py_solid
 
 __all__ = [
     '__version__',
@@ -30,4 +31,5 @@
     'calc_solid_earth_tides_point',
     'plot_solid_earth_tides_point',
     'plot_power_spectral_density4tides',
+    'py_solid',
 ]

src/pysolid/py_solid/__init__.py

@@ -0,0 +1,20 @@
+from .grid import (
+    calc_solid_earth_tides_grid,
+    plot_solid_earth_tides_grid,
+)
+
+from .point import (
+    TIDES,
+    calc_solid_earth_tides_point,
+    plot_solid_earth_tides_point,
+    plot_power_spectral_density4tides,
+)
+
+__all__ = [
+    'TIDES',
+    'calc_solid_earth_tides_grid',
+    'plot_solid_earth_tides_grid',
+    'calc_solid_earth_tides_point',
+    'plot_solid_earth_tides_point',
+    'plot_power_spectral_density4tides',
+]

src/pysolid/py_solid/grid.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python3
+#######################################################################
+# A Python wrapper for solid Earth tides calculation using solid.for.
+#   Fortran code is originally written by Dennis Milbert, 2018-06-01.
+#   Available at: http://geodesyworld.github.io/SOFTS/solid.htm.
+# Author: Zhang Yunjun, Simran Sangha, Sep 2020
+# Copyright 2020, by the California Institute of Technology.
+#######################################################################
+# Recommend usage:
+#   import pysolid
+#   pysolid.calc_solid_earth_tides_grid()
+
+
+import os
+
+import numpy as np
+from scipy import ndimage
+
+
+##################################  Earth tides - grid mode  ###################################
+def calc_solid_earth_tides_grid(dt_obj, atr, step_size=1e3, display=False, verbose=True):
+    """Calculate SET in east/north/up direction for a spatial grid at a given date/time.
+
+    Note that we use step_size to speedup >30 times, by feeding the Fortran code (SET calc and
+    text file writing) the coarse grid, then resize the output to the same shape as the original
+    input size. This uses the fact that SET varies slowly in space. Comparison w and w/o step_size
+    shows a difference in tide_u with max of 5e-8 m, thus negligible.
+
+    Parameters: dt_obj    - datetime.datetime object (with precision up to the second)
+                atr       - dict, metadata including the following keys:
+                                LENGTH/WIDTTH
+                                X/Y_FIRST
+                                X/Y_STEP
+                step_size - float, grid step feeded into the fortran code in meters
+                                to speedup the calculation
+                display   - bool, plot the calculated SET
+                verbose   - bool, print verbose message
+    Returns:    tide_e    - 2D np.ndarray, SET in east  direction in meters
+                tide_n    - 2D np.ndarray, SET in north direction in meters
+                tide_u    - 2D np.ndarray, SET in up    direction in meters
+    Examples:   atr = readfile.read_attribute('geo_velocity.h5')
+                tide_e, tide_n, tide_u = calc_solid_earth_tides_grid('20180219', atr)
+    """
+    try:
+        from .solid import solid_grid
+    except ImportError:
+        msg = "Cannot import name 'solid' from 'pysolid'!"
+        msg += '\n    Maybe solid.for is NOT compiled yet.'
+        msg += '\n    Check instruction at: https://github.com/insarlab/PySolid.'
+        raise ImportError(msg)
+
+    vprint = print if verbose else lambda *args, **kwargs: None
+
+    # location
+    lat0 = float(atr['Y_FIRST'])
+    lon0 = float(atr['X_FIRST'])
+    lat1 = lat0 + float(atr['Y_STEP']) * int(atr['LENGTH'])
+    lon1 = lon0 + float(atr['X_STEP']) * int(atr['WIDTH'])
+
+    vprint('PYSOLID: ----------------------------------------')
+    vprint('PYSOLID: datetime: {}'.format(dt_obj.isoformat()))
+    vprint('PYSOLID: SNWE: {}'.format((lat1, lat0, lon0, lon1)))
+
+    # step size
+    num_step = int(step_size / 108e3 / abs(float(atr['Y_STEP'])))
+    num_step = max(1, num_step)
+    length = np.rint(int(atr['LENGTH']) / num_step - 1e-4).astype(int)
+    width = np.rint(int(atr['WIDTH']) / num_step - 1e-4).astype(int)
+    lat_step = float(atr['Y_STEP']) * num_step
+    lon_step = float(atr['X_STEP']) * num_step
+    vprint('SOLID  : calculate solid Earth tides in east/north/up direction')
+    vprint('SOLID  : shape: {s}, step size: {la:.4f} by {lo:.4f} deg'.format(
+        s=(length, width), la=lat_step, lo=lon_step))
+
+    lats = lat0 + np.arange(length) * lat_step
+    lons = lon0 + np.arange(width) * lon_step
+
+    # calc solid Earth tides
+    tides = solid_grid(dt_obj, lats, lons)
+    tide_e = tides[0]
+    tide_n = tides[1]
+    tide_u = tides[2]
+
+    # resample to the input size
+    # via scipy.ndimage.zoom or skimage.transform.resize
+    if num_step > 1:
+        in_shape = tide_e.shape
+        out_shape = (int(atr['LENGTH']), int(atr['WIDTH']))
+        vprint(
+            f'PYSOLID: resize data to the shape of {out_shape} '
+            'using order-1 spline interpolation')
+
+        enu = np.stack([tide_e, tide_n, tide_u])
+        zoom_factors = [1, *np.divide(out_shape, in_shape)]
+        kwargs = dict(order=1, mode="nearest", grid_mode=True)
+        tide_e, tide_n, tide_u = ndimage.zoom(enu, zoom_factors, **kwargs)
+
+    # plot
+    if display:
+        plot_solid_earth_tides_grid(tide_e, tide_n, tide_u, dt_obj)
+
+    return tide_e, tide_n, tide_u
+
+
+#########################################  Plot  ###############################################
+def plot_solid_earth_tides_grid(tide_e, tide_n, tide_u, dt_obj=None,
+                                out_fig=None, save=False, display=True):
+    """Plot the solid Earth tides in ENU direction."""
+    from matplotlib import pyplot as plt, ticker
+
+    # plot
+    fig, axs = plt.subplots(nrows=1, ncols=3, figsize=[6, 3], sharex=True, sharey=True)
+    for ax, data, label in zip(axs.flatten(),
+                               [tide_e, tide_n, tide_u],
+                               ['East', 'North', 'Up']):
+        im = ax.imshow(data*100, cmap='RdBu')
+        ax.tick_params(which='both', direction='in', bottom=True, top=True, left=True, right=True)
+        fig.colorbar(im, ax=ax, orientation='horizontal', label=label+' [cm]', pad=0.1, ticks=ticker.MaxNLocator(3))
+    fig.tight_layout()
+
+    # super title
+    if dt_obj is not None:
+        axs[1].set_title('solid Earth tides at {}'.format(dt_obj.isoformat()), fontsize=12)
+
+    # output
+    if out_fig:
+        save = True
+
+    if save:
+        if not out_fig:
+            out_fig = os.path.abspath('SET_grid.png')
+        print('save figure to {}'.format(out_fig))
+        fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)
+
+    if display:
+        print('showing...')
+        plt.show()
+    else:
+        plt.close()
+
+    return

src/pysolid/py_solid/point.py

@@ -0,0 +1,295 @@
+#!/usr/bin/env python3
+#######################################################################
+# A Python wrapper for solid Earth tides calculation using solid.for.
+#   Fortran code is originally written by Dennis Milbert, 2018-06-01.
+#   Available at: http://geodesyworld.github.io/SOFTS/solid.htm.
+# Author: Zhang Yunjun, Jan 2021
+# Copyright 2020, by the California Institute of Technology.
+#######################################################################
+# Recommend usage:
+#   import pysolid
+#   pysolid.calc_solid_earth_tides_point()
+
+
+import collections
+import datetime as dt
+import os
+
+import numpy as np
+
+
+## Tidal constituents
+# https://en.wikipedia.org/wiki/Theory_of_tides#Tidal_constituents. Accessed on: 2022-03-07.
+# unit: period (hour), speed (deg per hour)
+Tag = collections.namedtuple('Tag', 'species symbol period speed doodson_num noaa_order')
+TIDES = (
+    # Semi-diurnal
+    Tag('Principal lunar semidiurnal'              , r'$M_2$'      , 12.4206012 , 28.9841042, 255.555, 1 ),
+    Tag('Principal solar semidiurnal'              , r'$S_2$'      , 12.0       , 30.0      , 273.555, 2 ),
+    Tag('Larger lunar elliptic semidiurnal'        , r'$N_2$'      , 12.65834751, 28.4397295, 245.655, 3 ),
+    Tag('Larger lunar evectional'                  , r'$v_2$'      , 12.62600509, 28.5125831, 247.455, 11),
+    Tag('Variational'                              , r'$\mu_2$'    , 12.8717576 , 27.9682084, 237.555, 13),
+    Tag('Lunar elliptical semidiurnal second-order', '2"N'+r'$_2$' , 12.90537297, 27.8953548, 235.755, 14),
+    Tag('Smaller lunar evectional'                 , r'$\lambda_2$', 12.22177348, 29.4556253, 263.655, 16),
+    Tag('Larger solar elliptic'                    , r'$T_2$'      , 12.01644934, 29.9589333, 272.555, 27),
+    Tag('Smaller solar elliptic'                   , r'$R_2$'      , 11.98359564, 30.0410667, 274.555, 28),
+    Tag('Shallow water semidiurnal'                , r'$2SM_2$'    , 11.60695157, 31.0158958, 291.555, 31),
+    Tag('Smaller lunar elliptic semidiurnal'       , r'$L_2$'      , 12.19162085, 29.5284789, 265.455, 33),
+    Tag('Lunisolar semidiurnal'                    , r'$K_2$'      , 11.96723606, 30.0821373, 275.555, 35),
+
+    # Diurnal
+    Tag('Lunar diurnal'                  , r'$K_1$' , 23.93447213, 15.0410686, 165.555, 4 ),
+    Tag('Lunar diurnal'                  , r'$O_1$' , 25.81933871, 13.9430356, 145.555, 6 ),
+    Tag('Lunar diurnal'                  , r'$OO_1$', 22.30608083, 16.1391017, 185.555, 15),
+    Tag('Solar diurnal'                  , r'$S_1$' , 24.0       , 15.0      , 164.555, 17),
+    Tag('Smaller lunar elliptic diurnal' , r'$M_1$' , 24.84120241, 14.4920521, 155.555, 18),
+    Tag('Smaller lunar elliptic diurnal' , r'$J_1$' , 23.09848146, 15.5854433, 175.455, 19),
+    Tag('Larger lunar evectional diurnal', r'$\rho$', 26.72305326, 13.4715145, 137.455, 25),
+    Tag('Larger lunar elliptic diurnal'  , r'$Q_1$' , 26.868350  , 13.3986609, 135.655, 26),
+    Tag('Larger elliptic diurnal'        , r'$2Q_1$', 28.00621204, 12.8542862, 125.755, 29),
+    Tag('Solar diurnal'                  , r'$P_1$' , 24.06588766, 14.9589314, 163.555, 30),
+
+    # Long period
+    Tag('Lunar monthly'                  , r'$M_m$'   ,  661.3111655, 0.5443747, 65.455, 20),  # period  27.554631896 days
+    Tag('Solar semiannual'               , r'$S_{sa}$', 4383.076325 , 0.0821373, 57.555, 21),  # period 182.628180208 days
+    Tag('Solar annual'                   , r'$S_a$'   , 8766.15265  , 0.0410686, 56.555, 22),  # period 365.256360417 days
+    Tag('Lunisolar synodic fortnightly'  , r'$MS_f$'  ,  354.3670666, 1.0158958, 73.555, 23),  # period  14.765294442 days
+    Tag('Lunisolar fortnightly'          , r'$M_f$'   ,  327.8599387, 1.0980331, 75.555, 24),  # period  13.660830779 days
+
+    # Short period
+    Tag('Shallow water overtides of principal lunar', r'$M_4$'      , 6.210300601,  57.9682084, 455.555, 5 ),
+    Tag('Shallow water overtides of principal lunar', r'$M_6$'      , 4.140200401,  86.9523127, 655.555, 7 ),
+    Tag('Shallow water terdiurnal'                  , r'$MK_3$'     , 8.177140247,  44.0251729, 365.555, 8 ),
+    Tag('Shallow water overtides of principal solar', r'$S_4$'      , 6.0        ,  60.0      , 491.555, 9 ),
+    Tag('Shallow water quarter diurnal'             , r'$MN_4$'     , 6.269173724,  57.4238337, 445.655, 10),
+    Tag('Shallow water overtides of principal solar', r'$S_6$'      , 4.0        ,  90.0      , np.nan , 12),
+    Tag('Lunar terdiurnal'                          , r'$M_3$'      , 8.280400802,  43.4761563, 355.555, 32),
+    Tag('Shallow water terdiurnal'                  , '2"MK'+r'$_3$', 8.38630265 ,  42.9271398, 345.555, 34),
+    Tag('Shallow water eighth diurnal'              , r'$M_8$'      , 3.105150301, 115.9364166, 855.555, 36),
+    Tag('Shallow water quarter diurnal'             , r'$MS_4$'     , 6.103339275,  58.9841042, 473.555, 37),
+)
+
+
+##################################  Earth tides - point mode  ##################################
+def calc_solid_earth_tides_point(lat, lon, dt0, dt1, step_sec=60, display=False, verbose=True):
+    """Calculate SET in east/north/up direction for the given time period at the given point (lat/lon).
+
+    Parameters: lat/lon  - float32, latitude/longitude of the point of interest
+                dt0/1    - datetime.datetime object, start/end date and time
+                step_sec - int16, time step in seconds
+                display  - bool, plot the calculated SET
+                verbose  - bool, print verbose message
+    Returns:    dt_out   - 1D np.ndarray in dt.datetime objects
+                tide_e   - 1D np.ndarray in float32, SET in east  direction in meters
+                tide_n   - 1D np.ndarray in float32, SET in north direction in meters
+                tide_u   - 1D np.ndarray in float32, SET in up    direction in meters
+    Examples:   dt0 = dt.datetime(2020,11,1,4,0,0)
+                dt1 = dt.datetime(2020,12,31,2,0,0)
+                (dt_out,
+                 tide_e,
+                 tide_n,
+                 tide_u) = calc_solid_earth_tides_point(34.0, -118.0, dt0, dt1)
+    """
+
+    print('PYSOLID: calculate solid Earth tides in east/north/up direction')
+    print(f'PYSOLID: lot/lon: {lat}/{lon} degree')
+    print(f'PYSOLID: start UTC: {dt0.isoformat()}')
+    print(f'PYSOLID: end   UTC: {dt1.isoformat()}')
+    print(f'PYSOLID: time step: {step_sec} seconds')
+
+    dt_out = []
+    tide_e = []
+    tide_n = []
+    tide_u = []
+
+    ndays = (dt1.date() - dt0.date()).days + 1
+    for i in range(ndays):
+        di = dt0.date() + dt.timedelta(days=i)
+        if verbose:
+            print(f'SOLID  : {di.isoformat()} {i+1}/{ndays} ...')
+
+        # calc tide_u/n/u for the whole day
+        (dt_outi,
+         tide_ei,
+         tide_ni,
+         tide_ui) = calc_solid_earth_tides_point_per_day(lat, lon,
+                                                         date_str=di.strftime('%Y%m%d'),
+                                                         step_sec=int(step_sec))
+
+        # flag to mark the first/last datetime
+        if i == 0:
+            flag = dt_outi >= dt0
+        elif i == ndays - 1:
+            flag = dt_outi <= dt1
+        else:
+            flag = np.ones(dt_outi.size, dtype=np.bool_)
+
+        # update
+        dt_out += dt_outi[flag].tolist()
+        tide_e += tide_ei[flag].tolist()
+        tide_n += tide_ni[flag].tolist()
+        tide_u += tide_ui[flag].tolist()
+
+    # list --> np.ndarray for output
+    dt_out = np.array(dt_out)
+    tide_e = np.array(tide_e)
+    tide_n = np.array(tide_n)
+    tide_u = np.array(tide_u)
+
+    # plot
+    if display:
+        plot_solid_earth_tides_point(dt_out, tide_e, tide_n, tide_u, lalo=[lat, lon])
+
+    return dt_out, tide_e, tide_n, tide_u
+
+
+def calc_solid_earth_tides_point_per_day(lat, lon, date_str, step_sec=60):
+    """Calculate solid Earth tides (SET) in east/north/up direction
+    for one day at the given point (lat/lon).
+
+    Parameters: lat/lon  - float32, latitude/longitude of the point of interest
+                date_str - str, date in YYYYMMDD
+                step_sec - int16, time step in seconds
+    Returns:    dt_out   - 1D np.ndarray in dt.datetime objects
+                tide_e   - 1D np.ndarray in float32, SET in east  direction in meters
+                tide_n   - 1D np.ndarray in float32, SET in north direction in meters
+                tide_u   - 1D np.ndarray in float32, SET in up    direction in meters
+    Examples:   (dt_out,
+                 tide_e,
+                 tide_n,
+                 tide_u) = calc_solid_earth_tides_point_per_day(34.0, -118.0, '20180219')
+    """
+    try:
+        from .solid import solid_point, LLH
+    except ImportError:
+        msg = "Cannot import name 'solid' from 'pysolid'!"
+        msg += '\n    Maybe solid.for is NOT compiled yet.'
+        msg += '\n    Check instruction at: https://github.com/insarlab/PySolid.'
+        raise ImportError(msg)
+
+    # calc solid Earth tides
+    t = dt.datetime.strptime(date_str, '%Y%m%d')
+    # secs, tide_e, tide_n, tide_u  = solid_point(
+    #     lat, lon, t.year, t.month, t.day, step_sec
+    # )
+    secs, tides = solid_point(
+        LLH(lat, lon, 0.), t.date(), step_sec
+    )
+    tide_e = tides[:, 0]
+    tide_n = tides[:, 1]
+    tide_u = tides[:, 2]
+
+    dt_out = [t + dt.timedelta(seconds=sec) for sec in secs]
+    dt_out = np.array(dt_out)
+
+    return dt_out, tide_e, tide_n, tide_u
+
+
+#########################################  Plot  ###############################################
+def plot_solid_earth_tides_point(dt_out, tide_e, tide_n, tide_u, lalo=None,
+                                 out_fig=None, save=False, display=True):
+    """Plot the solid Earth tides at one point."""
+    from matplotlib import pyplot as plt, dates as mdates
+
+    # plot
+    fig, axs = plt.subplots(nrows=3, ncols=1, figsize=[6, 4], sharex=True)
+    for ax, data, label in zip(axs.flatten(),
+                               [tide_e, tide_n, tide_u],
+                               ['East [cm]', 'North [cm]', 'Up [cm]']):
+        ax.plot(dt_out, data*100)
+        ax.set_ylabel(label, fontsize=12)
+
+    # axis format
+    for ax in axs:
+        ax.tick_params(which='both', direction='out', bottom=True, top=False, left=True, right=True)
+        ax.xaxis.set_major_locator(mdates.MonthLocator())
+        ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
+        ax.xaxis.set_minor_locator(mdates.DayLocator())
+
+    if lalo:
+        axs[0].set_title('solid Earth tides at (N{}, E{})'.format(lalo[0], lalo[1]), fontsize=12)
+    fig.tight_layout()
+    fig.align_ylabels()
+
+    # output
+    if out_fig:
+        save = True
+
+    if save:
+        if not out_fig:
+            out_fig = os.path.abspath('SET_point.png')
+        print('save figure to {}'.format(out_fig))
+        fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)
+
+    if display:
+        print('showing...')
+        plt.show()
+    else:
+        plt.close()
+
+    return
+
+
+def plot_power_spectral_density4tides(tide_ts, sample_spacing, out_fig=None, fig_dpi=300, min_psd=1500):
+    """Plot the power spectral density (PSD) of tides time-series.
+    Note: for accurate PSD analysis, a long time-series, e.g. one year, is recommended.
+    """
+    from matplotlib import pyplot as plt, ticker
+    from scipy import signal
+
+    ## calc PSD
+    freq, psd = signal.periodogram(tide_ts, fs=1/sample_spacing, scaling='density')
+    # get rid of zero in the first element
+    psd = psd[1:]
+    freq = freq[1:]
+    period = 1./3600./freq   # period (hour)
+
+    ## plot
+    fig, axs = plt.subplots(nrows=1, ncols=2, figsize=[8,3], sharey=True)
+    for ax in axs:
+        ax.plot(period, psd, '.', ms=16, mfc='none', mew=2)
+
+    # axis format
+    for ax in axs:
+        ax.tick_params(which='both', direction='out', bottom=True, top=True, left=True, right=True)
+        ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
+        ax.set_xlabel('Period [hour]')
+    axs[0].set_xlim(11.3, 13.2)
+    axs[1].set_xlim(22.0, 28.0)
+    ax = axs[0]
+    ax.set_ylabel('Power Spectral Density\n'+r'[$m^2/Hz$]')
+    ax.set_ylim(0, ymax=axs[0].get_ylim()[1] * 1.1)
+    #ax.set_ylim(1e-1,5e6);  plt.yscale('log')
+    #ax.yaxis.set_major_locator(ticker.FixedLocator([0,50e3,100e3,150e3]))
+    #ax.set_yticklabels(['0','50k','100k','150k'])
+    fig.tight_layout()
+
+    # Tidal constituents
+    for ax in axs:
+        add_tidal_constituents(ax, period, psd, min_psd=min_psd)
+    axs[0].annotate('semi-diurnal', xy=(0.05,0.85), xycoords='axes fraction')
+    axs[1].annotate('diurnal',      xy=(0.05,0.85), xycoords='axes fraction')
+
+    # output
+    if out_fig:
+        print('save figure to file:', os.path.abspath(out_fig))
+        fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=fig_dpi)
+    plt.show()
+
+    return
+
+
+def add_tidal_constituents(ax, period, psd, min_psd=1500, verbose=False):
+    """Mark tidal constituents covered in the axes."""
+    pmin, pmax = ax.get_xlim()
+    for tide in TIDES:
+        if pmin <= tide.period <= pmax:
+            tide_psd = psd[np.argmin(np.abs(period - tide.period))]
+            if tide_psd >= min_psd:
+                ymax = tide_psd / ax.get_ylim()[1]
+                ax.axvline(x=tide.period, ymax=ymax, color='k', ls='--', lw=1)
+                ax.annotate(tide.symbol, xy=(tide.period, tide_psd), xytext=(6, 4), textcoords='offset points')
+                if verbose:
+                    print('tide: speices={}, symbol={}, period={} hours, psd={} m^2/Hz'.format(
+                        tide.species, tide.symbol, tide.period, tide_psd))
+    return

tests/grid_pyimpl.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+# Author: Zhang Yunjun, Jan 2021
+# Copyright 2020, by the California Institute of Technology.
+
+
+import os
+import sys
+import datetime as dt
+
+import numpy as np
+
+from pysolid import py_solid
+
+
+if __name__ == '__main__':
+
+    # print the file/module path
+    print('-'*50)
+    print(os.path.abspath(__file__))
+
+    # prepare inputs
+    dt_obj = dt.datetime(2020, 12, 25, 14, 7, 44)
+    atr = {
+        'LENGTH'  : 400,
+        'WIDTH'   : 500,
+        'X_FIRST' : -118.2,
+        'Y_FIRST' : 33.8,
+        'X_STEP'  :  0.000833333,
+        'Y_STEP'  : -0.000833333,
+    }
+
+    # reference
+    # calculated based on version 0.3.2.post6 on Jun 24, 2024
+    # env: macOS with python-3.10, numpy-1.24
+    # install: manual compilation via f2py
+    tide_e_80_100 = np.array(
+        [[0.01628786, 0.01630887, 0.01633078, 0.01635247, 0.01637394],
+         [0.01633248, 0.01635348, 0.01637538, 0.01639706, 0.01641851],
+         [0.01638009, 0.01640107, 0.01642296, 0.01644462, 0.01646606],
+         [0.01642767, 0.01644864, 0.01647052, 0.01649217, 0.01651359],
+         [0.01647523, 0.01649619, 0.01651805, 0.01653968, 0.01656109]],
+    )
+    tide_n_80_100 = np.array(
+        [[-0.02406203, -0.02412341, -0.02418807, -0.02425273, -0.0243174 ],
+         [-0.02407558, -0.02413699, -0.02420168, -0.02426637, -0.02433107],
+         [-0.02408992, -0.02415136, -0.02421608, -0.02428081, -0.02434554],
+         [-0.02410413, -0.0241656 , -0.02423036, -0.02429511, -0.02435988],
+         [-0.02411821, -0.02417972, -0.0242445 , -0.02430929, -0.02437408]],
+    )
+    tide_u_80_100 = np.array(
+        [[-0.05548462, -0.05533455, -0.05517631, -0.05501789, -0.05485928],
+         [-0.05529561, -0.0551451 , -0.05498639, -0.0548275 , -0.05466843],
+         [-0.05509374, -0.05494276, -0.05478355, -0.05462417, -0.05446461],
+         [-0.05489176, -0.05474031, -0.05458061, -0.05442073, -0.05426067],
+         [-0.05468968, -0.05453776, -0.05437757, -0.05421719, -0.05405664]],
+    )
+
+    # calculate
+    (tide_e,
+     tide_n,
+     tide_u) = py_solid.calc_solid_earth_tides_grid(dt_obj, atr, verbose=True)
+
+    # compare
+    assert np.allclose(tide_e[::80, ::100], tide_e_80_100)
+    assert np.allclose(tide_n[::80, ::100], tide_n_80_100)
+    assert np.allclose(tide_u[::80, ::100], tide_u_80_100)
+
+    # plot
+    out_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'pic'))
+    os.makedirs(out_dir, exist_ok=True)
+
+    out_fig = os.path.join(out_dir, 'grid_pyimpl.png')
+    py_solid.plot_solid_earth_tides_grid(
+        tide_e, tide_n, tide_u, dt_obj,
+        out_fig=out_fig,
+        display=False)
+
+    # open the plotted figures
+    if sys.platform in ['linux']:
+        os.system(f'display {out_fig}')
+    elif sys.platform in ['darwin']:
+        os.system(f'open {out_fig}')
+    elif sys.platform.startswith('win'):
+        os.system(out_fig)
+    else:
+        print(f'Unknown OS system ({sys.platform}). Check results in file: {out_fig}.')

tests/point_pyimpl.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+# Author: Zhang Yunjun, Jan 2021
+# Copyright 2020, by the California Institute of Technology.
+
+
+import os
+import sys
+import datetime as dt
+
+import numpy as np
+
+from pysolid import py_solid
+
+
+if __name__ == '__main__':
+
+    # print the file/module path
+    print('-'*50)
+    print(os.path.abspath(__file__))
+
+    # prepare inputs
+    lat, lon = 34.0, -118.0     # Los Angles, CA
+    dt_obj0 = dt.datetime(2020, 11,  5, 12, 0, 0)
+    dt_obj1 = dt.datetime(2020, 12, 31,  0, 0, 0)
+
+    # reference
+    # calculated based on version 0.3.2.post6 on Jun 24, 2024
+    # env: macOS with python-3.10, numpy-1.24
+    # install: manual compilation via f2py
+    dt_out_8000 = np.array(
+        [dt.datetime(2020, 11,  5, 12,  0),
+         dt.datetime(2020, 11, 11,  1, 20),
+         dt.datetime(2020, 11, 16, 14, 40),
+         dt.datetime(2020, 11, 22,  4,  0),
+         dt.datetime(2020, 11, 27, 17, 20),
+         dt.datetime(2020, 12,  3,  6, 40),
+         dt.datetime(2020, 12,  8, 20,  0),
+         dt.datetime(2020, 12, 14,  9, 20),
+         dt.datetime(2020, 12, 19, 22, 40),
+         dt.datetime(2020, 12, 25, 12,  0)], dtype=object,
+    )
+    tide_e_8000 = np.array(
+        [-0.02975027,  0.04146837, -0.02690945, -0.00019223,  0.01624152,
+          0.0532655 , -0.02140918, -0.05554432,  0.01371739, -0.00516968],
+    )
+    tide_n_8000 = np.array(
+        [-0.01275229, -0.02834036,  0.00886857, -0.03247227, -0.05237735,
+         -0.00590791, -0.01990448, -0.01964124, -0.04439581, -0.00410378],
+    )
+    tide_u_8000 = np.array(
+        [ 0.16008235, -0.05721991, -0.15654693, -0.00041214,  0.03041098,
+          0.13082217, -0.1006462 ,  0.24870719, -0.02648802, -0.08420228],
+    )
+
+    # calculate
+    (dt_out,
+     tide_e,
+     tide_n,
+     tide_u) = py_solid.calc_solid_earth_tides_point(lat, lon, dt_obj0, dt_obj1, verbose=False)
+
+    # compare
+    assert all(dt_out[::8000] == dt_out_8000)
+    assert np.allclose(tide_e[::8000], tide_e_8000)
+    assert np.allclose(tide_n[::8000], tide_n_8000)
+    assert np.allclose(tide_u[::8000], tide_u_8000)
+
+    # plot
+    out_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'pic'))
+    os.makedirs(out_dir, exist_ok=True)
+
+    out_fig = os.path.join(out_dir, 'point_pyimpl.png')
+    py_solid.plot_solid_earth_tides_point(
+        dt_out, tide_e, tide_n, tide_u,
+        lalo=[lat, lon],
+        out_fig=out_fig,
+        display=False)
+
+    # open the saved figure
+    if sys.platform in ['linux']:
+        os.system(f'display {out_fig}')
+    elif sys.platform in ['darwin']:
+        os.system(f'open {out_fig}')
+    elif sys.platform.startswith('win'):
+        os.system(out_fig)
+    else:
+        print(f'Unknown OS system ({sys.platform}). Check results in file: {out_fig}.')

tests/test_grid.py

@@ -0,0 +1,107 @@
+"""Test suite for grid-based solid earth tide calculations."""
+
+import datetime as dt
+import numpy as np
+import pytest
+from pysolid import py_solid
+import pysolid
+
+
+@pytest.mark.parametrize(
+    "module",
+    [
+        pytest.param(py_solid, id="python_impl"),
+        pytest.param(pysolid, id="fortran_impl"),
+    ],
+)
+def test_grid_calculation(module):
+    """Test grid-based solid earth tide calculations."""
+    # Test inputs
+    dt_obj = dt.datetime(2020, 12, 25, 14, 7, 44)
+    atr = {
+        "LENGTH": 400,
+        "WIDTH": 500,
+        "X_FIRST": -118.2,
+        "Y_FIRST": 33.8,
+        "X_STEP": 0.000833333,
+        "Y_STEP": -0.000833333,
+    }
+
+    # Reference outputs
+    ref_tide_e = np.array(
+        [
+            [0.01628786, 0.01630887, 0.01633078, 0.01635247, 0.01637394],
+            [0.01633248, 0.01635348, 0.01637538, 0.01639706, 0.01641851],
+            [0.01638009, 0.01640107, 0.01642296, 0.01644462, 0.01646606],
+            [0.01642767, 0.01644864, 0.01647052, 0.01649217, 0.01651359],
+            [0.01647523, 0.01649619, 0.01651805, 0.01653968, 0.01656109],
+        ]
+    )
+
+    ref_tide_n = np.array(
+        [
+            [-0.02406203, -0.02412341, -0.02418807, -0.02425273, -0.0243174],
+            [-0.02407558, -0.02413699, -0.02420168, -0.02426637, -0.02433107],
+            [-0.02408992, -0.02415136, -0.02421608, -0.02428081, -0.02434554],
+            [-0.02410413, -0.0241656, -0.02423036, -0.02429511, -0.02435988],
+            [-0.02411821, -0.02417972, -0.0242445, -0.02430929, -0.02437408],
+        ]
+    )
+
+    ref_tide_u = np.array(
+        [
+            [-0.05548462, -0.05533455, -0.05517631, -0.05501789, -0.05485928],
+            [-0.05529561, -0.0551451, -0.05498639, -0.0548275, -0.05466843],
+            [-0.05509374, -0.05494276, -0.05478355, -0.05462417, -0.05446461],
+            [-0.05489176, -0.05474031, -0.05458061, -0.05442073, -0.05426067],
+            [-0.05468968, -0.05453776, -0.05437757, -0.05421719, -0.05405664],
+        ]
+    )
+
+    # Run calculation
+    tide_e, tide_n, tide_u = module.calc_solid_earth_tides_grid(
+        dt_obj, atr, verbose=False
+    )
+
+    # Check results (using subsampled grid)
+    np.testing.assert_array_almost_equal(tide_e[::80, ::100], ref_tide_e)
+    np.testing.assert_array_almost_equal(tide_n[::80, ::100], ref_tide_n)
+    np.testing.assert_array_almost_equal(tide_u[::80, ::100], ref_tide_u)
+
+
+def test_grid_input_validation():
+    """Test input validation for grid calculations."""
+    # Test with missing attributes
+    with pytest.raises(KeyError):
+        invalid_atr = {"LENGTH": 400}  # missing required attributes
+        py_solid.calc_solid_earth_tides_grid(
+            dt.datetime.now(), invalid_atr, verbose=False
+        )
+
+    # Test with invalid grid dimensions
+    with pytest.raises(ValueError):
+        # negative length
+        invalid_atr = {
+            "LENGTH": -1,
+            "WIDTH": 500,
+            "X_FIRST": -118.2,
+            "Y_FIRST": 33.8,
+            "X_STEP": 0.000833333,
+            "Y_STEP": -0.000833333,
+        }
+        py_solid.calc_solid_earth_tides_grid(
+            dt.datetime.now(), invalid_atr, verbose=False
+        )
+
+    # Test with invalid coordinates
+    # invalid longitude
+    with pytest.raises(ValueError):
+        invalid_atr = {
+            "LENGTH": 400,
+            "WIDTH": 500,
+            "X_FIRST": -500,
+            "Y_FIRST": 33.8,
+            "X_STEP": 0.000833333,
+            "Y_STEP": -0.000833333,
+        }
+        py_solid.calc_solid_earth_tides_grid(dt.datetime.now(), invalid_atr)

tests/test_point.py

@@ -0,0 +1,125 @@
+"""Test suite for point-based solid earth tide calculations."""
+
+import datetime as dt
+import numpy as np
+import pytest
+from pysolid import py_solid
+import pysolid
+
+
+@pytest.mark.parametrize(
+    "module",
+    [
+        pytest.param(py_solid, id="python_impl"),
+        pytest.param(pysolid, id="fortran_impl"),
+    ],
+)
+def test_point_calculation(module):
+    """Test point-based solid earth tide calculations."""
+    # Test inputs
+    lat, lon = 34.0, -118.0
+    dt_obj0 = dt.datetime(2020, 11, 5, 12, 0, 0)
+    dt_obj1 = dt.datetime(2020, 12, 31, 0, 0, 0)
+
+    # Reference outputs
+    ref_dates = np.array(
+        [
+            dt.datetime(2020, 11, 5, 12, 0),
+            dt.datetime(2020, 11, 11, 1, 20),
+            dt.datetime(2020, 11, 16, 14, 40),
+            dt.datetime(2020, 11, 22, 4, 0),
+            dt.datetime(2020, 11, 27, 17, 20),
+            dt.datetime(2020, 12, 3, 6, 40),
+            dt.datetime(2020, 12, 8, 20, 0),
+            dt.datetime(2020, 12, 14, 9, 20),
+            dt.datetime(2020, 12, 19, 22, 40),
+            dt.datetime(2020, 12, 25, 12, 0),
+        ],
+        dtype=object,
+    )
+
+    ref_tide_e = np.array(
+        [
+            -0.02975027,
+            0.04146837,
+            -0.02690945,
+            -0.00019223,
+            0.01624152,
+            0.0532655,
+            -0.02140918,
+            -0.05554432,
+            0.01371739,
+            -0.00516968,
+        ]
+    )
+
+    ref_tide_n = np.array(
+        [
+            -0.01275229,
+            -0.02834036,
+            0.00886857,
+            -0.03247227,
+            -0.05237735,
+            -0.00590791,
+            -0.01990448,
+            -0.01964124,
+            -0.04439581,
+            -0.00410378,
+        ]
+    )
+
+    ref_tide_u = np.array(
+        [
+            0.16008235,
+            -0.05721991,
+            -0.15654693,
+            -0.00041214,
+            0.03041098,
+            0.13082217,
+            -0.1006462,
+            0.24870719,
+            -0.02648802,
+            -0.08420228,
+        ]
+    )
+
+    # Run calculation
+    dt_out, tide_e, tide_n, tide_u = module.calc_solid_earth_tides_point(
+        lat, lon, dt_obj0, dt_obj1, verbose=False
+    )
+
+    # Check results
+    assert all(d1 == d2 for (d1, d2) in zip(dt_out[::8000], ref_dates))
+    np.testing.assert_array_almost_equal(tide_e[::8000], ref_tide_e)
+    np.testing.assert_array_almost_equal(tide_n[::8000], ref_tide_n)
+    np.testing.assert_array_almost_equal(tide_u[::8000], ref_tide_u)
+
+
+@pytest.mark.parametrize(
+    "module",
+    [
+        pytest.param(py_solid, id="python_impl"),
+        pytest.param(pysolid, id="fortran_impl"),
+    ],
+)
+def test_point_input_validation(module):
+    """Test input validation for point calculations."""
+    with pytest.raises(ValueError):
+        # Invalid latitude
+        module.calc_solid_earth_tides_point(
+            91.0,
+            -118.0,
+            dt.datetime.now(),
+            dt.datetime.now() + dt.timedelta(days=1),
+            verbose=False,
+        )
+
+    with pytest.raises(ValueError):
+        # invalid longitude
+        module.calc_solid_earth_tides_point(
+            34.0,
+            -500.0,
+            dt.datetime.now(),
+            dt.datetime.now() + dt.timedelta(days=1),
+            verbose=False,
+        )