feat: Add 3-hour AOD interval for RL06 #37 (#41)

feat: release monthly dealiasing (for CSR GAA etc)
fix: inputs to AOD-corrected SLR geocenter coefficients

Author: tsutterley

doc/source/getting_started/GRACE-Data-File-Formats.rst

@@ -5,8 +5,8 @@ Data File Formats
 Product Identifier
 ##################
 
-GRACE Level-2 products consist of spherical harmonic coefficients of the Earth's gravitational field.
-The data files are typically gzipped ascii files with names formatted as the following:
+GRACE and GRACE-FO Level-2 products consist of spherical harmonic coefficients of the Earth's gravitational field.
+The Spherical Harmonic Model (SHM) data files are typically gzipped ascii files with names formatted as the following:
 ``PID-2_YYYYDOY-yyyydoy_ndays_center_flag_rrrr`` or ``PID-2_YYYYDOY-yyyydoy_mssn_center_flag_rrrr``
 
  - ``PID`` is a product identification string (for standard products: GSM, GAD, GAC, GAA, GAB)

doc/source/index.rst

@@ -105,6 +105,7 @@ missions
     user_guide/grace_mean_harmonics.rst
     user_guide/grace_spatial_error.rst
     user_guide/grace_spatial_maps.rst
+    user_guide/make_monthly_dealiasing.md
     user_guide/mascon_reconstruct.rst
     user_guide/regress_grace_maps.rst
     user_guide/run_grace_date.md

doc/source/user_guide/aod1b_geocenter.md

@@ -6,7 +6,7 @@ aod1b_geocenter.py
     * `'ocn'`: oceanic loading from OMCT/MPIOM
     * `'glo'`: global atmospheric and oceanic loading
     * `'oba'`: ocean bottom pressure from OMCT/MPIOM
-- Creates monthly files of geocenter variations at 6-hour intervals
+- Creates monthly files of geocenter variations at 3 or 6-hour intervals
 
 #### Calling Sequence
 ```bash

doc/source/user_guide/aod1b_oblateness.md

@@ -6,7 +6,7 @@ aod1b_oblateness.py
     * `'ocn'`: oceanic loading from OMCT/MPIOM
     * `'glo'`: global atmospheric and oceanic loading
     * `'oba'`: ocean bottom pressure from OMCT/MPIOM
-- Creates monthly files of oblateness variations at 6-hour intervals
+- Creates monthly files of oblateness variations at 3 or 6-hour intervals
 
 #### Calling Sequence
 ```bash

doc/source/user_guide/make_monthly_dealiasing.md

