diff --git a/README.rst b/README.rst index 3ebcded3..5b4185a8 100644 --- a/README.rst +++ b/README.rst @@ -56,11 +56,11 @@ References T. C. Sutterley, I. Velicogna, and C.-W. Hsu, "Self-Consistent Ice Mass Balance and Regional Sea Level From Time-Variable Gravity", *Earth and Space Science*, 7, - (2020). `doi:10.1029/2019EA000860 `_ + (2020). `doi: 10.1029/2019EA000860 `_ T. C. Sutterley and I. Velicogna, "Improved estimates of geocenter variability from time-variable gravity and ocean model outputs", *Remote Sensing*, 11(18), - 2108, (2019). `doi:10.3390/rs11182108 `_ + 2108, (2019). `doi: 10.3390/rs11182108 `_ J. Wahr, S. C. Swenson, and I. Velicogna, "Accuracy of GRACE mass estimates", *Geophysical Research Letters*, 33(6), L06401, (2006). @@ -81,11 +81,11 @@ Data Repositories T. C. Sutterley, I. Velicogna, and C.-W. Hsu, "Ice Mass and Regional Sea Level Estimates from Time-Variable Gravity", (2020). - `doi:10.6084/m9.figshare.9702338 `_ + `doi: 10.6084/m9.figshare.9702338 `_ T. C. Sutterley and I. Velicogna, "Geocenter Estimates from Time-Variable Gravity and Ocean Model Outputs", (2019). - `doi:10.6084/m9.figshare.7388540 `_ + `doi: 10.6084/m9.figshare.7388540 `_ Download ######## diff --git a/doc/source/api_reference/utilities.rst b/doc/source/api_reference/utilities.rst index 1b5158ee..e1f41ee4 100644 --- a/doc/source/api_reference/utilities.rst +++ b/doc/source/api_reference/utilities.rst @@ -19,6 +19,9 @@ General Methods .. autofunction:: gravity_toolkit.utilities.get_data_path +.. autoclass:: gravity_toolkit.utilities.reify + :members: + .. autofunction:: gravity_toolkit.utilities.get_hash .. autofunction:: gravity_toolkit.utilities.get_git_revision_hash diff --git a/doc/source/getting_started/Citations.rst b/doc/source/getting_started/Citations.rst index 4bcfd22b..3637aef1 100644 --- a/doc/source/getting_started/Citations.rst +++ b/doc/source/getting_started/Citations.rst @@ -12,15 +12,15 @@ most recently to the following work: E. Rignot, T. C. Sutterley, M. van den Broeke, J. M. van Wessem, and D. Wiese, "Continuity of ice sheet mass loss in Greenland and Antarctica from the GRACE and GRACE Follow-On missions", *Geophysical Research Letters*, 47, - (2020). `doi:10.1029/2020GL087291 `_ + (2020). `doi: 10.1029/2020GL087291 `_ T. C. Sutterley, I. Velicogna, and C.-W. Hsu, "Self-Consistent Ice Mass Balance and Regional Sea Level From Time-Variable Gravity", *Earth and Space Science*, 7, - (2020). `doi:10.1029/2019EA000860 `_ + (2020). `doi: 10.1029/2019EA000860 `_ T. C. Sutterley and I. Velicogna, "Improved estimates of geocenter variability from time-variable gravity and ocean model outputs", *Remote Sensing*, 11(18), - 2108, (2019). `doi:10.3390/rs11182108 `_ + 2108, (2019). `doi: 10.3390/rs11182108 `_ Some of the text within this documentation are modifications from these original works, which is allowed under the diff --git a/gravity_toolkit/grace_find_months.py b/gravity_toolkit/grace_find_months.py index 17dbfed7..dd0ac642 100644 --- a/gravity_toolkit/grace_find_months.py +++ b/gravity_toolkit/grace_find_months.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" grace_find_months.py -Written by Tyler Sutterley (11/2022) +Written by Tyler Sutterley (05/2023) Parses date index file from grace_date program Finds the months available for a GRACE/GRACE-FO/Swarm product @@ -36,6 +36,7 @@ grace_date.py: reads GRACE index file and calculates dates for each month UPDATE HISTORY: + Updated 05/2023: use formatting for reading from date file Updated 11/2022: use f-strings for formatting verbose or ascii output Updated 04/2022: updated docstrings to numpy documentation format Updated 05/2021: define int/float precision to prevent deprecation warning @@ -109,22 +110,30 @@ def grace_find_months(base_dir, PROC, DREL, DSET='GSM'): if not os.access(date_file, os.F_OK): grace_date(base_dir, PROC=PROC, DREL=DREL, DSET=DSET, OUTPUT=True) - # read GRACE/GRACE-FO date ascii file from grace_date.py + # names and formats of GRACE/GRACE-FO date ascii file + names = ('t','mon','styr','stday','endyr','endday','total') + formats = ('f','i','i','i','i','i','i') + dtype = np.dtype({'names':names, 'formats':formats}) + # read GRACE/GRACE-FO date ascii file # skip the header row and extract dates (decimal format) and months - date_input = np.loadtxt(date_file, skiprows=1) - tdec = date_input[:,0] - months = date_input[:,1].astype(np.int64) + date_input = np.loadtxt(date_file, skiprows=1, dtype=dtype) + # date info dictionary + var_info = {} + var_info['time'] = date_input['t'] + var_info['months'] = date_input['mon'] + var_info['start'] = np.min(date_input['mon']) + var_info['end'] = np.max(date_input['mon']) # array of all possible months (or in case of CNES RL01/2: 10-day sets) - all_months = np.arange(1,months.max(),dtype=np.int64) + all_months = np.arange(1, var_info['end'], dtype=np.int64) # missing months (values in all_months but not in months) - missing = sorted(set(all_months)-set(months)) + var_info['missing'] = sorted(set(all_months) - set(date_input['mon'])) # If CNES RL01/2: simply convert into numpy array # else: remove months 1-3 and convert into numpy array if ((PROC == 'CNES') & (DREL in ('RL01','RL02'))): - missing = np.array(missing,dtype=np.int64) + var_info['missing'] = np.array(var_info['missing'], dtype=np.int64) else: - missing = np.array(missing[3:],dtype=np.int64) + var_info['missing'] = np.array(var_info['missing'][3:], dtype=np.int64) - return {'time':tdec, 'start':months[0], 'end':months[-1], 'months':months, - 'missing':missing} + # return the date information dictionary + return var_info diff --git a/gravity_toolkit/grace_months_index.py b/gravity_toolkit/grace_months_index.py index 209a06ad..c9888e6a 100644 --- a/gravity_toolkit/grace_months_index.py +++ b/gravity_toolkit/grace_months_index.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" grace_months_index.py -Written by Tyler Sutterley (11/2022) +Written by Tyler Sutterley (05/2023) Creates a file with the start and end days for each dataset Shows the range of each month for (CSR/GFZ/JPL) (RL04/RL05/RL06) @@ -40,6 +40,7 @@ time.py: utilities for calculating time operations UPDATE HISTORY: + Updated 05/2023: use formatting for reading from date file Updated 11/2022: use f-strings for formatting verbose or ascii output Updated 05/2022: use argparse descriptions within documentation use new GSFC release 6 version 2 mascons as the default @@ -108,34 +109,19 @@ def grace_months_index(base_dir, DREL=['RL06','rl06v2.0'], MODE=None): for pr in PROC: # Setting the data directory for processing center and release grace_dir = os.path.join(base_dir, pr, rl, DSET) - # read GRACE date ascii file - # file created in read_grace.py or grace_dates.py - grace_date_file = f'{pr}_{rl}_DATES.txt' - if os.access(os.path.join(grace_dir,grace_date_file), os.F_OK): - # skip the header line - date_input = np.loadtxt(os.path.join(grace_dir,grace_date_file), - skiprows=1) - # number of months - nmon = np.shape(date_input)[0] - + # read GRACE/GRACE-FO date ascii file + grace_date_file = os.path.join(grace_dir,f'{pr}_{rl}_DATES.txt') + # names and formats of GRACE/GRACE-FO date ascii file + names = ('t','mon','styr','stday','endyr','endday','total') + formats = ('f','i','i','i','i','i','i') + dtype = np.dtype({'names':names, 'formats':formats}) + # check that the GRACE/GRACE-FO date file exists + if os.access(grace_date_file, os.F_OK): # Setting the dictionary key e.g. 'CSR_RL04' var_name = f'{pr}_{rl}' - - # Creating a python dictionary for each dataset with parameters: - # month #, start year, start day, end year, end day - # Purpose is to get all of the dates loaded for each dataset - # Adding data to dictionary for data processing and release - var_info[var_name] = {} - # allocate for output variables - var_info[var_name]['mon'] = np.zeros((nmon),dtype=np.int64) - var_info[var_name]['styr'] = np.zeros((nmon),dtype=np.int64) - var_info[var_name]['stday'] = np.zeros((nmon),dtype=np.int64) - var_info[var_name]['endyr'] = np.zeros((nmon),dtype=np.int64) - var_info[var_name]['endday'] = np.zeros((nmon),dtype=np.int64) - # place output variables in dictionary - for i,key in enumerate(['mon','styr','stday','endyr','endday']): - # first column is date in decimal form (start at 1 not 0) - var_info[var_name][key] = date_input[:,i+1].astype(np.int64) + # skip the header line + var_info[var_name] = np.loadtxt(grace_date_file, + skiprows=1, dtype=dtype) # Finding the maximum month measured if (var_info[var_name]['mon'].max() > max_mon): # if the maximum month in this dataset is greater diff --git a/gravity_toolkit/harmonics.py b/gravity_toolkit/harmonics.py index b8a80f17..90f1f10b 100644 --- a/gravity_toolkit/harmonics.py +++ b/gravity_toolkit/harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" harmonics.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (05/2023) Contributions by Hugo Lecomte Spherical harmonic data class for processing GRACE/GRACE-FO Level-2 data @@ -25,6 +25,7 @@ destripe_harmonics.py: filters spherical harmonics for correlated errors UPDATE HISTORY: + Updated 05/2023: use reify decorators for complex form and amplitude Updated 03/2023: customizable file-level attributes to netCDF4 and HDF5 add attributes fetching to the from_dict and to_dict functions retrieve all root attributes from HDF5 and netCDF4 datasets @@ -101,6 +102,7 @@ from gravity_toolkit.destripe_harmonics import destripe_harmonics from gravity_toolkit.read_gfc_harmonics import read_gfc_harmonics from gravity_toolkit.read_GRACE_harmonics import read_GRACE_harmonics +from gravity_toolkit.utilities import reify # attempt imports try: @@ -1806,7 +1808,7 @@ def destripe(self, **kwargs): # return the destriped field return temp - @property + @reify def amplitude(self): """ Degree amplitude of the spherical harmonics @@ -1852,7 +1854,7 @@ def ndim(self): """ return np.ndim(self.clm) - @property + @reify def ilm(self): """ Complex form of the spherical harmonics diff --git a/gravity_toolkit/utilities.py b/gravity_toolkit/utilities.py index 58c5bb86..c5262c51 100644 --- a/gravity_toolkit/utilities.py +++ b/gravity_toolkit/utilities.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" utilities.py -Written by Tyler Sutterley (04/2023) +Written by Tyler Sutterley (05/2023) Download and management utilities for syncing time and auxiliary files PYTHON DEPENDENCIES: @@ -9,6 +9,7 @@ https://pypi.python.org/pypi/lxml UPDATE HISTORY: + Updated 05/2023: add reify decorator for evaluation of properties Updated 04/2023: use release-03 GFZ GravIS SLR and geocenter files Updated 03/2023: place boto3 import within try/except statement Updated 01/2023: add default ssl context attribute with protocol @@ -104,6 +105,21 @@ def get_data_path(relpath): elif isinstance(relpath,str): return os.path.join(filepath,relpath) +class reify(object): + """Class decorator that puts the result of the method it + decorates into the instance""" + def __init__(self, wrapped): + self.wrapped = wrapped + self.__name__ = wrapped.__name__ + self.__doc__ = wrapped.__doc__ + + def __get__(self, inst, objtype=None): + if inst is None: + return self + val = self.wrapped(inst) + setattr(inst, self.wrapped.__name__, val) + return val + # PURPOSE: get the hash value of a file def get_hash(local, algorithm='MD5'): """ diff --git a/requirements.txt b/requirements.txt index 68b83306..86591e4e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,7 @@ +boto3 future lxml +matplotlib numpy python-dateutil pyyaml diff --git a/scripts/dealiasing_monthly_mean.py b/scripts/dealiasing_monthly_mean.py index ac5c7751..35281870 100755 --- a/scripts/dealiasing_monthly_mean.py +++ b/scripts/dealiasing_monthly_mean.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" dealiasing_monthly_mean.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (05/2023) Reads GRACE/GRACE-FO AOD1B datafiles for a specific product and outputs monthly the mean for a specific GRACE/GRACE-FO processing center and data release @@ -48,6 +48,7 @@ time.py: utilities for calculating time operations UPDATE HISTORY: + Updated 05/2023: use formatting for reading from date file Updated 03/2023: read data into flattened harmonics objects debug-level logging of member names and header lines convert shape and ndim to harmonic class properties @@ -103,7 +104,7 @@ def dealiasing_monthly_mean(base_dir, PROC=None, DREL=None, DSET=None, aod1b_products = dict(GAA='atm',GAB='ocn',GAC='glo',GAD='oba') # compile regular expressions operator for the clm/slm headers # for the specific AOD1b product - hx = re.compile(r'^DATA.*SET.*{0}'.format(aod1b_products[DSET]),re.VERBOSE) + hx = re.compile(fr'^DATA.*SET.*{aod1b_products[DSET]}',re.VERBOSE) # compile regular expression operator to find numerical instances # will extract the data from the file regex_pattern = r'[-+]?(?:(?:\d*\.\d+)|(?:\d+\.?))(?:[Ee][+-]?\d+)?' @@ -163,16 +164,20 @@ def dealiasing_monthly_mean(base_dir, PROC=None, DREL=None, DSET=None, CENTER = default_center # read input DATE file from GSM data product - grace_datefile = '{0}_{1}_DATES.txt'.format(PROC, DREL) - date_input = np.loadtxt(os.path.join(grace_dir,'GSM',grace_datefile), - skiprows=1) - grace_month = date_input[:,1].astype(np.int64) - start_yr = date_input[:,2] - start_day = date_input[:,3].astype(np.int64) - end_yr = date_input[:,4] - end_day = date_input[:,5].astype(np.int64) + grace_date_file = f'{PROC}_{DREL}_DATES.txt' + # names and formats of GRACE/GRACE-FO date ascii file + names = ('t','mon','styr','stday','endyr','endday','total') + formats = ('f','i','i','i','i','i','i') + dtype = np.dtype({'names':names, 'formats':formats}) + date_input = np.loadtxt(os.path.join(grace_dir,'GSM',grace_date_file), + skiprows=1, dtype=dtype) + grace_month = date_input['mon'] + start_yr = date_input['styr'] + start_day = date_input['stday'] + end_yr = date_input['endyr'] + end_day = date_input['endday'] # output date file reduced to months with complete AOD - f_out = open(os.path.join(grace_dir,DSET,grace_datefile), + f_out = open(os.path.join(grace_dir,DSET,grace_date_file), mode='w', encoding='utf8') # date file header information args = ('Mid-date','Month','Start_Day','End_Day','Total_Days') diff --git a/scripts/regress_grace_maps.py b/scripts/regress_grace_maps.py index 51acc9a3..e2e27f1e 100755 --- a/scripts/regress_grace_maps.py +++ b/scripts/regress_grace_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" regress_grace_maps.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (05/2023) Reads in GRACE/GRACE-FO spatial files and fits a regression model at each grid point @@ -60,6 +60,7 @@ utilities.py: download and management utilities for files UPDATE HISTORY: + Updated 05/2023: split S2 tidal aliasing terms into GRACE and GRACE-FO eras Updated 03/2023: updated inputs to spatial from_ascii function use attributes from units class for writing to netCDF4/HDF5 files Updated 01/2023: refactored time series analysis functions @@ -176,38 +177,6 @@ def regress_grace_maps(LMAX, RAD, nlon = len(lon) nlat = len(lat) - # Setting output parameters for each fit type - coef_str = ['x{0:d}'.format(o) for o in range(ORDER+1)] - unit_suffix = [' yr^{0:d}'.format(-o) if o else '' for o in range(ORDER+1)] - if (ORDER == 0):# Mean - unit_longname = ['Mean'] - elif (ORDER == 1):# Trend - unit_longname = ['Constant','Trend'] - elif (ORDER == 2):# Quadratic - unit_longname = ['Constant','Linear','Quadratic'] - # filename strings for cyclical terms - cyclic_str = {} - cyclic_str['SEMI'] = ['SS','SC'] - cyclic_str['ANN'] = ['AS','AC'] - cyclic_str['S2'] = ['S2S','S2C'] - # unit longnames for cyclical terms - cyclic_longname = {} - cyclic_longname['SEMI'] = ['Semi-Annual Sine', 'Semi-Annual Cosine'] - cyclic_longname['ANN'] = ['Annual Sine', 'Annual Cosine'] - cyclic_longname['S2'] = ['S2 Tidal Alias Sine', 'S2 Tidal Alias Cosine'] - amp_str = [] - for i,c in enumerate(CYCLES): - if (c == 0.5): - flag = 'SEMI' - elif (c == 1.0): - flag = 'ANN' - elif (c == (161.0/365.25)): - flag = 'S2' - coef_str.extend(cyclic_str[flag]) - unit_longname.extend(cyclic_longname[flag]) - unit_suffix.extend(['','']) - amp_str.append(flag) - # input data spatial object spatial_list = [] for t,grace_month in enumerate(months): @@ -234,9 +203,62 @@ def regress_grace_maps(LMAX, RAD, grid = gravtk.spatial().from_list(spatial_list) spatial_list = None + # Setting output parameters for each fit type + coef_str = ['x{0:d}'.format(o) for o in range(ORDER+1)] + unit_suffix = [' yr^{0:d}'.format(-o) if o else '' for o in range(ORDER+1)] + if (ORDER == 0):# Mean + unit_longname = ['Mean'] + elif (ORDER == 1):# Trend + unit_longname = ['Constant','Trend'] + elif (ORDER == 2):# Quadratic + unit_longname = ['Constant','Linear','Quadratic'] + + # amplitude string for cyclical components + amp_str = [] + # extra terms for S2 tidal aliasing components or custom fits + TERMS = [] + for i,c in enumerate(CYCLES): + if (c == 0.5): + coef_str.extend(['SS','SC']) + amp_str.append('SEMI') + unit_longname.extend(['Semi-Annual Sine', 'Semi-Annual Cosine']) + unit_suffix.extend(['','']) + elif (c == 1.0): + coef_str.extend(['AS','AC']) + amp_str.append('ANN') + unit_longname.extend(['Annual Sine', 'Annual Cosine']) + unit_suffix.extend(['','']) + elif (c == (161.0/365.25)): + # create custom terms for S2 tidal aliasing during GRACE period + ii, = np.nonzero(grid.time < 2018.0) + S2sin = np.zeros_like(grid.time) + S2cos = np.zeros_like(grid.time) + S2sin[ii] = np.sin(2.0*np.pi*grid.time[ii]/np.float64(c)) + S2cos[ii] = np.cos(2.0*np.pi*grid.time[ii]/np.float64(c)) + TERMS.append(S2sin) + TERMS.append(S2cos) + coef_str.extend(['S2SGRC','S2CGRC']) + amp_str.append('S2GRC') + unit_longname.extend(['S2 Tidal Alias Sine', 'S2 Tidal Alias Cosine']) + unit_suffix.extend(['','']) + # create custom terms for S2 tidal aliasing during GRACE-FO period + ii, = np.nonzero(grid.time >= 2018.0) + S2sin = np.zeros_like(grid.time) + S2cos = np.zeros_like(grid.time) + S2sin[ii] = np.sin(2.0*np.pi*grid.time[ii]/np.float64(c)) + S2cos[ii] = np.cos(2.0*np.pi*grid.time[ii]/np.float64(c)) + TERMS.append(S2sin) + TERMS.append(S2cos) + coef_str.extend(['S2SGFO','S2CGFO']) + amp_str.append('S2GFO') + unit_longname.extend(['S2 Tidal Alias Sine', 'S2 Tidal Alias Cosine']) + unit_suffix.extend(['','']) + # remove the original S2 tidal aliasing term from CYCLES list + CYCLES.remove(c) + # Fitting seasonal components ncomp = len(coef_str) - ncycles = 2*len(CYCLES) + ncycles = 2*len(CYCLES) + len(TERMS) # Allocating memory for output variables out = dinput.zeros_like() @@ -257,7 +279,7 @@ def regress_grace_maps(LMAX, RAD, for j in range(nlon): # Calculating the regression coefficients tsbeta = gravtk.time_series.regress(grid.time, grid.data[i,j,:], - ORDER=ORDER, CYCLES=CYCLES, CONF=0.95) + ORDER=ORDER, CYCLES=CYCLES, TERMS=TERMS, CONF=0.95) # save regression components for k in range(0, ncomp): out.data[i,j,k] = tsbeta['beta'][k] @@ -302,10 +324,10 @@ def regress_grace_maps(LMAX, RAD, # if fitting coefficients with cyclical components if (ncycles > 0): # output spatial titles for amplitudes - amp_title = {'ANN':'Annual Amplitude','SEMI':'Semi-Annual Amplitude', - 'S2':'S2 Tidal Alias Amplitude'} - ph_title = {'ANN':'Annual Phase','SEMI':'Semi-Annual Phase', - 'S2':'S2 Tidal Alias Phase'} + amp_title = dict(ANN='Annual Amplitude',SEMI='Semi-Annual Amplitude', + S2GRC='S2 Tidal Alias Amplitude',S2GFO='S2 Tidal Alias Amplitude') + ph_title = dict(ANN='Annual Phase',SEMI='Semi-Annual Phase', + S2GRC='S2 Tidal Alias Phase',S2GFO='S2 Tidal Alias Phase') # output amplitude and phase of cyclical components for i,flag in enumerate(amp_str):