refactor: use prediction function to calculate cartesian tide displac…

…ements

refactor: use io function to extract ocean pole tide values
refactor: use rotation matrix to convert from cartesian to spherical

tides/compute_LPT_ICESat_GLA12.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_LPT_ICESat_GLA12.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (08/2024)
 Calculates radial load pole tide displacements for correcting
     ICESat/GLAS L2 GLA12 Antarctic and Greenland Ice Sheet
     elevation data following IERS Convention (2010) guidelines
@@ -45,6 +45,8 @@
         doi: 10.1007/s00190-015-0848-7
 
 UPDATE HISTORY:
+    Updated 08/2024: use prediction function for cartesian tidal displacements
+        use rotation matrix to convert from cartesian to spherical
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 08/2023: create s3 filesystem when using s3 urls as input
@@ -181,15 +183,28 @@ def compute_LPT_ICESat(INPUT_FILE,
     # calculate X, Y and Z from geodetic latitude and longitude
     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
-    # geocentric colatitude in radians
-    theta = dtr*(90.0 - latitude_geocentric)
+    # geocentric latitude (radians)
+    latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))
+    # geocentric colatitude (radians)
+    theta = (np.pi/2.0 - latitude_geocentric)
+    # calculate longitude (radians)
+    phi = np.arctan2(Y, X)
 
     # compute normal gravity at spatial location
     # p. 80, Eqn.(2-199)
     gamma_0 = wgs84.gamma_0(theta)
 
+    # rotation matrix for converting from cartesian coordinates
+    R = np.zeros((n_40HZ, 3, 3))
+    R[:,0,0] = np.cos(phi)*np.cos(theta)
+    R[:,1,0] = -np.sin(phi)
+    R[:,2,0] = np.cos(phi)*np.sin(theta)
+    R[:,0,1] = np.sin(phi)*np.cos(theta)
+    R[:,1,1] = np.cos(phi)
+    R[:,2,1] = np.sin(phi)*np.sin(theta)
+    R[:,0,2] = -np.sin(theta)
+    R[:,2,2] = np.cos(theta)
+
     # calculate load pole tides in cartesian coordinates
     XYZ = np.c_[X, Y, Z]
     dxi = pyTMD.predict.load_pole_tide(ts.tide, XYZ,
@@ -200,14 +215,12 @@ def compute_LPT_ICESat(INPUT_FILE,
         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)
+    # calculate components of load pole tides
+    S = np.einsum('ti...,tji...->tj...', dxi, R)
 
     # convert to masked array
     Srad = np.ma.zeros((n_40HZ),fill_value=fv)
-    Srad.data[:] = drad.copy()
+    Srad.data[:] = S[:,2].copy()
     # replace fill values
     Srad.mask = np.isnan(Srad.data)
     Srad.data[Srad.mask] = Srad.fill_value

tides/compute_LPT_icebridge_data.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_LPT_icebridge_data.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (08/2024)
 Calculates load pole tide displacements for correcting Operation IceBridge
     elevation data following IERS Convention (2010) guidelines
     http://maia.usno.navy.mil/conventions/2010officialinfo.php
@@ -41,6 +41,8 @@
     read_ATM1b_QFIT_binary.py: read ATM1b QFIT binary files (NSIDC version 1)
 
 UPDATE HISTORY:
+    Updated 08/2024: use prediction function for cartesian tidal displacements
+        use rotation matrix to convert from cartesian to spherical
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 05/2023: use timescale class for time conversion operations
@@ -212,15 +214,28 @@ def compute_LPT_icebridge_data(arg,
     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)
+    # geocentric latitude (radians)
+    latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))
+    # geocentric colatitude (radians)
+    theta = (np.pi/2.0 - latitude_geocentric)
+    # calculate longitude (radians)
+    phi = np.arctan2(Y, X)
 
     # compute normal gravity at spatial location
     # p. 80, Eqn.(2-199)
     gamma_0 = wgs84.gamma_0(theta)
 