@@ -0,0 +1,30 @@
+make_monthly_dealiasing.py
+==========================
+
+- Reads GRACE/GRACE-FO level-1b dealiasing data files for a specific product and outputs monthly the mean for a specific GRACE/GRACE-FO processing center and data release
+    * `'GAA'`: atmospheric loading from ECMWF
+    * `'GAB'`: oceanic loading from OMCT/MPIOM
+    * `'GAC'`: global atmospheric and oceanic loading
+    * `'GAD'`: ocean bottom pressure from OMCT/MPIOM
+- Creates monthly files of oblateness variations at 3 or 6-hour intervals
+
+#### Calling Sequence
+```bash
+python make_monthly_dealiasing.py --release RL06 --product GAD
+```
+[Source code](https://github.com/tsutterley/read-GRACE-harmonics/blob/main/scripts/make_monthly_dealiasing.py)
+
+#### Command Line Options
+- `-D X`, `--directory X`: Working Data Directory
+- `-c X`, `--center X`: GRACE/GRACE-FO Processing Center for dates
+- `-r X`, `--release X`: GRACE/GRACE-FO Data Release (RL05 or RL06)
+- `-p X`, `--product X`: GRACE/GRACE-FO dealiasing product (GAA, GAB, GAC, GAD)
+- `-L X`, `--lmax X`: Maximum spherical harmonic degree and order for output
+- `-F X`, `--format X`: Output data format
+    * `'ascii'`
+    * `'netCDF4'`
+    * `'HDF5'`
+    * `'SHM'`
+- `-C`, `--clobber`: Overwrite existing data
+- `-M X`, `--mode X`: Permission mode of directories and files
+- `-V`, `--verbose`: Output information for each output file

gravity_toolkit/grace_input_months.py

@@ -1,7 +1,8 @@
 #!/usr/bin/env python
 u"""
 grace_input_months.py
-Written by Tyler Sutterley (06/2021)
+Written by Tyler Sutterley (07/2021)
+Contributions by Hugo Lecomte and Yara Mohajerani
 
 Reads GRACE/GRACE-FO files for a specified spherical harmonic degree and order
     and for a specified date range
@@ -92,6 +93,7 @@
     read_GRACE_harmonics.py: reads an input GRACE data file and calculates date
 
 UPDATE HISTORY:
+    Updated 07/2021: fix inputs to AOD-corrected SLR geocenter coefficients
     Updated 06/2021: can use SLR figure axis harmonics produced by GSFC
         read GRACE/GRACE-FO fields before reading replacement values
     Updated 05/2021: can use SLR low-degree harmonic values produced by GFZ
@@ -373,7 +375,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon,
         #-- new file of degree-1 mass variations from Minkang Cheng
         #-- http://download.csr.utexas.edu/outgoing/cheng/gct2est.220_5s
         DEG1_file = os.path.join(base_dir,'geocenter','gct2est.220_5s')
-        DEG1_input = aod_corrected_SLR_geocenter(DEG1_file,HEADER=15,
+        DEG1_input = aod_corrected_SLR_geocenter(DEG1_file,DREL,HEADER=15,
             RADIUS=6.378136e9,COLUMNS=['MJD','time','X','Y','Z','XM','YM','ZM',
             'X_sigma','Y_sigma','Z_sigma','XM_sigma','YM_sigma','ZM_sigma'])
         FLAGS.append('_w{0}_DEG1'.format(DEG1))

scripts/aod1b_geocenter.py

@@ -1,15 +1,16 @@
 #!/usr/bin/env python
 u"""
 aod1b_geocenter.py
-Written by Tyler Sutterley (05/2021)
+Written by Tyler Sutterley (07/2021)
+Contributions by Hugo Lecomte (03/2021)
 
 Reads GRACE/GRACE-FO level-1b dealiasing data files for a specific product
     atm: atmospheric loading from ECMWF
     ocn: oceanic loading from OMCT/MPIOM
     glo: global atmospheric and oceanic loading
     oba: ocean bottom pressure from OMCT/MPIOM
 
-Creates monthly files of geocenter variations at 6-hour intervals
+Creates monthly files of geocenter variations at 3 or 6 hour intervals
 
 NOTE: this reads the GFZ AOD1B files downloaded from PO.DAAC
 https://podaac-uat.jpl.nasa.gov/drive/files/allData/grace/L1B/GFZ/AOD1B/RL06/
@@ -30,9 +31,12 @@
 
 PROGRAM DEPENDENCIES:
     geocenter.py: converts degree 1 spherical harmonics to geocenter variations
+    utilities.py: download and management utilities for files
 
 UPDATED HISTORY:
+    Updated 07/2021: can use default argument files to define options
     Updated 05/2021: define int/float precision to prevent deprecation warning
+    Updated 03/2021: Add 3-hour interval depending on Release
     Updated 10/2020: use argparse to set command line parameters
     Updated 07/2020: added function docstrings
     Updated 06/2019: using python3 compatible regular expression patterns
@@ -58,19 +62,13 @@
 import argparse
 import numpy as np
 from gravity_toolkit.geocenter import geocenter
-
-#-- aod1b data products
-product = {}
-product['atm'] = 'Atmospheric loading from ECMWF'
-product['ocn'] = 'Oceanic loading from OMCT'
-product['glo'] = 'Global atmospheric and oceanic loading'
-product['oba'] = 'Ocean bottom pressure from OMCT'
+import gravity_toolkit.utilities as utilities
 
 #-- program module to read the degree 1 coefficients of the AOD1b data
 def aod1b_geocenter(base_dir, DREL='', DSET='', CLOBBER=False, MODE=0o775,
     VERBOSE=False):
     """
-    Creates monthly files of geocenter variations at 6-hour intervals from
+    Creates monthly files of geocenter variations at 6-hour or 3-hour intervals from
     GRACE/GRACE-FO level-1b dealiasing data files
 
     Arguments
@@ -105,11 +103,32 @@ def aod1b_geocenter(base_dir, DREL='', DSET='', CLOBBER=False, MODE=0o775,
     #-- output formatting string
     fstr = '{0:4d}-{1:02d}-{2:02d}T{3:02d}:00:00 {4:12.8f} {5:12.8f} {6:12.8f}'
 
-    #-- Maximum spherical harmonic degree (LMAX)
-    LMAX = 100
+    #-- set number of hours in a file
+    #-- set the ocean model for a given release
+    if DREL in ('RL01','RL02','RL03','RL04','RL05'):
+        #-- for 00, 06, 12 and 18
+        n_time = 4
+        ATMOSPHERE = 'ECMWF'
+        OCEAN_MODEL = 'OMCT'
+        LMAX = 100
+    elif DREL in ('RL06',):
+        #-- for 00, 03, 06, 09, 12, 15, 18 and 21
+        n_time = 8
+        ATMOSPHERE = 'ECMWF'
+        OCEAN_MODEL = 'MPIOM'
+        LMAX = 180
+    else:
+        raise ValueError('Invalid data release')
     #-- Calculating the number of cos and sin harmonics up to LMAX
     n_harm = (LMAX**2 + 3*LMAX)//2 + 1
 
+    #-- AOD1B data products
+    product = {}
+    product['atm'] = 'Atmospheric loading from {0}'.format(ATMOSPHERE)
+    product['ocn'] = 'Oceanic loading from {0}'.format(OCEAN_MODEL)
+    product['glo'] = 'Global atmospheric and oceanic loading'
+    product['oba'] = 'Ocean bottom pressure from {0}'.format(OCEAN_MODEL)
+
     #-- AOD1B directory and output geocenter directory
     grace_dir = os.path.join(base_dir,'AOD1B',DREL)
     output_dir = os.path.join(grace_dir,'geocenter')
@@ -169,15 +188,15 @@ def aod1b_geocenter(base_dir, DREL='', DSET='', CLOBBER=False, MODE=0o775,
                 else:
                     fid = tar.extractfile(member)
                 #-- degree 1 spherical harmonics for day and hours
-                C10 = np.zeros((4))
-                C11 = np.zeros((4))
-                S11 = np.zeros((4))
-                hours = np.zeros((4),dtype=np.int64)
+                C10 = np.zeros((n_time))
+                C11 = np.zeros((n_time))
+                S11 = np.zeros((n_time))
+                hours = np.zeros((n_time),dtype=np.int64)
 
                 #-- create counter for hour in dataset
                 c = 0
                 #-- while loop ends when dataset is read
-                while (c < 4):
+                while (c < n_time):
                     #-- read line
                     file_contents = fid.readline().decode('ISO-8859-1')
                     #-- find file header for data product
@@ -220,9 +239,11 @@ def main():
     #-- Read the system arguments listed after the program
     parser = argparse.ArgumentParser(
         description="""Creates monthly files of geocenter variations
-            at 6-hour intervals
-            """
+            at 3 or 6-hour intervals
+            """,
+        fromfile_prefix_chars="@"
     )
