Core functionality

This module has both Spherical Harmonics Analysis and Spherical Harmonics Synthesis

Spherical Harmonic Analysis and Synthesis

gsha(f, method, grid=None, lmax=-9999)

GSHA - Global Spherical Harmonic Analysis

Parameters:

Name	Type	Description	Default
f	ndarray	global field of size \((l_{max} + 1) * 2 * l_{max}\) or \(l_{max} * 2 * l_{max}\)	required
method	str	method to be used	required
grid	str	choose between 'block' or 'cell'. Defaults to None.	None
lmax	int	maximum degree of development. Defaults to -9999.	-9999

Returns:

Type	Description
	np.ndarray: Clm, Slm in |C\S| format

Raises:

Type	Description
ValueError	grid argument can only be 'block' or 'cell'
ValueError	Grid type entered is not right
TypeError	Invalid size of matrix F
TypeError	GRID and METHOD must be strings
ValueError	2nd Neumann method ONLY on a ''neumann''/''gauss'' GRID'
ValueError	Block mean method ONLY on a ''block''/''cell'' GRID
ValueError	Maximum degree of development is higher than number of rows of input.

Uses

plm, neumann, iplm, sc2cs

Source code in pyshbundle/gsha.py

def gsha(f, method: str, grid: str = None, lmax: int = -9999):
    """ GSHA - Global Spherical Harmonic Analysis

    Args:
        f (np.ndarray): global field of size $(l_{max} + 1) * 2 * l_{max}$ or $l_{max} * 2 * l_{max}$
        method (str): method to be used
        grid (str, optional): choose between 'block' or 'cell'. Defaults to None.
        lmax (int, optional): maximum degree of development. Defaults to -9999.

    Returns:
        np.ndarray: Clm, Slm in |C\S| format

    Raises:
        ValueError: grid argument can only be 'block' or 'cell'
        ValueError: Grid type entered is not right
        TypeError: Invalid size of matrix F
        TypeError: GRID and METHOD must be strings
        ValueError: 2nd Neumann method ONLY on a ''neumann''/''gauss'' GRID'
        ValueError: Block mean method ONLY on a ''block''/''cell'' GRID
        ValueError: Maximum degree of development is higher than number of rows of input.

    Uses:
        `plm`, `neumann`, `iplm`, `sc2cs`
    """
    rows, cols = f.shape

    if cols == 2 * rows: #Check conditions
        if lmax == -9999:
            lmax = rows

        if grid == None:
            grid = 'block'

        if (grid != 'block') and (grid != 'cell'):
            raise ValueError("Your GRID variable should be either block or cell")

        n = rows
        dt = 180 / n
        theta = np.arange(dt/2, 180+(dt/4), dt)
        lam = np.arange(dt/2, 360+(dt/4), dt)

    elif cols == 2 * rows - 2:
        if lmax == -9999:
            lmax = rows - 1
        if grid == None:
            grid = 'pole'

        n = rows - 1
        dt = 180 / n

        if (grid == 'pole') or (grid == 'mesh'):                   
            theta = np.arange(0, 180+(dt/4), dt)
            lam = np.arange(0, 360+(dt/4) - dt, dt)
        elif (grid == 'neumann') or (grid == 'gauss'): 
        # gw, gx = neumann(n+1) #For some reason, grule does not work for even values
            gw, gx = neumann(n)
            theta = np.arccos(np.flipud(gx)) * 180 / np.pi
            lam = np.arange(0, 360+(dt/4)-dt, dt)

            if len(gw.shape) == 1:
                gw = gw.reshape(gw.shape[0],1)

            if len(gx.shape) == 1:
                gx = gx.reshape(gx.shape[0],1)
        else:
            raise ValueError("Grid type entered is not right")
    else:
        raise TypeError("Invalid size of matrix F")

    theRAD = theta * np.pi / 180
    # if len(theRAD.shape) == 1:
    # theRAD = theRAD.reshape(theRAD.shape[0],1)


    # further diagnostics

    if (type(grid) != str) or (type(method) != str):
        raise TypeError("GRID and METHOD must be strings.")

    if (method == 'snm') and ((grid != 'neumann') and (grid != 'gauss')):
        raise ValueError('2nd Neumann method ONLY on a ''neumann''/''gauss'' GRID')

    if (method == 'mean') and ((grid != 'block') and (grid != 'cell')):
        raise ValueError('Block mean method ONLY on a ''block''/''cell'' GRID')

    if lmax > n:
        raise ValueError('Maximum degree of development is higher than number of rows of input.')

    # Reshape variables as required

    if len(lam.shape) == 1:
        lam = lam.reshape(1,lam.shape[0])

    # Init

    L = n
    clm = np.zeros((L+1, L+1), dtype='longdouble')
    slm = np.zeros((L+1, L+1), dtype='longdouble')


    # First step of analysis

    m = np.arange(L+1).reshape(1,L+1)
    c = np.cos((lam.T @ m) * np.pi/180)
    s = np.sin((lam.T @ m) * np.pi/180)


    # % preserving the orthogonality (except for 'mean' case)
    # % we distinguish between 'block' and 'pole' type grids (in lambda)

    if (grid == 'block') or (grid == 'cell'):
        if method == 'mean':
            dl = dt
            c[:,0] = dl / 360
            m = np.arange(1, L+1)
            ms = 2 / m * np.sin(m * dl/2 * np.pi/180) / np.pi
            c[:,1:(L+1)+1] = c[:,1:(L+1)+1] * ms  
            s[:,1:(L+1)+1] = s[:,1:(L+1)+1] * ms

        else:
            c = c/L
            s = s/L
            c[:,0] = c[:,1]/2
            s[:,L] = s[:,L]/2
            c[:,L] = np.zeros(2*n)
            s[:,0] = np.zeros(2*n)
    else:
        c = c/L
        s = s/L
        c[:,[0, L]] = c[:,[0, L]]/2	
        s[:,[0, L]] = np.zeros((2*n,2))	  


    a = f @ c
    b = f @ s    

    # Second step of analysis: Clm and Slm

    if method == 'ls':
        for m in range(L+1):
