refactor: update pole tides to use pyTMD cartesian functions

.gitignore

@@ -9,6 +9,9 @@
 *.pyc
 *.m~
 *.m~~
+# Cython Files #
+################
+_*.c
 # Packages #
 ############
 *.7z
@@ -19,6 +22,12 @@
 *.rar
 *.tar
 *.zip
+# Data Files #
+##############
+*.nc
+*.h5
+*.hdf
+*.tif
 # Logs and databases #
 ######################
 *.log
@@ -37,8 +46,7 @@ wheels/
 *.egg-info/
 .installed.cfg
 *.egg
-.pytest_cache
-pythonenv*/
+setup-miniconda-patched-environment.yml
 # OS generated files #
 ######################
 .DS_Store
@@ -49,7 +57,7 @@ pythonenv*/
 ehthumbs.db
 Thumbs.db
 # LaTeX and Vim #
-######################
+#################
 *.aux
 *.bbl
 *.blg
@@ -80,6 +88,10 @@ Thumbs.db
 *.sw*
 *.hidden
 None*.png
+# JPL ephemerides #
+# #################
+*.bpc
+*.bsp
 # Jupyter Checkpoints #
 #######################
 .ipynb_checkpoints

test/test_tide_corrections.py

@@ -15,8 +15,8 @@
 timescale = gz.utilities.import_dependency('timescale')
 
 # path to an ATL03 file from NSIDC
-ATL03 = ['https://n5eil01u.ecs.nsidc.org','ATLAS','ATL03.005','2018.10.13',
-    'ATL03_20181013235645_02340114_005_01.h5']
+ATL03 = ['https://n5eil01u.ecs.nsidc.org','ATLAS','ATL03.006','2018.10.14',
+    'ATL03_20181014000347_02350101_006_02.h5']
 # PURPOSE: Download ATL03 granule from NSIDC
 @pytest.fixture(scope="module", autouse=True)
 def download_ATL03(username, password):
@@ -56,15 +56,10 @@ def test_ATL03_equilibrium_tides():
         fv = IS2_atl03_attrs[gtx]['geophys_corr']['tide_equilibrium']['_FillValue']
         tide_equilibrium = IS2_atl03_mds[gtx]['geophys_corr']['tide_equilibrium']
         # calculate tide time for beam
-        gps_seconds = atlas_sdp_gps_epoch + delta_time
-        leap_seconds = timescale.time.count_leap_seconds(gps_seconds)
-        tide_time = timescale.time.convert_delta_time(gps_seconds-leap_seconds,
-            epoch1=(1980,1,6,0,0,0), epoch2=(1992,1,1,0,0,0), scale=1.0/86400.0)
-        # interpolate delta times from calendar dates to tide time
-        delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
-        deltat = timescale.time.interpolate_delta_time(delta_file, tide_time)
+        ts = timescale.time.Timescale().from_deltatime(delta_time,
+            epoch=(2018,1,1), standard='GPS')
         # calculate long-period equilibrium tides
-        lpet = pyTMD.predict.equilibrium_tide(tide_time+deltat, latitude)
+        lpet = pyTMD.predict.equilibrium_tide(ts.tide + ts.tt_ut1, latitude)
         ii, = np.nonzero(tide_equilibrium != fv)
         # calculate differences between computed and data versions
         difference = np.ma.zeros((nref))
@@ -98,7 +93,7 @@ def test_ATL03_load_pole_tide():
         # calculate load pole tides from correction function
         Srad = pyTMD.compute.LPT_displacements(longitude, latitude,
             delta_time, EPSG=4326, EPOCH=(2018,1,1,0,0,0), TYPE='drift',
-            TIME='GPS', ELLIPSOID='IERS', CONVENTION='2010')
+            TIME='GPS', ELLIPSOID='IERS', CONVENTION='2003')
         # calculate differences between computed and data versions
         difference = np.ma.zeros((nref))
         difference.data[:] = Srad - tide_pole
@@ -131,7 +126,7 @@ def test_ATL03_ocean_pole_tide():
         # calculate ocean pole tides from correction function
         Urad = pyTMD.compute.OPT_displacements(longitude, latitude,
             delta_time, EPSG=4326, EPOCH=(2018,1,1,0,0,0), TYPE='drift',
-            TIME='GPS', ELLIPSOID='IERS', CONVENTION='2010')
+            TIME='GPS', ELLIPSOID='IERS', CONVENTION='2003')
         # calculate differences between computed and data versions
         difference = np.ma.zeros((nref))
         difference.data[:] = Urad - tide_oc_pole
