Auxillary Codes¶

The rest of the important functions have been bundled under the auxillary codes module.

Computing Legendre functions¶

The

`PLM(l, m, thetaRAD, nargin, nargout)` ¶

PLM Fully normalized associated Legendre functions for a selected order M

Parameters:

Name	Type	Description	Default
`l`	`array`	Degree, but not necessarily monotonic. For l < m a vector of zeros will be returned.	required
`m`	`int`	order (scalar). If absent, m = 0 is assumed.	required
`thetaRAD`	`array`	co-latitude [rad] (vector)	required
`nargin`	`int`	number of input argument	required
`nargout`	`int`	number of output argument	required

Returns: (np.array): PLM fully normalized

Source code in pyshbundle/plm.py

def PLM(l: np.array, m:int, thetaRAD, nargin, nargout): 
    """PLM Fully normalized associated Legendre functions for a selected order M

    Args:
        l (np.array): Degree, but not necessarily monotonic.
               For l < m a vector of zeros will be returned.
        m (int): order (scalar). If absent, m = 0 is assumed.
        thetaRAD (np.array): co-latitude [rad] (vector)
        nargin (int): number of input argument
        nargout (int): number of output argument
    Returns:
        (np.array): PLM fully normalized
    """

    if  min(l.shape) != 1:
        print('Degree l must be a vector (or scalar)') 
        sys.exit([])
    if  np.remainder(l,1).all() != 0:
        print('Vector l contains non-integers !!') 
        sys.exit([])
    # if 
    #     print('Order m must be a scalar') 
    #     sys.exit([])
    if  np.remainder(m,1) != 0:
        print('Order must be integer')
        sys.exit([])

# PRELIMINARIES
    lcol = len(l)
    trow = len(thetaRAD)
    lmax = int(max(l[0,:]))

    if lmax < m:
        p = np.zeros([len(thetaRAD), len(l)], dtype='longdouble')
        dp = np.zeros([len(thetaRAD), len(l)], dtype='longdouble')
        ddp = np.zeros([len(thetaRAD), len(l)], dtype='longdouble')
        sys.exit([])

    n = thetaRAD.size                                                               # number of latitudes
    t = thetaRAD[:]
    x = np.cos(t)
    y = np.sin(t)
    lvec = np.transpose(l)   
    lvec = np.intc(lvec)                                                       # l can now be used as running index

    if min(t).all() < 0 and max(t).all() > np.pi:
        print('Warning: Is co-latitude in radians ?')

    # Recursive computation of the temporary matrix ptmp, containing the Legendre
    # functions in its columns, with progressing degree l. The last column of
    # ptmp will contain zeros, which is useful for assignments when l < m.
    ptmp = np.zeros((n,lmax + 2 - m))
    if nargout >= 2:                                                                #  first derivative needs also P_{n,m+1} and P_{n,m-1}
        ptmp_m1 = np.zeros((n,lmax + 3 - m), dtype='longdouble')
        ptmp_p1 = np.zeros((n,lmax + 1 -m), dtype='longdouble')        
        dptmp = np.zeros((n,lmax + 2 - m), dtype='longdouble') 
    if nargout == 3:                                                                # second derivative needs also dP_{n,m+1} and dP_{n,m-1}
        dptmp_m1 = np.zeros((n,lmax + 3 -m), dtype='longdouble')
        dptmp_p1 = np.zeros((n,lmax + 1 -m), dtype='longdouble')
        ptmp_m2 = np.zeros((n,lmax + 4 -m), dtype='longdouble')                                         # but these first derivative need dP_{n,m+2} and dP_{n,m-2}
        ptmp_p2 = np.zeros((n,lmax - m), dtype='longdouble')
        ddptmp = np.zeros((n,lmax + 2 -m), dtype='longdouble')

    # sectorial recursion: PM (non-recursive, though)
    ptmp[:,0] = secrecur(m,y)
    if nargout >= 2:                                                                # frist derivative needs preceding and subsequent element
        if m > 0:    
            ptmp_m1[:,0] = secrecur(m-1,y)                                          # preceding elements
        if m < lmax: 
            ptmp_p1[:,0] = secrecur(m+1,y)                                          # subsequent elemtens
    if nargout == 3:                                                                # second derivative needs P_{n,m+2} and P_{n,m-2} as well
        if m > 1:           
            ptmp_m2[:,0] = secrecur(m-2,y)                                          # preceding elements
        if m < lmax-1: 
            ptmp_p2[:,0] = secrecur(m+2,y)                                          # subsequent elemtens

    # l-recursion: P
    ptmp = lrecur(ptmp,x,m,lmax);
    if nargout >= 2:                                                                # frist derivative needs preceding and subsequent element
        if m > 0:
            ptmp_m1 = lrecur(ptmp_m1,x,m-1,lmax)                                    # preceding elements
        if m < lmax:
            ptmp_p1 = lrecur(ptmp_p1,x,m+1,lmax)                                    # subsequent elemtens

    if nargout == 3:                                                                # second derivative needs P_{n,m+2} and P_{n,m-2} as well
        if m > 1:
            ptmp_m2 = lrecur(ptmp_m2,x,m-2,lmax)                                    # preceding elements
        if m < lmax-1:
            ptmp_p2 = lrecur(ptmp_p2,x,m+2,lmax)                                    # subsequent elemtens

    # now compute the derivatives 
    if nargout >= 2:                                                                # first derivative
        dptmp = derivALF(dptmp,ptmp_m1,ptmp_p1,m,lmax)
    if nargout == 3:                                                                # second derivative
        if m > 0:    
            dptmp_m1 = derivALF(dptmp_m1,ptmp,ptmp_m2,m-1,lmax)
        if m < lmax:
            dptmp_p1 = derivALF(dptmp_p1,ptmp,ptmp_p2,m+1,lmax)
        ddptmp = derivALF(ddptmp,dptmp_m1,dptmp_p1,m,lmax)


    # %--------------------------------------------------------------------
    # % The Legendre functions have been computed. What remains to be done, is to
    # % extract the proper columns from ptmp, corresponding to the vector lvec. 
    # % If l or thetaRAD is scalar the output matrix p reduces to a vector. It should
    # % have the shape of respectively thetaRAD or l in that case.
    # %--------------------------------------------------------------------
    lind       = (lvec < m)   	 # index into l < m
    pcol       = lvec - m + 0			                                            # index into columns of ptmp
    pcol[lind] = np.ndarray((lmax-m+2-6)*np.ones((sum(sum(lind)),1)))	            # Now l < m points to last col.
    p      = ptmp[:,pcol]			                                                # proper column extraction 
    if nargout >= 2:
        dp =  dptmp[:,pcol]                                                         # proper column extraction 
    if nargout == 3: 
        ddp = ddptmp[:,pcol]                                                        # proper column extraction  
    if max(lvec.shape)==1  and min(thetaRAD.shape)==1 and (trow == 1):
        p = p.T
        if nargout >= 2:
            dp  = np.transpose(dp)
        if nargout == 3:
            ddp = np.transpose(ddp)
    if max(thetaRAD.shape)==1 and min(lvec.shape)==1  and (lcol == 1):
        p = p.T
        if nargout >= 2:
            dp  = dp.T  
        if nargout == 3:
            ddp = ddp.T

    if nargout == 1: 
        return p
    if nargout == 2: 
        return p,dp
    if nargout == 3: 
        return p, dp, ddp