#            l = np.arange(m,L+1)
            l = np.arange(m,L+1).reshape(L+1-m, 1)
            l = l.T

            p = PLM(l,m,theRAD, 3, 1)
            p = p[:,:,0]
            ai = a[:, m]
            bi = b[:, m]

            clm[m+1:L+2, m+1] = linalg.lstsq(p, ai)
            slm[m+1:L+2, m+1] = linalg.lstsq(p, bi)


    elif method == 'aq': #Approximate Quadrature
        si = np.sin(theRAD)
        si = 2 * si / np.sum(si)

        for m in range(L+1):
            l = np.arange(m, L+1).reshape(L+1-m, 1)
            l = l.T

            p = PLM(l,m,theRAD, 3, 1)

            ai = a[:, m]
            bi = b[:, m]

            clm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ (si * ai)
            slm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ (si * bi)

    elif method == 'fnm': #1st Neumann method (exact upto L/2)
        w = neumann(np.cos(theRAD))

        for m in range(L+1):
            l = np.arange(m, L+1).reshape(L+1-m, 1)
            l = l.T

            p = PLM(l,m,theRAD, 3, 1)

            ai = a[:, m]
            bi = b[:, m]

            clm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ (w * ai)
            slm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ (w * bi)

    elif method == 'snm': #2nd Neumann method (exact)
        for m in range(L+1):
            l = np.arange(m, L+1).reshape(L+1-m, 1)
            l = l.T

            p = PLM(l,m,theRAD, 3, 1)

            ai = a[:, m]
            bi = b[:, m]

            clm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ (gw * ai)
            slm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ (gw * bi)

    elif method == 'mean':
        for m in range(L+1):
            print(m)
            #l = np.arange(m,L+1).reshape(L+1-m,1)
            #l = l.T


            l = np.array([np.arange(m,L+1, 1)])
        # l = np.array([[m]])

            p = iplm(l,m,theRAD)
        # p = p[:,-1]
            ai = a[:, m]
            bi = b[:, m]

            clm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ ai
            slm[m:L+1, m] = (1 + (m == 0))/ 4 * p.T @ bi

    # Write the coefficients Clm & Slm in |C\S| format

    slm = np.fliplr(slm)
    cs = sc2cs(np.concatenate((slm[:, np.arange(L)], clm), axis = 1))
    cs = cs[:int(lmax+1), :int(lmax+1)]


    # np.save('/path/csRb.npy',cs)

    return cs