@@ -142,8 +137,8 @@ def test_ATL03_ocean_pole_tide():
             assert np.all(np.abs(difference) < eps)
 
 # path to an ATL07 file from NSIDC
-ATL07 = ['https://n5eil01u.ecs.nsidc.org','ATLAS','ATL07.005','2018.10.14',
-    'ATL07-01_20181014000347_02350101_005_03.h5']
+ATL07 = ['https://n5eil01u.ecs.nsidc.org','ATLAS','ATL07.006','2018.10.14',
+    'ATL07-01_20181014000347_02350101_006_02.h5']
 # PURPOSE: Download ATL07 granule from NSIDC
 @pytest.fixture(scope="module", autouse=True)
 def download_ATL07(username, password):
@@ -186,15 +181,10 @@ def test_ATL07_equilibrium_tides():
         fv = attrs['geophysical']['height_segment_lpe']['_FillValue']
         height_segment_lpe = val['geophysical']['height_segment_lpe'][:]
         # calculate tide time for beam
-        gps_seconds = atlas_sdp_gps_epoch + delta_time
-        leap_seconds = timescale.time.count_leap_seconds(gps_seconds)
-        tide_time = timescale.time.convert_delta_time(gps_seconds-leap_seconds,
-            epoch1=(1980,1,6,0,0,0), epoch2=(1992,1,1,0,0,0), scale=1.0/86400.0)
-        # interpolate delta times from calendar dates to tide time
-        delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
-        deltat = timescale.time.interpolate_delta_time(delta_file, tide_time)
+        ts = timescale.time.Timescale().from_deltatime(delta_time,
+            epoch=(2018,1,1), standard='GPS')
         # calculate long-period equilibrium tides
-        lpet = pyTMD.predict.equilibrium_tide(tide_time+deltat, latitude)
+        lpet = pyTMD.predict.equilibrium_tide(ts.tide + ts.tt_ut1, latitude)
         # calculate differences between computed and data versions
         difference = np.ma.zeros((nseg))
         difference.data[:] = lpet - height_segment_lpe
@@ -204,7 +194,7 @@ def test_ATL07_equilibrium_tides():
         if not np.all(difference.mask):
             assert np.all(np.abs(difference) < eps)
         # calculate long-period equilibrium tides from correction function