`derivALF(inn, miin, plin, m, lmax)` ¶

HelpeFunction

Parameters:

Name	Type	Description	Default
`inn`	`_type_`	description	required
`miin`	`_type_`	description	required
`plin`	`_type_`	description	required
`m`	`_type_`	description	required
`lmax`	`_type_`	description	required

Returns:

Name	Type	Description
`_type_`		description

Source code in pyshbundle/plm.py

def derivALF(inn, miin, plin, m, lmax):
    """HelpeFunction

    Args:
        inn (_type_): _description_
        miin (_type_): _description_
        plin (_type_): _description_
        m (_type_): _description_
        lmax (_type_): _description_

    Returns:
        _type_: _description_
    """
    l = np.arange(m,lmax+2,1)
    #l=l.reshape(l.shape[0],1)
    if m == 0:
        inn[:,0] = 0
        if lmax > m:             
            inn[:,1:] = plin*(-np.sqrt(    ((l[1:]+1)*l[1:]   /2).real))            # (-ones(n,1)*realsqrt((l(2:end)+1).*l(2:end)./2)).*plin            
    elif m == 1:
        inn[:,0] = miin[:,1]
        if lmax > m: 
            inn[:,1:] =  miin[:,2:]*(np.sqrt((l[1:]+1)*l[1:]/2).real) -0.5*plin*(np.sqrt((l[1:]-1)*(l[1:]+2)).real)
    elif m == lmax:
        inn[:,0] = np.sqrt(m/2*miin[:,1:]).real
    else:
        inn[:,0] = np.sqrt((m/2)*miin[:,1:]).real
        if lmax > m: 
            inn[:,1:] = 0.5*miin[:,2:]*np.sqrt((l[:,1:]+m)*(l[:,1:]-m+1)).real - 0.5*plin*(np.sqrt((l[:,1:]-m)*(l[:,1:]+m+1)).real)
    return inn

`lrecur(inn, x, m, lmax)` ¶

[Helper Function]

Parameters:

Name	Type	Description	Default
`inn`	`_type_`	description	required
`x`	`_type_`	description	required
`m`	`_type_`	description	required
`lmax`	`_type_`	description	required

Returns:

Name	Type	Description
`_type_`		description

Source code in pyshbundle/plm.py

def lrecur(inn, x, m, lmax):
    """[Helper Function]  

    Args:
        inn (_type_): _description_
        x (_type_): _description_
        m (_type_): _description_
        lmax (_type_): _description_

    Returns:
        _type_: _description_
    """
    for ll in np.arange(int(m)+1,lmax+1,1):
       col   = ll - m+1			                                                # points to the next collumn of ptmp
       root1 = np.sqrt( (2*ll+1)*(2*ll-1)/((ll-m)*(ll+m)) ).real 
       root2 = np.sqrt( (2*ll+1)*(ll+m-1)*(ll-m-1) / ( (2*ll-3)*(ll-m)*(ll+m) ) ).real

       # % recursion 
       if ll == m+1:
           inn[:, col-1] = root1 *x*inn[:, col-2]
       else:
           inn[:, col-1] = root1 *x*inn[:, col-2] - root2 *inn[:, col-3] 
    return inn

`secrecur(m, y)` ¶

Helper Function:

Parameters:

Name	Type	Description	Default
`m`	`_type_`	description	required
`y`	`_type_`	description	required

Returns:

Name	Type	Description
`_type_`		description

Source code in pyshbundle/plm.py

def secrecur(m, y):
    """Helper Function: 

    Args:
        m (_type_): _description_
        y (_type_): _description_

    Returns:
        _type_: _description_
    """
    if m == 0:
       fac = 1
    else:
       mm  = np.array([2*x for x in range(1, m+1)])
       fac = np.sqrt(2*np.prod((mm+1)/mm))
    out = fac*np.power(y,m)                                                         # The 1st column of ptmp
    return out

`iplm(l, m, theRAD, dt=-9999)` ¶

IPLM Integrals of the fully normalized associated Legendre functions over blocks for a selected order M.

Parameters:

Name	Type	Description	Default
`l`	`_type_`	degree (vector). Integer, but not necessarily monotonic. For l < m a vector of zeros will be returned.	required
`m`	`int`	order (scalar)	required
`theRAD`	`array`	co-latitude [rad] (vector)	required
`dt`	`int`	integration block-size [rad] (scalar). Defaults to -9999.	`-9999`

Returns:

Type	Description
	np.ndarray: Matrix with integrated Legendre functions. Functions are integrated from theRAD(i)-dt/2 till theRAD(i)+dt/2. The matrix has length(TH) rows and length(L) columns, unless L or TH is scalar. Then the output vector follows the shape of respectively L or TH.