GSHS(field, quant='none', grd='mesh', n=-9999, h=0, jflag=1)

GSHS - Global Spherical Harmonic Synthesis

Parameters:

Name	Type	Description	Default
field	_type_	set of SH coefficients, either in SC-triangle or CS-square format	required
quant	str	defining the field quantity. Defaults to 'none'.	'none'
grd	str	defining the grid. Defaults to 'mesh'.	'mesh'
n	int	description. Defaults to -9999.	-9999
h	int	description. Defaults to 0.	0
jflag	int	description. Defaults to 1.	1

Returns:

Name	Type	Description
f	ndarray	the global field
		theRAD (): vector of co-latitudes [rad]
		lamRAD (): vector of longitudes [rad]

Raises:

Type	Description
Exception	Check format of the field
Exception	n must be scalar
Exception	n must be integer
Exception	Grid argument must be string
Exception	description

Uses

cs2sc, normalklm, plm, eigengrav, ispec

Todo

• Change general exceptions to specific and descriptive built-in ones
• using the not and then check is not always advisable
• Check how to document valid options

Source code in pyshbundle/gshs.py

def GSHS(field, quant = 'none', grd = 'mesh', n = -9999, h = 0, jflag = 1):
    """GSHS - Global Spherical Harmonic Synthesis

    Args:
        field (_type_): set of SH coefficients, either in SC-triangle or CS-square format
        quant (str, optional): defining the field quantity. Defaults to 'none'.
        grd (str, optional): defining the grid. Defaults to 'mesh'.
        n (int, optional): _description_. Defaults to -9999.
        h (int, optional): _description_. Defaults to 0.
        jflag (int, optional): _description_. Defaults to 1.

    Returns:
        f (np.ndarray): the global field
        theRAD (): vector of co-latitudes [rad]
        lamRAD (): vector of longitudes [rad]

    Raises:
        Exception: Check format of the field
        Exception: n must be scalar
        Exception: n must be integer
        Exception: Grid argument must be string
        Exception: _description_

    Uses:
        `cs2sc`, `normalklm`, `plm`, `eigengrav`, `ispec`

    Todo: 
        * Change general exceptions to specific and descriptive built-in ones
        + using the not and then check is not always advisable
        + Check how to document valid options
    """

    wd = getcwd()
    chdir(wd)

    rows, cols = field.shape

    if rows == cols:                    #field in CS-format 
        lmax = rows - 1
        field = cs2sc(field)
    elif cols - 2 * rows == -1:         #field in SC-format already
        lmax = rows - 1
    else:
        raise Exception("Check format of the field")

    if n == -9999:                      #(no value of n is input) -> set n = lmax
        n = lmax

    if not np.isscalar(n):
        raise Exception("n must be scalar")

    if not np.issubdtype(type(n), np.integer):
        raise Exception("n must be integer")

    if not type(grd) == str:
        raise Exception("Grid argument must be string")

    grd = grd.lower()



    #Grid Definition
    dt = np.pi/n

    if grd == 'pole' or grd == 'mesh':
        theRAD = np.arange(0, np.pi+dt*0.5, dt, dtype='longdouble')
        lamRAD = np.arange(0, 2*np.pi, dt, dtype='longdouble')
    elif grd == 'block' or grd == 'cell':
        theRAD = np.arange(dt/2, np.pi + dt*0.5, dt, dtype='longdouble')
        lamRAD = np.arange(dt/2, 2*np.pi + dt*0.5, dt, dtype='longdouble')
    else:
        raise Exception("Incorrect grid type input")

    nlat = len(theRAD)
    nlon = len(lamRAD)