+    # rotation matrix for converting from cartesian coordinates
+    R = np.zeros((file_lines, 3, 3))
+    R[:,0,0] = np.cos(phi)*np.cos(theta)
+    R[:,1,0] = -np.sin(phi)
+    R[:,2,0] = np.cos(phi)*np.sin(theta)
+    R[:,0,1] = np.sin(phi)*np.cos(theta)
+    R[:,1,1] = np.cos(phi)
+    R[:,2,1] = np.sin(phi)*np.sin(theta)
+    R[:,0,2] = -np.sin(theta)
+    R[:,2,2] = np.cos(theta)
+
     # calculate load pole tides in cartesian coordinates
     XYZ = np.c_[X, Y, Z]
     dxi = pyTMD.predict.load_pole_tide(ts.tide, XYZ,
@@ -231,10 +246,8 @@ def compute_LPT_icebridge_data(arg,
         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)
+    # calculate components of load pole tides
+    S = np.einsum('ti...,tji...->tj...', dxi, R)
 
     # output load pole tide HDF5 file
     # form: rg_NASA_LOAD_POLE_TIDE_WGS84_fl1yyyymmddjjjjj.H5
@@ -257,7 +270,7 @@ def compute_LPT_icebridge_data(arg,
 
     # convert to masked array
     Srad = np.ma.zeros((file_lines),fill_value=fill_value)
-    Srad.data[:] = drad.copy()
+    Srad.data[:] = S[:,2].copy()
     # replace fill values
     Srad.mask = np.isnan(Srad.data)
     Srad.data[Srad.mask] = Srad.fill_value

tides/compute_OPT_ICESat_GLA12.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_OPT_ICESat_GLA12.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (08/2024)
 Calculates radial ocean pole tide displacements for correcting
     ICESat/GLAS L2 GLA12 Antarctic and Greenland Ice Sheet
     elevation data following IERS Convention (2010) guidelines
@@ -50,6 +50,9 @@
         doi: 10.1007/s00190-015-0848-7
 
 UPDATE HISTORY:
+    Updated 08/2024: use io function to extract ocean pole tide values
+        use prediction function to calculate cartesian tide displacements
+        use rotation matrix to convert from cartesian to spherical
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 08/2023: create s3 filesystem when using s3 urls as input
@@ -198,29 +201,31 @@ def compute_OPT_ICESat(INPUT_FILE,
     # calculate X, Y and Z from geodetic latitude and longitude
     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
-    # geocentric colatitude and longitude in radians
-    theta = dtr*(90.0 - latitude_geocentric)
-    phi = dtr*lon_40HZ
+    # geocentric latitude (radians)
+    latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))
+    # geocentric colatitude (radians)
+    theta = (np.pi/2.0 - latitude_geocentric)
+    # calculate longitude (radians)
+    phi = np.arctan2(Y, X)
 
     # read ocean pole tide map from Desai (2002)
     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)
+    # rotation matrix for converting to/from cartesian coordinates
+    R = np.zeros((n_40HZ, 3, 3))
+    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)
+    Rinv = np.linalg.inv(R)
 
     # 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)
+    UXYZ = np.einsum('ti...,tji...->tj...', np.c_[un, ue, ur], R)
 
     # calculate ocean pole tides in cartesian coordinates
     XYZ = np.c_[X, Y, Z]
@@ -234,13 +239,11 @@ def compute_OPT_ICESat(INPUT_FILE,
         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)
+    # calculate components of ocean pole tides
+    U = np.einsum('ti...,tji...->tj...', dxi, Rinv)
     # convert to masked array
     Urad = np.ma.zeros((n_40HZ),fill_value=fv)
-    Urad.data[:] = drad.copy()
+    Urad.data[:] = U[:,2].copy()
     # replace fill values
     Urad.mask = np.isnan(Urad.data)
     Urad.data[Urad.mask] = Urad.fill_value

tides/compute_OPT_icebridge_data.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_OPT_icebridge_data.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (08/2024)
 Calculates radial ocean pole tide displacements for correcting Operation
     IceBridge elevation data following IERS Convention (2010) guidelines
     http://maia.usno.navy.mil/conventions/2010officialinfo.php
@@ -46,6 +46,9 @@
     read_ATM1b_QFIT_binary.py: read ATM1b QFIT binary files (NSIDC version 1)
 
 UPDATE HISTORY:
+    Updated 08/2024: use io function to extract ocean pole tide values
+        use prediction function to calculate cartesian tide displacements
+        use rotation matrix to convert from cartesian to spherical
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 05/2023: use timescale class for time conversion operations
@@ -234,28 +237,31 @@ def compute_OPT_icebridge_data(arg,
     # calculate X, Y and Z from geodetic latitude and longitude
     X,Y,Z = pyTMD.spatial.to_cartesian(lon, lat,
         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
-    theta = dtr*(90.0 - latitude_geocentric)
-    phi = dtr*lon
+    # geocentric latitude (radians)
+    latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))
+    # geocentric colatitude (radians)
+    theta = (np.pi/2.0 - latitude_geocentric)
+    # calculate longitude (radians)
+    phi = np.arctan2(Y, X)
 
     # read ocean pole tide map from Desai (2002)
     ur, un, ue = pyTMD.io.IERS.extract_coefficients(lon,
         latitude_geocentric, method=METHOD)