Notes

The blocks at the pole might become too large under circumstances. This is not treated separately, i.e. unwanted output may appear. In case TH is scalar, dt will be 1 (arbitrarily).

TO DO

Instead of using sys.exit() we could raise exceptions - that would be a better way of error handling

Uses

plm

Source code in pyshbundle/iplm.py

def iplm(l, m:int, theRAD, dt=-9999):
    """IPLM Integrals of the fully normalized associated Legendre functions
        over blocks for a selected order M. 

    Args:
        l (_type_): degree (vector). Integer, but not necessarily monotonic.
                For l < m a vector of zeros will be returned.
        m (int): order (scalar)
        theRAD (np.array): co-latitude [rad] (vector)
        dt (int, optional): integration block-size [rad] (scalar). Defaults to -9999.

    Returns:
        np.ndarray: Matrix with integrated Legendre functions.
                Functions are integrated from theRAD(i)-dt/2 till theRAD(i)+dt/2.
                The matrix has length(TH) rows and length(L) columns, unless L 
                or TH is scalar. Then the output vector follows the shape of 
                respectively L or TH.

    Notes:
        The blocks at the pole might become too large under circumstances.
        This is not treated separately, i.e. unwanted output may appear.
        In case TH is scalar, dt will be 1 (arbitrarily).

    TO DO:
        Instead of using sys.exit() we could raise exceptions - that would be a better way of error handling

    Uses:
        `plm`
    """
    if dt == -9999:
        if len(theRAD) == 1:
            dt = np.pi/180
        else:
            dt = theRAD[1] - theRAD[0]
    if  min(l.shape) != 1:
        print('Degree l must be a vector (or scalar)') 
        sys.exit([])
    if  np.remainder(l,1).all() != 0:
        print('Vector l contains non-integers !!') 
        sys.exit([])
    # if 
    #     print('Order m must be a scalar') 
    #     sys.exit([])
    if  np.remainder(m,1) != 0:
        print('Order must be integer')
        sys.exit([])
    # if min(dt.shape) !=1:
    #     print('Block size DT must be scalar.') 
        sys.exit([])
    if dt == 0: 
        print('DT cannot be zero') 
        sys.exit([])

    # init
    lcol = len(l)
    trow = len(theRAD)
    n = len(theRAD)
    theRAD.T
    if min(theRAD) < 0 or max(theRAD) > np.pi:
        print('Is the co-latitude ''theta'' given in radian?')
        sys.exit([])
    lmax = max(l[0])
    mfix = m
    lvec = np.transpose(l) 
    l = np.arange(mfix,lmax+1,1)

    # Initialization of cosine, sine and Plm functions
    stplus  = np.sin(theRAD+dt/2)
    stmin   = np.sin(theRAD-dt/2)
    ctplus  = np.cos(theRAD+dt/2)
    ctmin   = np.cos(theRAD-dt/2)
    plmplus = np.ones([n,lmax+1])
    plmmin = np.ones([n,lmax + 1])
    plmplus[:,l] = PLM(np.array([l]),mfix,(theRAD + dt/2),3,1)[:,:,0]                  # Tesserals
    plmmin[:,l] = PLM(np.array([l]),mfix,(theRAD - dt/2),3,1)[:,:,0] 
    if mfix > 0:
        m = np.arange(1,mfix + 1,1)
        mm = 2*m
        fac = np.sqrt(2*np.cumprod((mm+1)/mm))
        mgr, stp = np.meshgrid(m, stplus)
        fgr, stm = np.meshgrid(fac, stmin)
        plmplus[:, m] = fgr * np.power(stp, mgr)
        plmmin[:, m] = fgr * np.power(stm, mgr)
    ptmp = np.zeros([n, lmax +2 ])
    ptmp00 = np.cos(theRAD - dt/2) - ctplus
    ptmp11 = np.sqrt(3)/2 * (dt - ctplus* stplus + ctmin* stmin)
    ptmp10 = np.sqrt(3)/2 * (np.power(stplus,2) - np.power(stmin,2))
    ptmp[:,0] = ptmp00

    # Compute first the integrals of order m == 0
    if mfix == 0:
        ptmp[:,1] = ptmp10
        for l in range(2,lmax+1,1):              #loop over the degree l 
            rootnm = np.sqrt( (2*l+1)*(2*l-1)/np.power(l,2))
            root1nm = np.sqrt( (2*l-1)*(2*l-3)/np.power(l-1,2))
            ptmp[:,l] = rootnm/(l+1)*(((l-2)*ptmp[:,l-2]/root1nm).T + np.power(stplus,2)*plmplus[:,l-1].T - np.power(stmin,2)*plmmin[:,l-1].T )
    else:
        # Compute the integrals of order m > 0

        # First we compute the diagonal element IPmm (lmax == mfix)

        ptmp[:,1] = ptmp11
        for l in range(2,mfix+1,1):
            # print(l)
            rootmm = np.sqrt( (2*l+1)/(2*l) )
            root1mm = np.sqrt( (2*l-1)/(2*l-2))
            if l == 2:
                root1mm = np.sqrt(3)

            ptmp[:,l] = rootmm/(l+1)*( l*root1mm*ptmp[:,l-2].T - (ctplus*plmplus[:,l].T -ctmin*plmmin[:,l].T)/rootmm )
    #the arbitrary element IPlm ( computed only when lmax > mfix)        
        if lmax > mfix:
            l = mfix + 1
        #first we do the element IPlm, for which l - m = 1    
            rootnm = np.sqrt( (2*l+1)*(2*l-1)/(l+mfix)/(l-mfix))
            ptmp[:,l] = rootnm/(l+1)*(np.power(stplus,2)*plmplus[:,l-1].T - np.power(stmin,2)*plmmin[:,l-1].T)
        #now we do the rest
            for l in range(mfix+2,lmax+1,1):              #loop over the degree l
                rootnm = np.sqrt( (2*l+1) * (2*l-1) / (l+mfix) / (l-mfix) )
                root1nm = np.sqrt( (2*l-1) * (2*l-3) / (l-1+mfix) / (l-1-mfix) )
                ptmp[:,l] = rootnm/(l+1)*( (l-2)*ptmp[:,l-2].T/root1nm + np.power(stplus,2)*plmplus[:,l-1].T -np.power(stmin,2)*plmmin[:,l-1].T)