#    % -------------------------------------------------------------------------
#% Preprocessing on the coefficients: 
#%    - subtract reference field (if jflag is set)
#%    - specific transfer
#%    - upward continuation
#% -------------------------------------------------------------------------

    if jflag:
        field = field - cs2sc(normalklm(lmax+1))

    l = np.arange(0, lmax+1)
    transf = np.array([eigengrav(lmax, quant, h)])[0, :, :].T

    field = field * np.matmul(transf, np.ones((1, 2*lmax+1)), dtype='longdouble')


    '''
    % -------------------------------------------------------------------------
% Size declarations and start the waitbar:
% Note that the definition of dlam causes straight zero-padding in case N > L.
% When N < L, there will be zero-padding up to the smallest integer multiple
% of N larger than L. After the Fourier transformation (see below), the
% proper samples have to be picked out, with stepsize dlam.
% -------------------------------------------------------------------------
    '''

    dlam = int(np.ceil(lmax/n))             #longitude step size
    abcols = dlam*n + 1                     #columns required in A and B
    a = np.zeros((nlat, int(abcols)), dtype='longdouble')
    b = np.zeros((nlat, int(abcols)), dtype='longdouble')



    m = 0
    c = field[m:lmax+1, lmax+m] 
    l = np.array([np.arange(m,lmax+1)])
    p = PLM(l, m, theRAD, nargin = 3, nargout = 1)[:,:,0]
    a[:, m] = np.dot(p,c) 
    b[:, m] = np.zeros(nlat) 



    for m in range(1,lmax+1,1):
        c = field[m:lmax+1,lmax+m]
        s = field[m:lmax+1,lmax-m]

        l = np.array([np.arange(m,lmax+1)])
        p = PLM(l, m, theRAD, nargin = 3, nargout = 1)[:,:,0]
        a[:, m] = np.dot(p,c)
        b[:, m] = np.dot(p,s)

    del field
    '''
     -------------------------------------------------------------------------
 The second synthesis step consists of an inverse Fourier transformation
 over the rows of a and b. 
 In case of 'block', the spectrum has to be shifted first.
 When no zero-padding has been applied, the last b-coefficients must be set to
 zero. Otherwise the number of longitude samples will be 2N+1 instead of 2N.
 For N=L this corresponds to setting SLL=0!
 -------------------------------------------------------------------------
    '''

    if grd =='block' or grd == 'cell': 
      m      = np.arange(0,abcols,1)
      cshift = np.array([np.ones(nlat)], dtype='longdouble').T * np.array([np.cos(m*np.pi/2/n)], dtype='longdouble');	# cshift/sshift describe the 
      sshift = np.array([np.ones(nlat)], dtype='longdouble').T * np.array([np.sin(m*np.pi/2/n)], dtype='longdouble');	# half-blocksize lambda shift.
      atemp  =  cshift*a + sshift*b
      b      = -sshift*a + cshift*b
      a      = atemp



    if np.remainder(n,lmax) == 0:               #Case without zero-padding
        b[:,abcols-1] = np.zeros(nlat)

    #Code for ispec

    f = ispec(a.T, b.T).T
    if dlam > 1: 
        f = f[:,np.arange(1,dlam*nlon+1,dlam)]

    return f, theRAD, lamRAD

Intro to Grace Data Driven Correction

GRACE_Data_Driven_Correction_Vishwakarma(F, cf, GaussianR, basins)

summary

Parameters:

Name	Type	Description	Default
F	_type_	a cell matrix with one column containing SH coefficients	required
cf	_type_	the column in F that contains SH coefficients from GRACE	required
GaussianR	_type_	radius of the Gaussian filter (recommened = 400)	required
basins	_type_	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	required

Raises:

Type	Description
Exception	corrected data-driven time-series (Least Squares fit method)
Exception	corrected data-driven time-series (shift and amplify method)
Exception	gaussian filtered GRACE TWS time-series for all the basins.

Returns:

Name	Type	Description
_type_		description

Todo

• TypeError

Source code in pyshbundle/GRACE_Data_Driven_Correction_Vishwakarma.py

