refactor: add descriptive file-level attributes to output netCDF4/HDF…

…5 files (#115)

feat: include option to not compensate for elastic deformation for #37
feat: added regex formatting for CNES GRGS harmonics

gravity_toolkit/grace_date.py

@@ -47,6 +47,7 @@
 
 UPDATE HISTORY:
     Updated 03/2023: use f-strings for formatting output date lines
+        added regex formatting for CNES GRGS harmonics
     Updated 11/2022: use f-strings for formatting verbose output
     Updated 09/2022: raise exception if index file cannot be found
         use logging for debugging level verbose output
@@ -131,6 +132,7 @@ def grace_date(base_dir, PROC='', DREL='', DSET='', OUTPUT=True, MODE=0o775):
             - ``'GRAZ'``: Institute of Geodesy from GRAZ University of Technology
             - ``'COSTG'``: Combination Service for Time-variable Gravity Fields
             - ``'Swarm'``: Time-variable gravity data from Swarm satellites
+            - ``'GRGS'``: CNES Groupe de Recherche de Geodesie Spatiale
     DREL: str, default ''
         GRACE/GRACE-FO/Swarm data release
     DSET: str, default ''

gravity_toolkit/read_GRACE_harmonics.py

@@ -42,6 +42,7 @@
     time.py: utilities for calculating time operations
 
 UPDATE HISTORY:
+    Updated 03/2023: added regex formatting for CNES GRGS harmonics
     Updated 11/2022: use f-strings for formatting verbose or ascii output
     Updated 10/2022: make keyword arguments part of kwargs dictionary
     Updated 05/2022: updated comments
@@ -192,6 +193,8 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs):
         # clm and slm drift rates for RL04
         drift_c = np.zeros((LMAX+1,MMAX+1))
         drift_s = np.zeros((LMAX+1,MMAX+1))
+    # set default degree 0 harmonics for intercomparability between centers
+    grace_L2_input['clm'][0, 0] = 1.0
 
     # extract GRACE and GRACE-FO file headers
     # replace colons in header if within quotations
@@ -305,9 +308,11 @@ def parse_file(input_file):
     # JPLMSC: NASA Jet Propulsion Laboratory (mascon solutions)
     # GRGS: French Centre National D'Etudes Spatiales (CNES)
     # COSTG: International Combined Time-variable Gravity Fields
-    args = r'UTCSR|EIGEN|GFZOP|JPLEM|JPLMSC|GRGS|COSTG'
-    regex_pattern = (r'(.*?)-2_(\d{{4}})(\d{{3}})-(\d{{4}})(\d{{3}})_'
-        r'(.*?)_({0})_(.*?)_(\d+)(.*?)(\.gz|\.gfc)?$').format(args)
+    # GRGS: CNES Groupe de Recherche de Geodesie Spatiale
+    centers = r'UTCSR|EIGEN|GFZOP|JPLEM|JPLMSC|GRGS|COSTG|GRGS'
+    suffixes = r'\.gz|\.gfc|\.txt'
+    regex_pattern = (r'(.*?)-2_(\d{4})(\d{3})-(\d{4})(\d{3})_'
+        rf'(.*?)_({centers})_(.*?)_(\d+)(.*?)({suffixes})?$')
     rx = re.compile(regex_pattern, re.VERBOSE)
     # extract parameters from input filename
     if isinstance(input_file, io.IOBase):

gravity_toolkit/time.py

