From 081f8a8e3120164412ad96fbf54fe1564ad01f18 Mon Sep 17 00:00:00 2001 From: tsutterley Date: Fri, 30 Aug 2024 13:11:09 -0700 Subject: [PATCH] refactor: use prediction function to calculate cartesian tide displacements 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 | 33 ++++++++++++++------ tides/compute_LPT_icebridge_data.py | 33 ++++++++++++++------ tides/compute_OPT_ICESat_GLA12.py | 47 +++++++++++++++-------------- tides/compute_OPT_icebridge_data.py | 46 +++++++++++++++------------- tides/compute_SET_ICESat_GLA12.py | 30 ++++++++++++++---- tides/compute_SET_icebridge_data.py | 30 ++++++++++++++---- 6 files changed, 144 insertions(+), 75 deletions(-) diff --git a/tides/compute_LPT_ICESat_GLA12.py b/tides/compute_LPT_ICESat_GLA12.py index 980bcd6..24852de 100755 --- a/tides/compute_LPT_ICESat_GLA12.py +++ b/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 diff --git a/tides/compute_LPT_icebridge_data.py b/tides/compute_LPT_icebridge_data.py index 031d0d5..3536436 100755 --- a/tides/compute_LPT_icebridge_data.py +++ b/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 diff --git a/tides/compute_OPT_ICESat_GLA12.py b/tides/compute_OPT_ICESat_GLA12.py index d84e088..d2c257a 100755 --- a/tides/compute_OPT_ICESat_GLA12.py +++ b/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 diff --git a/tides/compute_OPT_icebridge_data.py b/tides/compute_OPT_icebridge_data.py index 73cea6b..af7ef6e 100755 --- a/tides/compute_OPT_icebridge_data.py +++ b/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 diff --git a/tides/compute_SET_ICESat_GLA12.py b/tides/compute_SET_ICESat_GLA12.py index c2e7c59..a1aeba1 100755 --- a/tides/compute_SET_ICESat_GLA12.py +++ b/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 diff --git a/tides/compute_SET_icebridge_data.py b/tides/compute_SET_icebridge_data.py index f987fcd..daabb8e 100755 --- a/tides/compute_SET_icebridge_data.py +++ b/tides/compute_SET_icebridge_data.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_SET_icebridge_data.py -Written by Tyler Sutterley (05/2024) +Written by Tyler Sutterley (08/2024) Calculates radial solid Earth tide displacements for correcting Operation IceBridge elevation data following IERS Convention (2010) guidelines http://maia.usno.navy.mil/conventions/2010officialinfo.php @@ -37,6 +37,7 @@ read_ATM1b_QFIT_binary.py: read ATM1b QFIT binary files (NSIDC version 1) 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 @@ -195,16 +196,33 @@ def compute_SET_icebridge_data(arg, 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((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) + # 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) # save solid earth tide displacements to output dictionary - dinput['tide_earth'] = drad.copy() + dinput['tide_earth'] = SE[:,2].copy() # calculate permanent tide offset (meters) dinput['tide_earth_free2mean'] = 0.06029 - \ 0.180873*np.sin(dinput['lat']*np.pi/180.0)**2