# The integrated functions have been computed. What remains to be done, is to
# extract the proper columns from ptmp, corresponding to the vector lvec. 
# If l or theta is scalar the output matrix p reduces to a vector. It should
# have the shape of respectively theta or l in that case.

# p     = zeros(n, length(lvec))
    lind = np.argwhere(lvec<mfix)[:,0]      #index into l < m
    pcol = lvec + 1                         #index into columns of ptmp
    pcol[lind] = (lmax + 2)*np.ones([len(lind),1])   #Now l < m points to last col
    p = ptmp[:,pcol[:,0]-1]                 #proper column extraction 

    if max(lvec.shape) == 1 and min(np.array([theRAD]).shape) == 1 and trow == 1:
        p = p.T
    if max(np.array([theRAD]).shape) == 1 and min(lvec.shape) == 1 and lcol == 1:
        p = p.T
    return p

GRACE Data Pre-Processing¶

`lovenr(lmax)` ¶

Created on Mon May 11 11:09:28 2022

Todo

Add type and input checking functionality

author: Dr. Bramha Dutt Vishwakarma, Interdisciplinary Center for Water Research (ICWaR), Indian Institute of Science (IISc)

Source code in pyshbundle/GRACEpy.py

def lovenr(lmax: int):
    """
    Created on Mon May 11 11:09:28 2022

    Todo:
        + Add type and input checking functionality

    _author_: Dr. Bramha Dutt Vishwakarma, Interdisciplinary Center for Water Research (ICWaR), Indian Institute of Science (IISc)
    """
    l  = [0,  1,    2,    3,    4,    5,   6,   7,   8,   9,  10,  12,  15,  20,  30,  40,  50,  70, 100, 150, 200]
    kl = numpy.divide([0,270,-3030,-1940,-1320,-1040,-890,-810,-760,-720,-690,-640,-580,-510,-400,-330,-270,-200,-140,-100, -700],1e4)
    n = range(0, lmax+1, 1)
    kn = numpy.interp(n,l,kl)
    return(kn)

`lovenrPREM(lmax, frame)` ¶

Created on Mon May 11 11:51:29 2022

@author: Dr. Bramha Dutt Vishwakarma, Interdisciplinary Center for Water Research (ICWaR), Indian Institute of Science (IISc)

Source code in pyshbundle/GRACEpy.py

def lovenrPREM(lmax:int, frame):
    """
    Created on Mon May 11 11:51:29 2022

    @author:  Dr. Bramha Dutt Vishwakarma, Interdisciplinary Center for Water Research (ICWaR), Indian Institute of Science (IISc)
    """

    data = numpy.array([[ 1,  -0.28476,   0.00000,   0.10462],
    [2,  -0.99297,  -0.61274,   0.04661] ,
    [3,  -1.05142,  -0.58897 ,  0.21048] ,
    [4,  -1.05378,  -0.53513 ,  0.23564] ,
    [5,  -1.08658,  -0.52382 ,  0.23186] ,
    [6,  -1.14404,  -0.54222 ,  0.23263] ,
    [7,  -1.21254,  -0.57464 ,  0.24058] ,
    [8,  -1.28403,  -0.61256 ,  0.25308] ,
    [9,  -1.35479,  -0.65203 ,  0.26799] ,
    [10,  -1.42330,  -0.69140,   0.28419] ,
    [11,  -1.48909,  -0.72998,   0.30121] ,
    [12,  -1.55204,  -0.76749,   0.31880] ,
    [13,  -1.61221,  -0.80381,   0.33684] ,
    [14,  -1.66968,  -0.83886,   0.35522] ,
    [15,  -1.72454,  -0.87260,   0.37382] ,
    [16,  -1.77684,  -0.90499,   0.39251] ,
    [17,  -1.82668,  -0.93599,   0.41119] ,
    [18,  -1.87414,  -0.96560,   0.42973] ,    
    [19,  -1.91928,  -0.99382,   0.44804] ,
    [20,  -1.96220,  -1.02066,   0.46603] ,
    [21,  -2.00297,  -1.04614,   0.48363] ,
    [22,  -2.04169,  -1.07029,   0.50078] ,
    [23,  -2.07844,  -1.09313,   0.51742] ,
    [24,  -2.11332,  -1.11472,   0.53355] ,
    [25,  -2.14642,  -1.13511,   0.54912] ,
    [30,  -2.28839,  -1.22067,   0.61848] ,
    [40,  -2.48641,  -1.33024,   0.71925] ,
    [50,  -2.61710,  -1.39016,   0.78410] ,
    [60,  -2.71254,  -1.42377,   0.82683] ,
    [70,  -2.78865,  -1.44313,   0.85550] ,
    [80,  -2.85368,  -1.45474,   0.87479] ,
    [90,  -2.91216,  -1.46226,   0.88764] ,
    [100,  -2.96672,  -1.46787,   0.89598] ,
    [120,  -3.06983,  -1.47811,   0.90421] ,
    [140,  -3.16950,  -1.49082,   0.90634] ,
    [160,  -3.26809,  -1.50771,   0.90603] ,
    [180,  -3.36633,  -1.52909,   0.90532] ,
    [200,  -3.48436,  -1.55473,   0.90547] ,
    [250,  -3.70773,  -1.63448,   0.91388] ,
    [300,  -3.94607,  -1.73053,   0.93714] ,
    [350,  -4.17591,  -1.83593,   0.97495] ,
    [400,  -4.39433,  -1.94515,   1.02467] ,
    [500,  -4.78872,  -2.15940,   1.14615] ,
    [600,  -5.12008,  -2.35243,   1.27714] ,
    [800,  -5.59959,  -2.64798,   1.50995] ,
    [1000,  -5.88447,  -2.83157,  1.67325] ,
    [1500,  -6.15106,  -3.00957,   1.84797] ,
    [2000,  -6.20058,  -3.04408,   1.88423] ,
    [3000,  -6.21044,  -3.05176,   1.89114] ,
    [5000,  -6.21155,  -3.05324,   1.89118] ,
    [10000,  -6.21226,  -3.05427,   1.89110]])

    l  =  data[:,0]
    hl =  data[:,1]
    kl =  numpy.divide(data[:,2], l)
    ll =  numpy.divide(data[:,3], l)

    if frame == 'CM':
        hl[0] = hl[0] - 1
        ll[0] = ll[0] - 1
        kl[0] = kl[0] - 1
        print('Love numbers are in center of mass frame')
    elif frame == 'CF':
        hlo = hl[0]
        llo = ll[0]
        hl[0] = (hlo - llo) * 2/3
        ll[0] = (hlo - llo) * (-1/3)
        kl[0] = ((-2/3)*llo) - ((-1/3)*hlo)
        print('Love numbers are in center of figure frame')
    elif frame == 'CE':
        print('Love numbers are in center of solid Earth frame')
    else:
        lovenrPREM.exit('Please choose a compatible frame of reference: one of CM, CF, or CE')


    n = range(0, lmax+1, 1)
    kn = numpy.interp(n,l,kl)
    hn = numpy.interp(n,l,hl)
    ln = numpy.interp(n,l,ll)
    kn[0] = 0 
    hn[0] = 0
    ln[0] = 0
    return(kn,hn,ln)

