Merge pull request #19 from tsutterley/perth3

add GOT4.7 test program

.github/workflows/python-package.yml

@@ -55,4 +55,7 @@ jobs:
         pip install --no-deps .
         git clone ${{ secrets.TMD_MATLAB_TOOLBOX }}
         pytest --username=${{ secrets.EARTHDATA_USERNAME }} \
-            --password=${{ secrets.EARTHDATA_PASSWORD }}
+            --password=${{ secrets.EARTHDATA_PASSWORD }} \
+            --aws-access=${{ secrets.AWS_ACCESS_KEY_ID }} \
+            --aws-secret=${{ secrets.AWS_SECRET_KEY }} \
+            --aws-region=${{ secrets.AWS_REGION_NAME }}

.github/workflows/python-request.yml

@@ -52,4 +52,7 @@ jobs:
         pip install --no-deps .
         git clone ${{ secrets.TMD_MATLAB_TOOLBOX }}
         pytest --username=${{ secrets.EARTHDATA_USERNAME }} \
-            --password=${{ secrets.EARTHDATA_PASSWORD }}
+            --password=${{ secrets.EARTHDATA_PASSWORD }} \
+            --aws-access=${{ secrets.AWS_ACCESS_KEY_ID }} \
+            --aws-secret=${{ secrets.AWS_SECRET_KEY }} \
+            --aws-region=${{ secrets.AWS_REGION_NAME }}

doc/source/conf.py

@@ -22,7 +22,7 @@
 author = 'Tyler C. Sutterley'
 
 # The full version, including alpha/beta/rc tags
-release = '1.0.2.8'
+release = '1.0.2.9'
 
 
 # -- General configuration ---------------------------------------------------

doc/source/index.rst

@@ -40,6 +40,8 @@ calculating radial pole tide displacements
     user_guide/calc_astrol_longitudes.md
     user_guide/calc_delta_time.md
     user_guide/calc_iers_mean_pole.md
+    user_guide/compute_equilibrium_tide.md
+    user_guide/compute_LPET_icebridge_data.md
     user_guide/compute_LPT_displacements.md
     user_guide/compute_LPT_icebridge_data.md
     user_guide/compute_OPT_displacements.md

doc/source/user_guide/calc_astrol_longitudes.md

