feat: added notebook for visualizing harmonic errors (#74)

Author: tsutterley

doc/source/user_guide/utilities.rst

@@ -55,6 +55,8 @@ General Methods
 
 .. autofunction:: gravity_toolkit.utilities.s3_client
 
+.. autofunction:: gravity_toolkit.utilities.s3_bucket
+
 .. autofunction:: gravity_toolkit.utilities.s3_key
 
 .. autofunction:: gravity_toolkit.utilities.check_credentials

gravity_toolkit/degree_amplitude.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 degree_amplitude.py
-Written Tyler Sutterley (04/2022)
+Written Tyler Sutterley (05/2022)
 
 Calculates the amplitude of each spherical harmonic degree
 
@@ -20,6 +20,7 @@
     numpy: Scientific Computing Tools For Python (https://numpy.org)
 
 UPDATE HISTORY:
+    Updated 05/2022: use numpy atleast_3d to add singleton dimensions
     Updated 04/2022: updated docstrings to numpy documentation format
     Updated 07/2020: added function docstrings
     Updated 05/2020: add singleton dimension to calculate time series amplitudes
@@ -49,9 +50,8 @@ def degree_amplitude(clm, slm, LMAX=None, MMAX=None):
         degree amplitude
     """
     #-- add a singleton dimension to input harmonics
-    if (np.ndim(clm) == 2):
-        clm = np.copy(clm[:,:,None])
-        slm = np.copy(slm[:,:,None])
+    clm = np.atleast_3d(clm)
+    slm = np.atleast_3d(slm)
     #-- check shape
     LMp1,MMp1,nt = np.shape(clm)
 
@@ -64,8 +64,8 @@ def degree_amplitude(clm, slm, LMAX=None, MMAX=None):
 
     #-- allocating for output array
     amp = np.zeros((LMAX+1,nt))
-    for l in range(LMAX+1):#-- m:LMAX
-        m = np.arange(l,MMAX+1)
+    for l in range(LMAX+1):
+        m = np.arange(0,MMAX+1)
         #-- degree amplitude of spherical harmonic degree
         amp[l,:] = np.sqrt(np.sum(clm[l,m,:]**2 + slm[l,m,:]**2,axis=0))
 

gravity_toolkit/grace_months_index.py

@@ -41,6 +41,7 @@
 
 UPDATE HISTORY:
     Updated 05/2022: use argparse descriptions within documentation
+        use new GSFC release 6 version 2 mascons as the default
     Updated 04/2022: updated docstrings to numpy documentation format
     Updated 09/2021: use functions for converting to and from GRACE months
     Updated 05/2021: define int/float precision to prevent deprecation warning
@@ -67,7 +68,7 @@
 import numpy as np
 from gravity_toolkit.time import grace_to_calendar
 
-def grace_months_index(base_dir, DREL=['RL06','v02.4'], MODE=None):
+def grace_months_index(base_dir, DREL=['RL06','rl06v2.0'], MODE=None):
     """
     Creates a file with the start and end days for each dataset
 
@@ -206,7 +207,7 @@ def arguments():
     #-- GRACE/GRACE-FO data release
     parser.add_argument('--release','-r',
         metavar='DREL', type=str, nargs='+',
-        default=['RL06','v02.4'],
+        default=['RL06','rl06v2.0'],
         help='GRACE/GRACE-FO Data Release')
     #-- permissions mode of the local directories and files (number in octal)
     parser.add_argument('--mode','-M',

gravity_toolkit/read_GRACE_harmonics.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 read_GRACE_harmonics.py
-Written by Tyler Sutterley (04/2022)
+Written by Tyler Sutterley (05/2022)
 Contributions by Hugo Lecomte
 
 Reads GRACE files and extracts spherical harmonic data and drift rates (RL04)
@@ -42,6 +42,7 @@
     time.py: utilities for calculating time operations
 
 UPDATE HISTORY:
+    Updated 05/2022: updated comments
     Updated 04/2022: updated docstrings to numpy documentation format
         include utf-8 encoding in reads to be windows compliant
         check if GRACE/GRACE-FO data file is present in file-system
@@ -105,7 +106,7 @@ def read_GRACE_harmonics(input_file, LMAX, MMAX=None, POLE_TIDE=False):
     eslm: float
         sine spherical harmonic uncalibrated standard deviations
     header: str
-        Header text from the GRACE file
+        Header text from the GRACE/GRACE-FO file
 
     References
     ----------
@@ -118,9 +119,12 @@ def read_GRACE_harmonics(input_file, LMAX, MMAX=None, POLE_TIDE=False):
 
     #-- parse filename
     PFX,SY,SD,EY,ED,N,PRC,F1,DRL,F2,SFX = parse_file(input_file)
-    file_contents = extract_file(input_file, (SFX=='.gz'))
+    #-- check if file is compressed
+    compressed = (SFX == '.gz')
+    #-- extract file contents
+    file_contents = extract_file(input_file, compressed)
 
-    #-- JPL Mascon solutions
+    #-- JPL mascon solutions in spherical harmonic form
     if PRC in ('JPLMSC',):
         DSET = 'GSM'
         DREL = np.int64(DRL)
@@ -134,18 +138,18 @@ def read_GRACE_harmonics(input_file, LMAX, MMAX=None, POLE_TIDE=False):
     #-- COST-G unfiltered combination solutions
     #-- https://doi.org/10.5880/ICGEM.COST-G.001
     elif PRC in ('COSTG',):
-        DSET, = re.findall('GSM|GAC',PFX)
+        DSET, = re.findall(r'GSM|GAC',PFX)
         DREL = np.int64(DRL)
         FLAG = r'gfc'
-    #-- Standard GRACE solutions
+    #-- Standard GRACE/GRACE-FO Level-2 solutions
     else:
         DSET = PFX
         DREL = np.int64(DRL)
         FLAG = r'GRCOF2'
 
-    #-- output python dictionary with GRACE data and date information
+    #-- output python dictionary with GRACE/GRACE-FO data and metadata
     grace_L2_input = {}
-    #-- extract GRACE date information from input file name
+    #-- extract GRACE/GRACE-FO date information from input file name
     start_yr = np.float64(SY)
     end_yr = np.float64(EY)
     start_day = np.float64(SD)
@@ -230,10 +234,9 @@ def read_GRACE_harmonics(input_file, LMAX, MMAX=None, POLE_TIDE=False):
     #-- if drift rates exist at any time, will add to harmonics
     #-- Will convert the secular rates into a stokes contribution
     #-- Currently removes 2003.3 to get the temporal average close to 0.
-    #-- note: += means grace_xlm = grace_xlm + drift_x
     if ((DREL == 4) and (DSET == 'GSM')):
         #-- time since 2003.3
-        dt = (grace_L2_input['time']-2003.3)
+        dt = (grace_L2_input['time'] - 2003.3)
         grace_L2_input['clm'][:,:] += dt*drift_c[:,:]
         grace_L2_input['slm'][:,:] += dt*drift_s[:,:]
 
@@ -261,8 +264,7 @@ def read_GRACE_harmonics(input_file, LMAX, MMAX=None, POLE_TIDE=False):
             #-- before computing their harmonic solutions
             C21_PT = -1.551e-9*(m1 - 0.62e-3*dt) - 0.012e-9*(m2 + 3.48e-3*dt)
             S21_PT = 0.021e-9*(m1 - 0.62e-3*dt) - 1.505e-9*(m2 + 3.48e-3*dt)
-            #-- correct GRACE spherical harmonics for pole tide
-            #-- note: -= means grace_xlm = grace_xlm - PT
+            #-- correct GRACE/GRACE-FO spherical harmonics for pole tide
             grace_L2_input['clm'][2,1] -= C21_PT
             grace_L2_input['slm'][2,1] -= S21_PT
         #-- GFZ Pole Tide Correction
@@ -271,13 +273,13 @@ def read_GRACE_harmonics(input_file, LMAX, MMAX=None, POLE_TIDE=False):
             #-- GFZ removes only a constant pole position
             C21_PT = -1.551e-9*(-0.62e-3*dt) - 0.012e-9*(3.48e-3*dt)
             S21_PT = 0.021e-9*(-0.62e-3*dt) - 1.505e-9*(3.48e-3*dt)
-            #-- correct GRACE spherical harmonics for pole tide
-            #-- note: -= means grace_xlm = grace_xlm - PT
+            #-- correct GRACE/GRACE-FO spherical harmonics for pole tide
             grace_L2_input['clm'][2,1] -= C21_PT
             grace_L2_input['slm'][2,1] -= S21_PT
 
-    #-- return the GRACE data, GRACE date (mid-month in decimal), and the
-    #-- start and end days as Julian dates
+    #-- return the header data, GRACE/GRACE-FO data
+    #-- GRACE/GRACE-FO date (mid-month in decimal)
+    #-- and the start and end days as Julian dates
     return grace_L2_input
 
 #-- PURPOSE: extract parameters from filename

gravity_toolkit/utilities.py

@@ -1,14 +1,15 @@
 #!/usr/bin/env python
 u"""
 utilities.py
-Written by Tyler Sutterley (04/2022)
+Written by Tyler Sutterley (05/2022)
 Download and management utilities for syncing time and auxiliary files
 
 PYTHON DEPENDENCIES:
     lxml: processing XML and HTML in Python
         https://pypi.python.org/pypi/lxml
 
 UPDATE HISTORY:
+    Updated 05/2022: function for extracting bucket name from presigned url
     Updated 04/2022: updated docstrings to numpy documentation format
         update CMR queries to prepare for version 1 of RL06
     Updated 03/2022: add NASA Common Metadata Repository (CMR) queries
@@ -789,6 +790,25 @@ def s3_client(HOST=None, timeout=None, region_name='us-west-2'):
     #-- return the AWS client for region
     return client
 
+#-- PURPOSE: get a s3 bucket name from a presigned url
+def s3_bucket(presigned_url):
+    """
+    Get a s3 bucket name from a presigned url
+
+    Parameters
+    ----------
+    presigned_url: str
+        s3 presigned url
+
+    Returns
+    -------
+    bucket: str
+        s3 bucket name
+    """
+    host = url_split(presigned_url)
+    bucket = re.sub(r's3:\/\/', r'', host[0], re.IGNORECASE)
+    return bucket
+
 #-- PURPOSE: get a s3 bucket key from a presigned url
 def s3_key(presigned_url):
     """
@@ -798,9 +818,15 @@ def s3_key(presigned_url):
     ----------
     presigned_url: str
         s3 presigned url
+
+    Returns
+    -------
+    key: str
+        s3 bucket key for object
     """
     host = url_split(presigned_url)
-    return posixpath.join(*host[1:])
+    key = posixpath.join(*host[1:])
+    return key
 
 #-- PURPOSE: check that entered NASA Earthdata credentials are valid
 def check_credentials(HOST='https://podaac-tools.jpl.nasa.gov/drive/files'):