`upwcon(degree, height)` ¶

Returns the upward continuation \((R/r)^l\)

Parameters:

Name	Type	Description	Default
`degree`	`int`	Spherical harmonic degree	required
`height`	`int`	Height above mean Earth radius [m] [scalar/vector]	required

Returns:

Name	Type	Description
`uc`	`_type_`	Upward continuation terms

Uses

GRACEconstants.GC

Todo

Add input checking functionality and raise exceptions
Add reference to formula

Source code in pyshbundle/GRACEpy.py

def upwcon(degree: int, height):
    """Returns the upward continuation $(R/r)^l$

    Args:
        degree (int): Spherical harmonic degree
        height (int): Height above mean Earth radius [m] [scalar/vector]

    Returns:
        uc (_type_): Upward continuation terms

    Uses:
        `GRACEconstants.GC`

    Todo:
        + Add input checking functionality and raise exceptions
        + Add reference to formula
    """
    # Created on Sat May  9 18:49:45 2022
    rr = numpy.divide(GC.ae, numpy.add(GC.ae,height))
    uc = numpy.power(rr, degree)

    return(uc)    

Filtering the GRACE Data¶

`Gaussian(L, cap)` ¶

The program delivers the spherical harmonic coefficients of a gaussian smoothing filter. The coefficients are calculated according to Wahr et. al.(1998) equation (34) and Swenson and Wahr equation (34)

Parameters:

Name	Type	Description	Default
`L`	`int`	maximum degree	required
`cap`	`int`	half width of Gaussian smoothing function [km]	required

Returns:

Type	Description
	np.ndarray: smoothing coefficients

Raises:

Type	Description
`TypeError`	Degree must be integer
`ValueError`	Maximum degree must be higher than 2
`TypeError`	Cap size must be an integer

References

Wahr et.al. (1998) equation (34) and Swenson and Wahr equation (34)

Source code in pyshbundle/gaussian.py

def Gaussian(L: int, cap: int):
    """The program delivers the spherical harmonic coefficients of a gaussian
    smoothing filter. The coefficients are calculated according to Wahr et. al.(1998)
    equation (34) and Swenson and Wahr equation (34)


    Args:
        L (int): maximum degree
        cap (int): half width of Gaussian smoothing function [km]

    Returns:
        np.ndarray: smoothing coefficients

    Raises:
        TypeError: Degree must be integer
        ValueError: Maximum degree must be higher than 2
        TypeError: Cap size must be an integer

    References:
        Wahr et.al. (1998) equation (34) and Swenson and Wahr equation (34)

    """

    #Check input
    if type(L) != int:
        raise TypeError('Degree must be integer')

    if L<2:
        raise ValueError('Maximum degree must be higher than 2')

    if type(cap) != int:
        raise TypeError('Cap size must be an integer')

    #Calculations
    W = np.zeros([L+1, 1])
    b = np.log(2)/(1 - np.cos(cap/6371))

    #Recursive calculation of the weighting coefficients
    W[0,0] = 1
    W[1,0] = np.power(10, np.log10( (1 + np.exp(-2*b))/(1-np.exp(-2*b)) - (1/b)))

    i = 1
    while i < L:        
        j = i + 1
        W[i+1][0] = W[i-1][0] - (2*(j-1) + 1)/b * W[i][0]
        if W[i+1, 0] > W[i] or W[i+1] < 0:
            W[i+1] = 0
        i = i + 1

    return W

Numerical Integration¶

`grule(n)` ¶

This function computes Gauss base points and weight factors using the algorithm-see Reference

Parameters:

Name	Type	Description	Default
`n`	`int`	number of base points required	required

Returns:

Type	Description
	np.array: cosine of the base points
	np.array: weight factors for computing integrals and such

References

'Methods of Numerical Integration' by Davis and Rabinowitz, page 365, Academic Press, 1975.

Source code in pyshbundle/grule.py

