From bd9d549fbb41d28fe7c66e58f17f057e40afe2e5 Mon Sep 17 00:00:00 2001 From: tsutterley Date: Tue, 27 Aug 2024 14:38:42 -0700 Subject: [PATCH] feat: add functions to calculate pole tides in cartesian coordinates for #323 refactor: renamed io for Desai ocean pole tide file to IERS --- .../io/{ocean_pole_tide.rst => IERS.rst} | 14 +- doc/source/api_reference/io/io.rst | 2 +- pyTMD/compute.py | 257 +++++---- pyTMD/io/IERS.py | 331 ++++++++++++ pyTMD/io/__init__.py | 2 +- pyTMD/io/ocean_pole_tide.py | 185 ------- pyTMD/predict.py | 235 ++++++++- scripts/compute_LPT_displacements.py | 91 ++-- scripts/compute_OPT_displacements.py | 117 +++-- scripts/compute_SET_displacements.py | 22 +- test/test_ocean_pole_tide.py | 126 ----- test/test_pole_tide.py | 494 ++++++++++++++++++ 12 files changed, 1365 insertions(+), 511 deletions(-) rename doc/source/api_reference/io/{ocean_pole_tide.rst => IERS.rst} (70%) create mode 100644 pyTMD/io/IERS.py delete mode 100644 pyTMD/io/ocean_pole_tide.py delete mode 100755 test/test_ocean_pole_tide.py create mode 100755 test/test_pole_tide.py diff --git a/doc/source/api_reference/io/ocean_pole_tide.rst b/doc/source/api_reference/io/IERS.rst similarity index 70% rename from doc/source/api_reference/io/ocean_pole_tide.rst rename to doc/source/api_reference/io/IERS.rst index 1e6991af..8916a3a7 100644 --- a/doc/source/api_reference/io/ocean_pole_tide.rst +++ b/doc/source/api_reference/io/IERS.rst @@ -1,6 +1,6 @@ -=============== -ocean_pole_tide -=============== +==== +IERS +==== - Reads ocean pole load tide coefficients provided by IERS as computed by `Desai et al. (2002) `_ and `Desai et al. (2015) `_ - See `materials from Chapter 7 of the IERS Conventions `_ @@ -14,10 +14,12 @@ Calling Sequence import pyTMD.io import pyTMD.utilities ocean_pole_tide_file = pyTMD.utilities.get_data_path(['data','opoleloadcoefcmcor.txt.gz']) - ur,un,ue,glon,glat = pyTMD.io.ocean_pole_tide(ocean_pole_tide_file) + ur,un,ue,glon,glat = pyTMD.io.IERS.read_binary_file(model_file=ocean_pole_tide_file) `Source code`__ -.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/io/ocean_pole_tide.py +.. __: https://github.com/tsutterley/pyTMD/blob/main/pyTMD/io/IERS.py -.. autofunction:: pyTMD.io.ocean_pole_tide +.. autofunction:: pyTMD.io.IERS.extract_coefficients + +.. autofunction:: pyTMD.io.IERS.read_binary_file diff --git a/doc/source/api_reference/io/io.rst b/doc/source/api_reference/io/io.rst index 17523d74..02041110 100644 --- a/doc/source/api_reference/io/io.rst +++ b/doc/source/api_reference/io/io.rst @@ -11,4 +11,4 @@ io OTIS.rst constituents.rst model.rst - ocean_pole_tide.rst + IERS.rst diff --git a/pyTMD/compute.py b/pyTMD/compute.py index e1018471..9dadf88a 100644 --- a/pyTMD/compute.py +++ b/pyTMD/compute.py @@ -69,6 +69,8 @@ renamed format for netcdf to ATLAS-netcdf renamed format for FES to FES-netcdf and added FES-ascii renamed format for GOT to GOT-ascii and added GOT-netcdf + drop use of heights when converting to cartesian coordinates + use prediction function to calculate cartesian tide displacements Updated 06/2024: use np.clongdouble instead of np.longcomplex Updated 04/2024: use wrapper to importlib for optional dependencies Updated 02/2024: changed class name for ellipsoid parameters to datum @@ -780,7 +782,6 @@ def LPET_elevations( # following IERS Convention (2010) guidelines def LPT_displacements( x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - h: float | np.ndarray = 0.0, EPSG: str | int = 3031, EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), TYPE: str | None = 'drift', @@ -802,8 +803,6 @@ def LPT_displacements( y-coordinates in projection EPSG delta_time: np.ndarray seconds since EPOCH or datetime array - h: float or np.ndarray, default 0.0 - height coordinates in meters EPSG: int, default: 3031 (Polar Stereographic South, WGS84) Input coordinate system EPOCH: tuple, default (2000,1,1,0,0,0) @@ -876,83 +875,106 @@ def LPT_displacements( ts = timescale.time.Timescale().from_deltatime(delta_time, epoch=EPOCH, standard=TIME) - # convert dynamic time to Modified Julian Days (MJD) - MJD = ts.tt - _jd_mjd - # 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) # number of time points - nt = len(time_decimal) + nt = len(ts) - # degrees and arcseconds to radians + # degrees to radians dtr = np.pi/180.0 - atr = np.pi/648000.0 # earth and physical parameters for ellipsoid units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') - # tidal love number appropriate for the load tide + # tidal love/shida numbers appropriate for the load tide hb2 = 0.6207 + lb2 = 0.0836 - # flatten heights - h = np.array(h).flatten() if np.ndim(h) else h # 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.flatten(), lat.flatten(), h=h, + X,Y,Z = pyTMD.spatial.to_cartesian(lon.flatten(), lat.flatten(), a_axis=units.a_axis, flat=units.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 - # geocentric colatitude and longitude in radians + # geocentric colatitude in radians theta = dtr*(90.0 - latitude_geocentric) - phi = dtr*lon.flatten() - # compute normal gravity at spatial location and elevation of points. - # Normal gravity at height h. p. 82, Eqn.(2-215) - gamma_h = units.gamma_h(theta, h) - dfactor = -hb2*atr*(units.omega**2*rr**2)/(2.0*gamma_h) - - # 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) + # compute normal gravity at spatial location + # p. 80, Eqn.(2-199) + gamma_0 = units.gamma_0(theta) # calculate radial displacement at time if (TYPE == 'grid'): ny,nx = np.shape(x) Srad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) Srad.mask = np.zeros((ny,nx,nt),dtype=bool) + XYZ = np.c_[X, Y, Z] for i in range(nt): - SRAD = dfactor*np.sin(2.0*theta)*(mx[i]*np.cos(phi)+my[i]*np.sin(phi)) - # reform grid - Srad.data[:,:,i] = np.reshape(SRAD, (ny,nx)) + # calculate load pole tides in cartesian coordinates + dxi = pyTMD.predict.load_pole_tide(ts.tide[i], XYZ, + deltat=ts.tt_ut1[i], + gamma_0=gamma_0, + omega=units.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=units.a_axis, flat=units.flat) + # reshape to output dimensions + Srad.data[:,:,i] = np.reshape(drad, (ny,nx)) Srad.mask[:,:,i] = np.isnan(Srad.data[:,:,i]) elif (TYPE == 'drift'): + # 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=units.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=units.a_axis, flat=units.flat) + # reshape to output dimensions Srad = np.ma.zeros((nt), fill_value=FILL_VALUE) - Srad.data[:] = dfactor*np.sin(2.0*theta)*(mx*np.cos(phi)+my*np.sin(phi)) + Srad.data[:] = drad.copy() Srad.mask = np.isnan(Srad.data) elif (TYPE == 'time series'): nstation = len(x) Srad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) Srad.mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): - SRAD = dfactor[s]*np.sin(2.0*theta[s])*(mx*np.cos(phi[s])+my*np.sin(phi[s])) - Srad.data[s,:] = np.copy(SRAD) + # convert coordinates to column arrays + XYZ = np.repeat(np.c_[X[s], Y[s], Z[s]], nt, axis=0) + # calculate load pole tides in cartesian coordinates + dxi = pyTMD.predict.load_pole_tide(ts.tide, XYZ, + deltat=ts.tt_ut1, + gamma_0=gamma_0[s], + omega=units.omega, + h2=hb2, + l2=lb2, + convention=CONVENTION + ) + # calculate radial component of load pole tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X[s] + dxi[:,0], Y[s] + dxi[:,1], Z[s] + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # reshape to output dimensions + Srad.data[s,:] = drad.copy() Srad.mask[s,:] = np.isnan(Srad.data[s,:]) + # replace invalid data with fill values Srad.data[Srad.mask] = Srad.fill_value - # return the solid earth pole tide displacements + # return the load pole tide displacements return Srad # PURPOSE: compute radial load pole tide displacements # following IERS Convention (2010) guidelines def OPT_displacements( x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - h: float | np.ndarray = 0.0, EPSG: str | int = 3031, EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), TYPE: str | None = 'drift', @@ -975,8 +997,6 @@ def OPT_displacements( y-coordinates in projection EPSG delta_time: np.ndarray seconds since EPOCH or datetime array - h: float or np.ndarray, default 0.0 - height coordinates in meters EPSG: int, default: 3031 (Polar Stereographic South, WGS84) Input coordinate system EPOCH: tuple, default (2000,1,1,0,0,0) @@ -1066,9 +1086,8 @@ def OPT_displacements( # number of time points nt = len(time_decimal) - # degrees and arcseconds to radians + # degrees to radians dtr = np.pi/180.0 - atr = np.pi/648000.0 # earth and physical parameters for ellipsoid units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') # mean equatorial gravitational acceleration [m/s^2] @@ -1078,72 +1097,108 @@ def OPT_displacements( # tidal love number differential (1 + kl - hl) for pole tide frequencies gamma = 0.6870 + 0.0036j - # flatten heights - h = np.array(h).flatten() if np.ndim(h) else h # 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.flatten(), lat.flatten(), h=h, + X,Y,Z = pyTMD.spatial.to_cartesian(lon.flatten(), lat.flatten(), a_axis=units.a_axis, flat=units.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 + npts = len(latitude_geocentric) + # geocentric colatitude and longitude in radians + theta = dtr*(90.0 - latitude_geocentric) + phi = dtr*lon.flatten() - # pole tide displacement scale factor - Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM - K = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis/(3.0*ge) - K1 = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis**3/(3.0*units.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) - - # 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((len(latitude_geocentric)), dtype=np.clongdouble) - UR.real = f1.ev(lon.flatten(), latitude_geocentric) - UR.imag = f2.ev(lon.flatten(), 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.flatten(), latitude_geocentric]) + # read and interpolate ocean pole tide map from Desai (2002) + ur, un, ue = pyTMD.io.IERS.extract_coefficients(lon.flatten(), + latitude_geocentric, method=METHOD) + # convert to cartesian coordinates + R = np.zeros((3, 3, npts)) + 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 radial displacement at time if (TYPE == 'grid'): ny,nx = np.shape(x) - Urad = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) + Urad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) Urad.mask = np.zeros((ny,nx,nt),dtype=bool) + XYZ = np.c_[X, Y, Z] for i in range(nt): - URAD = K*atr*np.real((mx[i]*gamma.real + my[i]*gamma.imag)*UR.real + - (my[i]*gamma.real - mx[i]*gamma.imag)*UR.imag) - # reform grid - Urad.data[:,:,i] = np.reshape(URAD, (ny,nx)) + # calculate ocean pole tides in cartesian coordinates + dxi = pyTMD.predict.ocean_pole_tide(ts.tide[i], XYZ, UXYZ, + deltat=ts.tt_ut1[i], + a_axis=units.a_axis, + gamma_0=ge, + GM=units.GM, + omega=units.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=units.a_axis, flat=units.flat) + # reshape to output dimensions + Urad.data[:,:,i] = np.reshape(drad, (ny,nx)) Urad.mask[:,:,i] = np.isnan(Urad.data[:,:,i]) elif (TYPE == 'drift'): - Urad = np.ma.zeros((nt),fill_value=FILL_VALUE) - Urad.data[:] = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real + - (my*gamma.real - mx*gamma.imag)*UR.imag) + # 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=units.a_axis, + gamma_0=ge, + GM=units.GM, + omega=units.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=units.a_axis, flat=units.flat) + # convert to masked array + Urad = np.ma.zeros((nt), fill_value=FILL_VALUE) + Urad.data[:] = drad.copy() Urad.mask = np.isnan(Urad.data) elif (TYPE == 'time series'): nstation = len(x) - Urad = np.ma.zeros((nstation,nt),fill_value=FILL_VALUE) + Urad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) Urad.mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): - URAD = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real[s] + - (my*gamma.real - mx*gamma.imag)*UR.imag[s]) - Urad.data[s,:] = np.copy(URAD) + # convert coordinates to column arrays + XYZ = np.repeat(np.c_[X[s], Y[s], Z[s]], nt, axis=0) + uxyz = np.repeat(np.atleast_2d(UXYZ[s,:]), nt, axis=0) + # calculate ocean pole tides in cartesian coordinates + dxi = pyTMD.predict.ocean_pole_tide(ts.tide, XYZ, uxyz, + deltat=ts.tt_ut1, + a_axis=units.a_axis, + gamma_0=ge, + GM=units.GM, + omega=units.omega, + rho_w=rho_w, + g2=gamma, + convention=CONVENTION + ) + # calculate radial component of ocean pole tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X[s] + dxi[:,0], Y[s] + dxi[:,1], Z[s] + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # reshape to output dimensions + Urad.data[s,:] = drad.copy() Urad.mask[s,:] = np.isnan(Urad.data[s,:]) + # replace invalid data with fill values Urad.data[Urad.mask] = Urad.fill_value @@ -1153,7 +1208,6 @@ def OPT_displacements( # PURPOSE: compute solid earth tidal elevations def SET_displacements( x: np.ndarray, y: np.ndarray, delta_time: np.ndarray, - h: float | np.ndarray = 0.0, EPSG: str | int = 3031, EPOCH: list | tuple = (2000, 1, 1, 0, 0, 0), TYPE: str | None = 'drift', @@ -1175,8 +1229,6 @@ def SET_displacements( y-coordinates in projection EPSG delta_time: np.ndarray seconds since EPOCH or datetime array - h: float or np.ndarray, default 0.0 - height coordinates in meters EPSG: int, default: 3031 (Polar Stereographic South, WGS84) Input coordinate system EPOCH: tuple, default (2000,1,1,0,0,0) @@ -1258,7 +1310,7 @@ def SET_displacements( units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') # convert input coordinates to cartesian - X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, h=h, + X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, a_axis=units.a_axis, flat=units.flat) # compute ephemerides for lunisolar coordinates SX, SY, SZ = pyTMD.astro.solar_ecef(ts.MJD, ephemerides=EPHEMERIDES) @@ -1284,9 +1336,8 @@ def SET_displacements( dln,dlt,drad = pyTMD.spatial.to_geodetic( X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se[:,:,i] = np.reshape(drad - h, (ny,nx)) + # reshape to output dimensions + tide_se[:,:,i] = np.reshape(drad, (ny,nx)) elif (TYPE == 'drift'): # convert coordinates to column arrays XYZ = np.c_[X, Y, Z] @@ -1300,9 +1351,8 @@ def SET_displacements( dln,dlt,drad = pyTMD.spatial.to_geodetic( X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se = drad - h + # reshape to output dimensions + tide_se = np.copy(drad) elif (TYPE == 'time series'): nstation = len(x) tide_se = np.zeros((nstation,nt)) @@ -1318,11 +1368,10 @@ def SET_displacements( 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], + X[s] + dxi[:,0], Y[s] + dxi[:,1], Z[s] + dxi[:,2], a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se[s,:] = drad - h + # reshape to output dimensions + tide_se[s,:] = np.copy(drad) # return the solid earth tide displacements return tide_se diff --git a/pyTMD/io/IERS.py b/pyTMD/io/IERS.py new file mode 100644 index 00000000..55da3cbd --- /dev/null +++ b/pyTMD/io/IERS.py @@ -0,0 +1,331 @@ +#!/usr/bin/env python +u""" +IERS.py +Written by Tyler Sutterley (08/2024) + +Reads ocean pole load tide coefficients provided by IERS +http://maia.usno.navy.mil/conventions/2010/2010_official/chapter7/tn36_c7.pdf +http://maia.usno.navy.mil/conventions/2010/2010_update/chapter7/icc7.pdf + +IERS 0.5x0.5 map of ocean pole tide coefficients: +ftp://maia.usno.navy.mil/conventions/2010/2010_update/chapter7/additional_info/ + opoleloadcoefcmcor.txt.gz + +OUTPUTS: + ur: radial ocean pole tide coefficients + un: north ocean pole tide coefficients + ue: east ocean pole tide coefficients + glon: ocean grid longitude + glat: ocean grid latitude + +PYTHON DEPENDENCIES: + numpy: Scientific Computing Tools For Python + https://numpy.org + https://numpy.org/doc/stable/user/numpy-for-matlab-users.html + +REFERENCES: + S. Desai, "Observing the pole tide with satellite altimetry", Journal of + Geophysical Research: Oceans, 107(C11), 2002. doi: 10.1029/2001JC001224 + S. Desai, J. Wahr and B. Beckley "Revisiting the pole tide for and from + satellite altimetry", Journal of Geodesy, 89(12), p1233-1243, 2015. + doi: 10.1007/s00190-015-0848-7 + +UPDATE HISTORY: + Updated 08/2024: convert outputs to be in -180:180 longitude convention + added function to interpolate ocean pole tide values to coordinates + renamed from ocean_pole_tide to IERS + Updated 06/2024: use np.clongdouble instead of np.longcomplex + Updated 05/2023: add default for ocean pole tide file + Updated 04/2023: using pathlib to define and expand paths + Updated 03/2023: add basic variable typing to function inputs + Updated 12/2022: refactor ocean pole tide read programs under io + Updated 04/2022: updated docstrings to numpy documentation format + use longcomplex data format to be windows compliant + Updated 07/2021: added check that ocean pole tide file is accessible + Updated 03/2021: replaced numpy bool/int to prevent deprecation warnings + Updated 08/2020: output north load and east load deformation components + Updated 07/2020: added function docstrings + Updated 12/2018: Compatibility updates for Python3 + Written 09/2017 +""" +from __future__ import annotations + +import re +import gzip +import pathlib +import warnings +import numpy as np +import scipy.interpolate +from pyTMD.utilities import get_data_path + +# ocean pole tide file from Desai (2002) and IERS conventions +_ocean_pole_tide_file = get_data_path(['data','opoleloadcoefcmcor.txt.gz']) + +# PURPOSE: extract ocean pole tide values from Desai (2002) at coordinates +def extract_coefficients( + lon: np.ndarray, lat: np.ndarray, + **kwargs + ): + """ + Reads ocean pole tide file from [1]_ [2]_ and spatially interpolates + to input coordinates + + Parameters + ---------- + lon: np.ndarray + longitude to interpolate + lat: np.ndarray + latitude to interpolate + model_file: str + IERS map of ocean pole tide coefficients + method: str, default 'spline' + Interpolation method + + - ``'spline'``: scipy bivariate spline interpolation + - ``'linear'``, ``'nearest'``: scipy regular grid interpolations + + Returns + ------- + ur: np.ndarray + radial ocean pole tide coefficients + un: np.ndarray + north ocean pole tide coefficients + ue: np.ndarray + east ocean pole tide coefficients + + References + ---------- + .. [1] S. Desai, "Observing the pole tide with satellite altimetry", *Journal of + Geophysical Research: Oceans*, 107(C11), (2002). + `doi: 10.1029/2001JC001224 `_ + .. [2] S. Desai, J. Wahr and B. Beckley "Revisiting the pole tide for and from + satellite altimetry", *Journal of Geodesy*, 89(12), p1233-1243, (2015). + `doi: 10.1007/s00190-015-0848-7 `_ + """ + # default keyword arguments + kwargs.setdefault('model_file', _ocean_pole_tide_file) + kwargs.setdefault('method', 'spline') + # number of points + npts = len(np.atleast_1d(lat)) + # read ocean pole tide map from Desai (2002) + Umap = {} + Umap['R'], Umap['N'], Umap['E'], ilon, ilat = read_binary_file(**kwargs) + # interpolate ocean pole tide map from Desai (2002) + Uint = {} + if (kwargs['method'] == 'spline'): + # use scipy bivariate splines to interpolate to output points + for key,val in Umap.items(): + Uint[key] = np.zeros((npts), dtype=np.clongdouble) + f1 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], + val[:,::-1].real, kx=1, ky=1) + f2 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], + val[:,::-1].imag, kx=1, ky=1) + Uint[key].real = f1.ev(lon, lat) + Uint[key].imag = f2.ev(lon, lat) + else: + # use scipy regular grid to interpolate values for a given method + for key,val in Umap.items(): + Uint[key] = np.zeros((npts), dtype=np.clongdouble) + r1 = scipy.interpolate.RegularGridInterpolator((ilon,ilat[::-1]), + val[:,::-1], bounds_error=False, method=kwargs['method']) + Uint[key][:] = r1.__call__(np.c_[lon, lat]) + # return the interpolated values + return (Uint['R'], Uint['N'], Uint['E']) + +# PURPOSE: read real and imaginary ocean pole tide coefficients +def read_binary_file(**kwargs): + """ + Read real and imaginary ocean pole tide coefficients from [1]_ [2]_ + + Parameters + ---------- + model_file: str or pathlib.Path + IERS map of ocean pole tide coefficients + + Returns + ------- + ur: np.ndarray + radial ocean pole tide coefficients + un: np.ndarray + north ocean pole tide coefficients + ue: np.ndarray + east ocean pole tide coefficients + glon: np.ndarray + ocean grid longitude + glat: np.ndarray + ocean grid latitude + + References + ---------- + .. [1] S. Desai, "Observing the pole tide with satellite altimetry", *Journal of + Geophysical Research: Oceans*, 107(C11), (2002). + `doi: 10.1029/2001JC001224 `_ + .. [2] S. Desai, J. Wahr and B. Beckley "Revisiting the pole tide for and from + satellite altimetry", *Journal of Geodesy*, 89(12), p1233-1243, (2015). + `doi: 10.1007/s00190-015-0848-7 `_ + """ + # default keyword arguments + kwargs.setdefault('model_file', _ocean_pole_tide_file) + # convert input file to tilde-expanded pathlib object + input_file = pathlib.Path(kwargs['model_file']).expanduser() + # check that ocean pole tide file is accessible + if not input_file.exists(): + raise FileNotFoundError(str(input_file)) + + # read GZIP ocean pole tide file + with gzip.open(input_file, 'rb') as f: + file_contents = f.read().splitlines() + + # counts the number of lines in the header + count = 0 + # Reading over header text + HEADER = True + while HEADER: + # file line at count + line = file_contents[count] + # find --------- within line to set HEADER flag to False when found + HEADER = not bool(re.match(rb'---------',line)) + # add 1 to counter + count += 1 + + # grid parameters and dimensions + dlon,dlat = (0.50,0.50) + glon = np.arange(dlon/2.0,360+dlon/2.0,dlon) + glat = np.arange(90.0-dlat/2.0,-90.0-dlat/2.0,-dlat) + nlon = len(glon) + nlat = len(glat) + # allocate for output grid maps + U = {} + U['R'] = np.zeros((nlon,nlat), dtype=np.clongdouble) + U['N'] = np.zeros((nlon,nlat), dtype=np.clongdouble) + U['E'] = np.zeros((nlon,nlat), dtype=np.clongdouble) + # read lines of file and add to output variables + for i,line in enumerate(file_contents[count:]): + ln,lt,urr,uri,unr,uni,uer,uei = np.array(line.split(), dtype='f8') + ilon = int(ln/dlon) + ilat = int((90.0-lt)/dlat) + U['R'][ilon,ilat] = urr + 1j*uri + U['N'][ilon,ilat] = unr + 1j*uni + U['E'][ilon,ilat] = uer + 1j*uei + # shift ocean pole tide grid to -180:180 + for key, val in U.items(): + U[key], glon = _shift(val, glon, lon0=180.0, + cyclic=360.0, direction='west') + # extend matrix for bilinear interpolation + glon = _extend_array(glon,dlon) + # pad ends of matrix for interpolation + for key, val in U.items(): + U[key] = _extend_matrix(val) + # return values + return (U['R'], U['N'], U['E'], glon, glat) + +# PURPOSE: deprecated function to read ocean pole tide coefficients +def ocean_pole_tide(input_file: str | pathlib.Path = _ocean_pole_tide_file): + warnings.warn("Deprecated. Please use pyTMD.io.IERS instead", + DeprecationWarning) + # pass the input file to the read_binary_file function + return read_binary_file(model_file=input_file) + +# PURPOSE: Extend a longitude array +def _extend_array(input_array: np.ndarray, step_size: float): + """ + Extends a longitude array + + Parameters + ---------- + input_array: np.ndarray + array to extend + step_size: float + step size between elements of array + + Returns + ------- + temp: np.ndarray + extended array + """ + n = len(input_array) + temp = np.zeros((n+2), dtype=input_array.dtype) + # extended array [x-1,x0,...,xN,xN+1] + temp[0] = input_array[0] - step_size + temp[1:-1] = input_array[:] + temp[-1] = input_array[-1] + step_size + return temp + +# PURPOSE: Extend a global matrix +def _extend_matrix(input_matrix: np.ndarray): + """ + Extends a global matrix + + Parameters + ---------- + input_matrix: np.ndarray + matrix to extend + + Returns + ------- + temp: np.ndarray + extended matrix + """ + nx,ny = np.shape(input_matrix) + temp = np.zeros((nx+2,ny), dtype=input_matrix.dtype) + temp[0,:] = input_matrix[-1,:] + temp[1:-1,:] = input_matrix[:,:] + temp[-1,:] = input_matrix[0,:] + return temp + +# PURPOSE: shift a grid east or west +def _shift( + input_matrix: np.ndarray, + ilon: np.ndarray, + lon0: int | float = 180, + cyclic: int | float = 360, + direction: str = 'west' + ): + """ + Shift global grid east or west to a new base longitude + + Parameters + ---------- + input_matrix: np.ndarray + input matrix to shift + ilon: np.ndarray + longitude of tidal model + lon0: int or float, default 180 + Starting longitude for shifted grid + cyclic: int or float, default 360 + width of periodic domain + direction: str, default 'west' + Direction to shift grid + + - ``'west'`` + - ``'east'`` + + Returns + ------- + temp: np.ndarray + shifted matrix + lon: np.ndarray + shifted longitude + """ + # find the starting index if cyclic + offset = 0 if (np.fabs(ilon[-1]-ilon[0]-cyclic) > 1e-4) else 1 + i0 = np.argmin(np.fabs(ilon - lon0)) + # shift longitudinal values + lon = np.zeros(ilon.shape, ilon.dtype) + lon[0:-i0] = ilon[i0:] + lon[-i0:] = ilon[offset: i0+offset] + # add or remove the cyclic + if (direction == 'east'): + lon[-i0:] += cyclic + elif (direction == 'west'): + lon[0:-i0] -= cyclic + # allocate for shifted data + if np.ma.isMA(input_matrix): + temp = np.ma.zeros(input_matrix.shape,input_matrix.dtype) + else: + temp = np.zeros(input_matrix.shape, input_matrix.dtype) + # shift data values + temp[:-i0,:] = input_matrix[i0:, :] + temp[-i0:,:] = input_matrix[offset: i0+offset, :] + # return the shifted values + return (temp, lon) diff --git a/pyTMD/io/__init__.py b/pyTMD/io/__init__.py index fafe7f62..2f792cb1 100644 --- a/pyTMD/io/__init__.py +++ b/pyTMD/io/__init__.py @@ -5,6 +5,6 @@ from .FES import * from .GOT import * from .OTIS import * +from .IERS import * from .constituents import constituents from .model import model -from .ocean_pole_tide import ocean_pole_tide diff --git a/pyTMD/io/ocean_pole_tide.py b/pyTMD/io/ocean_pole_tide.py deleted file mode 100644 index 3db2b448..00000000 --- a/pyTMD/io/ocean_pole_tide.py +++ /dev/null @@ -1,185 +0,0 @@ -#!/usr/bin/env python -u""" -ocean_pole_tide.py -Written by Tyler Sutterley (06/2024) - -Reads ocean pole load tide coefficients provided by IERS -http://maia.usno.navy.mil/conventions/2010/2010_official/chapter7/tn36_c7.pdf -http://maia.usno.navy.mil/conventions/2010/2010_update/chapter7/icc7.pdf - -IERS 0.5x0.5 map of ocean pole tide coefficients: -ftp://maia.usno.navy.mil/conventions/2010/2010_update/chapter7/additional_info/ - opoleloadcoefcmcor.txt.gz - -OUTPUTS: - ur: radial ocean pole tide coefficients - un: north ocean pole tide coefficients - ue: east ocean pole tide coefficients - glon: ocean grid longitude - glat: ocean grid latitude - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - -REFERENCES: - S. Desai, "Observing the pole tide with satellite altimetry", Journal of - Geophysical Research: Oceans, 107(C11), 2002. doi: 10.1029/2001JC001224 - S. Desai, J. Wahr and B. Beckley "Revisiting the pole tide for and from - satellite altimetry", Journal of Geodesy, 89(12), p1233-1243, 2015. - doi: 10.1007/s00190-015-0848-7 - -UPDATE HISTORY: - Updated 06/2024: use np.clongdouble instead of np.longcomplex - Updated 05/2023: add default for ocean pole tide file - Updated 04/2023: using pathlib to define and expand paths - Updated 03/2023: add basic variable typing to function inputs - Updated 12/2022: refactor ocean pole tide read programs under io - Updated 04/2022: updated docstrings to numpy documentation format - use longcomplex data format to be windows compliant - Updated 07/2021: added check that ocean pole tide file is accessible - Updated 03/2021: replaced numpy bool/int to prevent deprecation warnings - Updated 08/2020: output north load and east load deformation components - Updated 07/2020: added function docstrings - Updated 12/2018: Compatibility updates for Python3 - Written 09/2017 -""" -from __future__ import annotations - -import re -import gzip -import pathlib -import numpy as np -from pyTMD.utilities import get_data_path - -# ocean pole tide file from Desai (2002) and IERS conventions -_ocean_pole_tide_file = get_data_path(['data','opoleloadcoefcmcor.txt.gz']) - -# PURPOSE: read real and imaginary ocean pole tide coefficients -def ocean_pole_tide(input_file: str | pathlib.Path = _ocean_pole_tide_file): - """ - Read real and imaginary ocean pole tide coefficients - - Parameters - ---------- - input_file: str - IERS map of ocean pole tide coefficients - - Returns - ------- - ur: np.ndarray - radial ocean pole tide coefficients - un: np.ndarray - north ocean pole tide coefficients - ue: np.ndarray - east ocean pole tide coefficients - glon: np.ndarray - ocean grid longitude - glat: np.ndarray - ocean grid latitude - - References - ---------- - .. [1] S. Desai, "Observing the pole tide with satellite altimetry", *Journal of - Geophysical Research: Oceans*, 107(C11), (2002). - `doi: 10.1029/2001JC001224 `_ - .. [2] S. Desai, J. Wahr and B. Beckley "Revisiting the pole tide for and from - satellite altimetry", *Journal of Geodesy*, 89(12), p1233-1243, (2015). - `doi: 10.1007/s00190-015-0848-7 `_ - """ - # convert input file to tilde-expanded pathlib object - input_file = pathlib.Path(input_file).expanduser() - # check that ocean pole tide file is accessible - if not input_file.exists(): - raise FileNotFoundError(str(input_file)) - - # read GZIP ocean pole tide file - with gzip.open(input_file, 'rb') as f: - file_contents = f.read().splitlines() - - # counts the number of lines in the header - count = 0 - # Reading over header text - HEADER = True - while HEADER: - # file line at count - line = file_contents[count] - # find --------- within line to set HEADER flag to False when found - HEADER = not bool(re.match(rb'---------',line)) - # add 1 to counter - count += 1 - - # grid parameters and dimensions - dlon,dlat = (0.50,0.50) - glon = np.arange(dlon/2.0,360+dlon/2.0,dlon) - glat = np.arange(90.0-dlat/2.0,-90.0-dlat/2.0,-dlat) - nlon = len(glon) - nlat = len(glat) - # allocate for output grid maps - ur = np.zeros((nlon,nlat),dtype=np.clongdouble) - un = np.zeros((nlon,nlat),dtype=np.clongdouble) - ue = np.zeros((nlon,nlat),dtype=np.clongdouble) - # read lines of file and add to output variables - for i,line in enumerate(file_contents[count:]): - ln,lt,urr,uri,unr,uni,uer,uei = np.array(line.split(), dtype='f8') - ilon = int(ln/dlon) - ilat = int((90.0-lt)/dlat) - ur[ilon,ilat] = urr + 1j*uri - un[ilon,ilat] = unr + 1j*uni - ue[ilon,ilat] = uer + 1j*uei - - # extend matrix for bilinear interpolation - glon = _extend_array(glon,dlon) - ur = _extend_matrix(ur) - un = _extend_matrix(un) - ue = _extend_matrix(ue) - # return values - return (ur, un, ue, glon, glat) - -# PURPOSE: Extend a longitude array -def _extend_array(input_array: np.ndarray, step_size: float): - """ - Extends a longitude array - - Parameters - ---------- - input_array: np.ndarray - array to extend - step_size: float - step size between elements of array - - Returns - ------- - temp: np.ndarray - extended array - """ - n = len(input_array) - temp = np.zeros((n+2), dtype=input_array.dtype) - # extended array [x-1,x0,...,xN,xN+1] - temp[0] = input_array[0] - step_size - temp[1:-1] = input_array[:] - temp[-1] = input_array[-1] + step_size - return temp - -# PURPOSE: Extend a global matrix -def _extend_matrix(input_matrix: np.ndarray): - """ - Extends a global matrix - - Parameters - ---------- - input_matrix: np.ndarray - matrix to extend - - Returns - ------- - temp: np.ndarray - extended matrix - """ - nx,ny = np.shape(input_matrix) - temp = np.zeros((nx+2,ny), dtype=input_matrix.dtype) - temp[0,:] = input_matrix[-1,:] - temp[1:-1,:] = input_matrix[:,:] - temp[-1,:] = input_matrix[0,:] - return temp diff --git a/pyTMD/predict.py b/pyTMD/predict.py index d9d7f70f..e71a51ea 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -24,6 +24,8 @@ include inference of eps2 and eta2 when predicting from GOT models add keyword argument to allow inferring specific minor constituents use nodal arguments for all non-OTIS model type cases + add load pole tide function that exports in cartesian coordinates + add ocean pole tide function that exports in cartesian coordinates Updated 07/2024: use normalize_angle from pyTMD astro module make number of days to convert tide time to MJD a variable Updated 02/2024: changed class name for ellipsoid parameters to datum @@ -56,9 +58,14 @@ import pyTMD.arguments import pyTMD.astro from pyTMD.crs import datum +import timescale.time +# number of days between the Julian day epoch and MJD +_jd_mjd = 2400000.5 # number of days between MJD and the tide epoch (1992-01-01T00:00:00) _mjd_tide = 48622.0 +# number of days between the Julian day epoch and the tide epoch +_jd_tide = _jd_mjd + _mjd_tide # PURPOSE: Predict tides at single times def map(t: float | np.ndarray, @@ -407,7 +414,7 @@ def equilibrium_tide(t: np.ndarray, lat: np.ndarray): t: np.ndarray time (days relative to January 1, 1992) lat: np.ndarray - latitudes (degrees north) + latitude (degrees north) Returns ------- @@ -473,6 +480,218 @@ def equilibrium_tide(t: np.ndarray, lat: np.ndarray): # return the long-period equilibrium tides return lpet +# PURPOSE: estimate load pole tides in Cartesian coordinates +def load_pole_tide( + t: np.ndarray, + XYZ: np.ndarray, + deltat: float = 0.0, + gamma_0: float = 9.80665, + omega: float = 7.2921151467e-5, + h2: float = 0.6207, + l2: float = 0.0836, + convention: str = '2018' + ): + """ + Estimate load pole tide displacements in Cartesian coordinates + following IERS Convention (2010) guidelines + + Parameters + ---------- + t: np.ndarray + Time (days relative to January 1, 1992) + XYZ: np.ndarray + Cartesian coordinates of the prediction points (meters) + deltat: float or np.ndarray, default 0.0 + time correction for converting to Ephemeris Time (days) + gamma_0: float, default 9.80665 + Normal gravity (m/s^2) + omega: float, default 7.2921151467e-5 + Earth's rotation rate (radians/second) + h2: float, default 0.6207 + Degree-2 Love number of vertical displacement + l2: float, default 0.0836 + Degree-2 Love (Shida) number of horizontal displacement + convention: str, default '2018' + IERS Mean or Secular Pole Convention + + - ``'2003'`` + - ``'2010'`` + - ``'2015'`` + - ``'2018'`` + + Returns + ------- + dxt: np.ndarray + Load pole tide displacements in meters in Cartesian coordinates + """ + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # convert time to Terrestial Time (TT) + tt = t + _jd_tide + deltat + # convert time to Modified Julian Days (MJD) + MJD = tt - _jd_mjd + # convert Julian days to calendar dates + Y,M,D,h,m,s = timescale.time.convert_julian(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) + + # radius of the Earth + radius = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2 + XYZ[:,2]**2) + # geocentric latitude (radians) + latitude = np.arctan(XYZ[:,2] / np.sqrt(XYZ[:,0]**2.0 + XYZ[:,1]**2.0)) + # geocentric colatitude (radians) + theta = (np.pi/2.0 - latitude) + # calculate longitude (radians) + phi = np.arctan2(XYZ[:,1], XYZ[:,0]) + + # 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 + # using the latest definition from IERS Conventions (2010) + mx = px - mpx + my = -(py - mpy) + + # number of points + n = np.maximum(len(time_decimal), len(theta)) + # conversion factors in latitude, longitude, and radial directions + dfactor = np.zeros((n, 3)) + dfactor[:,0] = -l2*atr*(omega**2 * radius**2)/(gamma_0) + dfactor[:,1] = l2*atr*(omega**2 * radius**2)/(gamma_0) + dfactor[:,2] = -h2*atr*(omega**2 * radius**2)/(2.0*gamma_0) + + # calculate pole tide displacements (meters) + S = np.zeros((n, 3)) + # pole tide displacements in latitude, longitude, and radial directions + S[:,0] = dfactor[:,0]*np.cos(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi)) + S[:,1] = dfactor[:,1]*np.cos(theta)*(mx*np.sin(phi) - my*np.cos(phi)) + S[:,2] = dfactor[:,2]*np.sin(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi)) + + # rotation matrix + R = np.zeros((3, 3, n)) + 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) + # rotate displacements to ECEF coordinates + dxt = np.einsum('ti...,jit...->tj...', S, R) + + # return the pole tide displacements + # in Cartesian coordinates + return dxt + +# PURPOSE: estimate ocean pole tides in Cartesian coordinates +def ocean_pole_tide( + t: np.ndarray, + XYZ: np.ndarray, + UXYZ: np.ndarray, + deltat: float = 0.0, + gamma_0: float = 9.780325, + a_axis: float = 6378136.3, + GM: float = 3.986004418e14, + omega: float = 7.2921151467e-5, + rho_w: float = 1025.0, + g2: complex = 0.6870 + 0.0036j, + convention: str = '2018' + ): + """ + Estimate ocean pole tide displacements in Cartesian coordinates + following IERS Convention (2010) guidelines + + Parameters + ---------- + t: np.ndarray + Time (days relative to January 1, 1992) + XYZ: np.ndarray + Cartesian coordinates of the prediction points (meters) + UXYZ: np.ndarray + Ocean pole tide values from Desai (2002) + deltat: float or np.ndarray, default 0.0 + time correction for converting to Ephemeris Time (days) + a_axis: float, default 6378136.3 + Semi-major axis of the Earth (meters) + gamma_0: float, default 9.780325 + Normal gravity (m/s^2) + GM: float, default 3.986004418e14 + Earth's gravitational constant [m^3/s^2] + omega: float, default 7.2921151467e-5 + Earth's rotation rate (radians/second) + rho_w: float, default 1025.0 + Density of sea water [kg/m^3] + g2: complex, default 0.6870 + 0.0036j + Degree-2 Love number differential (1 + k2 - h2) + convention: str, default '2018' + IERS Mean or Secular Pole Convention + + - ``'2003'`` + - ``'2010'`` + - ``'2015'`` + - ``'2018'`` + + Returns + ------- + dxt: np.ndarray + Load pole tide displacements in meters in Cartesian coordinates + """ + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # convert time to Terrestial Time (TT) + tt = t + _jd_tide + deltat + # convert time to Modified Julian Days (MJD) + MJD = tt - _jd_mjd + # convert Julian days to calendar dates + Y,M,D,h,m,s = timescale.time.convert_julian(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) + + # radius of the Earth + radius = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2 + XYZ[:,2]**2) + # geocentric latitude (radians) + latitude = np.arctan(XYZ[:,2] / np.sqrt(XYZ[:,0]**2.0 + XYZ[:,1]**2.0)) + # geocentric colatitude (radians) + theta = (np.pi/2.0 - latitude) + # calculate longitude (radians) + phi = np.arctan2(XYZ[:,1], XYZ[:,0]) + # universal gravitational constant [N*m^2/kg^2] + G = 6.67430e-11 + + # 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 + # using the latest definition from IERS Conventions (2010) + mx = px - mpx + my = -(py - mpy) + + # pole tide displacement factors + Hp = np.sqrt(8.0*np.pi/15.0)*(omega**2 * a_axis**4)/GM + K = 4.0*np.pi*G*rho_w*Hp*a_axis/(3.0*gamma_0) + + # number of points + n = np.maximum(len(time_decimal), len(theta)) + # calculate ocean pole tide displacements (meters) + dxt = np.zeros((n, 3)) + for i in range(3): + dxt[:,i] = K*atr*np.real( + (mx*g2.real + my*g2.imag)*UXYZ.real[:,i] + + (my*g2.real - mx*g2.imag)*UXYZ.imag[:,i]) + + # return the ocean pole tide displacements + # in Cartesian coordinates + return dxt + # get IERS parameters _iers = datum(ellipsoid='IERS', units='MKS') @@ -495,7 +714,7 @@ def solid_earth_tide( t: np.ndarray Time (days relative to January 1, 1992) XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray @@ -616,7 +835,7 @@ def _out_of_phase_diurnal( Parameters ---------- XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray @@ -677,7 +896,7 @@ def _out_of_phase_semidiurnal( Parameters ---------- XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray @@ -745,7 +964,7 @@ def _latitude_dependence( Parameters ---------- XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) SXYZ: np.ndarray Earth-centered Earth-fixed coordinates of the sun (meters) LXYZ: np.ndarray @@ -812,7 +1031,7 @@ def _frequency_dependence_diurnal( Parameters ---------- XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) MJD: np.ndarray Modified Julian Day (MJD) """ @@ -894,7 +1113,7 @@ def _frequency_dependence_long_period( Parameters ---------- XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) MJD: np.ndarray Modified Julian Day (MJD) """ @@ -948,7 +1167,7 @@ def _free_to_mean( Parameters ---------- XYZ: np.ndarray - Cartesian coordinates of the station (meters) + Cartesian coordinates of the prediction points (meters) h2: float or np.ndarray Degree-2 Love number of vertical displacement l2: float or np.ndarray diff --git a/scripts/compute_LPT_displacements.py b/scripts/compute_LPT_displacements.py index 1104aaf3..4dedef65 100644 --- a/scripts/compute_LPT_displacements.py +++ b/scripts/compute_LPT_displacements.py @@ -80,6 +80,8 @@ UPDATE HISTORY: Updated 08/2024: changed from 'geotiff' to 'GTiff' and 'cog' formats + drop use of heights when converting to cartesian coordinates + use prediction function to calculate cartesian tide displacements Updated 07/2024: assert that data type is a known value Updated 06/2024: include attributes in output parquet files Updated 05/2024: use function to reading parquet files to allow @@ -246,64 +248,91 @@ def compute_LPT_displacements(input_file, output_file, # number of time points nt = len(time_decimal) - # degrees and arcseconds to radians + # degrees to radians dtr = np.pi/180.0 - atr = np.pi/648000.0 # earth and physical parameters for ellipsoid units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') - # tidal love number appropriate for the load tide + # tidal love/shida numbers appropriate for the load tide hb2 = 0.6207 + lb2 = 0.0836 - # flatten heights - h = np.ravel(dinput['data']) if ('data' in dinput.keys()) else 0.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(np.ravel(lon), np.ravel(lat), h=h, + X,Y,Z = pyTMD.spatial.to_cartesian(np.ravel(lon), np.ravel(lat), a_axis=units.a_axis, flat=units.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 - # geocentric colatitude and longitude in radians + # geocentric colatitude in radians theta = dtr*(90.0 - latitude_geocentric) - phi = dtr*np.ravel(lon) - # compute normal gravity at spatial location and elevation of points. - # Normal gravity at height h. p. 82, Eqn.(2-215) - gamma_h = units.gamma_h(theta, h) - dfactor = -hb2*atr*(units.omega**2*rr**2)/(2.0*gamma_h) - - # 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) + # compute normal gravity at spatial location + # p. 80, Eqn.(2-199) + gamma_0 = units.gamma_0(theta) # calculate radial displacement at time if (TYPE == 'grid'): Srad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) Srad.mask = np.zeros((ny,nx,nt),dtype=bool) + XYZ = np.c_[X, Y, Z] for i in range(nt): - SRAD = dfactor*np.sin(2.0*theta) * \ - (mx[i]*np.cos(phi) + my[i]*np.sin(phi)) - # reform grid - Srad.data[:,:,i] = np.reshape(SRAD, (ny,nx)) + # calculate load pole tides in cartesian coordinates + dxi = pyTMD.predict.load_pole_tide(ts.tide[i], XYZ, + deltat=ts.tt_ut1[i], + gamma_0=gamma_0, + omega=units.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=units.a_axis, flat=units.flat) + # reshape to output dimensions + Srad.data[:,:,i] = np.reshape(drad, (ny,nx)) Srad.mask[:,:,i] = np.isnan(Srad.data[:,:,i]) elif (TYPE == 'drift'): + # 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=units.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=units.a_axis, flat=units.flat) + # reshape to output dimensions Srad = np.ma.zeros((nt), fill_value=FILL_VALUE) - Srad.data[:] = dfactor*np.sin(2.0*theta) * \ - (mx*np.cos(phi) + my*np.sin(phi)) + Srad.data[:] = drad.copy() Srad.mask = np.isnan(Srad.data) elif (TYPE == 'time series'): Srad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) Srad.mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): - SRAD = dfactor[s]*np.sin(2.0*theta[s]) * \ - (mx*np.cos(phi[s]) + my*np.sin(phi[s])) - Srad.data[s,:] = np.copy(SRAD) + # convert coordinates to column arrays + XYZ = np.repeat(np.c_[X[s], Y[s], Z[s]], nt, axis=0) + # calculate load pole tides in cartesian coordinates + dxi = pyTMD.predict.load_pole_tide(ts.tide, XYZ, + deltat=ts.tt_ut1, + gamma_0=gamma_0[s], + omega=units.omega, + h2=hb2, + l2=lb2, + convention=CONVENTION + ) + # calculate radial component of load pole tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X[s] + dxi[:,0], Y[s] + dxi[:,1], Z[s] + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # reshape to output dimensions + Srad.data[s,:] = drad.copy() Srad.mask[s,:] = np.isnan(Srad.data[s,:]) + # replace invalid data with fill values Srad.data[Srad.mask] = Srad.fill_value diff --git a/scripts/compute_OPT_displacements.py b/scripts/compute_OPT_displacements.py index 788dc4e1..b26ef4bb 100644 --- a/scripts/compute_OPT_displacements.py +++ b/scripts/compute_OPT_displacements.py @@ -81,7 +81,7 @@ crs.py: Coordinate Reference System (CRS) routines spatial.py: utilities for reading and writing spatial data utilities.py: download and management utilities for syncing files - io/ocean_pole_tide.py: read ocean pole load tide map from IERS + io/IERS.py: read ocean pole load tide map from IERS REFERENCES: S. Desai, "Observing the pole tide with satellite altimetry", Journal of @@ -92,6 +92,9 @@ UPDATE HISTORY: Updated 08/2024: changed from 'geotiff' to 'GTiff' and 'cog' formats + drop use of heights when converting to cartesian coordinates + use io function to extract ocean pole tide values at coordinates + use prediction function to calculate cartesian tide displacements Updated 06/2024: include attributes in output parquet files use np.clongdouble instead of np.longcomplex Updated 05/2024: use function to reading parquet files to allow @@ -264,9 +267,8 @@ def compute_OPT_displacements(input_file, output_file, # number of time points nt = len(time_decimal) - # degrees and arcseconds to radians + # degrees to radians dtr = np.pi/180.0 - atr = np.pi/648000.0 # earth and physical parameters for ellipsoid units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') # mean equatorial gravitational acceleration [m/s^2] @@ -276,19 +278,20 @@ def compute_OPT_displacements(input_file, output_file, # tidal love number differential (1 + kl - hl) for pole tide frequencies gamma = 0.6870 + 0.0036j - # flatten heights - h = np.ravel(dinput['data']) if ('data' in dinput.keys()) else 0.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(np.ravel(lon), np.ravel(lat), h=h, + X,Y,Z = pyTMD.spatial.to_cartesian(np.ravel(lon), np.ravel(lat), a_axis=units.a_axis, flat=units.flat) # calculate geocentric latitude and convert to degrees latitude_geocentric = np.arctan(Z / np.sqrt(X**2.0 + Y**2.0))/dtr + npts = len(latitude_geocentric) + # geocentric colatitude and longitude in radians + theta = dtr*(90.0 - latitude_geocentric) + phi = dtr*lon.flatten() # pole tide displacement scale factor Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM K = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis/(3.0*ge) - K1 = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis**3/(3.0*units.GM) # calculate angular coordinates of mean/secular pole at time mpx, mpy, fl = timescale.eop.iers_mean_pole(time_decimal, @@ -300,52 +303,94 @@ def compute_OPT_displacements(input_file, output_file, my = -(py - mpy) # 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((len(latitude_geocentric)),dtype=np.clongdouble) - UR.real = f1.ev(np.ravel(lon), latitude_geocentric) - UR.imag = f2.ev(np.ravel(lon), 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_[np.ravel(lon), latitude_geocentric]) + ur, un, ue = pyTMD.io.IERS.extract_coefficients(lon.flatten(), + latitude_geocentric, method=METHOD) + # convert to cartesian coordinates + R = np.zeros((3, 3, npts)) + 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 radial displacement at time if (TYPE == 'grid'): Urad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) Urad.mask = np.zeros((ny,nx,nt),dtype=bool) + XYZ = np.c_[X, Y, Z] for i in range(nt): - URAD = K*atr*np.real( - (mx[i]*gamma.real + my[i]*gamma.imag)*UR.real + - (my[i]*gamma.real - mx[i]*gamma.imag)*UR.imag + # calculate ocean pole tides in cartesian coordinates + dxi = pyTMD.predict.ocean_pole_tide(ts.tide[i], XYZ, UXYZ, + deltat=ts.tt_ut1[i], + a_axis=units.a_axis, + gamma_0=ge, + GM=units.GM, + omega=units.omega, + rho_w=rho_w, + g2=gamma, + convention=CONVENTION ) - # reform grid - Urad.data[:,:,i] = np.reshape(URAD, (ny,nx)) + # 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=units.a_axis, flat=units.flat) + # reshape to output dimensions + Urad.data[:,:,i] = np.reshape(drad, (ny,nx)) Urad.mask[:,:,i] = np.isnan(Urad.data[:,:,i]) elif (TYPE == 'drift'): - Urad = np.ma.zeros((nt), fill_value=FILL_VALUE) - Urad.data[:] = K*atr*np.real( - (mx*gamma.real + my*gamma.imag)*UR.real + - (my*gamma.real - mx*gamma.imag)*UR.imag + # 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=units.a_axis, + gamma_0=ge, + GM=units.GM, + omega=units.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=units.a_axis, flat=units.flat) + # convert to masked array + Urad = np.ma.zeros((nt), fill_value=FILL_VALUE) + Urad.data[:] = drad.copy() Urad.mask = np.isnan(Urad.data) elif (TYPE == 'time series'): Urad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) Urad.mask = np.zeros((nstation,nt),dtype=bool) for s in range(nstation): - URAD = K*atr*np.real( - (mx*gamma.real + my*gamma.imag)*UR.real[s] + - (my*gamma.real - mx*gamma.imag)*UR.imag[s] + # convert coordinates to column arrays + XYZ = np.repeat(np.c_[X[s], Y[s], Z[s]], nt, axis=0) + uxyz = np.repeat(np.atleast_2d(UXYZ[s,:]), nt, axis=0) + # calculate ocean pole tides in cartesian coordinates + dxi = pyTMD.predict.ocean_pole_tide(ts.tide, XYZ, uxyz, + deltat=ts.tt_ut1, + a_axis=units.a_axis, + gamma_0=ge, + GM=units.GM, + omega=units.omega, + rho_w=rho_w, + g2=gamma, + convention=CONVENTION ) - Urad.data[s,:] = np.copy(URAD) + # calculate radial component of ocean pole tides + dln,dlt,drad = pyTMD.spatial.to_geodetic( + X[s] + dxi[:,0], Y[s] + dxi[:,1], Z[s] + dxi[:,2], + a_axis=units.a_axis, flat=units.flat) + # reshape to output dimensions + Urad.data[s,:] = drad.copy() Urad.mask[s,:] = np.isnan(Urad.data[s,:]) + # replace invalid data with fill values Urad.data[Urad.mask] = Urad.fill_value diff --git a/scripts/compute_SET_displacements.py b/scripts/compute_SET_displacements.py index fda2decf..6be75af3 100644 --- a/scripts/compute_SET_displacements.py +++ b/scripts/compute_SET_displacements.py @@ -100,6 +100,7 @@ UPDATE HISTORY: Updated 08/2024: changed from 'geotiff' to 'GTiff' and 'cog' formats + drop use of heights when converting to cartesian coordinates Updated 07/2024: assert that data type is a known value Updated 06/2024: include attributes in output parquet files Updated 05/2024: use function to reading parquet files to allow @@ -242,10 +243,8 @@ def compute_SET_displacements(input_file, output_file, # earth and physical parameters for ellipsoid units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') - # flatten heights - h = np.ravel(dinput['data']) if ('data' in dinput.keys()) else 0.0 # convert input coordinates to cartesian - X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, h=h, + X, Y, Z = pyTMD.spatial.to_cartesian(lon, lat, a_axis=units.a_axis, flat=units.flat) # compute ephemerides for lunisolar coordinates SX, SY, SZ = pyTMD.astro.solar_ecef(ts.MJD, ephemerides=EPHEMERIDES) @@ -270,9 +269,8 @@ def compute_SET_displacements(input_file, output_file, dln, dlt, drad = pyTMD.spatial.to_geodetic( X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se[:,:,i] = np.reshape(drad - h, (ny,nx)) + # reshape to output dimensions + tide_se[:,:,i] = np.reshape(drad, (ny,nx)) elif (TYPE == 'drift'): # convert coordinates to column arrays XYZ = np.c_[X, Y, Z] @@ -286,9 +284,8 @@ def compute_SET_displacements(input_file, output_file, dln, dlt, drad = pyTMD.spatial.to_geodetic( X + dxi[:,0], Y + dxi[:,1], Z + dxi[:,2], a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se = drad - h + # reshape to output dimensions + tide_se = drad.copy() elif (TYPE == 'time series'): tide_se = np.zeros((nstation,nt)) # convert coordinates to column arrays @@ -303,11 +300,10 @@ def compute_SET_displacements(input_file, output_file, 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], + X[s] + dxi[:,0], Y[s] + dxi[:,1], Z[s] + dxi[:,2], a_axis=units.a_axis, flat=units.flat) - # remove effects of original topography - # (if added when computing cartesian coordinates) - tide_se[s,:] = drad - h + # reshape to output dimensions + tide_se[s,:] = drad.copy() # output netCDF4 and HDF5 file attributes # will be added to YAML header in csv files diff --git a/test/test_ocean_pole_tide.py b/test/test_ocean_pole_tide.py deleted file mode 100755 index 435da295..00000000 --- a/test/test_ocean_pole_tide.py +++ /dev/null @@ -1,126 +0,0 @@ -#!/usr/bin/env python -u""" -test_ocean_pole_tide.py -Written by Tyler Sutterley (06/2024) - -UPDATE HISTORY: - Updated 06/2024: use np.clongdouble instead of np.longcomplex - Updated 02/2024: changed class name for ellipsoid parameters to datum - Updated 04/2023: using pathlib to define and expand paths - Updated 12/2022: single implicit import of pyTMD - use constants class for ellipsoidal parameters - read header from test file and compare more variables - Updated 07/2022: define variable formats of test inputs - Updated 05/2022: change to longcomplex for windows compatibility - Written 08/2020 -""" -import re -import inspect -import pathlib -import pytest -import numpy as np -import scipy.interpolate -import pyTMD - -# current file path -filename = inspect.getframeinfo(inspect.currentframe()).filename -filepath = pathlib.Path(filename).absolute().parent - -# parameterize interpolation method -@pytest.mark.parametrize("METHOD", ['spline','nearest','linear']) -# test the interpolation of ocean pole tide values -def test_ocean_pole_tide(METHOD): - # read ocean pole tide test file for header text - ocean_pole_test_file = filepath.joinpath('opoleloadcmcor.test') - with ocean_pole_test_file.open(mode='r', encoding='utf8') as f: - file_contents = f.read().splitlines() - - # extract header parameters - rx = re.compile(r'(.*?)\s+\=\s+([-+]?\d+\.\d+[e][-+]?\d+)', re.VERBOSE) - header_contents = [l for l in file_contents if rx.match(l)] - test_header = {} - for line in header_contents: - param,contents = rx.findall(line).pop() - test_header[param] = np.float64(contents) - - # extract longitude and latitude from header - lat = re.findall(r'latitude = ([-+]?\d+\.\d+) degrees', file_contents[1]).pop() - lon = re.findall(r'longitude = ([-+]?\d+\.\d+) degrees', file_contents[1]).pop() - test_header['latitude'] = np.float64(lat) - test_header['longitude'] = np.float64(lon) - - # degrees and arcseconds to radians - dtr = np.pi/180.0 - atr = np.pi/648000.0 - # earth and physical parameters for ellipsoid - units = pyTMD.datum(ellipsoid='IERS', units='MKS') - # universal constant of gravitation for test [m^3/(kg*s^2)] - G = 6.673e-11 - # mean equatorial gravitational acceleration [m/s^2] - ge = 9.7803278 - # density of sea water [kg/m^3] - rho_w = 1025.0 - # tidal love number differential (1 + kl - hl) for pole tide frequencies - gamma = 0.6870 + 0.0036j - - # pole tide displacement scale factor - Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM - K = 4.0*np.pi*G*rho_w*Hp*units.a_axis/(3.0*ge) - K1 = 4.0*np.pi*G*rho_w*Hp*units.a_axis**3/(3.0*units.GM) - - # determine differences with constants from test data - eps = np.finfo(np.float16).eps - assert np.isclose(units.a_axis, test_header['aE']) - assert np.isclose(units.omega, test_header['Omega']) - assert np.isclose(units.GM, test_header['GM']) - assert np.isclose(rho_w, test_header['rhow']) - assert np.isclose(gamma.real, test_header['gamma2(real)']) - assert np.isclose(gamma.imag, test_header['gamma2(imag)']) - assert np.isclose(G, test_header['G']) - assert np.isclose(Hp, test_header['Hp']) - assert np.isclose(K, test_header['K']) - - # read test file for values - names = ('MJD','xbar_p','ybar_p','x_p','y_p','m1','m2','u_radial','u_north','u_east') - formats = ('i4','f4','f4','f4','f4','f4','f4','f4','f4','f4') - validation = np.loadtxt(ocean_pole_test_file, skiprows=26, - dtype=dict(names=names, formats=formats)) - file_lines = len(validation) - # mean pole coordinates for test - xmean = np.array([test_header['xmean(t0)'], test_header['xmeandot']]) - ymean = np.array([test_header['ymean(t0)'], test_header['ymeandot']]) - - # read ocean pole tide map from Desai (2002) - ocean_pole_tide_file = pyTMD.utilities.get_data_path( - ['data','opoleloadcoefcmcor.txt.gz']) - iur, iun, iue, ilon, ilat = pyTMD.io.ocean_pole_tide(ocean_pole_tide_file) - - # 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((file_lines),dtype=np.clongdouble) - UR.real = f1.ev(test_header['longitude'], test_header['latitude']) - UR.imag = f2.ev(test_header['longitude'], test_header['latitude']) - 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_[test_header['longitude'], test_header['latitude']]) - - # calculate differences in time - dt = (validation['MJD'] - test_header['t0'])/test_header['1 year'] - # calculate angular coordinates of mean pole at time - mpx = xmean[0] + xmean[1]*dt - mpy = ymean[0] + ymean[1]*dt - - # calculate differentials from mean pole positions - mx = validation['x_p'] - mpx - my = -(validation['y_p'] - mpy) - # calculate radial displacement at time - u_radial = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real + - (my*gamma.real - mx*gamma.imag)*UR.imag) - assert np.all((u_radial - validation['u_radial']) < eps) diff --git a/test/test_pole_tide.py b/test/test_pole_tide.py new file mode 100755 index 00000000..c6662682 --- /dev/null +++ b/test/test_pole_tide.py @@ -0,0 +1,494 @@ +#!/usr/bin/env python +u""" +test_pole_tide.py +Written by Tyler Sutterley (08/2024) + +UPDATE HISTORY: + Updated 08/2024: add tests for new cartesian pole tides + Updated 06/2024: use np.clongdouble instead of np.longcomplex + Updated 02/2024: changed class name for ellipsoid parameters to datum + Updated 04/2023: using pathlib to define and expand paths + Updated 12/2022: single implicit import of pyTMD + use constants class for ellipsoidal parameters + read header from test file and compare more variables + Updated 07/2022: define variable formats of test inputs + Updated 05/2022: change to longcomplex for windows compatibility + Written 08/2020 +""" +import re +import inspect +import pathlib +import pytest +import pyproj +import numpy as np +import scipy.interpolate +import pyTMD.compute +import timescale.time + +# current file path +filename = inspect.getframeinfo(inspect.currentframe()).filename +filepath = pathlib.Path(filename).absolute().parent + +# number of days between the Julian day epoch and MJD +_jd_mjd = 2400000.5 +# number of days between MJD and the tide epoch (1992-01-01T00:00:00) +_mjd_tide = 48622.0 +# number of days between the Julian day epoch and the tide epoch +_jd_tide = _jd_mjd + _mjd_tide + +# PURPOSE: compute radial load pole tide displacements +# following IERS Convention (2010) guidelines +@pytest.mark.parametrize("TYPE", ['grid','drift','time series']) +def test_load_pole_tide_displacements(TYPE): + + # create a test dataset for data type + if (TYPE == 'drift'): + # number of data points + n_time = 3000 + y = np.random.randint(-90,90,size=n_time) + x = np.random.randint(-180,180,size=n_time) + delta_time = np.random.randint(0,31557600,size=n_time) + elif (TYPE == 'grid'): + # number of data points + n_lat,n_lon,n_time = (181,361,100) + y = np.linspace(-90,90,n_lat) + x = np.linspace(-180,180,n_lon) + delta_time = np.random.randint(0,31557600,size=n_time) + elif (TYPE == 'time series'): + n_station,n_time = (300,100) + y = np.random.randint(-90,90,size=n_station) + x = np.random.randint(-180,180,size=n_station) + delta_time = np.random.randint(0,31557600,size=n_time) + + # parameters for ocean pole tide + EPSG = 4326 + EPOCH = (2018, 1, 1, 0, 0, 0) + TIME = 'UTC' + ELLIPSOID = 'WGS84' + CONVENTION = '2018' + FILL_VALUE = np.nan + + # compute load pole tides using cartesian functions + test = pyTMD.compute.LPT_displacements(x, y, delta_time, + EPSG=EPSG, EPOCH=EPOCH, TYPE=TYPE, TIME=TIME, + ELLIPSOID=ELLIPSOID, CONVENTION=CONVENTION, + FILL_VALUE=FILL_VALUE) + + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon,lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + if (TIME.lower() == 'datetime'): + ts = timescale.time.Timescale().from_datetime( + delta_time.flatten()) + else: + ts = timescale.time.Timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + + # convert dynamic time to Modified Julian Days (MJD) + MJD = ts.tt - _jd_mjd + # 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) + # number of time points + nt = len(time_decimal) + + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # earth and physical parameters for ellipsoid + units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') + # tidal love number appropriate for the load tide + hb2 = 0.6207 + + # 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.flatten(), lat.flatten(), + a_axis=units.a_axis, flat=units.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 + # geocentric colatitude and longitude in radians + theta = dtr*(90.0 - latitude_geocentric) + phi = dtr*lon.flatten() + + # compute normal gravity at spatial location + gamma_0 = units.gamma_0(theta) + dfactor = -hb2*atr*(units.omega**2*rr**2)/(2.0*gamma_0) + + # 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 + if (TYPE == 'grid'): + ny,nx = np.shape(x) + Srad = np.ma.zeros((ny,nx,nt), fill_value=FILL_VALUE) + Srad.mask = np.zeros((ny,nx,nt),dtype=bool) + approx = np.zeros((ny,nx,nt)) + for i in range(nt): + SRAD = dfactor*np.sin(2.0*theta)*(mx[i]*np.cos(phi)+my[i]*np.sin(phi)) + # reform grid + Srad.data[:,:,i] = np.reshape(SRAD, (ny, nx)) + Srad.mask[:,:,i] = np.isnan(Srad.data[:,:,i]) + # approximate values from IERS (2010) conventions + S = -33.0*np.sin(2.0*theta)*(mx[i]*np.cos(phi) + my[i]*np.sin(phi)) + approx[:,:,i] = np.reshape(S, (ny, nx))/1e3 + elif (TYPE == 'drift'): + Srad = np.ma.zeros((nt), fill_value=FILL_VALUE) + Srad.data[:] = dfactor*np.sin(2.0*theta)*(mx*np.cos(phi)+my*np.sin(phi)) + Srad.mask = np.isnan(Srad.data) + # approximate values from IERS (2010) conventions + S = -33.0*np.sin(2.0*theta)*(mx*np.cos(phi) + my*np.sin(phi)) + approx = S/1e3 + elif (TYPE == 'time series'): + nstation = len(x) + Srad = np.ma.zeros((nstation,nt), fill_value=FILL_VALUE) + Srad.mask = np.zeros((nstation,nt),dtype=bool) + approx = np.zeros((nstation,nt)) + for s in range(nstation): + SRAD = dfactor[s]*np.sin(2.0*theta[s])*(mx*np.cos(phi[s])+my*np.sin(phi[s])) + Srad.data[s,:] = np.copy(SRAD) + Srad.mask[s,:] = np.isnan(Srad.data[s,:]) + # approximate values from IERS (2010) conventions + S = -33.0*np.sin(2.0*theta[s])*(mx*np.cos(phi[s]) + my*np.sin(phi[s])) + approx[s,:] = np.copy(S)/1e3 + # replace invalid data with fill values + Srad.data[Srad.mask] = Srad.fill_value + # compare with functional values + eps = np.finfo(np.float16).eps + assert np.all(np.abs(Srad - test) < eps) + assert np.all(np.abs(approx - test) < eps) + +# parameterize interpolation method +@pytest.mark.parametrize("METHOD", ['spline','nearest','linear']) +# test the interpolation of ocean pole tide values +def test_read_ocean_pole(METHOD): + # read ocean pole tide test file for header text + ocean_pole_test_file = filepath.joinpath('opoleloadcmcor.test') + with ocean_pole_test_file.open(mode='r', encoding='utf8') as f: + file_contents = f.read().splitlines() + + # extract header parameters + rx = re.compile(r'(.*?)\s+\=\s+([-+]?\d+\.\d+[e][-+]?\d+)', re.VERBOSE) + header_contents = [l for l in file_contents if rx.match(l)] + header = {} + for line in header_contents: + param,contents = rx.findall(line).pop() + header[param] = np.float64(contents) + + # extract longitude and latitude from header + lat = re.findall(r'latitude = ([-+]?\d+\.\d+) degrees', file_contents[1]).pop() + lon = re.findall(r'longitude = ([-+]?\d+\.\d+) degrees', file_contents[1]).pop() + header['latitude'] = np.float64(lat) + # convert longitude from 0:360 to -180:180 + header['longitude'] = np.float64(lon) - 360.0 + + # number of points + npts = len(np.atleast_1d(header['longitude'])) + # read ocean pole tide map from Desai (2002) + Umap = {} + Umap['R'], Umap['N'], Umap['E'], ilon, ilat = pyTMD.io.IERS.read_binary_file() + # interpolate ocean pole tide map from Desai (2002) + Uint = {} + if (METHOD == 'spline'): + # use scipy bivariate splines to interpolate to output points + for key,val in Umap.items(): + Uint[key] = np.zeros((npts), dtype=np.clongdouble) + f1 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], + val[:,::-1].real, kx=1, ky=1) + f2 = scipy.interpolate.RectBivariateSpline(ilon, ilat[::-1], + val[:,::-1].imag, kx=1, ky=1) + Uint[key].real = f1.ev(header['longitude'], header['latitude']) + Uint[key].imag = f2.ev(header['longitude'], header['latitude']) + else: + # use scipy regular grid to interpolate values for a given method + for key,val in Umap.items(): + Uint[key] = np.zeros((npts), dtype=np.clongdouble) + r1 = scipy.interpolate.RegularGridInterpolator((ilon,ilat[::-1]), + val[:,::-1], bounds_error=False, method=METHOD) + Uint[key][:] = r1.__call__(np.c_[header['longitude'], header['latitude']]) + + # extract coefficients from IERS pole tide map + U = pyTMD.io.IERS.extract_coefficients( + header['longitude'], header['latitude'], + method=METHOD) + + # compare with functional values + for i, key in enumerate(['R','N','E']): + assert np.isclose(Uint[key], U[i]) + +# parameterize interpolation method +@pytest.mark.parametrize("METHOD", ['spline','nearest','linear']) +# test the interpolation of ocean pole tide values +def test_ocean_pole_tide(METHOD): + # read ocean pole tide test file for header text + ocean_pole_test_file = filepath.joinpath('opoleloadcmcor.test') + with ocean_pole_test_file.open(mode='r', encoding='utf8') as f: + file_contents = f.read().splitlines() + + # extract header parameters + rx = re.compile(r'(.*?)\s+\=\s+([-+]?\d+\.\d+[e][-+]?\d+)', re.VERBOSE) + header_contents = [l for l in file_contents if rx.match(l)] + header = {} + for line in header_contents: + param,contents = rx.findall(line).pop() + header[param] = np.float64(contents) + + # extract longitude and latitude from header + lat = re.findall(r'latitude = ([-+]?\d+\.\d+) degrees', file_contents[1]).pop() + lon = re.findall(r'longitude = ([-+]?\d+\.\d+) degrees', file_contents[1]).pop() + header['latitude'] = np.float64(lat) + # convert longitude from 0:360 to -180:180 + header['longitude'] = np.float64(lon) - 360.0 + + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # earth and physical parameters for ellipsoid + units = pyTMD.datum(ellipsoid='IERS', units='MKS') + # universal constant of gravitation for test [m^3/(kg*s^2)] + G = 6.673e-11 + # mean equatorial gravitational acceleration [m/s^2] + ge = 9.7803278 + # density of sea water [kg/m^3] + rho_w = 1025.0 + # tidal love number differential (1 + kl - hl) for pole tide frequencies + gamma = 0.6870 + 0.0036j + + # pole tide displacement scale factor + Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM + K = 4.0*np.pi*G*rho_w*Hp*units.a_axis/(3.0*ge) + + # determine differences with constants from test data + eps = np.finfo(np.float16).eps + assert np.isclose(units.a_axis, header['aE']) + assert np.isclose(units.omega, header['Omega']) + assert np.isclose(units.GM, header['GM']) + assert np.isclose(rho_w, header['rhow']) + assert np.isclose(gamma.real, header['gamma2(real)']) + assert np.isclose(gamma.imag, header['gamma2(imag)']) + assert np.isclose(G, header['G']) + assert np.isclose(Hp, header['Hp']) + assert np.isclose(K, header['K']) + + # read test file for values + names = ('MJD','xbar_p','ybar_p','x_p','y_p','m1','m2','u_radial','u_north','u_east') + formats = ('i4','f4','f4','f4','f4','f4','f4','f4','f4','f4') + validation = np.loadtxt(ocean_pole_test_file, skiprows=26, + dtype=dict(names=names, formats=formats)) + file_lines = len(validation) + # mean pole coordinates for test + xmean = np.array([header['xmean(t0)'], header['xmeandot']]) + ymean = np.array([header['ymean(t0)'], header['ymeandot']]) + + # read ocean pole tide map from Desai (2002) + ocean_pole_tide_file = pyTMD.utilities.get_data_path( + ['data','opoleloadcoefcmcor.txt.gz']) + iur, iun, iue, ilon, ilat = pyTMD.io.ocean_pole_tide(ocean_pole_tide_file) + + # 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((file_lines),dtype=np.clongdouble) + UR.real = f1.ev(header['longitude'], header['latitude']) + UR.imag = f2.ev(header['longitude'], header['latitude']) + 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_[header['longitude'], header['latitude']]) + + # calculate differences in time + dt = (validation['MJD'] - header['t0'])/header['1 year'] + # calculate angular coordinates of mean pole at time + mpx = xmean[0] + xmean[1]*dt + mpy = ymean[0] + ymean[1]*dt + + # calculate differentials from mean pole positions + mx = validation['x_p'] - mpx + my = -(validation['y_p'] - mpy) + # calculate radial displacement at time + u_radial = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real + + (my*gamma.real - mx*gamma.imag)*UR.imag) + assert np.all(np.abs(u_radial - validation['u_radial']) < eps) + +# PURPOSE: compute radial load pole tide displacements +# following IERS Convention (2010) guidelines +@pytest.mark.parametrize("TYPE", ['grid','drift','time series']) +@pytest.mark.parametrize("METHOD", ['spline','nearest','linear']) +def test_ocean_pole_tide_displacements(TYPE, METHOD): + + # create a test dataset for data type + if (TYPE == 'drift'): + # number of data points + n_time = 3000 + y = np.random.randint(-90,90,size=n_time) + x = np.random.randint(-180,180,size=n_time) + delta_time = np.random.randint(0,31557600,size=n_time) + elif (TYPE == 'grid'): + # number of data points + n_lat,n_lon,n_time = (181,361,100) + y = np.linspace(-90,90,n_lat) + x = np.linspace(-180,180,n_lon) + delta_time = np.random.randint(0,31557600,size=n_time) + elif (TYPE == 'time series'): + n_station,n_time = (300,100) + y = np.random.randint(-90,90,size=n_station) + x = np.random.randint(-180,180,size=n_station) + delta_time = np.random.randint(0,31557600,size=n_time) + + # parameters for ocean pole tide + EPSG = 4326 + EPOCH = (2018, 1, 1, 0, 0, 0) + TIME = 'UTC' + ELLIPSOID = 'WGS84' + CONVENTION = '2018' + FILL_VALUE = np.nan + + # compute load pole tides using cartesian functions + test = pyTMD.compute.OPT_displacements(x, y, delta_time, + EPSG=EPSG, EPOCH=EPOCH, TYPE=TYPE, TIME=TIME, + ELLIPSOID=ELLIPSOID, CONVENTION=CONVENTION, + METHOD=METHOD, FILL_VALUE=FILL_VALUE) + + # reform coordinate dimensions for input grids + # or verify coordinate dimension shapes + if (TYPE.lower() == 'grid') and (np.size(x) != np.size(y)): + x,y = np.meshgrid(np.copy(x),np.copy(y)) + elif (TYPE.lower() == 'grid'): + x = np.atleast_2d(x) + y = np.atleast_2d(y) + elif TYPE.lower() in ('time series', 'drift'): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + # converting x,y from EPSG to latitude/longitude + crs1 = pyTMD.crs().from_input(EPSG) + crs2 = pyproj.CRS.from_epsg(4326) + transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) + lon,lat = transformer.transform(x.flatten(), y.flatten()) + + # assert delta time is an array + delta_time = np.atleast_1d(delta_time) + # convert delta times or datetimes objects to timescale + ts = timescale.time.Timescale().from_deltatime(delta_time, + epoch=EPOCH, standard=TIME) + + # convert dynamic time to Modified Julian Days (MJD) + MJD = ts.tt - _jd_mjd + # 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) + # number of time points + nt = len(time_decimal) + + # degrees and arcseconds to radians + dtr = np.pi/180.0 + atr = np.pi/648000.0 + # earth and physical parameters for ellipsoid + units = pyTMD.datum(ellipsoid=ELLIPSOID, units='MKS') + # mean equatorial gravitational acceleration [m/s^2] + ge = 9.7803278 + # density of sea water [kg/m^3] + rho_w = 1025.0 + # tidal love number differential (1 + kl - hl) for pole tide frequencies + gamma = 0.6870 + 0.0036j + + # 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.flatten(), lat.flatten(), + a_axis=units.a_axis, flat=units.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 + + # pole tide displacement scale factor + Hp = np.sqrt(8.0*np.pi/15.0)*(units.omega**2*units.a_axis**4)/units.GM + K = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis/(3.0*ge) + K1 = 4.0*np.pi*units.G*rho_w*Hp*units.a_axis**3/(3.0*units.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) + + # read ocean pole tide map from Desai (2002) + iur, iun, iue, ilon, ilat = pyTMD.io.IERS.read_binary_file() + # 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((len(latitude_geocentric)), dtype=np.clongdouble) + UR.real = f1.ev(lon.flatten(), latitude_geocentric) + UR.imag = f2.ev(lon.flatten(), latitude_geocentric) + else: + # use scipy regular grid to interpolate values for a given method + r1 = scipy.interpolate.RegularGridInterpolator((ilon,ilat[::-1]), + iur[:,::-1], bounds_error=False, method=METHOD) + UR = r1.__call__(np.c_[lon.flatten(), latitude_geocentric]) + + # calculate radial displacement at time + if (TYPE == 'grid'): + ny,nx = np.shape(x) + Urad = np.ma.zeros((ny,nx,nt),fill_value=FILL_VALUE) + Urad.mask = np.zeros((ny,nx,nt),dtype=bool) + for i in range(nt): + URAD = K*atr*np.real((mx[i]*gamma.real + my[i]*gamma.imag)*UR.real + + (my[i]*gamma.real - mx[i]*gamma.imag)*UR.imag) + # reform grid + Urad.data[:,:,i] = np.reshape(URAD, (ny,nx)) + Urad.mask[:,:,i] = np.isnan(Urad.data[:,:,i]) + elif (TYPE == 'drift'): + Urad = np.ma.zeros((nt),fill_value=FILL_VALUE) + Urad.data[:] = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real + + (my*gamma.real - mx*gamma.imag)*UR.imag) + Urad.mask = np.isnan(Urad.data) + elif (TYPE == 'time series'): + nstation = len(x) + Urad = np.ma.zeros((nstation,nt),fill_value=FILL_VALUE) + Urad.mask = np.zeros((nstation,nt),dtype=bool) + for s in range(nstation): + URAD = K*atr*np.real((mx*gamma.real + my*gamma.imag)*UR.real[s] + + (my*gamma.real - mx*gamma.imag)*UR.imag[s]) + Urad.data[s,:] = np.copy(URAD) + Urad.mask[s,:] = np.isnan(Urad.data[s,:]) + # replace invalid data with fill values + Urad.data[Urad.mask] = Urad.fill_value + # compare with functional values + eps = np.finfo(np.float16).eps + assert np.all(np.abs(Urad - test) < eps)