@@ -8,12 +8,12 @@ calc_astrol_longitudes.py
 #### Calling Sequence
 ```python
 from pyTMD.calc_astrol_longitudes import calc_astrol_longitudes
-s,h,p,N,PP = calc_astrol_longitudes(time, ASTRO5=True)
+s,h,p,N,PP = calc_astrol_longitudes(MJD, ASTRO5=True)
 ```
 [Source code](https://github.com/tsutterley/pyTMD/blob/master/pyTMD/calc_astrol_longitudes.py)
 
 #### Inputs
- 1. `time`: modified Julian day of input date
+ 1. `MJD`: Modified Julian Day of input date
 
 #### Options
  - `MEEUS`: use additional coefficients from Meeus Astronomical Algorithms

doc/source/user_guide/compute_LPET_icebridge_data.md

@@ -0,0 +1,18 @@
+compute_LPET_icebridge_data.py
+=============================
+
+ - Calculates long-period equilibrium tides for correcting Operation IceBridge elevation data
+ - Uses the summation of fifteen tidal spectral lines from [Cartwright and Edden, (1973)](https://doi.org/10.1111/j.1365-246X.1973.tb03420.x)
+
+#### Calling Sequence
+```bash
+python compute_LPET_icebridge_data.py --verbose input_file
+```
+[Source code](https://github.com/tsutterley/pyTMD/blob/master/compute_LPET_icebridge_data.py)
+
+#### Inputs
+ 1. `input_file`: input ATM1B, ATM icessn or LVIS file from NSIDC
+
+#### Command Line Options
+ - `-M X`, `--mode=X`: Permission mode of output file
+ - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_equilibrium_tide.md

@@ -0,0 +1,19 @@
+compute_equilibrium_tide.py
+===========================
+
+ - Calculates the long-period equilibrium ocean tides
+ - Can be used to calculate tidal corrections for imagery  
+
+#### Calling Sequence
+```python
+from pyTMD.compute_equilibrium_tide import compute_equilibrium_tide
+lpet = compute_equilibrium_tide(time,lat)
+```
+[Source code](https://github.com/tsutterley/pyTMD/blob/master/pyTMD/compute_equilibrium_tide.py)
+
+#### Inputs
+ 1. `t`: days relative to Jan 1, 1992 (48622mjd)
+ 2. `lat`: latitudes in degrees
+
+#### Outputs
+ - `lpet`: long-period equilibrium tide values

doc/source/user_guide/infer_minor_corrections.md

@@ -7,13 +7,13 @@ infer_minor_corrections.py
 #### Calling Sequence
 ```python
 from pyTMD.infer_minor_corrections import infer_minor_corrections
-dh = infer_minor_corrections(time, zmajor, constituents,
+dh = infer_minor_corrections(t, zmajor, constituents,
     DELTAT=DELTAT, CORRECTIONS=CORRECTIONS)
 ```
 [Source code](https://github.com/tsutterley/pyTMD/blob/master/pyTMD/infer_minor_corrections.py)
 
 #### Inputs
- 1. `time`: days relative to Jan 1, 1992 (48622mjd)
+ 1. `t`: days relative to Jan 1, 1992 (48622mjd)
  2. `zmajor`: Complex oscillations for given constituents/points
  3. `constituents`: tidal constituent IDs
 

doc/source/user_guide/load_nodal_corrections.md

@@ -7,12 +7,12 @@ load_nodal_corrections.py
 #### Calling Sequence
 ```python
 from pyTMD.load_nodal_corrections import load_nodal_corrections
-pu,pf,G = load_nodal_corrections(time,constituents)
+pu,pf,G = load_nodal_corrections(MJD,constituents)
 ```
 [Source code](https://github.com/tsutterley/pyTMD/blob/master/pyTMD/load_nodal_corrections.py)
 
 #### Inputs
- 1. `time`: modified julian day of input date
+ 1. `MJD`: Modified Julian Day of input date
  2. `zmajor`: Complex oscillations for given constituents/points
  3. `constituents`: tidal constituent IDs
 

Plot Antarctic Tidal Currents.ipynb → notebooks/Plot Antarctic Tidal Currents.ipynb

@@ -308,8 +308,10 @@
     "    elif (model_format == 'FES'):\n",
     "        amp,ph = extract_FES_constants(lon, lat, model_directory[TYPE], model_files,\n",
     "            TYPE=TYPE, VERSION=MODEL, METHOD=METHOD, SCALE=SCALE)\n",
-    "        DELTAT = np.zeros_like(tide_time)\n",
-    "\n",
+    "        #-- interpolate delta times from calendar dates to tide time\n",
+    "        delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])\n",
+    "        DELTAT = calc_delta_time(delta_file, tide_time)\n",
+    "        \n",
     "    #-- calculate complex phase in radians for Euler's\n",
     "    cph = -1j*ph*np.pi/180.0\n",
     "    #-- calculate constituent oscillation\n",

Plot Ross Ice Shelf Map.ipynb → notebooks/Plot Ross Ice Shelf Map.ipynb

@@ -310,8 +310,10 @@
     "elif (model_format == 'FES'):\n",
     "    amp,ph = extract_FES_constants(lon, lat, model_directory, model_files,\n",
     "        TYPE=TYPE, VERSION=MODEL, METHOD='spline', SCALE=SCALE)\n",
-    "    DELTAT = np.zeros_like(tide_time)\n",
-    "            \n",
+    "    #-- interpolate delta times from calendar dates to tide time\n",
+    "    delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])\n",
+    "    DELTAT = calc_delta_time(delta_file, tide_time)\n",
+    "    \n",
     "#-- calculate complex phase in radians for Euler's\n",
     "cph = -1j*ph*np.pi/180.0\n",
     "#-- calculate constituent oscillation\n",

Plot Tide Forecasts.ipynb → notebooks/Plot Tide Forecasts.ipynb

@@ -336,7 +336,9 @@
     "    amp,ph = extract_FES_constants(np.array([LON]), np.array([LAT]),\n",
     "        model_directory, model_files, TYPE=TYPE, VERSION=TIDE_MODEL,\n",
     "        METHOD='spline', SCALE=SCALE)\n",
-    "    deltat = np.zeros_like(tide_time)\n",
+    "    #-- interpolate delta times from calendar dates to tide time\n",
+    "    delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])\n",
+    "    deltat = calc_delta_time(delta_file, tide_time)\n",
     "    \n",
     "#-- calculate complex phase in radians for Euler's\n",
     "cph = -1j*ph*np.pi/180.0\n",

pyTMD/__init__.py

@@ -18,6 +18,7 @@
 from pyTMD.convert_calendar_decimal import convert_calendar_decimal
 from pyTMD.calc_delta_time import calc_delta_time
 from pyTMD.calc_astrol_longitudes import calc_astrol_longitudes
+from pyTMD.compute_equilibrium_tide import compute_equilibrium_tide
 from pyTMD.convert_ll_xy import convert_ll_xy
 from pyTMD.infer_minor_corrections import infer_minor_corrections
 from pyTMD.load_constituent import load_constituent

pyTMD/calc_astrol_longitudes.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 u"""
-calc_astrol_longitudes.py (07/2020)
+calc_astrol_longitudes.py (08/2020)
 Modification of ASTROL fortran subroutine by Richard Ray 03/1999
 
 Computes the basic astronomical mean longitudes: s, h, p, N and PP
@@ -10,10 +10,10 @@
 MEEUS and ASTRO5 formulae are from versions of Meeus's Astronomical Algorithms
 
 CALLING SEQUENCE:
-    s,h,p,N,PP = calc_astrol_longitudes(time, ASTRO5=True)
+    s,h,p,N,PP = calc_astrol_longitudes(MJD, ASTRO5=True)
 
 INPUTS:
-    time: modified julian day of input date
+    MJD: Modified Julian Day of input date
 
 OUTPUTS:
     s: mean longitude of moon (degrees)
@@ -35,6 +35,7 @@
     Jean Meeus, Astronomical Algorithms, 2nd edition, 1998.
 
 UPDATE HISTORY:
+    Updated 08/2020: change time variable names to not overwrite functions
     Updated 07/2020: added function docstrings
     Updated 07/2018: added option ASTRO5 to use coefficients from Richard Ray
         for use with the GSFC Global Ocean Tides (GOT) model
@@ -48,25 +49,27 @@
 import numpy as np
 
 #-- PURPOSE: calculate the sum of a polynomial function of time
-def polynomial_sum(coefficients, time):
+def polynomial_sum(coefficients, t):
     """
     Calculates the sum of a polynomial function of time
 
     Arguments
     ---------
     coefficients: leading coefficient of polynomials of increasing order
-    time: delta time in units for a given astrological longitudes calculation
+    t: delta time in units for a given astrological longitudes calculation
     """
-    return np.sum([c * (time ** i) for i,c in enumerate(coefficients)])
+    #-- convert time to array if importing a single value
+    t = np.array([t]) if (np.ndim(t) == 0) else np.copy(t)
+    return np.sum([c * (t ** i) for i,c in enumerate(coefficients)],axis=0)
 
 #-- PURPOSE: compute the basic astronomical mean longitudes
-def calc_astrol_longitudes(time, MEEUS=False, ASTRO5=False):
+def calc_astrol_longitudes(MJD, MEEUS=False, ASTRO5=False):
     """
     Computes the basic astronomical mean longitudes: s, h, p, N and PP
 
     Arguments
     ---------
-    time: Modified Julian Day (MJD) of input date
+    MJD: Modified Julian Day (MJD) of input date
 
     Keyword arguments
     -----------------
@@ -83,8 +86,8 @@ def calc_astrol_longitudes(time, MEEUS=False, ASTRO5=False):
     """
     circle = 360.0
     if MEEUS:
-        #-- convert from MJD to days relative to 2000-01-01
-        T = time - 51544.5
+        #-- convert from MJD to days relative to 2000-01-01T12:00:00
+        T = MJD - 51544.5
         #-- mean longitude of moon
         lunar_longitude = np.array([218.3164591, 13.17639647754579,
             -9.9454632e-13, 3.8086292e-20, -8.6184958e-27])
@@ -104,8 +107,8 @@ def calc_astrol_longitudes(time, MEEUS=False, ASTRO5=False):
         #-- mean longitude of solar perigee (Simon et al., 1994)
         PP = 282.94 + 1.7192 * T
     elif ASTRO5:
-        #-- convert from MJD to centuries relative to 2000-01-01
-        T = (time - 51544.5)/36525.0
+        #-- convert from MJD to centuries relative to 2000-01-01T12:00:00
+        T = (MJD - 51544.5)/36525.0
         #-- mean longitude of moon (p. 338)
         lunar_longitude = np.array([218.3164477, 481267.88123421, -1.5786e-3,
              1.855835e-6, -1.53388e-8])
@@ -124,8 +127,9 @@ def calc_astrol_longitudes(time, MEEUS=False, ASTRO5=False):
         #-- mean longitude of solar perigee (Simon et al., 1994)
         PP = 282.94 + 1.7192 * T
     else:
-        #-- convert from MJD to days relative to 2000-01-01
-        T = time - 51544.4993
+        #-- convert from MJD to days relative to 2000-01-01T12:00:00
+        #-- convert from Universal Time to Dynamic Time at 2000-01-01
+        T = MJD - 51544.4993
         #-- mean longitude of moon
         s = 218.3164 + 13.17639648 * T
         #-- mean longitude of sun

pyTMD/calc_delta_time.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 calc_delta_time.py
-Written by Tyler Sutterley (07/2020)
+Written by Tyler Sutterley (08/2020)
 Calculates the difference between dynamic time and universal time (TT - UT1)
     following Richard Ray's PERTH3 algorithms
 
@@ -12,7 +12,7 @@
     idays: input times to interpolate (days since 1992-01-01T00:00:00)
 
 OUTPUTS:
-    deltat: delta time estimates at the output times
+    deltat: delta time estimates at the output times in days
 
 PYTHON DEPENDENCIES:
     numpy: Scientific Computing Tools For Python
@@ -23,6 +23,7 @@
 
 UPDATE HISTORY:
     Updated 08/2020: using builtin time operations, interpolate with tide time
+        convert output units to be in days
     Updated 07/2020: added function docstrings. scipy interpolating splines
     Updated 11/2019: pad input time dimension if entering a single value
     Updated 07/2018: linearly extrapolate if using dates beyond the table
@@ -55,5 +56,5 @@ def calc_delta_time(delta_file,idays):
         epoch=(1992,1,1,0,0,0))
     #-- use scipy interpolating splines to interpolate delta times
     spl=scipy.interpolate.UnivariateSpline(days,dinput[:,3],k=1,s=0,ext=0)
-    #-- return the delta time for the input date
-    return spl(idays)
+    #-- return the delta time for the input date converted to days
+    return spl(idays)/86400.0

pyTMD/compute_equilibrium_tide.py

@@ -0,0 +1,56 @@
+"""
+compute_equilibrium_tide.py (08/2020)
+Calculates the long-period equilibrium ocean tides
+Fifteen tidal spectral lines from the Cartwright-Tayler-Edden
+    tables are summed over to compute the long-period tides
+
+INPUTS:
+    t: days relative to Jan 1, 1992 (48622 MJD)
+    lat: latitudes in degrees
+
+OUTPUTS:
+    lpet: long-period equilibrium tide in meters
+
+REFERENCES:
+    Cartwright & Tayler, Geophys. J. R.A.S., 23, 45, 1971.
+    Cartwright & Edden, Geophys. J. R.A.S., 33, 253, 1973.
+"""
+import numpy as np
+
+def compute_equilibrium_tide(t, lat):
+    #-- longitude of moon
+    #-- longitude of sun
+    #-- longitude of lunar perigee
+    #-- longitude of ascending lunar node
+    PHC = np.array([290.21,280.12,274.35,343.51])
+    DPD = np.array([13.1763965,0.9856473,0.1114041,0.0529539])
+
+    #-- convert time from days relative to 1992-01-01 to 1987-01-01
+    #-- Compute 4 principal mean longitudes in radians at delta time
+    npts = len(t)
+    SHPN = np.zeros((4,npts))
+    for N in range(4):
+        ANGLE = PHC[N] + (t + 1826.0)*DPD[N]
+        SHPN[N,:] = np.pi*np.mod(ANGLE, 360.0)/180.0
+
+    #-- Assemble long-period tide potential from 15 CTE terms > 1 mm.
+    #-- Nodal term is included but not the constant term.
+    PH = np.zeros((npts))
+    Z = np.zeros((npts))
+    Z += 2.79*np.cos(SHPN[3,:]) - 0.49*np.cos(SHPN[1,:] - \
+        283.0*np.pi/180.0) - 3.10*np.cos(2.0*SHPN[1,:])
+    PH += SHPN[0,:]
+    Z += -0.67*np.cos(PH - 2.0*SHPN[1,:] + SHPN[2,:]) - \
+        (3.52 - 0.46*np.cos(SHPN[3,:]))*np.cos(PH - SHPN[2,:])
+    PH += SHPN[0,:]
+    Z += - 6.66*np.cos(PH) - 2.76*np.cos(PH + SHPN[3,:]) - \
+        0.26 * np.cos(PH + 2.*SHPN[3,:]) - 0.58 * np.cos(PH - 2.*SHPN[1,:]) - \
+        0.29 * np.cos(PH - 2.*SHPN[2,:])
+    PH += SHPN[0,:]
+    Z += - 1.27*np.cos(PH - SHPN[2,:]) - \
+        0.53*np.cos(PH - SHPN[2,:] + SHPN[3,:]) - \
+        0.24*np.cos(PH - 2.0*SHPN[1,:] + SHPN[2,:])
+
+    #-- Multiply by gamma_2 * sqrt(5/4 pi) * P20(lat)
+    lpet = 0.437*Z*(1.5*np.sin(lat*np.pi/180.0)**2 - 0.5)/100.0
+    return lpet

pyTMD/compute_tide_corrections.py

@@ -345,7 +345,9 @@ def compute_tide_corrections(x, y, delta_time, DIRECTORY=None, MODEL=None,
     elif (model_format == 'FES'):
         amp,ph = extract_FES_constants(lon, lat, model_directory, model_files,
             TYPE=model_type, VERSION=MODEL, METHOD=METHOD, SCALE=SCALE)
-        deltat = np.zeros_like(t)
+        #-- interpolate delta times from calendar dates to tide time
+        delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
+        deltat = calc_delta_time(delta_file, t)
 
     #-- calculate complex phase in radians for Euler's
     cph = -1j*ph*np.pi/180.0