def GRACE_Data_Driven_Correction_Vishwakarma(F, cf, GaussianR, basins):
    """_summary_

    Args:
        F (_type_): a cell matrix with one column containing SH coefficients
        cf (_type_): the column in F that contains SH coefficients from GRACE
        GaussianR (_type_): radius of the Gaussian filter (recommened = 400)
        basins (_type_): 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

    Raises:
        Exception: corrected data-driven time-series (Least Squares fit method)
        Exception: corrected data-driven time-series (shift and amplify method)
        Exception: gaussian filtered GRACE TWS time-series for all the basins.

    Returns:
        _type_: _description_

    Todo:
        + TypeError
    """
    deg = 0.5
    deg_rad = deg_to_rad(deg)

    x = np.linspace(0, 360-deg, int(360/deg))
    y = np.linspace(0, 180-deg, int(180/deg))
    x1 = np.linspace(deg, 360, int(360/deg))
    y1 = np.linspace(deg, 180, int(180/deg))
    lambdd,theta = np.meshgrid(x,y)  
    lambdd1,theta1 = np.meshgrid(x1,y1)

    theta_rad = deg_to_rad(theta)
    theta1_rad = deg_to_rad(theta1)

    #Areahalfdeg = (6378.137**2)*np.power(10,6)*np.pi/180*(np.multiply(a,b)) #Area matrix
    Areahalfdeg = (6378.137**2)*(((np.pi/180)*lambdd1) - ((np.pi/180)*lambdd))*(np.sin((np.pi/2) - theta_rad) - np.sin((np.pi/2) - theta1_rad))

    qty = 'water'

    if type(F) != np.ndarray:
        raise Exception("input GRACE field should be in Numpy Ndarray format, please check guidelines")


    if type(basins) != np.ndarray:
        raise Exception("input basin field should be in Numpy NdArray format, please check guidelines")


    r = F.shape[0] #No of entries in F numpy ndarrray

    cid = 1 #number of river catchments

    f = F[:,cf-1:cf]
    l = f[0][0].shape[0]
    cfield = f[0][0].shape[1]
    if cfield == l:
        flag_cs = 0
    else:
        flag_cs = 1

    Weights = Gaussian(l-1, GaussianR) 
    #gaussian returns weights as a list #gaussian is np.array()

    try: #Broadcase Weights into dimensions
        filter_ = np.ones([1,(2*(l-1))+1]) * Weights
    except:
        w0 = Weights.shape[0]
        Weights = Weights.reshape(w0,1)
        filter_ = np.ones([1,(2*(l-1))+1]) * Weights


    #SH Synthesis
    if l == cfield:
        for m in range(r):
            if flag_cs == 0:
                Ft = cs2sc(f[m][0]).astype('longdouble') 
            else:
                Ft = f[m][0].astype('longdouble') 


            fFld__, _, _ = GSHS(Ft * filter_, qty, 'cell', int(180/deg), 0, 0) 
            ffFld__, _, _ = GSHS((Ft * filter_ * filter_), qty, 'cell', int(180/deg), 0, 0)

            if m == 0:
                fFld = np.zeros((r,fFld__.shape[0],fFld__.shape[1]), dtype = 'longdouble') 
                ffFld = np.zeros((r, ffFld__.shape[0], ffFld__.shape[1]), dtype = 'longdouble')

            fFld[m] = fFld__
            ffFld[m] = ffFld__

        long = 360/deg
        Area = Areahalfdeg
    else:
        raise Exception("enter CS coefficients")



    #Declaration of size of the vectors:
    cid = len(basins) #Here basins is a dictionary with each element storing nd array
    tsleaktotalf = np.zeros([r, cid], dtype = 'longdouble')
    tsleaktotalff = np.zeros([r, cid], dtype = 'longdouble')

    ftsleaktotal = np.zeros([r, cid], dtype = 'longdouble')
    fftsleaktotal = np.zeros([r, cid], dtype = 'longdouble')

    lhat = np.zeros([r, cid], dtype = 'longdouble')

    bfDevRegAv = np.zeros([r, cid], dtype = 'longdouble')
    bbfDevRegAv = np.zeros([r, cid], dtype = 'longdouble')

    FilteredTS = np.zeros([r, cid], dtype = 'longdouble')
    filfilts = np.zeros([r, cid], dtype = 'longdouble')

    leakage = np.zeros([r, cid], dtype = 'longdouble')
    leakager = np.zeros([r, cid], dtype = 'longdouble')   



    for rbasin in range(0, cid):
        #Get the basin functions ready

        #Basin functions, filtered basin function and transfer function Kappa
        Rb = basins[rbasin][0] 
        csRb = gsha(Rb, 'mean', 'block', long/2) 
        csF = cs2sc(csRb[0:l, 0:l]) 
        filRb_ = GSHS(csF * filter_, 'none', 'cell', int(long/2), 0, 0) 
        filRb = filRb_[0]
        kappa = (1-Rb) * filRb



        fF = np.zeros((fFld__.shape[0],fFld__.shape[1]), dtype = 'longdouble')
        ffF = np.zeros((fFld__.shape[0],fFld__.shape[1]), dtype = 'longdouble')
        for m in range(0,r):


            fF = np.concatenate((fFld[m,:,int(fF.shape[1]/2):], fFld[m,:,:int(fF.shape[1]/2)]), axis=1)
            ffF = np.concatenate((ffFld[m,:,int(ffF.shape[1]/2):], ffFld[m,:,:int(ffF.shape[1]/2)]), axis=1)
            #if False:    
            if np.isnan(fF[:20,:20]).any(): #if there is a gap in time series, fill it with NaNs


                tsleaktotalf[m][rbasin] = np.nan
                tsleaktotalff[m][rbasin] = np.nan
                FilteredTS[m][rbasin] = np.nan
                filfilts[m][0:rbasin] = np.nan
                bfDevRegAv[m][rbasin] = np.nan
                bbfDevRegAv[m][0:rbasin] = np.nan

            else:
                #leakage time series from filtered and twice filtered fields
                tsleaktotalf[m][rbasin] = np.sum(fF * kappa * Area) / np.sum(Rb * Area)
                tsleaktotalff[m][rbasin] = np.sum(ffF * kappa * Area) / np.sum(Rb * Area)

                #time series from filtered fields
                FilteredTS[m][rbasin] = np.sum(fF * Rb * Area) / np.sum(Rb * Area)
                filfilts[m][rbasin] = np.sum(ffF * Rb * Area) / np.sum(Rb * Area)

                #Deviation integral timeseries
                bfDevRegAv[m][rbasin] = np.sum((fF * Rb - FilteredTS[m][rbasin]) * filRb * Area) / np.sum(Rb * Area) #working 2022-10-20
                bbfDevRegAv[m][rbasin] = np.sum((ffF * Rb - filfilts[m][rbasin]) * filRb * Area) / np.sum(Rb * Area)
                print(m)





    b = list()
    bl = list()
    for i in range(0, cid):

        A = np.ones([60,2])
        A[:,1] = naninterp(bbfDevRegAv[:, i]) #Pchip interpolate should contain atleast two elements

        lssol_ = sc.linalg.lstsq(A, naninterp(bfDevRegAv[:, i])) #returns a tuple of solution "x", residue and rank of matrix A; for A x = B
        lssol = lssol_[0] 

        b.append(lssol[2-1])


        A = np.ones([60,2])
        A[:,1] = naninterp(tsleaktotalff[:, i])
        lssol_ = sc.linalg.lstsq(A, naninterp(tsleaktotalf[:, i])) #returns a tuple of solution "x", residue and rank of matrix A; for A x = B
        lssol = lssol_[0]
        bl.append(lssol[2-1])
        #Working till here 2022-10-21 1530pm

    multp = npm.repmat(b, r, 1) 
    devint = bfDevRegAv * multp
    multp = npm.repmat(bl, r, 1)
    leakLS = tsleaktotalf * multp


    ps = PhaseCalc(tsleaktotalf,tsleaktotalff)



    #Compute the near true leakage

    for i in range(0, cid):   
        ftsleaktotal[:,i] = naninterp(tsleaktotalf[:,i]) #Replaces gaps (NaN values) with an itnerpolated value in the leakage time series from once filtered fields
        fftsleaktotal[:,i] = naninterp(tsleaktotalff[:,i]) #replace the gaps (NaN values) with an interpolated value in leakage time series from twice filtered fields

        X = sc.fft.fft(ftsleaktotal[:,i]) #take fast Fourier transform #check shape of X 2022-10-21
        p = -ps[0,i] / r #compute the fraction of the time period by which the time series is to be shiftes
        Y = np.exp(1j * np.pi * p * ((np.arange(r)) - r/2) / r) #compute the Conjugate-Symmetric shift 
        Z = X * Y #Apply the shift

        a = sc.fft.ifft(Z) #apply inverse fft

        con = np.conj(a)

        s = a + con

        z = s/2

        leakage[:,i] = z #shifted timeseries


        #Shift timeseriecs from once filtered fields in the direction of the time series from twice filtered fields, to later compute the amplitude ratio
        p = ps[0,i] / r #Fraction of a time period to shift data
        Y = np.exp(1j * np.pi * p * ((np.arange(r)) - r/2) / r) #compute the Conjugate-Symmetric shift
        Z = X * Y

        a = sc.fft.ifft(Z) #apply inverse fft

        con = np.conj(a)

        s = a + con

        z = s/2

        leakager[:,i] = z #shifted timeseries



    #compute the ratio between the amplitude of the shifted leakage from once filtered fields and leakage from twice filtered fields
    rfn = leakage/fftsleaktotal
    rfn[(rfn) >= 2] = 1
    rfn[(rfn) <= -2] = -1
    rfn = np.sum(np.abs(rfn), axis = 0)
    rfn=rfn/r # amplitude ratio


    lhat = leakager * rfn #apply the amplitude ratio to the shifted leakage timeseries from the once filtered fields to get the near true leakage
    lhat[np.isnan(FilteredTS)] = np.nan #reintroduce nan for data gaps
    leakLS[np.isnan(FilteredTS)] = np.nan
    RecoveredTWS = FilteredTS - leakLS - devint
    RecoveredTWS2 = FilteredTS - lhat - devint

    return RecoveredTWS, RecoveredTWS2, FilteredTS