def grule(n: int):
    """This function computes Gauss base points and weight factors
    using the algorithm-see Reference

    Args:
        n (int): number of base points required

    Returns:
        np.array: cosine of the base points
        np.array: weight factors for computing integrals and such

    References:
        1. 'Methods of Numerical Integration' by Davis and Rabinowitz, page 365, Academic Press, 1975.

    """
    bp = np.zeros((n,1))
    wf = bp
    iter = 2
    m = np.floor((n+1)/2)
    e1 = n * (n+1)


    mm = 4*m - 1
    t = (np.pi / (4*n + 2)) * np.arange(3,mm+4,4)
    nn = (1 - (1 - 1/n)/(8*n*n))
    x0 = nn * np.cos(t)


    for i in range(iter):
        pkm1 = 1
        pk = x0

        for kk in range(n-1):
            k = kk + 2
            t1 = x0 * pk
            pkp1 = t1 - pkm1 - (t1-pkm1)/k  + t1
            pkm1=pk
            pk=pkp1

        den = 1 - x0*x0
        d1 = n * (pkm1 - x0*pk)
        dpn = d1/den


        d2pn = (2*x0*dpn - e1*pk) / den
        d3pn = (4*x0*d2pn + (2-e1)*dpn)/den
        d4pn = (6*x0*d3pn + (6-e1)*d2pn)/den
        u = pk/dpn
        v = d2pn/dpn
        h = -u * (1+(.5*u)*(v+u*(v*v - u*d3pn/(3*dpn))))
        p = pk + h*(dpn+(0.5*h)*(d2pn+(h/3)*(d3pn + 0.25*h*d4pn)))
        dp = dpn + h*(d2pn+(0.5*h)*(d3pn+h*d4pn/3))
        h = h-p/dp
        x0 = x0+h

    bp = -x0-h
    fx = d1 - h*e1*(pk+(h/2)*(dpn+(h/3)*(d2pn+(h/4)*(d3pn+(0.2*h)*d4pn))))
    wf = [2 * (1 - np.power(bp,2))]/(fx*fx)


    for i in range(len(bp),n):
        bp = np.append(bp,[0])
        wf = np.append(wf,[0])

    if ((m)+(m)) != (n):
        m = m-1

    for i in range(1,int(m+1)):
        bp[-i] = -bp[i-1]
        wf[-i] = wf[i-1] 
    return bp, wf

`naninterp(X)` ¶

This function uses cubic interpolation to replace NaNs

Parameters:

Name	Type	Description	Default
`X`	`_type_`	description	required

Returns:

Name	Type	Description
`_type_`		description

Source code in pyshbundle/naninterp.py

def naninterp(X):
    """This function uses cubic interpolation to replace NaNs

    Args:
        X (_type_): _description_

    Returns:
        _type_: _description_
    """

    nan = np.nan

    ok = ~np.isnan(X)
    xp = ok.ravel().nonzero()[0] #Indices of xs with values
    fp = X[~np.isnan(X)]

    x  = np.isnan(X).ravel().nonzero()[0] #Indices of xs without values

    pchip = PchipInterpolator(xp,fp) #Initialize scipy PHCIP cubic interpolation
    X[np.isnan(X)] = pchip(x) #Interpolate Nan values in X

    return X

`neumann(inn)` ¶

Returns the weights and nodes for Neumann's numerical integration

Parameters:

Name	Type	Description	Default
`inn`	`(int, array)`	base points (nodes) in the interval [-1;1]	required

Raises:

Type	Description
`TypeError`	Integer input argument required
`ValueError`	Error in input dimensions

Returns:

Name	Type	Description
`_type_`		quadrature weights
`_type_`		base points (nodes) in the interval [-1;1]

Remarks

1st N.-method: see Sneeuw (1994) GJI 118, pp 707-716, eq. 19.5
2nd N.-method: see uberall/GRULE

Todo

TypeError is more relavant and shape error from numpy

Uses

grule, plm

Examples:

>>> TO DO: write example how to use the function

Source code in pyshbundle/neumann.py

def neumann(inn):
    """Returns the weights and nodes for Neumann's numerical integration

    Args:
        inn (int, np.array): base points (nodes) in the interval [-1;1]

    Raises:
        TypeError: Integer input argument required
        ValueError: Error in input dimensions

    Returns:
        _type_: quadrature weights
        _type_: base points (nodes) in the interval [-1;1]

    Remarks:
        * 1st N.-method: see Sneeuw (1994) GJI 118, pp 707-716, eq. 19.5
        * 2nd N.-method: see uberall/GRULE

    Todo: 
        + TypeError is more relavant and shape error from numpy

    Uses:
        `grule`, `plm`

    Examples:
        >>> TO DO: write example how to use the function
    """

    try: #if input is an integer
        x, w = grule(inn)
    except: #if input is an array
        if(len(inn)==1): #2nd Neumann method
            x, w = grule(inn)
            if(np.not_equal(np.mod(x, 1), 0)): #Not integer
                raise TypeError("Integer input argument required")



        elif min(inn.shape) == 1: #1st Neumann method #Size gives 2 outputs for 2d array in matlab; for row and column
            x = inn
            theRAD = np.arccos(x) #x in radian
            l = np.array(list(range(len(x))))
            pp = PLM(l, theRAD)

            rr = list([2])
            for i in len(x-1):
                rr.append(0)
            r = np.asarray(rr)

            w,resid,rank,s = np.linalg.lstsq(pp,r) #Solve system of equations; Double check this operation
            if(x.shape != w.shape):
                w = w.T

        else:
            raise ValueError("Error in input dimensions")
            # TO DO: Write more descriptive exception messages

    return w, x

Important for Spherical Harmonic Synthesis¶

`ispec(a, b=-9999)` ¶

Returns the function F from the spectra A and B

Parameters:

Name	Type	Description	Default
`a`	`_type_`	cosine coefficients	required
`b`	`int`	sine coefficients. Defaults to -9999.	`-9999`

Returns:

Name	Type	Description
`f`	`_type_`	description

`normalklm(lmax, typ='wgs84')` ¶

NORMALKLM returns an ellipsoidal normal field consisting of normalized -Jn, n=0,2,4,6,8

Parameters:

Name	Type	Description	Default
`lmax`	`int`	maximum degree	required
`typ`	`str`	Ellipsoids can be either 'wgs84' - World Geodetic System 84, 'grs80' - , 'he' - hydrostatic equilibrium ellipsoid	`'wgs84'`