@@ -1,7 +1,9 @@
 #!/usr/bin/env python
 u"""
 time.py
-Written by Tyler Sutterley (11/2022)
+Written by Tyler Sutterley (03/2023)
+Contributions by Hugo Lecomte
+
 Utilities for calculating time operations
 
 PYTHON DEPENDENCIES:
@@ -11,6 +13,7 @@
         https://dateutil.readthedocs.io/en/stable/
 
 UPDATE HISTORY:
+    Updated 03/2023: added regex formatting for CNES GRGS harmonics
     Updated 11/2022: use f-strings for formatting verbose or ascii output
     Updated 10/2022: added more time parsing for longer periods
     Updated 08/2022: added file parsing functions from GRACE date utilities
@@ -147,9 +150,11 @@ def parse_grace_file(granule):
     # JPLMSC: NASA Jet Propulsion Laboratory (mascon solutions)
     # GRGS: French Centre National D'Etudes Spatiales (CNES)
     # COSTG: International Combined Time-variable Gravity Fields
-    args = r'UTCSR|EIGEN|GFZOP|JPLEM|JPLMSC|GRGS|COSTG'
-    regex_pattern = (r'(.*?)-2_(\d{{4}})(\d{{3}})-(\d{{4}})(\d{{3}})_'
-        r'(.*?)_({0})_(.*?)_(\d+)(.*?)(\.gz|\.gfc)?$').format(args)
+    # GRGS: CNES Groupe de Recherche de Geodesie Spatiale
+    centers = r'UTCSR|EIGEN|GFZOP|JPLEM|JPLMSC|GRGS|COSTG|GRGS'
+    suffixes = r'\.gz|\.gfc|\.txt'
+    regex_pattern = (r'(.*?)-2_(\d{4})(\d{3})-(\d{4})(\d{3})_'
+        rf'(.*?)_({centers})_(.*?)_(\d+)(.*?)({suffixes})?$')
     rx = re.compile(regex_pattern, re.VERBOSE)
     # extract parameters from input filename
     PFX,SY,SD,EY,ED,AUX,PRC,F1,DRL,F2,SFX = rx.findall(file_basename).pop()

gravity_toolkit/units.py

@@ -1,14 +1,17 @@
 #!/usr/bin/env python
 u"""
 units.py