-    # convert to cartesian coordinates
-    R = np.zeros((3, 3, file_lines))
-    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)
+    # rotation matrix for converting to/from cartesian coordinates
+    R = np.zeros((file_lines, 3, 3))
+    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)
+    Rinv = np.linalg.inv(R)
 
     # 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)
+    UXYZ = np.einsum('ti...,tji...->tj...', np.c_[un, ue, ur], R)
 
     # calculate ocean pole tides in cartesian coordinates
     XYZ = np.c_[X, Y, Z]
@@ -269,10 +275,8 @@ def compute_OPT_icebridge_data(arg,
         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)
+    # calculate components of ocean pole tides
+    U = np.einsum('ti...,tji...->tj...', dxi, Rinv)
 
     # output ocean pole tide HDF5 file
     # form: rg_NASA_OCEAN_POLE_TIDE_WGS84_fl1yyyymmddjjjjj.H5
@@ -295,7 +299,7 @@ def compute_OPT_icebridge_data(arg,
 
     # convert to masked array
     Urad = np.ma.zeros((file_lines),fill_value=fill_value)
-    Urad.data[:] = drad.copy()
+    Urad.data[:] = U[:,2].copy()
     # replace fill values
     Urad.mask = np.isnan(Urad.data)
     Urad.data[Urad.mask] = Urad.fill_value

tides/compute_SET_ICESat_GLA12.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_SET_ICESat_GLA12.py
-Written by Tyler Sutterley (05/2024)
+Written by Tyler Sutterley (08/2024)
 Calculates radial solid Earth tide displacements for correcting
     ICESat/GLAS L2 GLA12 Antarctic and Greenland Ice Sheet
     elevation data following IERS Convention (2010) guidelines
@@ -34,6 +34,7 @@
     predict.py: calculates solid Earth tides
 
 UPDATE HISTORY:
+    Updated 08/2024: use rotation matrix to convert from cartesian to spherical
     Updated 05/2024: use wrapper to importlib for optional dependencies
     Updated 04/2024: use timescale for temporal operations
     Updated 01/2024: refactored lunisolar ephemerides functions
@@ -161,17 +162,34 @@ def compute_SET_ICESat(INPUT_FILE,
     XYZ = np.c_[X, Y, Z]
     SXYZ = np.c_[SX, SY, SZ]
     LXYZ = np.c_[LX, LY, LZ]
+
+    # geocentric latitude (radians)
+    latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))
+    # geocentric colatitude (radians)
+    theta = (np.pi/2.0 - latitude_geocentric)
+    # calculate longitude (radians)
+    phi = np.arctan2(Y, X)
+
+    # rotation matrix for converting from cartesian coordinates
+    R = np.zeros((n_40HZ, 3, 3))
+    R[:,0,0] = np.cos(phi)*np.cos(theta)
+    R[:,1,0] = -np.sin(phi)
+    R[:,2,0] = np.cos(phi)*np.sin(theta)
+    R[:,0,1] = np.sin(phi)*np.cos(theta)
+    R[:,1,1] = np.cos(phi)
+    R[:,2,1] = np.sin(phi)*np.sin(theta)
+    R[:,0,2] = -np.sin(theta)
+    R[:,2,2] = np.cos(theta)
+
     # predict solid earth tides (cartesian)
     dxi = pyTMD.predict.solid_earth_tide(tide_time,
         XYZ, SXYZ, LXYZ, a_axis=wgs84.a_axis,
         tide_system=TIDE_SYSTEM)
-    # calculate radial component of solid earth 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)
+    # calculate components of solid earth tides
+    SE = np.einsum('ti...,tji...->tj...', dxi, R)
     # create masked array of solid earth tide displacements
     tide_se = np.ma.zeros((n_40HZ),fill_value=fv)
-    tide_se.data[:] = drad.copy()
+    tide_se.data[:] = SE[:,2].copy()
     # replace fill values
     tide_se.mask = np.isnan(tide_se.data) | (elev_40HZ == fv)
     tide_se.data[tide_se.mask] = tide_se.fill_value