deg_to_rad(deg)

Converts angle from degree to radian

Parameters:

Name	Type	Description	Default
deg	float	Angle in degree	required

Returns:

Name	Type	Description
float		Angle in Radian

Todo

• Inbuilt function available in numpy module

Source code in pyshbundle/GRACE_Data_Driven_Correction_Vishwakarma.py

def deg_to_rad(deg: float):
    """Converts angle from degree to radian

    Args:
        deg (float): Angle in degree

    Returns:
        float: Angle in Radian

    Todo:
        + Inbuilt function available in numpy module
    """
    return deg * np.pi/180

Hydrological Applications with GRACE

BasinAvg(data, path, c_rs, m, gs)

Computes the TWSA time-series for a given basin shape file, using the SH data.

Parameters:

Name	Type	Description	Default
data	Dataset	xarray dataset with format - {coordinates: [time, lat, lon], Data variables: [tws]}	required
path	str	valid path to the basin shape file with extension (.shp)	required
c_rs	crs	the crs into which the dataframe must be transformed (related to salem module)	required
m	int	number of files read	required
gs	float	grid size	required

Returns:

Name	Type	Description
		xarray.Dataset: basin averaged values of TWS
_type_		description

Source code in pyshbundle/basin_avg.py

def BasinAvg(data, path: str, c_rs, m, gs):
    """Computes the TWSA time-series for a given basin shape file, using the SH data.

    Args:
        data (xarray.Dataset): xarray dataset with format - {coordinates: [time, lat, lon], Data variables: [tws]}
        path (str): valid path to the basin shape file with extension (.shp)
        c_rs (crs): the crs into which the dataframe must be transformed (related to salem module)
        m (int): number of files read
        gs (float): grid size 

    Returns:
        xarray.Dataset: basin averaged values of TWS
        _type_: _description_
    """
    shdf = salem.read_shapefile(path)
    shdf.crs
    shdf.plot()
    shdf_area = sum(shdf.to_crs(c_rs).area)
    print('Area of basin in km2:',shdf_area/1e6)
    if shdf_area < 63*1e9:
        print('Warning basin too small for GRACE data application')

    tws_val = data.tws.values
    dates = data.time
    lat,lon = data.lat, data.lon
    lat_shape, lon_shape = data.tws.shape[1],data.tws.shape[2]

    # Calculation of area of each corresponding to  the latitudes and longitudes
    # not sure if ';' is proper syntax may be the octave residu

    deg = gs
    x = np.linspace(0, 359+(1-deg), int(360/deg), dtype='double')
    y = np.linspace(0, 179+(1-deg), int(180/deg), dtype='double')
    x1 = np.linspace(1*deg, 360, int(360/deg), dtype='double')
    y1 = np.linspace(1*deg, 180, int(180/deg), dtype='double')
    lambd,theta = np.meshgrid(x,y)
    lambd1,theta1 = np.meshgrid(x1,y1)
    a = np.sin(np.deg2rad(90-theta))-np.sin(np.deg2rad(90-theta1))
    b = (lambd1 - lambd)*np.pi/180


    # Area of each grid (360*720)
    area = (6378.137**2)*pow(10, 6)*(np.multiply(a, b))        # units m^2
    tot_area = np.sum(np.sum(area))
    tws_m = np.zeros([m, lat_shape, lon_shape])
    for i in range(0,m,1):
        tws_m[i, :, :] = np.multiply(tws_val[i, :, :],area)
    ds_area_w = xr.Dataset(
    data_vars=dict(
        tws=(["time","lat", "lon"], tws_m)
    ),
    coords = {
        "time":dates,
        "lat":lat,
        "lon":lon },
    attrs=dict(description="TWS Anomaly corresponding to long term (2004-2010) mean \n lmax=96 and half radius of Gaussian filter = 500Km"),
    )

    ds_area_w_clp= ds_area_w.salem.roi(shape=shdf)
    # Time series for the whole basin(shapefile) in user defined range
    alpha = ds_area_w_clp.tws.sum(dim=('lon','lat'))/shdf_area
    fig,ax = plt.subplots(figsize=(15,5))
    alpha.plot(ax=ax, color='b')
    ax.set_box_aspect(0.33)
    ax.set_title('Time series for the basin', size=15)
    ax.set_ylabel('TWS anomaly in mm ', size=15)
    plt.tight_layout()

    return alpha, ds_area_w