Returns:

Name	Type	Description
`nklm`	`array`	normal field in CS-format (sparse array - [1, -J2, -J4, -J6, -J8])

TODO

Find type of nklm; I think raising TypeError, VlueError or NameError instad of general Exception

Raises:

Type	Description
`TypeError`	lmax should be an integer
`ValueError`	lmax should be positive
`ValueError`	Unknown type of ellipsoid, supports 'wgs84', `GRS80` and 'he'

References

J2,J4 values for hydrostatic equilibrium ellipsoid from Lambeck (1988) "Geophysical Geodesy", p.18

Source code in pyshbundle/normalklm.py

def normalklm(lmax: int, typ: str = 'wgs84'):
    """ NORMALKLM returns an ellipsoidal normal field
    consisting of normalized -Jn, n=0,2,4,6,8

    Args:
        lmax (int): maximum degree
        typ (str): Ellipsoids can be either 
                    'wgs84' - World Geodetic System 84, 
                    'grs80' - , 
                    'he' - hydrostatic equilibrium ellipsoid

    Returns:
        nklm (np.array): normal field in CS-format (sparse array - [1, -J2, -J4, -J6, -J8])

    TODO: 
        Find type of nklm; I think raising TypeError, VlueError or NameError instad of general Exception

    Raises:
        TypeError: lmax should be an integer
        ValueError: lmax should be positive
        ValueError: Unknown type of ellipsoid, supports 'wgs84', `GRS80` and 'he'

    References:
        1. J2,J4 values for hydrostatic equilibrium ellipsoid from Lambeck (1988)
        "Geophysical Geodesy", p.18    
    """

    if type(lmax) != int:
        raise TypeError("lmax should be integer")

    if lmax < 0:
        raise ValueError("lmax should be positive")


    typ_ = typ.lower()
    if (typ_ == 'wgs84'):
        J2     =  1.08262982131e-3     #% earth's dyn. form factor (= -C20 unnormalized)
        J4     = -2.37091120053e-6    #% -C40 unnormalized
        J6     =  6.08346498882e-9     #% -C60 unnormalized
        J8     = -1.42681087920e-11    #% -C80 unnormalized
        jcoefs = np.array([1, -J2, -J4, -J6, -J8]).T.reshape(5,1)
        # as lmax + 2 is requires 
        l      = np.arange(0,min(lmax + 2,8 + 2), 2).T
        l.reshape(l.shape[0],1)

    elif (typ_ == 'grs80'):
        J2     =  1.08263e-3         # % earth's dyn. form factor (= -C20 unnormalized)
        J4     = -2.37091222e-6     #% -C40 unnormalized
        J6     =  6.08347e-9        #% -C60 unnormalized
        J8     = -1.427e-11         #% -C80 unnormalized
        jcoefs = np.array([1, -J2, -J4, -J6, -J8]).reshape(5,1)
        l      = np.arange(0,min(lmax + 2,8 + 2), 2).T
        l.reshape(l.shape[0],1)

    elif ((typ_ == 'he') or (typ_ == 'hydro')):
        J2     = 1.072618e-3		#% earth's dyn. form factor (= -C20 unnormalized)
        J4     = 0.2992e-5     	#% -C40 unnormalized
        jcoefs = np.array([1, -J2, -J4]).T.reshape(5,1)
        # adding (2) beacuse of arange function is only uptp last integer and not including last
        l      = np.arange(0,min(lmax + 2,4 + 2), 2).T
        l.reshape(l.shape[0],1)

    else:
        raise ValueError("Unknown type of ellipsoid:   ", typ)

    coefs = jcoefs[:len(l)].T / np.sqrt(2*l + 1)
#    coefs.reshape(coefs.shape[0],1)


    data = np.array(coefs)[0]
    row = np.array(l)
    col = np.zeros(len(l))
    # lmax = 96 then shape=(97, 97) -> consisitent with everything else
    nklm = sparse.coo_matrix((data,(row,col)),shape=(lmax+1,lmax+1)).toarray()
    return nklm

`eigengrav(lmax, fstr, h)` ¶

summary

Parameters:

Name	Type	Description	Default
`lmax`	`int`	Maximum degree of Spherical Coefficients	required
`fstr`	`str`	gravity quantity, options: 'None', 'geoid', 'potential', 'gravity', 'tr', 'trr', 'slope' 'water', 'smd', 'height'	required
`h`	`float`	description	required

Returns:

Name	Type	Description
`_type_`		description

Raises:

Type	Description
`TypeError`	Enter a valid lmax value

TO DO

Can we think about the raising a ValueError instead of instantly terminating the function Adding comments as variable names are not much descriptive

Source code in pyshbundle/eigengrav.py