+    parser.convert_arg_line_to_args = utilities.convert_arg_line_to_args
     #-- command line parameters
     #-- working data directory
     parser.add_argument('--directory','-D',

scripts/aod1b_oblateness.py

@@ -1,15 +1,16 @@
 #!/usr/bin/env python
 u"""
 aod1b_oblateness.py
-Written by Tyler Sutterley (05/2021)
+Written by Tyler Sutterley (07/2021)
+Contributions by Hugo Lecomte (03/2021)
 
 Reads GRACE/GRACE-FO level-1b dealiasing data files for a specific product
     atm: atmospheric loading from ECMWF
     ocn: oceanic loading from OMCT
     glo: global atmospheric and oceanic loading
     oba: ocean bottom pressure from OMCT
 
-Creates monthly files of oblateness (C20) variations at 6-hour intervals
+Creates monthly files of oblateness (C20) variations at 3 or 6 hour intervals
 
 NOTE: this reads the GFZ AOD1B files downloaded from PO.DAAC
 https://podaac-uat.jpl.nasa.gov/drive/files/allData/grace/L1B/GFZ/AOD1B/RL05/
@@ -31,7 +32,12 @@
 PYTHON DEPENDENCIES:
     numpy: Scientific Computing Tools For Python (https://numpy.org)
 
+PROGRAM DEPENDENCIES:
+    utilities.py: download and management utilities for files
+
 UPDATED HISTORY:
+    Updated 07/2021: can use default argument files to define options
+        Add 3-hour interval depending on Release
     Updated 05/2021: define int/float precision to prevent deprecation warning
     Updated 10/2020: use argparse to set command line parameters
     Updated 07/2020: added function docstrings
@@ -56,13 +62,7 @@
 import tarfile
 import argparse
 import numpy as np
-
-#-- aod1b data products
-product = {}
-product['atm'] = 'Atmospheric loading from ECMWF'
-product['ocn'] = 'Oceanic loading from OMCT'
-product['glo'] = 'Global atmospheric and oceanic loading'
-product['oba'] = 'Ocean bottom pressure from OMCT'
+import gravity_toolkit.utilities as utilities
 
 #-- program module to read the C20 coefficients of the AOD1b data
 def aod1b_oblateness(base_dir, DREL='', DSET='', CLOBBER=False, MODE=0o775,
@@ -103,11 +103,32 @@ def aod1b_oblateness(base_dir, DREL='', DSET='', CLOBBER=False, MODE=0o775,
     #-- output formatting string
     fstr = '{0:4d}-{1:02d}-{2:02d}T{3:02d}:00:00 {4:+16.8E}'
 
-    #-- Maximum spherical harmonic degree (LMAX)
-    LMAX = 100
+    #-- set number of hours in a file
+    #-- set the ocean model for a given release
+    if DREL in ('RL01','RL02','RL03','RL04','RL05'):
+        #-- for 00, 06, 12 and 18
+        n_time = 4
+        ATMOSPHERE = 'ECMWF'
+        OCEAN_MODEL = 'OMCT'
+        LMAX = 100
+    elif DREL in ('RL06',):
+        #-- for 00, 03, 06, 09, 12, 15, 18 and 21
+        n_time = 8
+        ATMOSPHERE = 'ECMWF'
+        OCEAN_MODEL = 'MPIOM'
+        LMAX = 180
+    else:
+        raise ValueError('Invalid data release')
     #-- Calculating the number of cos and sin harmonics up to LMAX
     n_harm = (LMAX**2 + 3*LMAX)//2 + 1
 
+    #-- AOD1B data products
+    product = {}
+    product['atm'] = 'Atmospheric loading from {0}'.format(ATMOSPHERE)
+    product['ocn'] = 'Oceanic loading from {0}'.format(OCEAN_MODEL)
+    product['glo'] = 'Global atmospheric and oceanic loading'
+    product['oba'] = 'Ocean bottom pressure from {0}'.format(OCEAN_MODEL)
+
     #-- AOD1B directory and output oblateness directory
     grace_dir = os.path.join(base_dir,'AOD1B',DREL)
     output_dir = os.path.join(grace_dir,'oblateness')
@@ -166,13 +187,13 @@ def aod1b_oblateness(base_dir, DREL='', DSET='', CLOBBER=False, MODE=0o775,
                 else:
                     fid = tar.extractfile(member)
                 #-- C20 spherical harmonics for day and hours
-                C20 = np.zeros((4))
-                hours = np.zeros((4),dtype=np.int64)
+                C20 = np.zeros((n_time))
+                hours = np.zeros((n_time),dtype=np.int64)
 
                 #-- create counter for hour in dataset
                 c = 0
                 #-- while loop ends when dataset is read
-                while (c < 4):
+                while (c < n_time):
                     #-- read line
                     file_contents = fid.readline().decode('ISO-8859-1')
                     #-- find file header for data product
@@ -210,9 +231,11 @@ def main():
     #-- Read the system arguments listed after the program
     parser = argparse.ArgumentParser(
         description="""Creates monthly files of oblateness (C20)
-            variations at 6-hour intervals
-            """
+            variations at 3 or 6-hour intervals
+            """,
+        fromfile_prefix_chars="@"
     )
+    parser.convert_arg_line_to_args = utilities.convert_arg_line_to_args
     #-- command line parameters
     #-- working data directory
     parser.add_argument('--directory','-D',