Syntax and codes Amin
the function uses data-driven approach published and explained in Vishwakarma et al. 2017, AGU WRR paper. doi:10.1002/2017WR021150.
F a cell matrix with one column containing SH coefficients
cf the column in F that contains SH coefficients from GRACE
GaussianR radius of the Gaussian filter (recommened = 400)
basins mask functions of basin, a cell data format with one column and each entry is a 360 x 720 matrix with 1 inside the catchment and 0 outside
every output has a size (number of months x basins)
RecoveredTWS corrected data-driven time-series (Least Squares fit method)
RecoveredTWS2 corrected data-driven time-series (shift and amplify method)
FilteredTS gaussian filtered GRACE TWS time-series for all the basins.
gaussian, cs2sc, gshs, gsha, naninterp, Phase_calc
Global Spherical Harmonic Synthesis
field set of SH coefficients, either in SC-triangle or CS-square format
quant [‘none’, ‘water’]
grd [‘pole’, ‘mesh’, ‘block’, ‘cell’]
n grid size parameter n. (default: n = lmax, determined from field)
longitude samples: 2*n
latitude samples n ('blocks') or n+1.
h (default: 0), height above Earth mean radius [m]
jflag (default: true), determines whether to subtract GRS80.
gshs Global spheric harmonic synthesis
cs2sc, normalklm, plm, eigengrav, ispec
calculates the spherical harmonic coefficients for a gridded global field by:
• least squares
• weighted least squares
• approximate quadrature
• first Neumann method
• second Neumann method
• block mean values
f set of SH coefficients, either in SC-triangle or CS-square format
quant global field of size (lmax+1)2lmax or lmax2lmax
method string argument, defining the analysis method:
ls least squares
wls weighted least squares
aq approximate quadrature
fnm first Neumann method
snm second Neumann method
mean block mean values
grd optional string argument, defining the grid:
pole, mesh (default if lmax+1), equi-angular (lmax+1)2lmax, including poles and Greenwich meridian.
block, cell (default if lmax), equi-angular block midpoints lmax2lmax
neumann, gauss Gauss-Neumann grid (lmax+1)2lmax; where lmax is the maximum degree of development
longitude samples: 2n
latitude samples n ('blocks') or n+1.
h (default: 0), height above Earth mean radius [m]
jflag (default: true), determines whether to subtract GRS80.
cs Clm, Slm in |C\S| format
plm, iplm, neumann, sc2cs
returns the weights and nodes for Neumann's numerical integration
NB1 In 1st N-method, length(x) should not become too high, since a linear system of this size is solved. Eg: 500.
NB2 No use is made of potential symmetries of nodes.
inn input; can either be n -> number of weights and number of base points; or x -> base points (nodes) in the interval [-1;1]
If the input is n, we apply 1st Neumann method
If the input is x, we apply 2nd Neumann method
w quadrature weights
x base points (nodes) in the interval [-1;1]
plm, grule
NORMALKLM returns an ellipsoidal normal field consisting of normalized -Jn, n=0,2,4,6,8
Reference: J2,J4 values for hydrostatic equilibrium ellipsoid from Lambeck (1988) "Geophysical Geodesy", p.18
lmax maximum degree
typ (default: wgs84) either wgs84 (equipotential ellipsoid), grs80,or he (hydrostatic equilibrium ellipsoid)
nklm normal field in CS-format (sparse)
The function converts the rectangular (L+1)x(2L+1) matrix field, containing spherical harmonics coefficients in /S|C\ storage format into a square (L+1)x(L+1) matrix in |C\S| format. (inverse of cs2sc)
field the rectangular (L+1)x(2L+1) matrix FIELD, containing spherical harmonics coefficients in /S|C\ storage format
cs square (L+1)x(L+1) matrix in |C\S| format
The function converts the storage format into a square (L+1)x(L+1) matrix in |C\S| format to rectangular (L+1)x(2L+1) matrix field, containing spherical harmonics coefficients in /S|C\ . (inverse of sc2cs)
field the square (L+1)x(L+1) matrix FIELD , containing spherical harmonics coefficients in |C\S| storage format
cs rectangular (L+1)x(2L+1) matrix in /S|C\format
The function computes Gauss base points and weight factors using the algorithm given by Davis and Rabinowitz in 'Methods of Numerical Integration', page 365, Academic Press, 1975.
n number of base points required
bp A numpy array containing cosines of the base points
wf A numpy array containing weight factor for computing integrals, etc.
Computes the phase difference between two timeseries based on the Hilbert transform method by Philip et al. (2012) doi:10.1029/2012GL052495. The function is called by the Grace_Data_Driven_Correction function based on the paper Vishwakarma et al. 2017, AGU WRR paper. doi:10.1002/2017WR021150.
fts filtered monthly fields
ffts twice-filtered monthly fields
ps phase difference between filtered monthly fields and twice-filtered monthly fields
ispec(a,b) returns the function F from the spectra a to b
a cosine coefficients
b (default: b = -9999) sine coefficients
a and b are defined by:
f(t) = a_0 + SUM_(i=1)^n2 a_i*cos(iwt) + b_i*sin(iwt)
with w = ground-frequency and n2 = half the number of samples (+1).
Note that no factor 2 appears in front of the sum.
F = ispec(a,b) considers A the cosine- and B the sine-spectrum.
F = ispec(s) assumes s is a list of [a,b]
If A and B are matrices, Fourier operations are columnwise.
Delivers the spherical harmonic coefficients of a gaussian smoothing filter. The coefficients are calculates according to Wahr et. al. (1998) equation (34) and Swenson and Wahr equation (34)
L (integer) maximum degree
cap (integer) half width of Gaussian smoothing function [km]
Wl [L+1 x 1] smoothing coefficients
Interpolates NaNs in the input timeseries using scipy PchipInterpolator (cubic interpolation)
X numpy timeseries with some NaN values
X input numpy timeseries with NaN values filled using scipy PchipInterpolator (cubic interpolation)
A python file consisting of frequently used constants for GRACE signal processing.
clight = 2.99792458e8 speed of light
G = 6.67259e-11 gravitational constant [m^3 /(kg s^2)]
au = 149.597870691e9 astronomical unit [m]
GRS80 defining constants
ae = 6378137 semi-major axis of ellipsoid [m]
GM = 3.986005e14 geocentric grav. constant [m^3 / s^2]
J2 = 1.08263e-3 earth's dyn. form factor (= -C20 unnormalized)
Omega = 7.292115e-5 mean ang. velocity [rad/s]
GRS80 derived constants
flat = 1/298.257222101 flattening
J4 = -0.237091222e-5 -C40 unnormalized
J6 = 0.608347e-8 -C60 unnormalized
J8 = -0.1427e-10 -C80 unnormalized