-        lpet = pyTMD.compute.LPET_displacements(longitude, latitude,
+        lpet = pyTMD.compute.LPET_elevations(longitude, latitude,
             delta_time, EPSG=4326, EPOCH=(2018,1,1,0,0,0), TYPE='drift',
             TIME='GPS')
         # calculate differences between computed and data versions

tides/compute_LPT_ICESat_GLA12.py

@@ -160,13 +160,6 @@ def compute_LPT_ICESat(INPUT_FILE,
     # create timescale from J2000: seconds since 2000-01-01 12:00:00 UTC
     ts = timescale.time.Timescale().from_deltatime(DS_UTCTime_40HZ[:],
         epoch=timescale.time._j2000_epoch, standard='UTC')
-    # convert dynamic time to Modified Julian Days (MJD)
-    MJD = ts.tt - 2400000.5
-    # convert Julian days to calendar dates
-    Y,M,D,h,m,s = timescale.time.convert_julian(ts.tt, format='tuple')
-    # calculate time in year-decimal format
-    time_decimal = timescale.time.convert_calendar_decimal(Y,M,day=D,
-        hour=h,minute=m,second=s)
 
     # parameters for Topex/Poseidon and WGS84 ellipsoids
     topex = pyTMD.datum(ellipsoid='TOPEX', units='MKS')
@@ -180,38 +173,41 @@ def compute_LPT_ICESat(INPUT_FILE,
 
     # degrees to radians
     dtr = np.pi/180.0
-    atr = np.pi/648000.0
-    # tidal love number appropriate for the load tide
+    # tidal love/shida numbers appropriate for the load tide
     hb2 = 0.6207
+    lb2 = 0.0847
 
     # 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_40HZ, lat_40HZ, h=elev_40HZ,
+    X,Y,Z = pyTMD.spatial.to_cartesian(lon_40HZ, lat_40HZ,
         a_axis=wgs84.a_axis, flat=wgs84.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
-    # colatitude and longitude in radians
+    # geocentric colatitude in radians
     theta = dtr*(90.0 - latitude_geocentric)
-    phi = lon_40HZ*dtr
-
-    # compute normal gravity at spatial location and elevation of points.
-    # Normal gravity at height h. p. 82, Eqn.(2-215)
-    gamma_h = wgs84.gamma_h(theta, elev_40HZ)
-
-    # calculate angular coordinates of mean/secular pole at time
-    mpx, mpy, fl = timescale.eop.iers_mean_pole(time_decimal,
-        convention=CONVENTION)
-    # read and interpolate IERS daily polar motion values
-    px, py = timescale.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
-    dfactor = -hb2*atr*(wgs84.omega**2*rr**2)/(2.0*gamma_h)
+
+    # compute normal gravity at spatial location
+    # p. 80, Eqn.(2-199)
+    gamma_0 = wgs84.gamma_0(theta)
+
+    # calculate load pole tides in cartesian coordinates
+    XYZ = np.c_[X, Y, Z]
+    dxi = pyTMD.predict.load_pole_tide(ts.tide, XYZ,
+        deltat=ts.tt_ut1,
+        gamma_0=gamma_0,
+        omega=wgs84.omega,
+        h2=hb2,
+        l2=lb2,
+        convention=CONVENTION
+    )
+    # calculate radial component of load pole tides
+    dln,dlt,drad = pyTMD.spatial.to_geodetic(
+        X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2],
+        a_axis=wgs84.a_axis, flat=wgs84.flat)
+
+    # convert to masked array
     Srad = np.ma.zeros((n_40HZ),fill_value=fv)
-    Srad.data[:] = dfactor*np.sin(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi))
+    Srad.data[:] = drad.copy()
     # replace fill values
     Srad.mask = np.isnan(Srad.data)
     Srad.data[Srad.mask] = Srad.fill_value

tides/compute_LPT_icebridge_data.py

@@ -196,40 +196,45 @@ def compute_LPT_icebridge_data(arg,
     # create timescale from J2000: seconds since 2000-01-01 12:00:00 UTC
     ts = timescale.time.Timescale().from_deltatime(dinput['time'],
         epoch=timescale.time._j2000_epoch, standard='UTC')
-    # convert dynamic time to Modified Julian Days (MJD)
-    MJD = ts.tt - 2400000.5
-    # convert Julian days to calendar dates
-    Y,M,D,h,m,s = timescale.time.convert_julian(ts.tt, format='tuple')
-    # calculate time in year-decimal format
-    time_decimal = timescale.time.convert_calendar_decimal(Y,M,day=D,
-        hour=h,minute=m,second=s)
-    # elevation
-    h1 = np.copy(dinput['data'][:])
 
     # degrees to radians
     dtr = np.pi/180.0
-    atr = np.pi/648000.0
     # earth and physical parameters for ellipsoid
     wgs84 = pyTMD.datum(ellipsoid='WGS84', units='MKS')
-    # tidal love number appropriate for the load tide
+    # tidal love/shida numbers appropriate for the load tide
     hb2 = 0.6207
+    lb2 = 0.0847
     # bad value
     fill_value = -9999.0
 
     # 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, lat, h=h1,
+    X,Y,Z = pyTMD.spatial.to_cartesian(lon, lat,
         a_axis=wgs84.a_axis, flat=wgs84.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
     # colatitude and longitude in radians
     theta = dtr*(90.0 - latitude_geocentric)
-    phi = lon*dtr
 
-    # compute normal gravity at spatial location and elevation of points.
-    # Normal gravity at height h. p. 82, Eqn.(2-215)
-    gamma_h = wgs84.gamma_h(theta, h1)
+    # compute normal gravity at spatial location
+    # p. 80, Eqn.(2-199)
+    gamma_0 = wgs84.gamma_0(theta)
+
+    # calculate load pole tides in cartesian coordinates
+    XYZ = np.c_[X, Y, Z]
+    dxi = pyTMD.predict.load_pole_tide(ts.tide, XYZ,
+        deltat=ts.tt_ut1,
+        gamma_0=gamma_0,
+        omega=wgs84.omega,
+        h2=hb2,
+        l2=lb2,
+        convention=CONVENTION
+    )
+    # calculate radial component of load pole tides
+    dln,dlt,drad = pyTMD.spatial.to_geodetic(
+        X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2],
+        a_axis=wgs84.a_axis, flat=wgs84.flat)
 
     # output load pole tide HDF5 file
     # form: rg_NASA_LOAD_POLE_TIDE_WGS84_fl1yyyymmddjjjjj.H5
@@ -250,18 +255,9 @@ def compute_LPT_icebridge_data(arg,
     # open output HDF5 file
     fid = h5py.File(output_file, mode='w')
 
-    # calculate angular coordinates of mean/secular pole at time
-    mpx, mpy, fl = timescale.eop.iers_mean_pole(time_decimal,
-        convention=CONVENTION)
-    # read and interpolate IERS daily polar motion values
-    px, py = timescale.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
-    dfactor = -hb2*atr*(wgs84.omega**2*rr**2)/(2.0*gamma_h)
+    # convert to masked array
     Srad = np.ma.zeros((file_lines),fill_value=fill_value)
-    Srad.data[:] = dfactor*np.sin(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi))
+    Srad.data[:] = drad.copy()
     # replace fill values
     Srad.mask = np.isnan(Srad.data)
     Srad.data[Srad.mask] = Srad.fill_value

tides/compute_OPT_ICESat_GLA12.py

@@ -185,9 +185,8 @@ def compute_OPT_ICESat(INPUT_FILE,
         wgs84.a_axis, wgs84.flat,
         eps=1e-12, itmax=10)
 
-    # degrees to radians and arcseconds to radians
+    # degrees to radians
     dtr = np.pi/180.0
-    atr = np.pi/648000.0
     # mean equatorial gravitational acceleration [m/s^2]
     ge = 9.7803278
     # density of sea water [kg/m^3]
@@ -197,47 +196,51 @@ def compute_OPT_ICESat(INPUT_FILE,
 
     # 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_40HZ, lat_40HZ, h=elev_40HZ,
+    X,Y,Z = pyTMD.spatial.to_cartesian(lon_40HZ, lat_40HZ,
         a_axis=wgs84.a_axis, flat=wgs84.flat)
     # 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)*(wgs84.omega**2*wgs84.a_axis**4)/wgs84.GM
-    K = 4.0*np.pi*wgs84.G*rho_w*Hp*wgs84.a_axis/(3.0*ge)
-    K1 = 4.0*np.pi*wgs84.G*rho_w*Hp*wgs84.a_axis**3/(3.0*wgs84.GM)
-
-    # calculate angular coordinates of mean/secular pole at time
-    mpx, mpy, fl = timescale.eop.iers_mean_pole(time_decimal,
-        convention=CONVENTION)
-    # read and interpolate IERS daily polar motion values
-    px, py = timescale.eop.iers_polar_motion(MJD, k=3, s=0)
-    # calculate differentials from mean/secular pole positions
-    mx = px - mpx
-    my = -(py - mpy)
+    # geocentric colatitude and longitude in radians
+    theta = dtr*(90.0 - latitude_geocentric)
+    phi = dtr*lon_40HZ
 
     # 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((n_40HZ),dtype=np.clongdouble)
-        UR.real = f1.ev(lon_40HZ,latitude_geocentric)
-        UR.imag = f2.ev(lon_40HZ,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_40HZ,latitude_geocentric])
-
-    # calculate radial displacement at time
+    ur, un, ue = pyTMD.io.IERS.extract_coefficients(lon_40HZ,
+        latitude_geocentric, method=METHOD)
+    # convert to cartesian coordinates
+    R = np.zeros((3, 3, n_40HZ))
+    R[0,0,:] = np.cos(phi)*np.cos(theta)
+    R[0,1,:] = -np.sin(phi)
+    R[0,2,:] = np.cos(phi)*np.sin(theta)
+    R[1,0,:] = np.sin(phi)*np.cos(theta)
+    R[1,1,:] = np.cos(phi)
+    R[1,2,:] = np.sin(phi)*np.sin(theta)
+    R[2,0,:] = -np.sin(theta)
+    R[2,2,:] = np.cos(theta)
+
+    # calculate pole tide displacements in Cartesian coordinates
+    # coefficients reordered to N, E, R to match IERS rotation matrix
+    UXYZ = np.einsum('ti...,jit...->tj...', np.c_[un, ue, ur], R)
+
+    # calculate ocean pole tides in cartesian coordinates
+    XYZ = np.c_[X, Y, Z]
+    dxi = pyTMD.predict.ocean_pole_tide(ts.tide, XYZ, UXYZ,
+        deltat=ts.tt_ut1,
+        a_axis=wgs84.a_axis,
+        gamma_0=ge,
+        GM=wgs84.GM,
+        omega=wgs84.omega,
+        rho_w=rho_w,
+        g2=gamma,
+        convention=CONVENTION
+    )
+    # calculate radial component of ocean pole tides
+    dln,dlt,drad = pyTMD.spatial.to_geodetic(
+        X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2],
+        a_axis=wgs84.a_axis, flat=wgs84.flat)
+    # convert to masked array
     Urad = np.ma.zeros((n_40HZ),fill_value=fv)
-    Urad.data[:] = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real +
-        (my*gamma.real - mx*gamma.imag)*UR.imag)
+    Urad.data[:] = drad.copy()
     # replace fill values
     Urad.mask = np.isnan(Urad.data)
     Urad.data[Urad.mask] = Urad.fill_value