def eigengrav(lmax: int, fstr: str, h: float):
    """_summary_

    Args:
        lmax (int): Maximum degree of Spherical Coefficients
        fstr (str): gravity quantity, options: 'None', 'geoid', 'potential', 'gravity', 'tr', 'trr', 'slope'
                    'water', 'smd', 'height'
        h (float): _description_

    Returns:
        _type_: _description_

    Raises:
        TypeError: Enter a valid lmax value

    TO DO:
        Can we think about the raising a ValueError instead of instantly terminating the function
        Adding comments as variable names are not much descriptive
    """

    if type(lmax) == int:
        rows = 1
    else:
        rows = len(lmax)
    # rows = len(l)

    if rows > 1 or lmax < 0:
        raise TypeError("Enter a valid lmax value")

    r = GC.ae + h

    # lmax issue - using lmax as per used in shbundle
    # no reference for height

    if fstr == 'none':
        tf = numpy.ones((1, lmax+1))
    elif fstr == 'geoid':
        tf = numpy.ones((1, lmax+1)) * r
    elif fstr == 'potential':
        tf = numpy.ones((1, lmax+1)) * (GC.GM/r)
    elif fstr == 'gravity' or fstr == 'dg':
        tf = numpy.multiply(range(-1, lmax, 1), ((GC.GM/r/r) * 1e5))
        tf = tf.reshape((1, lmax+1))
    elif fstr == 'tr':
        tf = numpy.multiply(range(-1, -(lmax+2), -1), ((GC.GM/r/r) * 1e5))
        tf = tf.reshape((1, lmax+1))
    elif fstr == 'trr':
        tf = numpy.multiply(range(1, (lmax+2), 1),
                            range(2, (lmax + 3), 1))*((GC.GM/r/r) * 1e9)
        tf = tf.reshape((1, lmax+1))
    elif fstr == 'slope':
        tf = numpy.sqrt(numpy.multiply(
            range(0, lmax+1, 1), range(1, lmax+2, 1)))
        tf = tf.reshape((1, lmax+1))
    elif fstr == 'water':
        ln = GB.lovenr(lmax)
        tf = numpy.divide(numpy.multiply(
            5.517*r, numpy.add(range(0, 2*lmax + 1, 2), 1)), numpy.multiply(3, (1+ln)))
        tf = tf.reshape((1, lmax+1))
    elif fstr == 'smd':
        ln = GB.lovenr(lmax)
        tf = numpy.divide(numpy.multiply(
            5517*r, numpy.add(range(0, 2*lmax + 1, 2), 1)), numpy.multiply(3, (1+ln)))
        tf = tf.reshape((1, lmax+1))
    elif fstr == 'height':
        # not sure aobut heights - kept it unchanged ... abhishek
        kl, hl, ll = GB.lovenrPREM(90, 'CF')
        tf = numpy.divide(numpy.multiply(hl, (GC.ae*1000)), numpy.add(kl, 1))
    else:
        ValueError('Please choose a valid quantity for fstr')

    if h > 0:
        upConTerm = GB.upwcon(lmax, h)
        tf = numpy.multiply(tf, upConTerm)

    return(tf)

Time Series¶

`PhaseCalc(fts, ffts)` ¶

calculates the phase difference between two time series based on the Hilbert transform method explained by Phillip et al.

Parameters:

Name	Type	Description	Default
`fts`	`ndarray`	description	required
`ffts`	`ndarray`	description	required

Returns:

Name	Type	Description
`_type_`		description

References

Phillips, T., R. S. Nerem, B. Fox-Kemper, J. S. Famiglietti, and B. Rajagopalan (2012), The influence of ENSO on global terrestrial water storage using GRACE, Geophysical Research Letters, 39 (16), L16,705, doi:10.1029/2012GL052495.

Source code in pyshbundle/Phase_calc.py

def PhaseCalc(fts, ffts):
    """calculates the phase difference between two time series based on the
    Hilbert transform method explained by Phillip et al.

    Args:
        fts (np.ndarray): _description_
        ffts (np.ndarray): _description_

    Returns:
        _type_: _description_

    References:
        1. Phillips, T., R. S. Nerem, B. Fox-Kemper, J. S. Famiglietti, and B. Rajagopalan (2012),
        The influence of ENSO on global terrestrial water storage using GRACE, Geophysical
        Research Letters, 39 (16), L16,705, doi:10.1029/2012GL052495.
    """
    c = fts.shape[1]

    ps = np.zeros((1, c))

    filter_ = ~np.isnan(fts)
    filter__ = ~np.isnan(ffts)

    fts_ = fts[filter_] #Extract values and leave Nan
    ffts_ = ffts[filter__] #Extract values and leave Nan

    fts = fts_.reshape(int(fts_.shape[0]/c),c)
    ffts = ffts_.reshape(int(ffts_.shape[0]/c),c)

    rn = fts.shape[0]

    for i in range(c):
        # A = np.concatenate(np.ones((rn,1)), np.real(signal.hilbert(ffts[:, i])), np.imag(signal.hilbert(ffts[:, i]))) #design matrix

        A = np.array((np.ones((rn)), np.real(signal.hilbert(ffts[:, i])), np.imag(signal.hilbert(ffts[:, i])))).T

        A = A.astype('double')
        B = fts[:,i]
        B = B.astype('double')
        abc = np.linalg.lstsq(A.T @ A, A.T @ B)[0]

        ps[0,i] = np.arctan2(abc[3-1],abc[2-1])*(180/np.pi) #check indices and degree/radian
    return ps

Auxillary Codes¶

Computing Legendre functions¶

PLM(l, m, thetaRAD, nargin, nargout) ¶

derivALF(inn, miin, plin, m, lmax) ¶

lrecur(inn, x, m, lmax) ¶

secrecur(m, y) ¶

iplm(l, m, theRAD, dt=-9999) ¶

GRACE Data Pre-Processing¶

lovenr(lmax) ¶

lovenrPREM(lmax, frame) ¶

upwcon(degree, height) ¶

Filtering the GRACE Data¶

Gaussian(L, cap) ¶

Numerical Integration¶

grule(n) ¶

naninterp(X) ¶

neumann(inn) ¶

Important for Spherical Harmonic Synthesis¶

ispec(a, b=-9999) ¶

normalklm(lmax, typ='wgs84') ¶

eigengrav(lmax, fstr, h) ¶

Time Series¶

PhaseCalc(fts, ffts) ¶

`PLM(l, m, thetaRAD, nargin, nargout)` ¶

`derivALF(inn, miin, plin, m, lmax)` ¶

`lrecur(inn, x, m, lmax)` ¶

`secrecur(m, y)` ¶

`iplm(l, m, theRAD, dt=-9999)` ¶

`lovenr(lmax)` ¶

`lovenrPREM(lmax, frame)` ¶

`upwcon(degree, height)` ¶

`Gaussian(L, cap)` ¶

`grule(n)` ¶

`naninterp(X)` ¶

`neumann(inn)` ¶

`ispec(a, b=-9999)` ¶

`normalklm(lmax, typ='wgs84')` ¶

`eigengrav(lmax, fstr, h)` ¶

`PhaseCalc(fts, ffts)` ¶