TWSCalc(data, lmax, gs, r, m)

summary

Parameters:

Name	Type	Description	Default
data	ndarray	SC coefficients	required
lmax	int	maximum degree	required
gs	float	grid size	required
r	_type_	description	required
m	_type_	description	required

Source code in pyshbundle/tws_cal.py

def TWSCalc(data, lmax: int, gs: float, r, m):
    """_summary_

    Args:
        data (np.ndarray): SC coefficients
        lmax (int): maximum degree
        gs (float): grid size
        r (_type_): _description_
        m (_type_): _description_
    """
    SC = data

    gfilter = Gaussian(lmax,r)
    grid_y = int(180/gs)
    grid_x = int(360/gs)
    tws_f = np.zeros([m,grid_y,grid_x], dtype ='longdouble')
    for i in tqdm(range(0,m,1)):
        field = SC[i,0:lmax+1,96-lmax:96+lmax+1]
        shfil = np.zeros([lmax+1,2*lmax+1])

        for j in range(0,2*lmax+1,1):
            shfil[:,j] = gfilter[:,0] * field[:,j]

        quant = 'water' 
        grd = 'cell'
        n = int(180/gs) 
        h = 0 
        jflag = 0

        ff = GSHS(shfil, quant, grd, n, h, jflag)[0]

        ff = ff*1000
        tws_f[i,:,0:int(grid_x/2)] = ff[:,int(grid_x/2):]
        tws_f[i,:,int(grid_x/2):] = ff[:,0:int(grid_x/2)]   

    plt.imshow(tws_f[0,:,:])
    return(tws_f)