fix: version check for mission in GFZ ftp sync (#87)

feat: add spatial gradient programs

doc/source/api_reference/fourier_legendre.rst

@@ -17,3 +17,5 @@ Calling Sequence
 .. __: https://github.com/tsutterley/read-GRACE-harmonics/blob/main/gravity_toolkit/fourier_legendre.py
 
 .. autofunction:: gravity_toolkit.fourier_legendre
+
+.. autofunction:: gravity_toolkit.legendre_gradient

doc/source/api_reference/harmonic_gradients.rst

@@ -0,0 +1,19 @@
+=====================
+harmonic_gradients.py
+=====================
+
+- Calculates the zonal and meridional gradients of a scalar field from a series of spherical harmonics
+
+Calling Sequence
+################
+
+.. code-block:: python
+
+    from gravity_toolkit.harmonic_gradients import harmonic_gradients
+    gradients = harmonic_gradients(clm,slm,lon,lat,LMAX=60)
+
+`Source code`__
+
+.. __: https://github.com/tsutterley/read-GRACE-harmonics/blob/main/gravity_toolkit/harmonic_gradients.py
+
+.. autofunction:: gravity_toolkit.harmonic_gradients

doc/source/index.rst

@@ -51,6 +51,7 @@ missions
     api_reference/grace_find_months.rst
     api_reference/grace_input_months.rst
     api_reference/grace_months_index.rst
+    api_reference/harmonic_gradients.rst
     api_reference/harmonic_summation.rst
     api_reference/harmonics.rst
     api_reference/legendre.rst

gravity_toolkit/__init__.py

@@ -23,7 +23,7 @@
 from gravity_toolkit.clenshaw_summation import clenshaw_summation
 from gravity_toolkit.degree_amplitude import degree_amplitude
 from gravity_toolkit.destripe_harmonics import destripe_harmonics
-from gravity_toolkit.fourier_legendre import fourier_legendre
+from gravity_toolkit.fourier_legendre import fourier_legendre, legendre_gradient
 from gravity_toolkit.gauss_weights import gauss_weights
 from gravity_toolkit.gen_averaging_kernel import gen_averaging_kernel
 from gravity_toolkit.gen_disc_load import gen_disc_load
@@ -37,6 +37,7 @@
 from gravity_toolkit.grace_input_months import grace_input_months, read_ecmwf_corrections
 from gravity_toolkit.grace_months_index import grace_months_index
 from gravity_toolkit.harmonics import harmonics
+from gravity_toolkit.harmonic_gradients import harmonic_gradients
 from gravity_toolkit.harmonic_summation import harmonic_summation
 from gravity_toolkit.legendre_polynomials import legendre_polynomials
 from gravity_toolkit.legendre import legendre

gravity_toolkit/fourier_legendre.py

@@ -2,7 +2,7 @@
 u"""
 fourier_legendre.py
 Original IDL code gen_plms.pro written by Sean Swenson
-Adapted by Tyler Sutterley (04/2022)
+Adapted by Tyler Sutterley (10/2022)
 
 Computes Fourier coefficients of the associated Legendre functions
 
@@ -20,6 +20,7 @@
     numpy: Scientific Computing Tools For Python (https://numpy.org)
 
 UPDATE HISTORY:
+    Updated 10/2022: add polynomial function for calculating gradients
     Updated 04/2022: updated docstrings to numpy documentation format
     Updated 09/2021: cleaned up program for public release
     Updated 07/2020: added function docstrings
@@ -42,7 +43,7 @@ def fourier_legendre(lmax, mmax):
 
     Returns
     -------
-    plms: float
+    plm: float
         Fourier coefficients
     """
 
@@ -208,3 +209,82 @@ def fourier_legendre(lmax, mmax):
 
     #-- return the fourier coefficients
     return plm
+
+def legendre_gradient(lmax, mmax):
+    """
+    Calculates functions for evaluating the integral of a
+        product of two fourier series
+
+    Parameters
+    ----------
+    lmax: int
+        Upper bound of Spherical Harmonic Degrees
+    mmax: int
+        Upper bound of Spherical Harmonic Orders
+
+    Returns
+    -------
+    vlm: float
+        Fourier coefficients for meridional gradients
+    wlm: float
+        Fourier coefficients for zonal gradients
+    """
+
+    plm = fourier_legendre(lmax, mmax)
+    vlm = np.zeros((lmax+1,lmax+1,lmax+1))
+    wlm = np.zeros((lmax+1,lmax+1,lmax+1))
+
+    #-- l=0 zero by definition
+    lind = np.arange(1,lmax+1)
+    #-- m=0 special case
+    #-- terms with m=0, m=1 have different coefficients
+    vlm[lind,0,:] = 2.0*np.dot(np.diag(np.sqrt((lind+1)*lind/2.0)), plm[lind,1,:])
+
+    #-- m+1 terms
+    for l in range(2,lmax+1):#-- from 2 to lmax
+        m = np.arange(1,l)#-- from 1 to l-1
+        lplus = np.arange(l+2,2*l+1)#-- from l+2 to 2*l
+        lminus = np.arange(l-1,0,-1)#-- from l-1 to 1
+        vlm[l,m,:] = np.dot(np.diag(np.sqrt(lplus*lminus/4.0)), plm[l,m+1,:])
+
+    #-- m-1 terms, m-1=0 has different coefficients
+    vlm[lind,1,:] -= np.dot(np.diag(np.sqrt((lind+1)*lind/2.0)), plm[lind,0,:])
+
+    for l in range(2,lmax+1):
+        m = np.arange(2,l+1)#-- from 2 to l
+        lplus = np.arange(l+2,2*l+1)#-- from l+2 to 2*l
+        lminus = np.arange(l-1,0,-1)#-- from l-1 to 1
+        vlm[l,m,:] -= np.dot(np.diag(np.sqrt(lplus*lminus/4.0)), plm[l,m-1,:])
+    #-- normalizations
+    for l in range(1,lmax+1):
+        vlm[l,:,:] /= np.sqrt((l+1)*l)
+
+    #-- m+1 terms
+    for l in range(2, lmax+1):
+        m = np.arange(1,l)#-- from 1 to l-1
+        dfactor = (2.0*l+1.0)/(2.0*l-1.0)
+        lminus2 = np.arange(l-2,-1,-1)#-- from l-2 to 0
+        lminus1 = np.arange(l-1,0,-1)#-- from l-1 to 1
+        wlm[l,m,:] = np.sqrt(dfactor) * \
+            np.dot(np.diag(np.sqrt(lminus2*lminus1/4.0)), plm[l,m+1,:])
+
+    #-- m-1 terms, m-1=0 has different coefficients
+    #-- m=1 term
+    for l in range(1, lmax+1):
+        dfactor = (2.0*l+1.0)/(2.0*l-1.0)
+        wlm[l,1,:] += np.sqrt(dfactor)*np.sqrt(l*(l+1)/2.0)*plm[l-1,0,:]
+
+    for l in range(2,lmax+1):
+        m = np.arange(2,l+1)
+        dfactor = (2.0*l+1.0)/(2.0*l-1.0)
+        lplus2 = np.arange(l+2,2*l+1)#-- from l+2 to 2*l
+        lplus1 = np.arange(l+1,2*l)#-- from l+1 to (2*l-1)
+        wlm[l,m,:] += np.sqrt(dfactor) * \
+            np.dot(np.diag(np.sqrt(lplus2*lplus1)/4.0), plm[l-1,m-1,:])
+    #-- normalizations
+    for l in range(1, lmax+1):
+        wlm[l,:,:] /= np.sqrt((l+1)*l)
+    #-- normalize vlm
+    vlm[:,0,:] /= 2.0
+
+    return (vlm, wlm)

gravity_toolkit/gen_harmonics.py

@@ -396,3 +396,4 @@ def fourier(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs):
 
     #-- return the output spherical harmonics object
     return Ylms
+

gravity_toolkit/harmonic_gradients.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+u"""
+harmonic_gradients.py
+Original IDL code calc_grad.pro written by Sean Swenson
+Adapted by Tyler Sutterley (10/2022)
+
+Calculates the zonal and meridional gradients of a scalar field
+from a series of spherical harmonics
+
+CALLING SEQUENCE:
+    gradient = harmonic_gradients(clm, slm, lon, lat)
+
+INPUTS:
+    clm1: cosine spherical harmonic coefficients in output units
+    slm1: sine spherical harmonic coefficients in output units
+    lon: longitude array for output spatial field
+    lat: latitude array for output spatial field
+
+OPTIONS:
+    LMIN: Lower bound of Spherical Harmonic Degrees
+    LMAX: Upper bound of Spherical Harmonic Degrees
+    MMAX: Upper bound of Spherical Harmonic Orders (default = LMAX)
+
+PYTHON DEPENDENCIES:
+    numpy: Scientific Computing Tools For Python (https://numpy.org)
+
+PROGRAM DEPENDENCIES:
+    fourier_legendre.py: Computes the Fourier coefficients of the associated
+        Legendre functions
+
+UPDATE HISTORY:
+    Updated 10/2022: cleaned up program for public release
+    Updated 07/2020: added function docstrings
+    Updated 06/2019: using Python3 compatible division
+    Updated 05/2015: code updates
+    Written 05/2013
+"""
+from __future__ import division
+import numpy as np
+from gravity_toolkit.fourier_legendre import legendre_gradient
+
+def harmonic_gradients(clm1, slm1, lon, lat,
+    LMIN=0, LMAX=60, MMAX=None):
+    """
+    Calculates the gradient of a scalar field from a series of
+    spherical harmonics
+
+    Parameters
+    ----------
+    clm1: float
+        cosine spherical harmonic coefficients in output units
+    slm1: float
+        sine spherical harmonic coefficients in output units
+    lon: float
+        longitude array
+    lat: float
+        latitude array
+    LMIN: int, default 0
+        Lower bound of Spherical Harmonic Degrees
+    LMAX: int, default 60
+        Upper bound of Spherical Harmonic Degrees
+    MMAX: int or NoneType, default None
+        Upper bound of Spherical Harmonic Orders
+
+    Returns
+    -------
+    gradients: float
+        zonal and meridional gradient fields
+    """
+
+    #-- if LMAX is not specified, will use the size of the input harmonics
+    if (LMAX == 0):
+        LMAX = np.shape(clm1)[0]-1
+    #-- upper bound of spherical harmonic orders (default = LMAX)
+    if MMAX is None:
+        MMAX = np.copy(LMAX)
+
+    #-- Longitude in radians
+    phi = (np.squeeze(lon)*np.pi/180.0)[np.newaxis,:]
+    #-- Colatitude in radians
+    th = (90.0 - np.squeeze(lat))*np.pi/180.0
+    thmax = len(np.squeeze(lat))
+    phimax = len(np.squeeze(lon))
+
+    #-- Truncating harmonics to degree and order LMAX
+    #-- removing coefficients below LMIN and above MMAX
+    mm = np.arange(0,MMAX+1)
+    clm = np.zeros((LMAX+1,MMAX+1))
+    slm = np.zeros((LMAX+1,MMAX+1))
+    clm[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm]
+    slm[LMIN:LMAX+1,mm] = slm1[LMIN:LMAX+1,mm]
+    #-- spherical harmonic degree and order
+    ll = np.arange(0,LMAX+1)[np.newaxis, :]#-- lmax+1 to include lmax
+    mm = np.arange(0,MMAX+1)[:, np.newaxis]#-- mmax+1 to include mmax
+
+    #-- generate Vlm coefficients (vlm and wlm)
+    vlm, wlm = legendre_gradient(LMAX, MMAX)
+
+    dlm = np.zeros((LMAX+1,LMAX+1,2))
+    #--  minus sign is because lat and theta change with opposite sign
+    for l in range(0,LMAX+1):
+        dlm[l,:,0] = -clm[l,:]*np.sqrt((l+1.0)*l)
+        dlm[l,:,1] = -slm[l,:]*np.sqrt((l+1.0)*l)
+
+    m_even = np.arange(0,MMAX+2,2)
+    m_odd = np.arange(1,MMAX,2)
+
+    #--  Calculate fourier coefficients from legendre coefficients
+    d_cos = np.zeros((LMAX+1,thmax,2))
+    d_sin = np.zeros((LMAX+1,thmax,2))
+    cnk = np.cos(np.dot(th[:,np.newaxis],ll))
+    snk = np.sin(np.dot(th[:,np.newaxis],ll))
+
+    wtmp = np.zeros((len(m_even),LMAX+1,2))
+    vtmp = np.zeros((len(m_even),LMAX+1,2))
+    #-- m = even terms (vlm,wlm sine series)
+    for n in range(0,LMAX+1):
+        wtmp[:,n,0] = np.sum(wlm[:,m_even,n]*dlm[:,m_even,0],axis=0)
+        wtmp[:,n,1] = np.sum(wlm[:,m_even,n]*dlm[:,m_even,1],axis=0)
+        vtmp[:,n,0] = np.sum(vlm[:,m_even,n]*dlm[:,m_even,0],axis=0)
+        vtmp[:,n,1] = np.sum(vlm[:,m_even,n]*dlm[:,m_even,1],axis=0)
+
+    d_cos[m_even,:,0] = np.dot(wtmp[:,:,1],np.transpose(snk))
+    d_sin[m_even,:,0] = np.dot(-wtmp[:,:,0],np.transpose(snk))
+    d_cos[m_even,:,1] = np.dot(vtmp[:,:,1],np.transpose(snk))
+    d_sin[m_even,:,1] = np.dot(-vtmp[:,:,0],np.transpose(snk))
+
+    #-- m = odd terms (vlm,wlm cosine series)
+    wtmp = np.zeros((len(m_odd),LMAX+1,2))
+    vtmp = np.zeros((len(m_odd),LMAX+1,2))
+    for n in range(0,LMAX+1):
+        wtmp[:,n,0] = np.sum(wlm[:,m_odd,n]*dlm[:,m_odd,0],axis=0)
+        wtmp[:,n,1] = np.sum(wlm[:,m_odd,n]*dlm[:,m_odd,1],axis=0)
+        vtmp[:,n,0] = np.sum(vlm[:,m_odd,n]*dlm[:,m_odd,0],axis=0)
+        vtmp[:,n,1] = np.sum(vlm[:,m_odd,n]*dlm[:,m_odd,1],axis=0)
+
+    d_cos[m_odd,:,0] = np.dot(wtmp[:,:,1],np.transpose(cnk))
+    d_sin[m_odd,:,0] = np.dot(-wtmp[:,:,0],np.transpose(cnk))
+    d_cos[m_odd,:,1] = np.dot(vtmp[:,:,1],np.transpose(cnk))
+    d_sin[m_odd,:,1] = np.dot(-vtmp[:,:,0],np.transpose(cnk))
+
+    #-- Calculating cos(m*phi) and sin(m*phi)
+    ccos = np.cos(np.dot(mm,phi))
+    ssin = np.sin(np.dot(mm,phi))
+    #-- Final signal recovery from fourier coefficients
+    gradients = np.zeros((phimax,thmax,2))
+    gradients[:,:,0] = np.dot(np.transpose(ccos), d_cos[:,:,0]) + \
+        np.dot(np.transpose(ssin), d_sin[:,:,0])
+    gradients[:,:,1] = np.dot(np.transpose(ccos), d_cos[:,:,1]) + \
+        np.dot(np.transpose(ssin), d_sin[:,:,1])
+    #-- return the gradient fields
+    return gradients

gravity_toolkit/harmonic_summation.py

@@ -78,7 +78,7 @@ def harmonic_summation(clm1, slm1, lon, lat,
     th = (90.0 - np.squeeze(lat))*np.pi/180.0
     thmax = len(th)
 
-    #--  Calculate fourier coefficients from legendre coefficients
+    #-- Calculate fourier coefficients from legendre coefficients
     d_cos = np.zeros((MMAX+1,thmax))#-- [m,th]
     d_sin = np.zeros((MMAX+1,thmax))#-- [m,th]
     if PLM is None:

scripts/gfz_isdc_grace_ftp.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 gfz_isdc_grace_ftp.py
-Written by Tyler Sutterley (08/2022)
+Written by Tyler Sutterley (10/2022)
 Syncs GRACE/GRACE-FO data from the GFZ Information System and Data Center (ISDC)
 Syncs CSR/GFZ/JPL files for RL06 GAA/GAB/GAC/GAD/GSM
     GAA and GAB are GFZ/JPL only
@@ -40,6 +40,7 @@
     utilities.py: download and management utilities for syncing files
 
 UPDATE HISTORY:
+    Updated 10/2022: fix version check for mission
     Updated 08/2022: moved regular expression function to utilities
         Dynamically select newest version of granules for index
     Updated 04/2022: added option for GRACE/GRACE-FO Level-2 data version
@@ -236,7 +237,7 @@ def gfz_isdc_grace_ftp(DIRECTORY, PROC=[], DREL=[], VERSION=[],
                 #-- for each satellite mission (grace, grace-fo)
                 for i,mi in enumerate(['grace','grace-fo']):
                     #-- modifiers for intermediate data releases
-                    if (int(VERSION) > 0):
+                    if (int(VERSION[i]) > 0):
                         drel_str = '{0}.{1}'.format(rl,VERSION[i])
                     else:
                         drel_str = copy.copy(rl)
@@ -364,7 +365,7 @@ def arguments():
         help='GRACE/GRACE-FO data release')
     #-- GRACE/GRACE-FO data version
     parser.add_argument('--version','-v',
-        metavar='VERSION', type=str, nargs='+',
+        metavar='VERSION', type=str, nargs=2,
         default=['0','1'], choices=['0','1','2','3'],
         help='GRACE/GRACE-FO Level-2 data version')
     #-- GRACE/GRACE-FO newsletters