-Written by Tyler Sutterley (02/2023)
+Written by Tyler Sutterley (03/2023)
+Contributions by Hugo Lecomte
 
 Class for converting GRACE/GRACE-FO Level-2 data to specific units
 
 PYTHON DEPENDENCIES:
     numpy: Scientific Computing Tools For Python (https://numpy.org)
 
 UPDATE HISTORY:
+    Updated 03/2023: include option to not compensate for elastic deformation
+        include option to include effects for Earth's oblateness
     Updated 02/2023: fixed case where maximum spherical harmonic degree is 0
     Updated 01/2023: added function to retrieve named units
     Updated 12/2022: set average Earth's density and radius as class properties
@@ -93,7 +96,7 @@ def rho_e(self):
         """
         return 0.75*self.GM/(self.G*np.pi*self.rad_e**3)
 
-    def harmonic(self, hl, kl, ll):
+    def harmonic(self, hl, kl, ll, **kwargs):
         """
         Calculates degree dependent factors for converting spherical
         harmonic units from [Wahr1998]_
@@ -106,6 +109,10 @@ def harmonic(self, hl, kl, ll):
             Gravitational Potential load Love numbers
         ll: float
             Horizontal Displacement load Love numbers
+        include_elastic: bool, default True
+            Include compensation for elastic response
+        include_ellipsoidal: bool, default False
+            Include consideration for Earth's oblateness
 
         Attributes
         ----------
@@ -132,58 +139,89 @@ def harmonic(self, hl, kl, ll):
         Pa: float
             pascals equivalent surface pressure
         """
+        # set default keyword arguments
+        kwargs.setdefault('include_elastic', True)
+        kwargs.setdefault('include_ellipsoidal', False)
+        fraction = np.ones((self.lmax+1))
+        # compensate for elastic deformation within the solid earth
+        if kwargs['include_elastic']:
+            fraction += kl[self.l]
+        # include effects for Earth's oblateness
+        if kwargs['include_ellipsoidal']:
+            fraction /= (1.0 - self.flat)
         # degree dependent coefficients
         # norm, fully normalized spherical harmonics
         self.norm = np.ones((self.lmax+1))
         # cmwe, centimeters water equivalent [g/cm^2]
-        self.cmwe = self.rho_e*self.rad_e*(2.0*self.l+1.0)/(1.0+kl[self.l])/3.0
+        self.cmwe = self.rho_e*self.rad_e*(2.0*self.l+1.0)/fraction/3.0
         # mmwe, millimeters water equivalent [kg/m^2]
-        self.mmwe = 10.0*self.rho_e*self.rad_e*(2.0*self.l+1.0)/(1.0+kl[self.l])/3.0
+        self.mmwe = 10.0*self.rho_e*self.rad_e*(2.0*self.l+1.0)/fraction/3.0
         # mmGH, millimeters geoid height
         self.mmGH = np.ones((self.lmax+1))*(10.0*self.rad_e)
         # mmCU, millimeters elastic crustal deformation (uplift)
-        self.mmCU = 10.0*self.rad_e*hl[self.l]/(1.0+kl[self.l])
+        self.mmCU = 10.0*self.rad_e*hl[self.l]/fraction
         # mmCH, millimeters elastic crustal deformation (horizontal)
-        self.mmCH = 10.0*self.rad_e*ll[self.l]/(1.0+kl[self.l])
+        self.mmCH = 10.0*self.rad_e*ll[self.l]/fraction
         # cmVCU, centimeters viscoelastic crustal uplift
         self.cmVCU = self.rad_e*(2.0*self.l+1.0)/2.0
         # mVCU, meters viscoelastic crustal uplift
         self.mVCU = self.rad_e*(2.0*self.l+1.0)/200.0
         # microGal, microGal gravity perturbations
         self.microGal = 1.e6*self.GM*(self.l+1.0)/(self.rad_e**2.0)
         # mbar, millibar equivalent surface pressure
-        self.mbar = self.g_wmo*self.rho_e*self.rad_e*(2.0*self.l+1.0)/(1.0+kl[self.l])/3e3
+        self.mbar = self.g_wmo*self.rho_e*self.rad_e*(2.0*self.l+1.0)/fraction/3e3
         # Pa, pascals equivalent surface pressure
-        self.Pa = self.g_wmo*self.rho_e*self.rad_e*(2.0*self.l+1.0)/(1.0+kl[self.l])/30.0
+        self.Pa = self.g_wmo*self.rho_e*self.rad_e*(2.0*self.l+1.0)/fraction/30.0
         # return the degree dependent unit conversions
         return self
 
-    def spatial(self, hl, kl, ll):
+    def spatial(self, hl, kl, ll, **kwargs):
         """
         Calculates degree dependent factors for converting spatial units
         from [Wahr1998]_
 
         Parameters
         ----------
+
         hl: float
             Vertical Displacement load Love numbers
         kl: float
             Gravitational Potential load Love numbers
         ll: float
             Horizontal Displacement load Love numbers
+        include_elastic: bool, default True
+            Include compensation for elastic response
 
         Attributes
         ----------
+        norm: float
+            fully normalized spherical harmonics
         cmwe: float
             centimeters water equivalent
         mmwe: float
             millimeters water equivalent
+        mmGH: float
+            millimeters geoid height
+        microGal: float
+            microGal gravity perturbations
         """
+        # set default keyword arguments
+        kwargs.setdefault('include_elastic', True)
+        fraction = np.ones((self.lmax+1))
+        # compensate for elastic deformation within the solid earth
+        if kwargs['include_elastic']:
+            fraction += kl[self.l]
         # degree dependent coefficients
+        # norm, fully normalized spherical harmonics
+        self.norm = np.ones((self.lmax+1))
         # cmwe, centimeters water equivalent [g/cm^2]
-        self.cmwe = 3.0*(1.0+kl[self.l])/(1.0+2.0*self.l)/(4.0*np.pi*self.rad_e*self.rho_e)
+        self.cmwe = 3.0*fraction/(1.0+2.0*self.l)/(4.0*np.pi*self.rad_e*self.rho_e)
         # mmwe, millimeters water equivalent [kg/m^2]
-        self.mmwe = 3.0*(1.0+kl[self.l])/(1.0+2.0*self.l)/(40.0*np.pi*self.rad_e*self.rho_e)
+        self.mmwe = 3.0*fraction/(1.0+2.0*self.l)/(40.0*np.pi*self.rad_e*self.rho_e)
+        # mmGH, millimeters geoid height
+        self.mmGH = np.ones((self.lmax+1))/(4.0*np.pi*self.rad_e)
+        # microGal, microGal gravity perturbations
+        self.microGal = (self.rad_e**2.0)/(4.0*np.pi*1.e6*self.GM)/(self.l+1.0)
         # return the degree dependent unit conversions
         return self
 

scripts/combine_harmonics.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 combine_harmonics.py
-Written by Tyler Sutterley (02/2023)
+Written by Tyler Sutterley (03/2023)
 Converts a file from the spherical harmonic domain into the spatial domain
 
 CALLING SEQUENCE:
@@ -71,6 +71,8 @@
     utilities.py: download and management utilities for files
 
 UPDATE HISTORY:
+    Updated 03/2023: add index ascii/netCDF4/HDF5 datatypes as possible inputs
+        add descriptive file-level attributes to output netCDF4/HDF5 files
     Updated 02/2023: use get function to retrieve specific units
         use love numbers class with additional attributes
     Updated 01/2023: refactored associated legendre polynomials
@@ -143,13 +145,27 @@ def combine_harmonics(INPUT_FILE, OUTPUT_FILE,
     DIRECTORY = os.path.abspath(os.path.dirname(OUTPUT_FILE))
     if not os.access(DIRECTORY, os.F_OK):
         os.makedirs(DIRECTORY,MODE,exist_ok=True)
+    # attributes for output files
+    attributes = dict(ROOT={})
+    attributes['ROOT']['product_type'] = 'gravity_field'
+    attributes['ROOT']['reference'] = f'Output from {os.path.basename(sys.argv[0])}'
 
     # upper bound of spherical harmonic orders (default = LMAX)
     if MMAX is None:
         MMAX = np.copy(LMAX)
 
     # read input spherical harmonic coefficients from file
-    input_Ylms = gravtk.harmonics().from_file(INPUT_FILE, format=DATAFORM)
+    if DATAFORM in ('ascii', 'netCDF4', 'HDF5'):
+        dataform = copy.copy(DATAFORM)
+        attributes['ROOT']['lineage'] = os.path.basename(INPUT_FILE)
+        input_Ylms = gravtk.harmonics().from_file(INPUT_FILE,
+            format=DATAFORM)
+    elif DATAFORM in ('index-ascii', 'index-netCDF4', 'index-HDF5'):
+        # read from index file
+        _,dataform = DATAFORM.split('-')
+        attributes['ROOT']['lineage'] = [os.path.basename(f) for f in INPUT_FILE]
+        input_Ylms = gravtk.harmonics().from_index(INPUT_FILE,
+            format=dataform)
     # reform harmonic dimensions to be l,m,t
     # truncate to degree and order LMAX, MMAX
     input_Ylms = input_Ylms.truncate(lmax=LMAX, mmax=MMAX).expand_dims()
@@ -163,6 +179,13 @@ def combine_harmonics(INPUT_FILE, OUTPUT_FILE,
     # read arrays of kl, hl, and ll Love Numbers
     LOVE = gravtk.load_love_numbers(LMAX, LOVE_NUMBERS=LOVE_NUMBERS,
         REFERENCE=REFERENCE, FORMAT='class')
+    # add attributes for earth parameters
+    attributes['ROOT']['earth_model'] = LOVE.model
+    attributes['ROOT']['earth_love_numbers'] = LOVE.citation
+    attributes['ROOT']['reference_frame'] = LOVE.reference
+    # add attributes for maximum degree and order
+    attributes['ROOT']['max_degree'] = LMAX
+    attributes['ROOT']['max_order'] = MMAX
 
     # distribute total mass uniformly over the ocean
     if REDISTRIBUTE:
@@ -225,23 +248,21 @@ def combine_harmonics(INPUT_FILE, OUTPUT_FILE,
     # Setting units factor for output
     # dfactor computes the degree dependent coefficients
     factors = gravtk.units(lmax=LMAX).harmonic(*LOVE)
-    if (UNITS == 1):
-        # 1: cmwe, centimeters water equivalent
-        dfactor = factors.get('cmwe')
-    elif (UNITS == 2):
-        # 2: mmGH, millimeters geoid height
-        dfactor = factors.get('mmGH')
-    elif (UNITS == 3):
-        # 3: mmCU, millimeters elastic crustal deformation
-        dfactor = factors.get('mmCU')
-    elif (UNITS == 4):
-        # 4: micGal, microGal gravity perturbations
-        dfactor = factors.get('microGal')
-    elif (UNITS == 5):
-        # 5: mbar, millibars equivalent surface pressure
-        dfactor = factors.get('mbar')
-    else:
-        raise ValueError(f'Invalid units code {UNITS:d}')
+    # output units and units longname
+    # 1: cmwe, centimeters water equivalent
+    # 2: mmGH, millimeters geoid height
+    # 3: mmCU, millimeters elastic crustal deformation
+    # 4: micGal, microGal gravity perturbations
+    # 5: mbar, millibars equivalent surface pressure
+    unit_short = ['cmwe', 'mmGH', 'mmCU', 'microGal', 'mbar']
+    unit_name = ['Equivalent_Water_Thickness', 'Geoid_Height',
+        'Elastic_Crustal_Uplift', 'Gravitational_Undulation',
+        'Equivalent_Surface_Pressure']
+    dfactor = factors.get(unit_short[UNITS-1])
+    # add attributes for earth parameters
+    attributes['ROOT']['earth_radius'] = f'{factors.rad_e:0.3f} cm'
+    attributes['ROOT']['earth_density'] = f'{factors.rho_e:0.3f} g/cm'
+    attributes['ROOT']['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2'
 
     # Computing plms for converting to spatial domain
     theta = (90.0 - grid.lat)*np.pi/180.0
@@ -256,33 +277,12 @@ def combine_harmonics(INPUT_FILE, OUTPUT_FILE,
             grid.lon, grid.lat, LMAX=LMAX, PLM=PLM).T
 
     # outputting data to file
-    output_data(grid.squeeze(), FILENAME=OUTPUT_FILE,
-        DATAFORM=DATAFORM, UNITS=UNITS, MODE=MODE)
+    grid.squeeze().to_file(filename=OUTPUT_FILE, format=dataform,
+        units=unit_short[UNITS-1], longname=unit_name[UNITS-1],
+        attributes=attributes)
 
-# PURPOSE: wrapper function for outputting data to file
-def output_data(data, FILENAME=None, DATAFORM=None, UNITS=None, MODE=None):
-    # output units and units longname
-    unit_short = ['cmwe', 'mmGH', 'mmCU', 'microGal', 'mbar']
-    unit_name = ['Equivalent_Water_Thickness', 'Geoid_Height',
-        'Elastic_Crustal_Uplift', 'Gravitational_Undulation',
-        'Equivalent_Surface_Pressure']
-    # attributes for output files
-    attributes = {}
-    attributes['units'] = copy.copy(unit_short[UNITS-1])
-    attributes['longname'] = copy.copy(unit_name[UNITS-1])
-    attributes['reference'] = f'Output from {os.path.basename(sys.argv[0])}'
-    # output to file
-    if (DATAFORM == 'ascii'):
-        # ascii (.txt)
-        data.to_ascii(FILENAME)
-    elif (DATAFORM == 'netCDF4'):
-        # netcdf (.nc)
-        data.to_netCDF4(FILENAME, **attributes)
-    elif (DATAFORM == 'HDF5'):
-        # HDF5 (.H5)
-        data.to_HDF5(FILENAME, **attributes)
     # change output permissions level to MODE
-    os.chmod(FILENAME, MODE)
+    os.chmod(OUTPUT_FILE, MODE)
 
 # PURPOSE: create argument parser
 def arguments():
@@ -359,8 +359,11 @@ def arguments():
         type=lambda p: os.path.abspath(os.path.expanduser(p)),
         help='Mean file to remove from the harmonic data')
     # input and output data format (ascii, netCDF4, HDF5)
+    choices = []
+    choices.extend(['ascii','netCDF4','HDF5'])
+    choices.extend(['index-ascii','index-netCDF4','index-HDF5'])
     parser.add_argument('--format','-F',
-        type=str, default='netCDF4', choices=['ascii','netCDF4','HDF5'],
+        type=str, default='netCDF4', choices=choices,
         help='Input and output data format')
     # print information about each input and output file
     parser.add_argument('--verbose','-V',