diff --git a/gravity_toolkit/__init__.py b/gravity_toolkit/__init__.py index 74db266f..e7d8816e 100644 --- a/gravity_toolkit/__init__.py +++ b/gravity_toolkit/__init__.py @@ -49,7 +49,6 @@ from gravity_toolkit.read_gfc_harmonics import read_gfc_harmonics from gravity_toolkit.read_GIA_model import read_GIA_model, gia from gravity_toolkit.read_GRACE_harmonics import read_GRACE_harmonics -from gravity_toolkit.read_ICGEM_harmonics import read_ICGEM_harmonics from gravity_toolkit.read_love_numbers import read_love_numbers, load_love_numbers, love_numbers from gravity_toolkit.read_SLR_harmonics import read_SLR_harmonics, convert_weekly from gravity_toolkit.sea_level_equation import sea_level_equation diff --git a/gravity_toolkit/hdf5_read.py b/gravity_toolkit/hdf5_read.py deleted file mode 100755 index 537fc398..00000000 --- a/gravity_toolkit/hdf5_read.py +++ /dev/null @@ -1,206 +0,0 @@ -#!/usr/bin/env python -u""" -hdf5_read.py -Written by Tyler Sutterley (11/2021) - -Reads spatial data from HDF5 files - -CALLING SEQUENCE: - dinput = hdf5_read(filename, DATE=False, VERBOSE=False) - -INPUTS: - filename: HDF5 file to be opened and read - -OUTPUTS: - data: z value of dataset - lon: longitudinal array - lat: latitudinal array - time: time value of dataset (if specified by DATE) - attributes: HDF5 attributes (for variables and title) - -OPTIONS: - DATE: HDF5 file has date information - VARNAME: z variable name in HDF5 file - LONNAME: longitude variable name in HDF5 file - LATNAME: latitude variable name in HDF5 file - TIMENAME: time variable name in HDF5 file - COMPRESSION: HDF5 file is compressed or streaming as bytes - gzip - zip - bytes - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - h5py: Python interface for Hierarchal Data Format 5 (HDF5) - (https://www.h5py.org) - -UPDATE HISTORY: - Updated 11/2021: try to get more global attributes. use kwargs - Updated 10/2021: using python logging for handling verbose output - Updated 02/2021: prevent warnings with python3 compatible regex strings - Updated 12/2020: try/except for getting variable unit attributes - attempt to get a standard set of attributes from each variable - add fallback for finding HDF5 file within from zip files - added bytes option for COMPRESSION if streaming from memory - Updated 08/2020: add options to read from gzip or zip compressed files - Updated 07/2020: added function docstrings - Updated 06/2020: output data as lat/lon following spatial module - attempt to read fill value attribute and set to None if not present - Updated 10/2019: changing Y/N flags to True/False - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 06/2018: extract fill_value and title without variable attributes - Updated 06/2016: using __future__ print function - Updated 05/2016: will only transpose if data is 2 dimensional (not 3) - added parameter to read the TITLE variable - Updated 05/2015: added parameter TIMENAME for time variable name - Updated 11/2014: got back to writing this - in working condition with updated attributes as in netcdf equivalent - Updated 12/2013: converted ncdf code to HDF5 code (alternative data type) - Updated 07/2013: switched from Scientific Python to Scipy - Updated 05/2013: converted to Python - Updated 03/2013: converted to Octave - Updated 01/2013: adding time variable - Written 07/2012 for GMT and for archiving datasets - Motivation for archival: netCDF files are much smaller than ascii - files and more portable/transferable than IDL .sav files - (possible to connect with geostatistics packages in R?) -""" -from __future__ import print_function - -import os -import io -import re -import gzip -import logging -import zipfile -import numpy as np -import warnings - -# attempt imports -try: - import h5py -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("h5py not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def hdf5_read(filename, **kwargs): - """ - Reads spatial data from HDF5 files - - Parameters - ---------- - filename: HDF5 file to be opened and read - - DATE: HDF5 file has date information - VARNAME: z variable name in HDF5 file - LONNAME: longitude variable name in HDF5 file - LATNAME: latitude variable name in HDF5 file - TIMENAME: time variable name in HDF5 file - COMPRESSION: HDF5 file is compressed or streaming as bytes - gzip - zip - bytes - - Returns - ------- - data: z value of dataset - lon: longitudinal array - lat: latitudinal array - time: time value of dataset - attributes: HDF5 attributes - """ - # set default keyword arguments - kwargs.setdefault('DATE',False) - kwargs.setdefault('VARNAME','z') - kwargs.setdefault('LONNAME','lon') - kwargs.setdefault('LATNAME','lat') - kwargs.setdefault('TIMENAME','time') - kwargs.setdefault('COMPRESSION',None) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use spatial.from_HDF5", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # Open the HDF5 file for reading - if (kwargs['COMPRESSION'] == 'gzip'): - # read gzip compressed file and extract into in-memory file object - with gzip.open(os.path.expanduser(filename),'r') as f: - fid = io.BytesIO(f.read()) - # set filename of BytesIO object - fid.filename = os.path.basename(filename) - # rewind to start of file - fid.seek(0) - # read as in-memory (diskless) HDF5 dataset from BytesIO object - fileID = h5py.File(fid, 'r') - elif (kwargs['COMPRESSION'] == 'zip'): - # read zipped file and extract file into in-memory file object - fileBasename,_ = os.path.splitext(os.path.basename(filename)) - with zipfile.ZipFile(os.path.expanduser(filename)) as z: - # first try finding a HDF5 file with same base filename - # if none found simply try searching for a HDF5 file - try: - f,=[f for f in z.namelist() if re.match(fileBasename,f,re.I)] - except: - f,=[f for f in z.namelist() if re.search(r'\.H(DF)?5$',f,re.I)] - # read bytes from zipfile into in-memory BytesIO object - fid = io.BytesIO(z.read(f)) - # set filename of BytesIO object - fid.filename = os.path.basename(filename) - # rewind to start of file - fid.seek(0) - # read as in-memory (diskless) HDF5 dataset from BytesIO object - fileID = h5py.File(fid, 'r') - elif (kwargs['COMPRESSION'] == 'bytes'): - # read as in-memory (diskless) HDF5 dataset - fileID = h5py.File(filename, 'r') - else: - # read HDF5 dataset - fileID = h5py.File(os.path.expanduser(filename), 'r') - # allocate python dictionary for output variables - dinput = {} - dinput['attributes'] = {} - - # Output HDF5 file information - logging.info(fileID.filename) - logging.info(list(fileID.keys())) - - # mapping between output keys and HDF5 variable names - keys = ['lon','lat','data'] - h5keys = [kwargs['LONNAME'],kwargs['LATNAME'],kwargs['VARNAME']] - if kwargs['DATE']: - keys.append('time') - h5keys.append(kwargs['TIMENAME']) - - # list of variable attributes - attributes_list = ['description','units','long_name','calendar', - 'standard_name','_FillValue','missing_value'] - # for each variable - for key,h5key in zip(keys,h5keys): - # Getting the data from each HDF5 variable - dinput[key] = np.squeeze(fileID[h5key][:]) - # Getting attributes of included variables - dinput['attributes'][key] = {} - for attr in attributes_list: - try: - dinput['attributes'][key][attr] = fileID[h5key].attrs[attr] - except (KeyError, AttributeError): - pass - - # switching data array to lat/lon if lon/lat - sz = dinput['data'].shape - if (dinput['data'].ndim == 2) and (len(dinput['lon']) == sz[0]): - dinput['data'] = dinput['data'].T - - # Global attributes - for att_name in ['title','description','reference']: - try: - dinput['attributes'][att_name] = fileID.attrs[att_name] - except (ValueError, KeyError, AttributeError): - pass - - # Closing the HDF5 file - fileID.close() - return dinput diff --git a/gravity_toolkit/hdf5_read_stokes.py b/gravity_toolkit/hdf5_read_stokes.py deleted file mode 100755 index 2773a219..00000000 --- a/gravity_toolkit/hdf5_read_stokes.py +++ /dev/null @@ -1,222 +0,0 @@ -#!/usr/bin/env python -u""" -hdf5_read_stokes.py -Written by Tyler Sutterley (11/2021) - -Reads spherical harmonic data from HDF5 files - -CALLING SEQUENCE: - Ylms = hdf5_read_stokes(filename, DATE=True, VERBOSE=False) - -INPUTS: - filename: HDF5 file to be opened and read - -OUTPUTS: - clm: cosine spherical harmonic coefficients - slm: sine spherical harmonic coefficients - l: degree (l) - m: order (m) - time: time of measurement (if specified by DATE) - month: GRACE/GRACE-FO month (if specified by DATE) - attributes: HDF5 attributes for: - spherical harmonics (clm,slm) - variables (l,m,time,month) - file description - -OPTIONS: - DATE: HDF5 file has date information - COMPRESSION: HDF5 file is compressed or streaming as bytes - gzip - zip - bytes - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - h5py: Python interface for Hierarchal Data Format 5 (HDF5) - (https://www.h5py.org) - -UPDATE HISTORY: - Updated 11/2021: try to get more global attributes. use kwargs - Updated 10/2021: using python logging for handling verbose output - Updated 02/2021: prevent warnings with python3 compatible regex strings - Updated 12/2020: try/except for getting variable unit attributes - add fallback for finding HDF5 file within from zip files - added bytes option for COMPRESSION if streaming from memory - Updated 09/2020: use try/except for reading attributes - Updated 08/2020: add options to read from gzip or zip compressed files - Updated 07/2020: added function docstrings - Updated 03/2020: added ATTRIBUTES option to check if file has attributes - Updated 10/2019: changing Y/N flags to True/False. check if time is array - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 06/2016: using __future__ print function - Updated 02/2016: capitalized LMAX and MMAX variables to match other programs - Updated 05/2015: minor change for MMAX != LMAX - Updated 11/2014: got back to writing this - in working condition with updated attributes as in netcdf equivalent - Updated 12/2013: converted ncdf code to HDF5 code (alternative data type) - Updated 07/2013: switched from Scientific Python to Scipy - Updated 03/2013: switched I/O to column arrays instead of matrix - Written 07/2012 -""" -from __future__ import print_function - -import os -import io -import re -import gzip -import logging -import zipfile -import numpy as np -import warnings - -# attempt imports -try: - import h5py -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("h5py not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def hdf5_read_stokes(filename, **kwargs): - """ - Reads spherical harmonic data from HDF5 files - - Parameters - ---------- - filename: HDF5 file to be opened and read - - DATE: HDF5 file has date information - VERBOSE: will print to screen the HDF5 structure parameters - COMPRESSION: HDF5 file is compressed or streaming as bytes - gzip - zip - bytes - - Returns - ------- - clm: cosine spherical harmonic coefficients - slm: sine spherical harmonic coefficients - l: degree - m: order - time: time of measurement - month: GRACE/GRACE-FO month - attributes: HDF5 attributes for variables and file - """ - # set default keyword arguments - kwargs.setdefault('DATE',True) - kwargs.setdefault('COMPRESSION',None) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use harmonics.from_HDF5", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # Open the HDF5 file for reading - if (kwargs['COMPRESSION'] == 'gzip'): - # read gzip compressed file and extract into in-memory file object - with gzip.open(os.path.expanduser(filename),'r') as f: - fid = io.BytesIO(f.read()) - # set filename of BytesIO object - fid.filename = os.path.basename(filename) - # rewind to start of file - fid.seek(0) - # read as in-memory (diskless) HDF5 dataset from BytesIO object - fileID = h5py.File(fid, 'r') - elif (kwargs['COMPRESSION'] == 'zip'): - # read zipped file and extract file into in-memory file object - fileBasename,_ = os.path.splitext(os.path.basename(filename)) - with zipfile.ZipFile(os.path.expanduser(filename)) as z: - # first try finding a HDF5 file with same base filename - # if none found simply try searching for a HDF5 file - try: - f,=[f for f in z.namelist() if re.match(fileBasename,f,re.I)] - except: - f,=[f for f in z.namelist() if re.search(r'\.H(DF)?5$',f,re.I)] - # read bytes from zipfile into in-memory BytesIO object - fid = io.BytesIO(z.read(f)) - # set filename of BytesIO object - fid.filename = os.path.basename(filename) - # rewind to start of file - fid.seek(0) - # read as in-memory (diskless) HDF5 dataset from BytesIO object - fileID = h5py.File(fid, 'r') - elif (kwargs['COMPRESSION'] == 'bytes'): - # read as in-memory (diskless) HDF5 dataset - fileID = h5py.File(filename, 'r') - else: - # read HDF5 dataset - fileID = h5py.File(os.path.expanduser(filename), 'r') - # allocate python dictionary for output variables - dinput = {} - dinput['attributes'] = {} - - # Output HDF5 file information - logging.info(fileID.filename) - logging.info(list(fileID.keys())) - - # output variable keys - h5keys = ['l','m','clm','slm'] - # Getting the data from each HDF5 variable - # converting HDF5 objects into numpy arrays - ll = np.array(fileID['l'][:]) - mm = np.array(fileID['m'][:]) - # Spherical harmonic files have date information - if kwargs['DATE']: - h5keys.extend(['time','month']) - dinput['time'] = fileID['time'][:].copy() - dinput['month'] = fileID['month'][:].copy() - n_time = len(dinput['time']) - else: - n_time = 0 - - # Restructuring input array back into matrix format - LMAX = np.max(ll) - MMAX = np.max(mm) - - # LMAX+1 to include LMAX (LMAX+1 elements) - dinput['l'] = np.arange(0,LMAX+1) - dinput['m'] = np.arange(0,MMAX+1) - # convert input clm/slm to numpy arrays - CLM = np.array(fileID['clm'][:]) - SLM = np.array(fileID['slm'][:]) - # size of the input grids - n_harm, = fileID['l'].shape - # import spherical harmonic data - if (kwargs['DATE'] and (n_time > 1)): - # contains multiple dates - dinput['clm'] = np.zeros((LMAX+1,MMAX+1,n_time)) - dinput['slm'] = np.zeros((LMAX+1,MMAX+1,n_time)) - for lm in range(n_harm): - dinput['clm'][ll[lm],mm[lm],:] = CLM[lm,:] - dinput['slm'][ll[lm],mm[lm],:] = SLM[lm,:] - else: - # contains either no dates or a single date - dinput['clm'] = np.zeros((LMAX+1,MMAX+1)) - dinput['slm'] = np.zeros((LMAX+1,MMAX+1)) - for lm in range(n_harm): - dinput['clm'][ll[lm],mm[lm]] = CLM[lm] - dinput['slm'][ll[lm],mm[lm]] = SLM[lm] - - # Getting attributes of clm/slm and included variables - # get attributes for the included variables - for key in h5keys: - try: - dinput['attributes'][key] = [ - fileID[key].attrs['units'], - fileID[key].attrs['long_name'] - ] - except (KeyError, AttributeError): - pass - # Global attributes - for att_name in ['title','description','reference']: - try: - dinput['attributes'][att_name] = fileID.attrs[att_name] - except (ValueError, KeyError, AttributeError): - pass - - # Closing the HDF5 file - fileID.close() - - # return the output variable - return dinput diff --git a/gravity_toolkit/hdf5_stokes.py b/gravity_toolkit/hdf5_stokes.py deleted file mode 100755 index 4dce7735..00000000 --- a/gravity_toolkit/hdf5_stokes.py +++ /dev/null @@ -1,214 +0,0 @@ -#!/usr/bin/env python -u""" -hdf5_stokes.py -Written by Tyler Sutterley (11/2021) - -Writes spherical harmonic coefficients to HDF5 files - -CALLING SEQUENCE: - hdf5_stokes(clm1,slm1,linp,minp,tinp,month,FILENAME=output_HDF5_file) - -INPUTS: - clm1: cosine spherical harmonic coefficients - slm1: sine spherical harmonic coefficients - linp: spherical harmonic degree (l) - minp: spherical harmonic order (m) - tinp: date of measurement - month: GRACE/GRACE-FO month - -OPTIONS: - FILENAME: output filename HDF5 - UNITS: spherical harmonic units - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - MONTHS_NAME: name of months variable within HDF5 file - MONTHS_UNITS: months variable units - MONTHS_LONGNAME: months variable description - TITLE: description attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing HDF5 file - DATE: harmonics have date information - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - h5py: Python interface for Hierarchal Data Format 5 (HDF5) - (https://www.h5py.org) - -UPDATE HISTORY: - Updated 11/2021: use remapped dictionary for filling HDF5 variables - Updated 10/2021: using python logging for handling verbose output - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 12/2020: added REFERENCE option to set file attribute - Updated 07/2020: added function docstrings - Updated 03/2020: only include title if not None - Updated 10/2019: changing Y/N flags to True/False - Updated 08/2019: don't include time (HH:MM:SS) in creation date - Updated 07/2019: added creation date as a global attribute - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 12/2018: using python dictionaries to improve readability - Updated 10/2018: using future division for python3 Compatibility - Updated 02/2017: added MONTHS_UNITS, MONTHS_LONGNAME, MONTHS_NAME parameters - aligned TIME_LONGNAME and TIME_UNITS with attributes - can output a HDF5 file with multiple dates similar to the netcdf program - Updated 06/2016: using __future__ print function - Updated 03/2016: direct calculation of number of harmonics n_harm - Updated 05/2015: minor change for MMAX != LMAX - Updated 11/2014: got back to writing this - in working condition with updated attributes as in netcdf equivalent - Updated 12/2013: converted ncdf code to HDF5 code (alternative data type) - Updated 07/2013: switched from Scientific Python to Scipy - Updated 05/2013 made UNITS an option in case converting the units to - mass harmonics or other harmonic variant - Updated 03/2013: added units to clm and slm as 'Geodesy Normalization' - switched I/O to column arrays for smaller file sizes and compatibility - between languages - made date an option for datasets that have no date - Updated 01/2013 to add time and GRACE/GRACE-FO month number - Written 07/2012 -""" -from __future__ import print_function, division - -import time -import logging -import numpy as np -import warnings - -# attempt imports -try: - import h5py -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("h5py not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def hdf5_stokes(clm1, slm1, linp, minp, tinp, month, **kwargs): - """ - Writes spherical harmonic coefficients to HDF5 files - - Parameters - ---------- - clm1: cosine spherical harmonic coefficients - slm1: sine spherical harmonic coefficients - linp: spherical harmonic degree (l) - minp: spherical harmonic order (m) - tinp: date of measurement - month: GRACE/GRACE-FO month - - FILENAME: HDF5 filename - UNITS: spherical harmonic units - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - MONTHS_NAME: name of months variable within HDF5 file - MONTHS_UNITS: months variable units - MONTHS_LONGNAME: months variable description - TITLE: description attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing HDF5 file - DATE: harmonics have date information - """ - # set default keyword arguments - kwargs.setdefault('FILENAME',None) - kwargs.setdefault('UNITS','Geodesy_Normalization') - kwargs.setdefault('TIME_UNITS',None) - kwargs.setdefault('TIME_LONGNAME',None) - kwargs.setdefault('MONTHS_NAME','month') - kwargs.setdefault('MONTHS_UNITS','number') - kwargs.setdefault('MONTHS_LONGNAME','GRACE_month') - kwargs.setdefault('TITLE',None) - kwargs.setdefault('REFERENCE',None) - kwargs.setdefault('DATE',True) - kwargs.setdefault('CLOBBER',True) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use harmonics.to_HDF5", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # setting HDF5 clobber attribute - clobber = 'w' if kwargs['CLOBBER'] else 'w-' - # opening HDF5 file for writing - fileID = h5py.File(kwargs['FILENAME'], clobber) - - # Maximum spherical harmonic degree (LMAX) and order (MMAX) - LMAX = np.max(linp) - MMAX = np.max(minp) - # Calculating the number of cos and sin harmonics up to LMAX - # taking into account MMAX (if MMAX == LMAX then LMAX-MMAX=0) - n_harm = (LMAX**2 + 3*LMAX - (LMAX-MMAX)**2 - (LMAX-MMAX))//2 + 1 - - # dictionary with output variables - output = {} - # restructured degree and order - output['l'] = np.zeros((n_harm,), dtype=np.int32) - output['m'] = np.zeros((n_harm,), dtype=np.int32) - # Restructuring output matrix to array format - # will reduce matrix size and insure compatibility between platforms - if kwargs['DATE']: - n_time = len(np.atleast_1d(tinp)) - output['time'] = np.copy(tinp) - output['month'] = np.copy(month) - if (n_time == 1): - output['clm'] = np.zeros((n_harm)) - output['slm'] = np.zeros((n_harm)) - else: - output['clm'] = np.zeros((n_harm,n_time)) - output['slm'] = np.zeros((n_harm,n_time)) - else: - n_time = 0 - output['clm'] = np.zeros((n_harm)) - output['slm'] = np.zeros((n_harm)) - - # create counter variable lm - lm = 0 - for m in range(0,MMAX+1):# MMAX+1 to include MMAX - for l in range(m,LMAX+1):# LMAX+1 to include LMAX - output['l'][lm] = np.int64(l) - output['m'][lm] = np.int64(m) - if (kwargs['DATE'] and (n_time > 1)): - output['clm'][lm,:] = clm1[l,m,:] - output['slm'][lm,:] = slm1[l,m,:] - else: - output['clm'][lm] = clm1[l,m] - output['slm'][lm] = slm1[l,m] - # add 1 to lm counter variable - lm += 1 - - # Defining the HDF5 dataset variables - h5 = {} - for key,val in output.items(): - h5[key] = fileID.create_dataset(key, val.shape, - data=val, dtype=val.dtype, compression='gzip') - - # filling HDF5 dataset attributes - # Defining attributes for degree and order - h5['l'].attrs['long_name'] = 'spherical_harmonic_degree'# degree long name - h5['l'].attrs['units'] = 'Wavenumber'# SH degree units - h5['m'].attrs['long_name'] = 'spherical_harmonic_order'# order long name - h5['m'].attrs['units'] = 'Wavenumber'# SH order units - # Defining attributes for dataset - h5['clm'].attrs['long_name'] = 'cosine_spherical_harmonics' - h5['clm'].attrs['units'] = kwargs['UNITS'] - h5['slm'].attrs['long_name'] = 'sine_spherical_harmonics' - h5['slm'].attrs['units'] = kwargs['UNITS'] - if kwargs['DATE']: - # Defining attributes for date and month (or integer date) - h5['time'].attrs['long_name'] = kwargs['TIME_LONGNAME'] - h5['time'].attrs['units'] = kwargs['TIME_UNITS'] - h5['month'].attrs['long_name'] = kwargs['MONTHS_LONGNAME'] - h5['month'].attrs['units'] = kwargs['MONTHS_UNITS'] - # description of file - if kwargs['TITLE']: - fileID.attrs['description'] = kwargs['TITLE'] - # reference of file - if kwargs['REFERENCE']: - fileID.attrs['reference'] = kwargs['REFERENCE'] - # date created - fileID.attrs['date_created'] = time.strftime('%Y-%m-%d',time.localtime()) - - # Output HDF5 structure information - logging.info(kwargs['FILENAME']) - logging.info(list(fileID.keys())) - - # Closing the HDF5 file - fileID.close() diff --git a/gravity_toolkit/hdf5_write.py b/gravity_toolkit/hdf5_write.py deleted file mode 100755 index 83bf18cf..00000000 --- a/gravity_toolkit/hdf5_write.py +++ /dev/null @@ -1,185 +0,0 @@ -#!/usr/bin/env python -u""" -hdf5_write.py -Written by Tyler Sutterley (11/2021) - -Writes spatial data to HDF5 files - -CALLING SEQUENCE: - hdf5_write(data, lon, lat, tim, FILENAME=output_HDF5_file) - -INPUTS: - data: z data - lon: longitude array - lat: latitude array - tim: time array - -OPTIONS: - FILENAME: output filename HDF5 - VARNAME: z variable name in HDF5 file - LONNAME: longitude variable name in HDF5 file - LATNAME: latitude variable name in HDF5 file - UNITS: z variable units - LONGNAME: z variable description - FILL_VALUE: missing value for z variable - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - TITLE: description attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing HDF5 file - DATE: data has date information - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - h5py: Python interface for Hierarchal Data Format 5 (HDF5) - (https://www.h5py.org) - -UPDATE HISTORY: - Updated 11/2021: use remapped dictionary for filling HDF5 variables - Updated 10/2021: using python logging for handling verbose output - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 12/2020: added REFERENCE option to set file attribute - Updated 07/2020: added function docstrings - Updated 04/2020: added option DATE if including time data - Updated 03/2020: only include title if not None - Updated 10/2019: changing Y/N flags to True/False - Updated 09/2019 for public release - Updated 08/2019: don't include time (HH:MM:SS) in creation date - Updated 07/2019: added creation date as a global attribute - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 08/2018: use n_time variable for output HDF5 dimensions - Updated 03/2018: added option TIMENAME to specify the variable name of time - Updated 02/2017: TIME_LONGNAME and TIME_UNITS with attributes, - updated TIME_LONGNAME to Date_in_Decimal_Years - Updated 06/2016: using __future__ print function - Updated 05/2016: output data types same as input data types - Updated 11/2014: got back to writing this - in working condition with updated attributes as in netcdf equivalent - Updated 12/2013: converted ncdf code to HDF5 code (alternative data type) - Updated 07/2013: switched from Scientific Python to Scipy - Updated 01/2013: adding time as a variable - Updated 10/2012: changed from variable names x and y to lon and lat. - Written 07/2012 -""" -from __future__ import print_function - -import time -import logging -import numpy as np -import warnings - -# attempt imports -try: - import h5py -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("h5py not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def hdf5_write(data, lon, lat, tim, **kwargs): - """ - Writes spatial data to HDF5 files - - Parameters - ---------- - data: z data - lon: longitude array - lat: latitude array - tim: time array - - FILENAME: HDF5 filename - VARNAME: z variable name in HDF5 file - LONNAME: longitude variable name in HDF5 file - LATNAME: latitude variable name in HDF5 file - UNITS: z variable units - LONGNAME: z variable description - FILL_VALUE: missing value for z variable - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - TITLE: description attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing HDF5 file - DATE: data has date information - """ - # set default keyword arguments - kwargs.setdefault('FILENAME',None) - kwargs.setdefault('VARNAME','z') - kwargs.setdefault('LONNAME','lon') - kwargs.setdefault('LATNAME','lat') - kwargs.setdefault('TIMENAME','time') - kwargs.setdefault('UNITS',None) - kwargs.setdefault('LONGNAME',None) - kwargs.setdefault('FILL_VALUE',None) - kwargs.setdefault('TIME_UNITS',None) - kwargs.setdefault('TIME_LONGNAME',None) - kwargs.setdefault('TITLE',None) - kwargs.setdefault('REFERENCE',None) - kwargs.setdefault('DATE',True) - kwargs.setdefault('CLOBBER',True) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use spatial.to_HDF5", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # setting HDF5 clobber attribute - clobber = 'w' if kwargs['CLOBBER'] else 'w-' - - # opening HDF5 file for writing - fileID = h5py.File(kwargs['FILENAME'], clobber) - - # create output dictionary with key mapping - output = {} - output[kwargs['LONNAME']] = np.copy(lon) - output[kwargs['LATNAME']] = np.copy(lat) - dimensions = [kwargs['LATNAME'],kwargs['LONNAME']] - # extend with date variables - if kwargs['DATE']: - output[kwargs['TIMENAME']] = np.array(tim,dtype='f') - output[kwargs['VARNAME']] = np.atleast_3d(data) - dimensions.append(kwargs['TIMENAME']) - else: - output[kwargs['VARNAME']] = np.copy(data) - - # Defining the HDF5 dataset variables - h5 = {} - for key,val in output.items(): - h5[key] = fileID.create_dataset(key, val.shape, data=val, - dtype=val.dtype, compression='gzip') - # add dimensions - for i,dim in enumerate(dimensions): - h5[kwargs['VARNAME']].dims[i].label = dim - h5[kwargs['VARNAME']].dims[i].attach_scale(h5[dim]) - - # filling HDF5 dataset attributes - # Defining attributes for longitude and latitude - h5[kwargs['LONNAME']].attrs['long_name'] = 'longitude' - h5[kwargs['LONNAME']].attrs['units'] = 'degrees_east' - h5[kwargs['LATNAME']].attrs['long_name'] = 'latitude' - h5[kwargs['LATNAME']].attrs['units'] = 'degrees_north' - # Defining attributes for dataset - h5[kwargs['VARNAME']].attrs['long_name'] = kwargs['LONGNAME'] - h5[kwargs['VARNAME']].attrs['units'] = kwargs['UNITS'] - # Dataset contains missing values - if (kwargs['FILL_VALUE'] is not None): - h5[kwargs['VARNAME']].attrs['_FillValue'] = kwargs['FILL_VALUE'] - # Defining attributes for date - if kwargs['DATE']: - h5[kwargs['TIMENAME']].attrs['long_name'] = kwargs['TIME_LONGNAME'] - h5[kwargs['TIMENAME']].attrs['units'] = kwargs['TIME_UNITS'] - # description of file - if kwargs['TITLE']: - fileID.attrs['description'] = kwargs['TITLE'] - # reference of file - if kwargs['REFERENCE']: - fileID.attrs['reference'] = kwargs['REFERENCE'] - # date created - fileID.attrs['date_created'] = time.strftime('%Y-%m-%d',time.localtime()) - - # Output HDF5 structure information - logging.info(kwargs['FILENAME']) - logging.info(list(fileID.keys())) - - # Closing the HDF5 file - fileID.close() diff --git a/gravity_toolkit/ncdf_read.py b/gravity_toolkit/ncdf_read.py deleted file mode 100755 index 4aaa5ffc..00000000 --- a/gravity_toolkit/ncdf_read.py +++ /dev/null @@ -1,204 +0,0 @@ -#!/usr/bin/env python -u""" -ncdf_read.py -Written by Tyler Sutterley (11/2021) - -Reads spatial data from COARDS-compliant netCDF4 files - -CALLING SEQUENCE: - dinput = ncdf_read(filename, DATE=False, VERBOSE=False) - -INPUTS: - filename: netCDF4 file to be opened and read - -OUTPUTS: - data: z value of dataset - lon: longitudinal array - lat: latitudinal array - time: time value of dataset (if specified by DATE) - attributes: netCDF4 attributes (for variables and title) - -OPTIONS: - DATE: netCDF4 file has date information - VERBOSE: will print to screen the netCDF4 structure parameters - VARNAME: z variable name in netCDF4 file - LONNAME: longitude variable name in netCDF4 file - LATNAME: latitude variable name in netCDF4 file - TIMENAME: time variable name in netCDF4 file - COMPRESSION: netCDF4 file is compressed or streaming as bytes - gzip - zip - bytes - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - netCDF4: Python interface to the netCDF C library - (https://unidata.github.io/netcdf4-python/netCDF4/index.html) - -UPDATE HISTORY: - Updated 11/2021: try to get more global attributes. use kwargs - Updated 10/2021: using python logging for handling verbose output - Updated 02/2021: prevent warnings with python3 compatible regex strings - Updated 12/2020: try/except for getting variable unit attributes - attempt to get a standard set of attributes from each variable - add fallback for finding netCDF4 file within from zip files - added bytes option for COMPRESSION if streaming from memory - Updated 08/2020: flake8 compatible regular expression strings - add options to read from gzip or zip compressed files - Updated 07/2020: added function docstrings - Updated 06/2020: output data as lat/lon following spatial module - attempt to read fill value attribute and set to None if not present - Updated 10/2019: changing Y/N flags to True/False - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 06/2018: extract fill_value and title without variable attributes - Updated 07-09/2016: using netCDF4-python - Updated 06/2016: using __future__ print, output filename if VERBOSE - Updated 05/2016: will only transpose if data is 2 dimensional (not 3) - added parameter to read the TITLE variable - Updated 07/2015: updated read title for different cases with regex - Updated 05/2015: added parameter TIMENAME for time variable name - Updated 04/2015: fix attribute outputs (forgot to copy to new dictionary) - Updated 02/2015: added copy for variable outputs - fixes new error flag from mmap=True - Updated 11/2014: new parameters for variable names and attributes - all variables in a single python dictionary - Updated 05/2014: new parameter for missing value - new outputs: all attributes, fill value - added try for TITLE attribute - converting time to numpy array - Updated 02/2014: minor update to if statements - Updated 07/2013: switched from Scientific Python to Scipy - Updated 01/2013: adding time variable - Written 07/2012 -""" -from __future__ import print_function - -import os -import re -import uuid -import gzip -import logging -import zipfile -import numpy as np -import warnings - -# attempt imports -try: - import netCDF4 -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("netCDF4 not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def ncdf_read(filename, **kwargs): - """ - Reads spatial data from COARDS-compliant netCDF4 files - - Parameters - ---------- - filename: netCDF4 file to be opened and read - - DATE: netCDF4 file has date information - VERBOSE: will print to screen the netCDF4 structure parameters - VARNAME: z variable name in netCDF4 file - LONNAME: longitude variable name in netCDF4 file - LATNAME: latitude variable name in netCDF4 file - TIMENAME: time variable name in netCDF4 file - COMPRESSION: netCDF4 file is compressed or streaming as bytes - gzip - zip - bytes - - Returns - ------- - data: z value of dataset - lon: longitudinal array - lat: latitudinal array - time: time value of dataset - attributes: netCDF4 attributes - """ - # set default keyword arguments - kwargs.setdefault('DATE',False) - kwargs.setdefault('VARNAME','z') - kwargs.setdefault('LONNAME','lon') - kwargs.setdefault('LATNAME','lat') - kwargs.setdefault('TIMENAME','time') - kwargs.setdefault('COMPRESSION',None) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use spatial.from_netCDF4", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # Open the NetCDF4 file for reading - if (kwargs['COMPRESSION'] == 'gzip'): - # read as in-memory (diskless) netCDF4 dataset - with gzip.open(os.path.expanduser(filename),'r') as f: - fileID = netCDF4.Dataset(os.path.basename(filename),memory=f.read()) - elif (kwargs['COMPRESSION'] == 'zip'): - # read zipped file and extract file into in-memory file object - fileBasename,_ = os.path.splitext(os.path.basename(filename)) - with zipfile.ZipFile(os.path.expanduser(filename)) as z: - # first try finding a netCDF4 file with same base filename - # if none found simply try searching for a netCDF4 file - try: - f,=[f for f in z.namelist() if re.match(fileBasename,f,re.I)] - except: - f,=[f for f in z.namelist() if re.search(r'\.nc(4)?$',f)] - # read bytes from zipfile as in-memory (diskless) netCDF4 dataset - fileID = netCDF4.Dataset(uuid.uuid4().hex, memory=z.read(f)) - elif (kwargs['COMPRESSION'] == 'bytes'): - # read as in-memory (diskless) netCDF4 dataset - fileID = netCDF4.Dataset(uuid.uuid4().hex, memory=filename.read()) - else: - # read netCDF4 dataset - fileID = netCDF4.Dataset(os.path.expanduser(filename), 'r') - # create python dictionary for output variables - dinput = {} - dinput['attributes'] = {} - - # Output NetCDF file information - logging.info(fileID.filepath()) - logging.info(list(fileID.variables.keys())) - - # mapping between output keys and netCDF4 variable names - keys = ['lon','lat','data'] - nckeys = [kwargs['LONNAME'],kwargs['LATNAME'],kwargs['VARNAME']] - if kwargs['DATE']: - keys.append('time') - nckeys.append(kwargs['TIMENAME']) - # list of variable attributes - attributes_list = ['description','units','long_name','calendar', - 'standard_name','_FillValue','missing_value'] - # for each variable - for key,nckey in zip(keys,nckeys): - # Getting the data from each NetCDF variable - dinput[key] = np.squeeze(fileID.variables[nckey][:].data) - # Getting attributes of included variables - dinput['attributes'][key] = {} - for attr in attributes_list: - # try getting the attribute - try: - dinput['attributes'][key][attr] = \ - fileID.variables[nckey].getncattr(attr) - except (KeyError,ValueError,AttributeError): - pass - - # switching data array to lat/lon if lon/lat - sz = dinput['data'].shape - if (dinput['data'].ndim == 2) and (len(dinput['lon']) == sz[0]): - dinput['data'] = dinput['data'].T - - # Global attributes - for att_name in ['title','description','reference']: - try: - ncattr, = [s for s in fileID.ncattrs() - if re.match(att_name,s,re.I)] - dinput['attributes'][att_name] = fileID.getncattr(ncattr) - except (ValueError, KeyError, AttributeError): - pass - # Closing the NetCDF file - fileID.close() - # return the output variable - return dinput diff --git a/gravity_toolkit/ncdf_read_stokes.py b/gravity_toolkit/ncdf_read_stokes.py deleted file mode 100755 index a5994881..00000000 --- a/gravity_toolkit/ncdf_read_stokes.py +++ /dev/null @@ -1,220 +0,0 @@ -#!/usr/bin/env python -u""" -ncdf_read_stokes.py -Written by Tyler Sutterley (11/2021) - -Reads spherical harmonic data from netCDF4 files - -CALLING SEQUENCE: - Ylms = ncdf_read_stokes(filename, DATE=True, VERBOSE=False) - -INPUTS: - filename: netCDF4 file to be opened and read - -OUTPUTS: - clm: Cosine Stokes Coefficient - slm: Sine Stokes Coefficient - l: degree (l) - m: order (m) - time: time of measurement (if specified by DATE) - month: GRACE/GRACE-FO month (if specified by DATE) - attributes: netCDF4 attributes for: - spherical harmonics (clm,slm) - variables (l,m,time,month) - file title - -OPTIONS: - DATE: netCDF4 file has date information - COMPRESSION: netCDF4 file is compressed or streaming as bytes - gzip - zip - bytes - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - netCDF4: Python interface to the netCDF C library - (https://unidata.github.io/netcdf4-python/netCDF4/index.html) - -UPDATE HISTORY: - Updated 11/2021: try to get more global attributes. use kwargs - Updated 10/2021: using python logging for handling verbose output - Updated 02/2021: prevent warnings with python3 compatible regex strings - Updated 12/2020: try/except for getting variable unit attributes - add fallback for finding netCDF4 file within from zip files - added bytes option for COMPRESSION if streaming from memory - Updated 09/2020: use try/except for reading attributes - Updated 08/2020: flake8 compatible regular expression strings - add options to read from gzip or zip compressed files - Updated 07/2020: added function docstrings - Updated 03/2020: added ATTRIBUTES option to check if file has attributes - Updated 10/2019: changing Y/N flags to True/False - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 09/2016: slicing of clm and slm on numpy arrays not netcdf variables - Updated 07/2016: using netCDF4-python - Updated 06/2016: using __future__ print, output filename if VERBOSE - Updated 02/2016: capitalized LMAX and MMAX variables to match other programs - Updated 07/2015: updated read title for different cases with regex - Updated 06/2015: can input a single netcdf with multiple dates - Updated 05/2015: minor change for MMAX != LMAX - Updated 02/2015: simplified attributes with for loop - Updated 01/2015: added copy for variable outputs - fixes new error flag from mmap=True - Updated 11/2014: all variables in a single python dictionary - Updated 05/2014: converted time and month to numpy arrays - Updated 05/2014: output all attributes under single variable - added try for TITLE attribute - Updated 02/2014: minor update to if statements - Updated 07/2013: switched from Scientific Python to Scipy - Updated 07/2013: switched from Scientific Python to Scipy - Updated 03/2013: switched I/O to column arrays instead of matrix - Written 07/2012 -""" -from __future__ import print_function - -import os -import re -import uuid -import gzip -import logging -import zipfile -import numpy as np -import warnings - -# attempt imports -try: - import netCDF4 -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("netCDF4 not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def ncdf_read_stokes(filename, **kwargs): - """ - Reads spherical harmonic data from netCDF4 files - - Parameters - ---------- - filename: netCDF4 file to be opened and read - - DATE: netCDF4 file has date information - COMPRESSION: netCDF4 file is compressed or streaming as bytes - gzip - zip - bytes - - Returns - ------- - clm: cosine spherical harmonic coefficients - slm: sine spherical harmonic coefficients - l: degree - m: order - time: time of measurement - month: GRACE/GRACE-FO month - attributes: netCDF4 attributes for variables and file - """ - # set default keyword arguments - kwargs.setdefault('DATE',True) - kwargs.setdefault('COMPRESSION',None) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use harmonics.from_netCDF4", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # Open the NetCDF4 file for reading - if (kwargs['COMPRESSION'] == 'gzip'): - # read as in-memory (diskless) netCDF4 dataset - with gzip.open(os.path.expanduser(filename),'r') as f: - fileID = netCDF4.Dataset(os.path.basename(filename),memory=f.read()) - elif (kwargs['COMPRESSION'] == 'zip'): - # read zipped file and extract file into in-memory file object - fileBasename,_ = os.path.splitext(os.path.basename(filename)) - with zipfile.ZipFile(os.path.expanduser(filename)) as z: - # first try finding a netCDF4 file with same base filename - # if none found simply try searching for a netCDF4 file - try: - f,=[f for f in z.namelist() if re.match(fileBasename,f,re.I)] - except: - f,=[f for f in z.namelist() if re.search(r'\.nc(4)?$',f)] - # read bytes from zipfile as in-memory (diskless) netCDF4 dataset - fileID = netCDF4.Dataset(uuid.uuid4().hex, memory=z.read(f)) - elif (kwargs['COMPRESSION'] == 'bytes'): - # read as in-memory (diskless) netCDF4 dataset - fileID = netCDF4.Dataset(uuid.uuid4().hex, memory=filename.read()) - else: - # read netCDF4 dataset - fileID = netCDF4.Dataset(os.path.expanduser(filename), 'r') - # create python dictionary for output variables - dinput = {} - dinput['attributes'] = {} - - # Output NetCDF file information - logging.info(fileID.filepath()) - logging.info(list(fileID.variables.keys())) - - # Getting the data from each NetCDF variable - # converting NetCDF objects into numpy arrays - nckeys = ['l','m','clm','slm'] - ll = fileID.variables['l'][:].copy() - mm = fileID.variables['m'][:].copy() - clm = fileID.variables['clm'][:].copy() - slm = fileID.variables['slm'][:].copy() - # read date variables if specified - if kwargs['DATE']: - nckeys.extend(['time','month']) - dinput['time'] = fileID.variables['time'][:].copy() - dinput['month'] = fileID.variables['month'][:].copy() - n_time = len(dinput['time']) - else: - n_time = 0 - - # Restructuring input array back into matrix format - LMAX = np.max(ll) - MMAX = np.max(mm) - # output spherical harmonic degree and order - # LMAX+1 to include LMAX (LMAX+1 elements) - dinput['l'] = np.arange(0,LMAX+1) - dinput['m'] = np.arange(0,MMAX+1) - # number of harmonics - n_harm, = fileID.variables['l'].shape - # import spherical harmonic data - if (kwargs['DATE'] and (n_time > 1)): - # contains multiple dates - dinput['clm'] = np.zeros((LMAX+1,MMAX+1,n_time)) - dinput['slm'] = np.zeros((LMAX+1,MMAX+1,n_time)) - for lm in range(n_harm): - dinput['clm'][ll[lm],mm[lm],:] = clm[lm,:] - dinput['slm'][ll[lm],mm[lm],:] = slm[lm,:] - else: - # contains either no dates or a single date - dinput['clm'] = np.zeros((LMAX+1,MMAX+1)) - dinput['slm'] = np.zeros((LMAX+1,MMAX+1)) - for lm in range(n_harm): - dinput['clm'][ll[lm],mm[lm]] = clm[lm] - dinput['slm'][ll[lm],mm[lm]] = slm[lm] - - # Getting attributes of clm/slm and included variables - # get attributes for the included variables - for key in nckeys: - try: - dinput['attributes'][key] = [ - fileID.variables[key].units, - fileID.variables[key].long_name - ] - except (KeyError,ValueError,AttributeError): - pass - # Global attributes - for att_name in ['title','description','reference']: - try: - ncattr, = [s for s in fileID.ncattrs() - if re.match(att_name,s,re.I)] - dinput['attributes'][att_name] = fileID.getncattr(ncattr) - except (ValueError, KeyError, AttributeError): - pass - - # Closing the NetCDF file - fileID.close() - - # return output variable - return dinput diff --git a/gravity_toolkit/ncdf_stokes.py b/gravity_toolkit/ncdf_stokes.py deleted file mode 100755 index e46ec6bf..00000000 --- a/gravity_toolkit/ncdf_stokes.py +++ /dev/null @@ -1,235 +0,0 @@ -#!/usr/bin/env python -u""" -ncdf_stokes.py -Written by Tyler Sutterley (11/2021) - -Writes spherical harmonic coefficients to netCDF4 files - -CALLING SEQUENCE: - ncdf_stokes(clm1,slm1,linp,minp,tinp,month,FILENAME=output_netcdf4_file) - -INPUTS: - clm1: cosine spherical harmonic coefficients - slm1: sine spherical harmonic coefficients - linp: spherical harmonic degree (l) - minp: spherical harmonic order (m) - tinp: date of measurement - month: GRACE/GRACE-FO month - -OPTIONS: - FILENAME: output netCDF4 filename - UNITS: spherical harmonic units - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - MONTHS_NAME: name of months variable within netCDF4 file - MONTHS_UNITS: months variable units - MONTHS_LONGNAME: months variable description - TITLE: title attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing netCDF4 file - DATE: harmonics have date information - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - netCDF4: Python interface to the netCDF C library - (https://unidata.github.io/netcdf4-python/netCDF4/index.html) - -UPDATE HISTORY: - Updated 11/2021: use remapped dictionary for filling HDF5 variables - Updated 10/2021: using python logging for handling verbose output - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 12/2020: added REFERENCE option to set file attribute - Updated 07/2020: added function docstrings - Updated 03/2020: only include title if not None - Updated 10/2019: changing Y/N flags to True/False - Updated 08/2019: don't include time (HH:MM:SS) in creation date - Updated 07/2019: added creation date as a global attribute - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 12/2018: using python dictionaries to improve readability - Updated 10/2018: using future division for python3 Compatibility - Updated 02/2017: added MONTHS_UNITS, MONTHS_LONGNAME, MONTHS_NAME parameters - aligned TIME_LONGNAME and TIME_UNITS with attributes - Updated 07/2016: using netCDF4-python - Updated 06/2016: using __future__ print function - Updated 03/2016: direct calculation of number of harmonics n_harm - Updated 07/2015: forgot to add case for DATE=False - Updated 06/2015: can output single netcdf with multiple dates - Updated 05/2015: minor change for MMAX != LMAX - Updated 05/2014: new parameters for time attributes - Updated 02/2014: minor update to if statements - Updated 07/2013: switched from Scientific Python to Scipy - Updated 05/2013 made UNITS an option in case converting the units to - mass harmonics or other harmonic variant - Updated 03/2013: added units to clm and slm as 'Geodesy Normalization' - switched I/O to column arrays for smaller file sizes and compatibility - between languages - made date an option for datasets that have no date (e.g. GIA) - Updated 01/2013 to add time and GRACE/GRACE-FO month number - Written 07/2012 -""" -from __future__ import print_function, division - -import time -import logging -import numpy as np -import warnings - -# attempt imports -try: - import netCDF4 -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("netCDF4 not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def ncdf_stokes(clm1, slm1, linp, minp, tinp, month, **kwargs): - """ - Writes spherical harmonic coefficients to netCDF4 files - - Parameters - ---------- - clm1: cosine spherical harmonic coefficients - slm1: sine spherical harmonic coefficients - linp: spherical harmonic degree (l) - minp: spherical harmonic order (m) - tinp: date of measurement - month: GRACE/GRACE-FO month - - FILENAME: netCDF4 filename - UNITS: spherical harmonic units - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - MONTHS_NAME: name of months variable within netCDF4 file - MONTHS_UNITS: months variable units - MONTHS_LONGNAME: months variable description - TITLE: title attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing netCDF4 file - DATE: harmonics have date information - """ - # set default keyword arguments - kwargs.setdefault('FILENAME',None) - kwargs.setdefault('UNITS','Geodesy_Normalization') - kwargs.setdefault('TIME_UNITS',None) - kwargs.setdefault('TIME_LONGNAME',None) - kwargs.setdefault('MONTHS_NAME','month') - kwargs.setdefault('MONTHS_UNITS','number') - kwargs.setdefault('MONTHS_LONGNAME','GRACE_month') - kwargs.setdefault('TITLE',None) - kwargs.setdefault('REFERENCE',None) - kwargs.setdefault('DATE',True) - kwargs.setdefault('CLOBBER',True) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use harmonics.to_netCDF4", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # setting NetCDF clobber attribute - clobber = 'w' if kwargs['CLOBBER'] else 'a' - # opening netCDF file for writing - fileID = netCDF4.Dataset(kwargs['FILENAME'], clobber, format="NETCDF4") - - # Maximum spherical harmonic degree (LMAX) and order (MMAX) - LMAX = np.max(linp) - MMAX = np.max(minp) - # Calculating the number of cos and sin harmonics up to LMAX - # taking into account MMAX (if MMAX == LMAX then LMAX-MMAX=0) - n_harm = (LMAX**2 + 3*LMAX - (LMAX-MMAX)**2 - (LMAX-MMAX))//2 + 1 - - # dictionary with output variables - output = {} - # restructured degree and order - output['l'] = np.zeros((n_harm,), dtype=np.int32) - output['m'] = np.zeros((n_harm,), dtype=np.int32) - # Restructuring output matrix to array format - # will reduce matrix size and insure compatibility between platforms - if kwargs['DATE']: - n_time = len(np.atleast_1d(tinp)) - output['time'] = np.copy(tinp) - output['month'] = np.copy(month) - if (n_time == 1): - output['clm'] = np.zeros((n_harm)) - output['slm'] = np.zeros((n_harm)) - else: - output['clm'] = np.zeros((n_harm,n_time)) - output['slm'] = np.zeros((n_harm,n_time)) - else: - n_time = 0 - output['clm'] = np.zeros((n_harm)) - output['slm'] = np.zeros((n_harm)) - - # create counter variable lm - lm = 0 - for m in range(0,MMAX+1):# MMAX+1 to include MMAX - for l in range(m,LMAX+1):# LMAX+1 to include LMAX - output['l'][lm] = np.int64(l) - output['m'][lm] = np.int64(m) - if (kwargs['DATE'] and (n_time > 1)): - output['clm'][lm,:] = clm1[l,m,:] - output['slm'][lm,:] = slm1[l,m,:] - else: - output['clm'][lm] = clm1[l,m] - output['slm'][lm] = slm1[l,m] - # add 1 to lm counter variable - lm += 1 - - # Defining the netCDF dimensions - fileID.createDimension('lm', n_harm) - if kwargs['DATE']: - fileID.createDimension('time', n_time) - - # defining the netCDF variables - nc = {} - # degree and order - nc['l'] = fileID.createVariable('l', 'i', ('lm',)) - nc['m'] = fileID.createVariable('m', 'i', ('lm',)) - # spherical harmonics - if (kwargs['DATE'] and (n_time > 1)): - nc['clm'] = fileID.createVariable('clm', 'd', ('lm','time',)) - nc['slm'] = fileID.createVariable('slm', 'd', ('lm','time',)) - else: - nc['clm'] = fileID.createVariable('clm', 'd', ('lm',)) - nc['slm'] = fileID.createVariable('slm', 'd', ('lm',)) - if kwargs['DATE']: - # time (in decimal form) - nc['time'] = fileID.createVariable('time', 'd', ('time',)) - # GRACE/GRACE-FO month (or integer date) - nc['month'] = fileID.createVariable(kwargs['MONTHS_NAME'], - 'i', ('time',)) - - # filling netCDF variables - for key,val in output.items(): - nc[key][:] = val.copy() - - # Defining attributes for degree and order - nc['l'].long_name = 'spherical_harmonic_degree'# SH degree long name - nc['l'].units = 'Wavenumber'# SH degree units - nc['m'].long_name = 'spherical_harmonic_order'# SH order long name - nc['m'].units = 'Wavenumber'# SH order units - # Defining attributes for harmonics - nc['clm'].long_name = 'cosine_spherical_harmonics' - nc['clm'].units = kwargs['UNITS'] - nc['slm'].long_name = 'sine_spherical_harmonics' - nc['slm'].units = kwargs['UNITS'] - if kwargs['DATE']: - # Defining attributes for date and month - nc['time'].long_name = kwargs['TIME_LONGNAME'] - nc['time'].units = kwargs['TIME_UNITS'] - nc['month'].long_name = kwargs['MONTHS_LONGNAME'] - nc['month'].units = kwargs['MONTHS_UNITS'] - # global variables of NetCDF file - if kwargs['TITLE']: - fileID.title = kwargs['TITLE'] - if kwargs['REFERENCE']: - fileID.reference = kwargs['REFERENCE'] - # date created - fileID.date_created = time.strftime('%Y-%m-%d',time.localtime()) - - # Output netCDF structure information - logging.info(kwargs['FILENAME']) - logging.info(list(fileID.variables.keys())) - - # Closing the netCDF file - fileID.close() diff --git a/gravity_toolkit/ncdf_write.py b/gravity_toolkit/ncdf_write.py deleted file mode 100755 index c2397188..00000000 --- a/gravity_toolkit/ncdf_write.py +++ /dev/null @@ -1,184 +0,0 @@ -#!/usr/bin/env python -u""" -ncdf_write.py -Written by Tyler Sutterley (11/2021) - -Writes spatial data to COARDS-compliant netCDF4 files - -CALLING SEQUENCE: - ncdf_write(data, lon, lat, tim, FILENAME=output_netcdf4_file) - -INPUTS: - data: z data - lon: longitude array - lat: latitude array - tim: time array - -OPTIONS: - FILENAME: output netCDF4 filename - VARNAME: z variable name in netCDF4 file - LONNAME: longitude variable name in netCDF4 file - LATNAME: latitude variable name in netCDF4 file - UNITS: z variable units - LONGNAME: z variable description - FILL_VALUE: missing value for z variable - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - TITLE: title attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing netCDF4 file - DATE: data has date information - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - netCDF4: Python interface to the netCDF C library - (https://unidata.github.io/netcdf4-python/netCDF4/index.html) - -UPDATE HISTORY: - Updated 11/2021: use remapped dictionary for filling netCDF4 variables - Updated 10/2021: using python logging for handling verbose output - Updated 12/2020: added REFERENCE option to set file attribute - Updated 07/2020: added function docstrings - Updated 04/2020: added option DATE if including time data - Updated 03/2020: only include title if not None - Updated 10/2019: changing Y/N flags to True/False - Updated 09/2019 for public release - Updated 08/2019: don't include time (HH:MM:SS) in creation date - Updated 07/2019: added creation date as a global attribute - Updated 03/2019: print variables keys in list for Python3 compatibility - Updated 03/2018: added option TIMENAME to specify the variable name of time - Updated 02/2017: TIME_LONGNAME and TIME_UNITS with attributes, - updated TIME_LONGNAME to Date_in_Decimal_Years - Updated 07/2016: using netCDF4-python with zlib compression - Updated 06/2016: using __future__ print function - Updated 05/2016: output data types same as input data types - Updated 11/2014: new parameters for variable names and attributes - Updated 05/2014: new parameters for time attributes, and missing values - Updated 02/2014: minor update to if statements - Updated 07/2013: switched from Scientific Python to Scipy - Updated 01/2013: adding time as a variable - Updated 10/2012: changed from variable names x and y to lon and lat. - Written 07/2012 -""" -from __future__ import print_function - -import time -import logging -import numpy as np -import warnings - -# attempt imports -try: - import netCDF4 -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("netCDF4 not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -def ncdf_write(data, lon, lat, tim, **kwargs): - """ - Writes spatial data to COARDS-compliant netCDF4 files - - Parameters - ---------- - data: z data - lon: longitude array - lat: latitude array - tim: time array - - FILENAME: netCDF4 filename - VARNAME: z variable name in netCDF4 file - LONNAME: longitude variable name in netCDF4 file - LATNAME: latitude variable name in netCDF4 file - UNITS: z variable units - LONGNAME: z variable description - FILL_VALUE: missing value for z variable - TIME_UNITS: time variable units - TIME_LONGNAME: time variable description - TITLE: title attribute of dataset - REFERENCE: reference attribute of dataset - CLOBBER: will overwrite an existing netCDF4 file - DATE: data has date information - """ - # set default keyword arguments - kwargs.setdefault('FILENAME',None) - kwargs.setdefault('VARNAME','z') - kwargs.setdefault('LONNAME','lon') - kwargs.setdefault('LATNAME','lat') - kwargs.setdefault('TIMENAME','time') - kwargs.setdefault('UNITS',None) - kwargs.setdefault('LONGNAME',None) - kwargs.setdefault('FILL_VALUE',None) - kwargs.setdefault('TIME_UNITS',None) - kwargs.setdefault('TIME_LONGNAME',None) - kwargs.setdefault('TITLE',None) - kwargs.setdefault('REFERENCE',None) - kwargs.setdefault('DATE',True) - kwargs.setdefault('CLOBBER',True) - # set deprecation warning - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use spatial.to_netCDF4", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # setting NetCDF clobber attribute - clobber = 'w' if kwargs['CLOBBER'] else 'a' - # opening NetCDF file for writing - # Create the NetCDF file - fileID = netCDF4.Dataset(kwargs['FILENAME'], clobber, format="NETCDF4") - - # create output dictionary with key mapping - output = {} - output[kwargs['LONNAME']] = np.copy(lon) - output[kwargs['LATNAME']] = np.copy(lat) - dimensions = [kwargs['LATNAME'],kwargs['LONNAME']] - # extend with date variables - if kwargs['DATE']: - output[kwargs['TIMENAME']] = np.atleast_1d(tim).astype('f') - output[kwargs['VARNAME']] = np.atleast_3d(data) - dimensions.append(kwargs['TIMENAME']) - else: - output[kwargs['VARNAME']] = np.copy(data) - - # defining the NetCDF dimensions and variables - nc = {} - # NetCDF dimensions - for i,dim in enumerate(dimensions): - fileID.createDimension(dim, len(output[dim])) - nc[dim] = fileID.createVariable(dim, output[dim].dtype, (dim,)) - # NetCDF spatial data - for key in [kwargs['VARNAME']]: - nc[key] = fileID.createVariable(key, output[key].dtype, - tuple(dimensions), fill_value=kwargs['FILL_VALUE'], - zlib=True) - # filling NetCDF variables - for key,val in output.items(): - nc[key][:] = val.copy() - - # Defining attributes for longitude and latitude - nc[kwargs['LONNAME']].long_name = 'longitude' - nc[kwargs['LONNAME']].units = 'degrees_east' - nc[kwargs['LATNAME']].long_name = 'latitude' - nc[kwargs['LATNAME']].units = 'degrees_north' - # Defining attributes for dataset - nc[kwargs['VARNAME']].long_name = kwargs['LONGNAME'] - nc[kwargs['VARNAME']].units = kwargs['UNITS'] - # Defining attributes for date if applicable - if kwargs['DATE']: - nc[kwargs['TIMENAME']].long_name = kwargs['TIME_LONGNAME'] - nc[kwargs['TIMENAME']].units = kwargs['TIME_UNITS'] - # global variables of NetCDF file - if kwargs['TITLE']: - fileID.title = kwargs['TITLE'] - if kwargs['REFERENCE']: - fileID.reference = kwargs['REFERENCE'] - # date created - fileID.date_created = time.strftime('%Y-%m-%d',time.localtime()) - - # Output NetCDF structure information - logging.info(kwargs['FILENAME']) - logging.info(list(fileID.variables.keys())) - - # Closing the NetCDF file - fileID.close() diff --git a/gravity_toolkit/piecewise_regress.py b/gravity_toolkit/piecewise_regress.py deleted file mode 100755 index cfbb7ab8..00000000 --- a/gravity_toolkit/piecewise_regress.py +++ /dev/null @@ -1,179 +0,0 @@ -#!/usr/bin/env python -u""" -piecewise_regress.py -Written by Tyler Sutterley (04/2022) - -Fits a synthetic signal to data over a time period by ordinary or weighted - least-squares for breakpoint analysis - -Derivation of Sharp Breakpoint Piecewise Regression: - https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/02-0472 - y = beta_0 + beta_1*t + e (for x <= alpha) - y = beta_0 + beta_1*t + beta_2*(t-alpha) + e (for x > alpha) - -Fit significance derivations are based on Burnham and Anderson (2002) - Model Selection and Multimodel Inference - -CALLING SEQUENCE: - pcwbeta = piecewise_regress(tdec,data,CYCLES=[0.5,1.0],BREAKPOINT=ind) - -INPUTS: - t_in: input time array - d_in: input data array - -OUTPUTS: - beta: regressed coefficients array - error: regression fit error for each coefficient for an input deviation - STDEV: standard deviation of output error - CONF: confidence interval of output error - std_err: standard error for each coefficient - R2: coefficient of determination (r**2). - Proportion of variability accounted by the model - R2Adj: adjusted r**2. adjusts the r**2 for the number of terms in the model - MSE: mean square error - WSSE: Weighted sum of squares error - NRMSE: normalized root mean square error - AIC: Akaike information criterion (Second-Order, AICc) - BIC: Bayesian information criterion (Schwarz criterion) - model: modeled timeseries - simple: modeled timeseries without oscillating components - residual: model residual - DOF: degrees of freedom - N: number of terms used in fit - cov_mat: covariance matrix - -OPTIONS: - BREAK_TIME: breakpoint time for piecewise regression - BREAKPOINT: breakpoint indice of piecewise regression - DATA_ERR: data precision - single value if equal - array if unequal for weighted least squares - WEIGHT: Set if measurement errors for use in weighted least squares - CYCLES: list of cyclical terms (0.5=semi-annual, 1=annual) - STDEV: standard deviation of output error - CONF: confidence interval of output error - AICc: use second order AIC - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - scipy: Scientific Tools for Python (https://docs.scipy.org/doc/) - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 01/2021: added function docstrings - Updated 10/2019: changing Y/N flags to True/False - Updated 01/2019: added option S2 to include 161-day tidal aliasing terms - Updated 12/2018: put transpose of design matrix within FIT_TYPE if statement - Updated 08/2018: import packages before function definition - Updated 09/2017: use rcond=-1 in numpy least-squares algorithms - Updated 08/2015: changed sys.exit to raise ValueError - Updated 11/2014: added simple output for model without climate oscillations - Updated 07/2014: import scipy.stats and scipy.special - Updated 06/2014: changed message to sys.exit - Updated 02/2014: minor update to if statements - Updated 10/2013: updated Log-likelihood (converted from Least-Squares (LS) - log-likelihood to maximum likelihood (ML) log-likelihood) - Added calculation for AICc (corrected for small sample size - Updated 09/2013: updated weighted least-squares and added AIC, BIC and LOGLIK - options for parameter evaluation - Added cases with known standard deviations and weighted least-squares - Minor change: P_cons, P_lin1 and P_lin2 changed to P_x0, P_x1a and P_x1b - Updated 07/2013: updated errors for beta2 - Updated 05/2013: converted to python - Updated 06/2012: added options for MSE and NRMSE. adjusted rsq_adj - Updated 06/2012: changed design matrix creation to 'FIT_TYPE' - added r_square to calcuating goodness of fit compared to others - Updated 03/2012: combined polynomial and harmonic regression functions - Updated 01/2012: added std weighting for a error weighted least-squares - Written 10/2011 -""" -import warnings -import gravity_toolkit.time_series - -def piecewise_regress(*args, **kwargs): - """ - Fits a synthetic signal to data over a time period by ordinary or - weighted least-squares for breakpoint analysis [Toms2003]_ - - Parameters - ---------- - t_in: float - input time array - d_in: float - input data array - BREAK_TIME: float or NoneType, default None - breakpoint time for piecewise regression - BREAKPOINT: int or NoneType, default None - breakpoint indice of piecewise regression - CYCLES: list, default [0.5, 1.0] - list of cyclical terms in fractions of year - DATA_ERR: float or list - data precision - - - single value if equal - - array if unequal for weighted least squares - WEIGHT: bool, default False - Use weighted least squares with measurement errors - STDEV: float, default 0 - Standard deviation of output error - CONF: float, default 0 - Confidence interval of output error - AICc: bool, default False - Use second order AIC for small sample sizes [Burnham2002]_ - - Returns - ------- - beta: float - regressed coefficients array - error: float - regression fit error for each coefficient for an input deviation - - - ``STDEV``: standard deviation of output error - - ``CONF``: confidence interval of output error - std_err: float - standard error for each coefficient - R2: float - coefficient of determination (r\ :sup:`2`) - R2Adj: float - r\ :sup:`2` adjusted for the number of terms in the model - MSE: float - mean square error - WSSE: float - Weighted sum of squares error - NRMSE: float - normalized root mean square error - AIC: float - Akaike information criterion - BIC: float - Bayesian information criterion (Schwarz criterion) - model: float - modeled timeseries - simple: float - modeled timeseries without oscillating components - residual: float - model residual - DOF: int - degrees of freedom - N: int - number of terms used in fit - cov_mat: float - covariance matrix - - References - ---------- - .. [Toms2003] J. D. Toms and M. L. Lesperance, - "Piecewise Regression: A Tool For Identifying Ecological - Thresholds", *Ecology*, 84, 2034-2041, (2003). - `doi: 10.1890/02-0472 `_ - .. [Burnham2002] K. P. Burnham and D. R. Anderson, - *Model Selection and Multimodel Inference*, - 2nd Edition, 488 pp., (2002). - `doi: 10.1007/b97636 `_ - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.time_series instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.time_series.piecewise(*args,**kwargs) diff --git a/gravity_toolkit/plm_colombo.py b/gravity_toolkit/plm_colombo.py deleted file mode 100755 index f0617d35..00000000 --- a/gravity_toolkit/plm_colombo.py +++ /dev/null @@ -1,138 +0,0 @@ -#!/usr/bin/env python -u""" -plm_colombo.py -Written by Tyler Sutterley (04/2022) - -Computes fully-normalized associated Legendre Polynomials - for a vector of x values (can also be singular) -Uses the Colombo (1981) recursion relation - Listed in the Geoid Cookbook and Holmes-Featherstone (2002) - as the most popular recursive algorithm used for computing - fully-normalized Legendre Polynomials in Geodesy -This is a Standard forward column method - -CALLING SEQUENCE: - plm, dplm = plm_colombo(LMAX, np.cos(theta)) - -INPUTS: - LMAX: Upper bound of Spherical Harmonic Degrees - x: elements ranging from -1 to 1 - typically cos(theta), where theta is the colatitude in radians - -OUTPUT: - plms: Legendre polynomials of x (geodesy normalization) - dplms: first differentials of Legendre polynomials of x - -OPTIONS: - ASTYPE: output variable type. Default is np.float64 - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - -REFERENCES: - Geoid Cookbook: http://mitgcm.org/~mlosch/geoidcookbook.pdf - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 09/2020: verify dimensions of input x variable - Updated 08/2020: prevent zero divisions by changing u==0 to eps of data type - Updated 07/2020: added function docstrings - Updated 07/2017: output first differential of legendre polynomials - Updated 09/2013: new format for file headers - Written 03/2013 -""" -import warnings -import numpy as np - -def plm_colombo(LMAX, x, ASTYPE=np.float64): - """ - Computes fully-normalized associated Legendre Polynomials and their - first derivative using a Standard forward column method [Colombo1981]_ - - Parameters - ---------- - LMAX: int - maximum degree of Legrendre polynomials - x: float - elements ranging from -1 to 1 - - Typically ``cos(theta)``, where ``theta`` is the colatitude in radians - ASTYPE: obj, default np.float64 - output variable data type - - Returns - ------- - plms: float - fully-normalized Legendre polynomials - dplms: float - first derivative of Legendre polynomials - - References - ---------- - .. [Colombo1981] O. L. Colombo, - "Numerical Methods for Harmonic Analysis on the Sphere", - Air Force Contract No. F19628-79-C-0027, - *OSURF Proj. No. 711664*, 140 pp., (1981). - .. [Losch2003] M. Losch and V. Seufer, - "How to Compute Geoid Undulations (Geoid Height Relative - to a Given Reference Ellipsoid) from Spherical Harmonic - Coefficients for Satellite Altimetry Applications", (2003). - `eprint ID: 11802 `_ - .. [Holmes2002] S. A. Holmes and W. E. Featherstone, - "A unified approach to the Clenshaw summation and the - recursive computation of very high degree and order - normalised associated Legendre functions", - *Journal of Geodesy*, 76, 279--299, (2002). - `doi: 10.1007/s00190-002-0216-2 `_ - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.associated_legendre instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # removing singleton dimensions of x - x = np.atleast_1d(x).flatten().astype(ASTYPE) - # length of the x array - jm = len(x) - # verify data type of spherical harmonic truncation - LMAX = np.int64(LMAX) - - # allocating for the plm matrix and differentials - plm = np.zeros((LMAX+1,LMAX+1,jm)) - dplm = np.zeros((LMAX+1,LMAX+1,jm)) - - # u is sine of colatitude (cosine of latitude) so that 0 <= s <= 1 - # for x=cos(th): u=sin(th) - u = np.sqrt(1.0 - x**2) - # update where u==0 to eps of data type to prevent invalid divisions - u[u == 0] = np.finfo(u.dtype).eps - - # Calculating the initial polynomials for the recursion - plm[0,0,:] = 1.0 - plm[1,0,:] = np.sqrt(3.0)*x - plm[1,1,:] = np.sqrt(3.0)*u - # calculating first derivatives for harmonics of degree 1 - dplm[1,0,:] = (1.0/u)*(x*plm[1,0,:] - np.sqrt(3)*plm[0,0,:]) - dplm[1,1,:] = (x/u)*plm[1,1,:] - for l in range(2, LMAX+1): - for m in range(0, l):# Zonal and Tesseral harmonics (non-sectorial) - # Computes the non-sectorial terms from previously computed - # sectorial terms. - alm = np.sqrt(((2.0*l-1.0)*(2.0*l+1.0))/((l-m)*(l+m))) - blm = np.sqrt(((2.0*l+1.0)*(l+m-1.0)*(l-m-1.0))/((l-m)*(l+m)*(2.0*l-3.0))) - # if (m == l-1): plm[l-2,m,:] will be 0 - plm[l,m,:] = alm*x*plm[l-1,m,:] - blm*plm[l-2,m,:] - # calculate first derivatives - flm = np.sqrt(((l**2.0 - m**2.0)*(2.0*l + 1.0))/(2.0*l - 1.0)) - dplm[l,m,:] = (1.0/u)*(l*x*plm[l,m,:] - flm*plm[l-1,m,:]) - - # Sectorial harmonics - # The sectorial harmonics serve as seed values for the recursion - # starting with P00 and P11 (outside the loop) - plm[l,l,:] = u*np.sqrt((2.0*l+1.0)/(2.0*l))*np.squeeze(plm[l-1,l-1,:]) - # calculate first derivatives for sectorial harmonics - dplm[l,l,:] = np.longdouble(l)*(x/u)*plm[l,l,:] - - # return the legendre polynomials and their first derivative - return plm, dplm diff --git a/gravity_toolkit/plm_holmes.py b/gravity_toolkit/plm_holmes.py deleted file mode 100755 index b9cc1979..00000000 --- a/gravity_toolkit/plm_holmes.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/env python -u""" -plm_holmes.py -Written by Tyler Sutterley (04/2022) - -Computes fully-normalized associated Legendre Polynomials - for a vector of x values (can also be singular) - -Uses Holmes and Featherstone (2002) recursion relation - -This recursion relation is better conditioned for high - degree and order than the Martin Mohlenkamp relation -It is stable up to very high degree and order (at least 3000). - -This is achieved by individually computing the sectorials to P_m,m -and then iterating up to P_l,m divided by P_m,m and the scale factor 1e280. -Eventually, the result is multiplied again with these to terms. - -CALLING SEQUENCE: - plm, dplm = plm_holmes(LMAX, np.cos(theta)) - -INPUTS: - LMAX: Upper bound of Spherical Harmonic Degrees - x: elements ranging from -1 to 1 - typically cos(theta), where theta is the colatitude in radians - -OUTPUT: - plms: Legendre polynomials of x (geodesy normalization) - dplms: first differentials of Legendre polynomials of x - -OPTIONS: - ASTYPE: output variable type (e.g. np.longdouble). Default is np.float64 - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - -REFERENCES: - S. A. Holmes and W. E. Featherstone, "A unified approach to the Clenshaw - summation and the recursive computation of very high degree and order - normalised associated Legendre functions" Journal of Geodesy, - 76: 279-299, 2002. https://doi.org/10.1007/s00190-002-0216-2 - - Geoid Cookbook: http://mitgcm.org/~mlosch/geoidcookbook.pdf - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 09/2020: verify dimensions of input x variable - Updated 08/2020: prevent zero divisions by changing u==0 to eps of data type - Updated 07/2020: added function docstrings - Updated 10/2018: using future division for python3 Compatibility - Updated 07/2017: output first differential of legendre polynomials - Written 05/2015 -""" -from __future__ import division -import warnings -import numpy as np - -def plm_holmes(LMAX, x, ASTYPE=np.float64): - """ - Computes fully-normalized associated Legendre Polynomials and their - first derivative using Holmes and Featherstone relation [Holmes2002]_ - - Parameters - ---------- - LMAX: int - maximum degree of Legrendre polynomials - x: float - elements ranging from -1 to 1 - - Typically ``cos(theta)``, where ``theta`` is the colatitude in radians - ASTYPE: obj, default np.float64 - output variable data type - - Returns - ------- - plms: float - fully-normalized Legendre polynomials - dplms: float - first derivative of Legendre polynomials - - References - ---------- - .. [Losch2003] M. Losch and V. Seufer, - "How to Compute Geoid Undulations (Geoid Height Relative - to a Given Reference Ellipsoid) from Spherical Harmonic - Coefficients for Satellite Altimetry Applications", (2003). - `eprint ID: 11802 `_ - .. [Holmes2002] S. A. Holmes and W. E. Featherstone, - "A unified approach to the Clenshaw summation and the - recursive computation of very high degree and order - normalised associated Legendre functions", - *Journal of Geodesy*, 76, 279--299, (2002). - `doi: 10.1007/s00190-002-0216-2 `_ - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.associated_legendre instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # removing singleton dimensions of x - x = np.atleast_1d(x).flatten().astype(ASTYPE) - # length of the x array - jm = len(x) - # verify data type of spherical harmonic truncation - LMAX = np.int64(LMAX) - # scaling factor - scalef = 1.0e-280 - - # allocate for multiplicative factors, and plms - f1 = np.zeros(((LMAX+1)*(LMAX+2)//2),dtype=ASTYPE) - f2 = np.zeros(((LMAX+1)*(LMAX+2)//2),dtype=ASTYPE) - p = np.zeros(((LMAX+1)*(LMAX+2)//2,jm),dtype=ASTYPE) - plm = np.zeros((LMAX+1,LMAX+1,jm),dtype=ASTYPE) - dplm = np.zeros((LMAX+1,LMAX+1,jm),dtype=ASTYPE) - - # Precompute multiplicative factors used in recursion relationships - # Note that prefactors are not used for the case when m=l and m=l-1, - # as a different recursion is used for these two values. - k = 2# k = l*(l+1)/2 + m - for l in range(2, LMAX+1): - k += 1 - f1[k] = np.sqrt(2.0*l-1.0)*np.sqrt(2.0*l+1.0)/np.longdouble(l) - f2[k] = np.longdouble(l-1.0)*np.sqrt(2.0*l+1.0)/(np.sqrt(2.0*l-3.0)*np.longdouble(l)) - for m in range(1, l-1): - k += 1 - f1[k] = np.sqrt(2.0*l+1.0)*np.sqrt(2.0*l-1.0)/(np.sqrt(l+m)*np.sqrt(l-m)) - f2[k] = np.sqrt(2.0*l+1.0)*np.sqrt(l-m-1.0)*np.sqrt(l+m-1.0)/ \ - (np.sqrt(2.0*l-3.0)*np.sqrt(l+m)*np.sqrt(l-m)) - k += 2 - - # u is sine of colatitude (cosine of latitude) so that 0 <= s <= 1 - # for x=cos(th): u=sin(th) - u = np.sqrt(1.0 - x**2) - # update where u==0 to eps of data type to prevent invalid divisions - u[u == 0] = np.finfo(u.dtype).eps - - # Calculate P(l,0). These are not scaled. - p[0,:] = 1.0 - p[1,:] = np.sqrt(3.0)*x - k = 1 - for l in range(2, LMAX+1): - k += l - p[k,:] = f1[k]*x*p[k-l,:] - f2[k]*p[k-2*l+1,:] - - # Calculate P(m,m), P(m+1,m), and P(l,m) - pmm = np.sqrt(2.0)*scalef - rescalem = 1.0/scalef - kstart = 0 - - for m in range(1, LMAX): - rescalem = rescalem * u - # Calculate P(m,m) - kstart += m+1 - pmm = pmm * np.sqrt(2*m+1)/np.sqrt(2*m) - p[kstart,:] = pmm - # Calculate P(m+1,m) - k = kstart+m+1 - p[k,:] = x*np.sqrt(2*m+3)*pmm - # Calculate P(l,m) - for l in range(m+2, LMAX+1): - k += l - p[k,:] = x*f1[k]*p[k-l,:] - f2[k]*p[k-2*l+1,:] - p[k-2*l+1,:] = p[k-2*l+1,:] * rescalem - # rescale - p[k,:] = p[k,:] * rescalem - p[k-LMAX,:] = p[k-LMAX,:] * rescalem - - # Calculate P(LMAX,LMAX) - rescalem = rescalem * u - kstart += m+2 - p[kstart,:] = pmm * np.sqrt(2*LMAX+1) / np.sqrt(2*LMAX) * rescalem - # reshape Legendre polynomials to output dimensions - for m in range(LMAX+1): - for l in range(m,LMAX+1): - lm = (l*(l+1))//2 + m - plm[l,m,:] = p[lm,:] - # calculate first derivatives - if (l == m): - dplm[l,m,:] = np.longdouble(m)*(x/u)*plm[l,m,:] - else: - flm = np.sqrt(((l**2.0 - m**2.0)*(2.0*l + 1.0))/(2.0*l - 1.0)) - dplm[l,m,:]= (1.0/u)*(l*x*plm[l,m,:] - flm*plm[l-1,m,:]) - - # return the legendre polynomials and their first derivative - return plm, dplm diff --git a/gravity_toolkit/plm_mohlenkamp.py b/gravity_toolkit/plm_mohlenkamp.py deleted file mode 100755 index f58623fd..00000000 --- a/gravity_toolkit/plm_mohlenkamp.py +++ /dev/null @@ -1,148 +0,0 @@ -#!/usr/bin/env python -u""" -plm_mohlenkamp.py -Written by Tyler Sutterley (04/2022) - -Computes fully-normalized associated Legendre Polynomials - for an array of x values -Uses Martin Mohlenkamp's recursion relation derived from the - Szego (1939) Recurrence formula for Jacobi Polynomials (Pg 71) - -With this algorithm, the associated Legendre Functions are - constructed as an amplitude times a Jacobi Polynomial - P[l,m](cos(theta)) = (sin(theta)^2)*J[l-m,m,m](cos(theta)) - -CALLING SEQUENCE: - plm = plm_mohlenkamp(LMAX, np.cos(theta)) - -INPUTS: - LMAX: Upper bound of Spherical Harmonic Degrees - x: elements ranging from -1 to 1 - typically cos(theta), where theta is the colatitude in radians - -OUTPUT: - plm: Legendre polynomials (geodesy normalization) - -OPTIONS: - MMAX: Upper bound of Spherical Harmonic Orders (default = LMAX) - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - -NOTES: - Modified and updated from IDL plm_x.pro coded by Sean Swenson - Difference from martin function in geoid_mk.mac.f: - plms from plm_mohlenkamp are normalized inside the function - plms from martin are normalized outside the function - For large spherical harmonic degrees this recurrence relation - is poorly conditioned - For spherical harmonic orders above ~1000 can cause overflows - -REFERENCES: - Martin Mohlenkamp, "A User's Guide to Spherical Harmonics" - http://www.ohiouniversityfaculty.com/mohlenka/research/uguide.pdf - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 09/2020: verify dimensions of input x variable - Updated 07/2020: added function docstrings - Updated 05/2015: added parameter MMAX for MMAX != LMAX - Written 09/2013 -""" -import warnings -import numpy as np - -def plm_mohlenkamp(LMAX, x, MMAX=None): - """ - Computes fully-normalized associated Legendre Polynomials - using Martin Mohlenkamp's recursion relation [Mohlenkamp2016]_ - - Derived from [Szego1939]_ recurrence formula for Jacobi Polynomials - - Parameters - ---------- - LMAX: int - maximum degree of Legrendre polynomials - x: float - elements ranging from -1 to 1 - - Typically ``cos(theta)``, where ``theta`` is the colatitude in radians - MMAX: int or NoneType, default None - maximum order of Associated Legrendre polynomials - - Returns - ------- - plms: float - fully-normalized Legendre polynomials - dplms: float - first derivative of Legendre polynomials - - References - ---------- - .. [Mohlenkamp2016] M. J. Mohlenkamp, - "A User's Guide to Spherical Harmonics", (2016). - `[pdf] `_ - .. [Szego1939] Gabor Szeg\ |ouml|\ , "Orthogonal Polynomials", 440 pp., (1939). - `[pdf] `_ - - .. |ouml| unicode:: U+00F6 .. LATIN SMALL LETTER O WITH DIAERESIS - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.associated_legendre instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - - # Verify LMAX as integer - LMAX = np.int64(LMAX) - # upper bound of spherical harmonic orders (default = LMAX) - if MMAX is None: - MMAX = np.copy(LMAX) - - # removing singleton dimensions of x - x = np.atleast_1d(x).flatten() - # length of the x array - sx = len(x) - - # Initialize the output Legendre polynomials - plm=np.zeros((LMAX+1,MMAX+1,sx)) - # Jacobi polynomial for the recurrence relation - jlmm=np.zeros((LMAX+1,MMAX+1,sx)) - # for x=cos(th): rsin= sin(th) - rsin=np.sqrt(1.0 - x**2) - - # for all spherical harmonic orders of interest - for mm in range(0,MMAX+1):# equivalent to 0:MMAX - # Initialize the recurrence relation - # J-1,m,m Term == 0 - # J0,m,m Term - if (mm > 0): - # j ranges from 1 to mm for the product - j = np.arange(0,mm)+1.0 - jlmm[0,mm,:] = np.prod(np.sqrt(1.0 + 1.0/(2.0*j)))/np.sqrt(2.0) - else: # if mm == 0: jlmm = 1/sqrt(2) - jlmm[0,mm,:] = 1.0/np.sqrt(2.0) - # Jk,m,m Terms - for k in range(1, LMAX+1):# computation for SH degrees - # Initialization begins at -1 - # this is to make the formula parallel the function written in - # Martin Mohlenkamp's Guide to Spherical Harmonics - # Jacobi General Terms - if (k == 1):# for degree 1 terms - jlmm[k,mm,:] = 2.0*x * jlmm[k-1,mm,:] * \ - np.sqrt(1.0 + (mm - 0.5)/k) * \ - np.sqrt(1.0 - (mm - 0.5)/(k + 2.0*mm)) - else:# for all other spherical harmonic degrees - jlmm[k,mm,:] = 2.0*x * jlmm[k-1,mm,:] * \ - np.sqrt(1.0 + (mm - 0.5)/k) * \ - np.sqrt(1.0 - (mm - 0.5)/(k + 2.0*mm)) - \ - jlmm[k-2,mm,:] * np.sqrt(1.0 + 4.0/(2.0*k + 2.0*mm - 3.0)) * \ - np.sqrt(1.0 - (1.0/k)) * np.sqrt(1.0 - 1.0/(k + 2.0*mm)) - # Normalization is geodesy convention - for l in range(mm,LMAX+1): # equivalent to mm:LMAX - if (mm == 0):# Geodesy normalization (m=0) == sqrt(2)*sin(th)^0 - # rsin^mm term is dropped as rsin^0 = 1 - plm[l,mm,:] = np.sqrt(2.0)*jlmm[l-mm,mm,:] - else:# Geodesy normalization all others == 2*sin(th)^mm - plm[l,mm,:] = 2.0*(rsin**mm)*jlmm[l-mm,mm,:] - return plm diff --git a/gravity_toolkit/read_ICGEM_harmonics.py b/gravity_toolkit/read_ICGEM_harmonics.py deleted file mode 100644 index adbf1742..00000000 --- a/gravity_toolkit/read_ICGEM_harmonics.py +++ /dev/null @@ -1,61 +0,0 @@ -#!/usr/bin/env python -u""" -read_ICGEM_harmonics.py -Written by Tyler Sutterley (07/2020) - -Read gfc files and extract gravity model spherical harmonics from the GFZ ICGEM - -GFZ International Centre for Global Earth Models (ICGEM) - http://icgem.gfz-potsdam.de/ - -INPUTS: - model_file: GFZ ICGEM gfc spherical harmonic data file - -OPTIONS: - FLAG: string denoting data lines (default gfc) - -OUTPUTS: - clm: cosine spherical harmonics of input data - slm: sine spherical harmonics of input data - eclm: cosine spherical harmonic standard deviations of type errors - eslm: sine spherical harmonic standard deviations of type errors - modelname: name of the gravity model - earth_gravity_constant: GM constant of the Earth for the gravity model - radius: semi-major axis of the Earth for the gravity model - max_degree: maximum degree and order for the gravity model - errors: error type of the gravity model - norm: normalization of the spherical harmonics - tide_system: tide system of gravity model (mean_tide, zero_tide, tide_free) - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - -PROGRAM DEPENDENCIES: - read_ICGEM_harmonics.py Reads the coefficients for a given gravity model file - calculate_tidal_offset.py: calculates the C20 offset for a tidal system - -UPDATE HISTORY: - Updated 03/2021: convert to a simple wrapper function to geoid toolkit - Updated 07/2020: added function docstrings - Updated 07/2017: include parameters to change the tide system - Written 12/2015 -""" -import warnings - -# attempt imports -try: - import geoid_toolkit.read_ICGEM_harmonics -except (ImportError, ModuleNotFoundError) as exc: - warnings.filterwarnings("module") - warnings.warn("geoid_toolkit not available", ImportWarning) -# ignore warnings -warnings.filterwarnings("ignore") - -# PURPOSE: read spherical harmonic coefficients of a gravity model -def read_ICGEM_harmonics(*args,**kwargs): - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use geoid-toolkit instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return geoid_toolkit.read_ICGEM_harmonics(*args,**kwargs) diff --git a/gravity_toolkit/read_SLR_C20.py b/gravity_toolkit/read_SLR_C20.py deleted file mode 100644 index 3b29473a..00000000 --- a/gravity_toolkit/read_SLR_C20.py +++ /dev/null @@ -1,143 +0,0 @@ -#!/usr/bin/env python -u""" -read_SLR_C20.py -Written by Tyler Sutterley (04/2022) - -Reads in C20 spherical harmonic coefficients derived from SLR measurements - -Dataset distributed by NASA PO.DAAC - https://podaac-tools.jpl.nasa.gov/drive/files/GeodeticsGravity/grace/docs - TN-05_C20_SLR.txt - TN-07_C20_SLR.txt - TN-11_C20_SLR.txt - TN-14_C30_C30_GSFC_SLR.txt -Dataset distributed by UTCSR - ftp://ftp.csr.utexas.edu/pub/slr/degree_2/C20_RL05.txt -Datasets distributed by GFZ - ftp://isdcftp.gfz-potsdam.de/grace/Level-2/GFZ/RL06_SLR_C20/ - GFZ_RL06_C20_SLR.dat - ftp://isdcftp.gfz-potsdam.de/grace/GravIS/GFZ/Level-2B/aux_data/ - GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat - -CALLING SEQUENCE: - SLR_C20 = read_SLR_C20(SLR_file) - -INPUTS: - SLR_file: - RL04: TN-05_C20_SLR.txt - RL05: TN-07_C20_SLR.txt - RL06: TN-11_C20_SLR.txt - CSR: C20_RL05.txt - GFZ: GFZ_RL06_C20_SLR.dat - GFZ (combined): GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat - -OUTPUTS: - data: SLR degree 2 order 0 cosine stokes coefficients (C20) - error: SLR degree 2 order 0 cosine stokes coefficient error (eC20) - month: GRACE/GRACE-FO month of measurement (April 2002 = 004) - date: date of SLR measurement - -OPTIONS: - AOD: remove background De-aliasing product from the SLR solution (for CSR) - HEADER: file contains header text to be skipped (default: True) - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - -REFERENCES: - Cheng, M. and Tapley, B. D., "Variations in the Earth's oblateness during - the past 28 years", Journal of Geophysical Research: Solid Earth, - 109(B9), B09402, 2004. 10.1029/2004JB003028 - Loomis, B. D., Rachlin, K. E., and Luthcke, S. B., "Improved Earth - Oblateness Rate Reveals Increased Ice Sheet Losses and Mass-Driven Sea - Level Rise", Geophysical Research Letters, 46(12), 6910-6917, 2019. - https://doi.org/10.1029/2019GL082929 - Koenig, R., Schreiner, P, and Dahle, C. "Monthly estimates of C(2,0) - generated by GFZ from SLR satellites based on GFZ GRACE/GRACE-FO - RL06 background models." V. 1.0. GFZ Data Services, (2019). - https://doi.org/10.5880/GFZ.GRAVIS_06_C20_SLR - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - include utf-8 encoding in reads to be windows compliant - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: added GFZ SLR and GravIS oblateness solutions - define int/float precision to prevent deprecation warning - Updated 02/2021: use adjust_months function to fix special months cases - replaced numpy bool to prevent deprecation warning - Updated 12/2020: using utilities from time module - Updated 08/2020: flake8 compatible regular expression strings - Updated 07/2020: added function docstrings - Updated 08/2019: add catch to verify input SLR file exists - Updated 07/2019: added tilde-expansion of input SLR file - Updated 06/2019: added new GRACE-FO special month (October 2018) - Updated 11/2018: new TN-11 files only list GRACE months available - Updated 06/2016: added option HEADER for files that do not have header text - Updated 05/2016: added option AOD to not remove the AOD correction - Updated 03/2016: minor update to read PO.DAAC - Updated 05/2015: minor change to file determination (only regular expressions) - Updated 02/2015: updated UT/CSR portion and comments - Updated 09/2014: rewrite of the TN-07 read program - using regular expressions and convert_calendar_decimal - Updated 01/2014: updated to use UT/CSR monthly time-series - as an alternative to PO.DAAC as it is updated more regularly - Updated 05/2013: adapted for python - Updated 09/2012: Changed month scheme to output. - Used to remove the GRACE missing months in this program by feeding in the GRACE months - BUT, as the new SLR files start with an earlier date, decided to parallel - the degree-1 read program, and remove the missing months in the read_grace program - Updated 06/2012: OVERHAUL of dating and modification for 'special' GRACE months - Initiated from an incorrect date tag in the SLR data file - New dating will convert from the MJD file into date fraction - Some GRACE 'months' have the accelerometer turned off - for half the month to preserve battery power - These months use half of the prior month in the GRACE global gravity solution - For these months the SLR file has a second dataline for the modified period - Will use these marked (*) data to replace the GRACE C2,0 - ALSO converted the mon and slrdate inputs into options - Updated 01/2012: Updated to feed in SLR file from outside - Will accommodate upcoming GRACE RL05, which will use different SLR files - Written 12/2011 -""" -import warnings -import gravity_toolkit.SLR - -# PURPOSE: read oblateness data from Satellite Laser Ranging (SLR) -def read_SLR_C20(*args, **kwargs): - """ - Reads C20 spherical harmonic coefficients from SLR measurements - - Parameters - ---------- - SLR_file: str - Satellite Laser Ranging file - AOD: bool, default True - Remove background De-aliasing product from the - Center for Space Research (CSR) SLR solutions - HEADER: bool, default True - File contains header text to be skipped - - Returns - ------- - data: float - SLR degree 2 order 0 cosine stokes coefficients - error: float - SLR degree 2 order 0 cosine stokes coefficient error - month: int - GRACE/GRACE-FO month of measurement - date: float - date of SLR measurement - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.SLR instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.SLR.C20(*args,**kwargs) diff --git a/gravity_toolkit/read_SLR_C30.py b/gravity_toolkit/read_SLR_C30.py deleted file mode 100644 index dcd538d9..00000000 --- a/gravity_toolkit/read_SLR_C30.py +++ /dev/null @@ -1,114 +0,0 @@ -#!/usr/bin/env python -u""" -read_SLR_C30.py -Written by Yara Mohajerani and Tyler Sutterley (04/2022) - -Reads monthly degree 3 zonal spherical harmonic data files from SLR - -Dataset distributed by NASA PO.DAAC - https://podaac-tools.jpl.nasa.gov/drive/files/GeodeticsGravity/gracefo/docs - TN-14_C30_C30_GSFC_SLR.txt - ftp://ftp.csr.utexas.edu/pub/slr/degree_5/ - CSR_Monthly_5x5_Gravity_Harmonics.txt -Dataset distributed by GFZ - ftp://isdcftp.gfz-potsdam.de/grace/GravIS/GFZ/Level-2B/aux_data/ - GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat - -CALLING SEQUENCE: - SLR_C30 = read_SLR_C30(SLR_file) - -INPUTS: - SLR_file: - CSR: CSR_Monthly_5x5_Gravity_Harmonics.txt - GFZ: GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat - GSFC: TN-14_C30_C30_GSFC_SLR.txt - LARES: C30_LARES_filtered.txt - -OUTPUTS: - data: SLR degree 3 order 0 cosine stokes coefficients (C30) - error: SLR degree 3 order 0 cosine stokes coefficient error (eC30) - month: GRACE/GRACE-FO month of measurement (April 2002 = 004) - time: date of SLR measurement - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - C30_MEAN: mean C30 to add to LARES C30 anomalies - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - read_SLR_harmonics.py: low-degree spherical harmonic coefficients from SLR - -REFERENCES: - Loomis, Rachlin, and Luthcke, "Improved Earth Oblateness Rate Reveals - Increased Ice Sheet Losses and Mass-Driven Sea Level Rise", - Geophysical Research Letters, 46(12), 6910-6917, (2019). - https://doi.org/10.1029/2019GL082929 - Loomis, Rachlin, Wiese, Landerer, and Luthcke, "Replacing GRACE/GRACE-FO - C30 with satellite laser ranging: Impacts on Antarctic Ice Sheet - mass change". Geophysical Research Letters, 47, (2020). - https://doi.org/10.1029/2019GL085488 - Dahle and Murboeck, "Post-processed GRACE/GRACE-FO Geopotential - GSM Coefficients GFZ RL06 (Level-2B Product)." - V. 0002. GFZ Data Services, (2019). - https://doi.org/10.5880/GFZ.GRAVIS_06_L2B - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - include utf-8 encoding in reads to be windows compliant - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: added GFZ GravIS GRACE/SLR low degree solutions - define int/float precision to prevent deprecation warning - Updated 04/2021: renamed SLR monthly 5x5 function from CSR - Updated 02/2021: use adjust_months function to fix special months cases - Updated 12/2020: using utilities from time module - Updated 08/2020: flake8 compatible regular expression strings - Updated 07/2020: added function docstrings - Updated 08/2019: new GSFC format with more columns - add catch to verify input SLR file exists - added LARES filtered C30 files from John Ries (C30_LARES_filtered.txt) - add C30 mean (9.5717395773300e-07) to LARES solutions - Updated 07/2019: added SLR C3,0 files from PO.DAAC (GSFC) - read CSR monthly 5x5 file and extract C3,0 coefficients - Written 05/2019 -""" -import warnings -import gravity_toolkit.SLR - -# PURPOSE: read Degree 3 zonal data from Satellite Laser Ranging (SLR) -def read_SLR_C30(*args, **kwargs): - """ - Reads C30 spherical harmonic coefficients from SLR measurements - - Parameters - ---------- - SLR_file: str - Satellite Laser Ranging file - C30_MEAN: float, default 9.5717395773300e-07 - Mean C30 to add to LARES C30 anomalies - HEADER: bool, default True - File contains header text to be skipped - - Returns - ------- - data: float - SLR degree 3 order 0 cosine stokes coefficients - error: float - SLR degree 3 order 0 cosine stokes coefficient error - month: int - GRACE/GRACE-FO month of measurement - time: float - date of SLR measurement - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.SLR instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.SLR.C30(*args,**kwargs) diff --git a/gravity_toolkit/read_SLR_C40.py b/gravity_toolkit/read_SLR_C40.py deleted file mode 100644 index 7cb03d62..00000000 --- a/gravity_toolkit/read_SLR_C40.py +++ /dev/null @@ -1,81 +0,0 @@ -#!/usr/bin/env python -u""" -read_SLR_C40.py -Written by Tyler Sutterley (09/2022) - -Reads monthly degree 4 zonal spherical harmonic data files from SLR - -Dataset distributed by CSR - ftp://ftp.csr.utexas.edu/pub/slr/degree_5/ - CSR_Monthly_5x5_Gravity_Harmonics.txt -Dataset distributed by GSFC - https://earth.gsfc.nasa.gov/geo/data/slr - gsfc_slr_5x5c61s61.txt - -CALLING SEQUENCE: - SLR_C40 = read_SLR_C40(SLR_file) - -INPUTS: - SLR_file: - GSFC: gsfc_slr_5x5c61s61.txt - CSR: CSR_Monthly_5x5_Gravity_Harmonics.txt - -OUTPUTS: - data: SLR degree 4 order 0 cosine stokes coefficients (C40) - error: SLR degree 4 order 0 cosine stokes coefficient error (eC40) - month: GRACE/GRACE-FO month of measurement (April 2002 = 004) - time: date of SLR measurement - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - C40_MEAN: mean C40 to add to LARES C40 anomalies - DATE: mid-point of monthly solution for calculating 28-day arc averages - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - read_SLR_harmonics.py: low-degree spherical harmonic coefficients from SLR - -UPDATE HISTORY: - Written 09/2022 -""" -import warnings -import gravity_toolkit.SLR - -# PURPOSE: read Degree 4 zonal data from Satellite Laser Ranging (SLR) -def read_SLR_C40(*args, **kwargs): - """ - Reads C40 spherical harmonic coefficients from SLR measurements - - Parameters - ---------- - SLR_file: str - Satellite Laser Ranging file - C40_MEAN: float, default 0.0 - Mean C40 to add to LARES C40 anomalies - DATE: float or NoneType, default None - Mid-point of monthly solution for calculating 28-day arc averages - - Returns - ------- - data: float - SLR degree 4 order 0 cosine stokes coefficients - error: float - SLR degree 4 order 0 cosine stokes coefficient error - month: int - GRACE/GRACE-FO month of measurement - time: float - date of SLR measurement - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.SLR instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.SLR.C40(*args,**kwargs) diff --git a/gravity_toolkit/read_SLR_C50.py b/gravity_toolkit/read_SLR_C50.py deleted file mode 100644 index cf605234..00000000 --- a/gravity_toolkit/read_SLR_C50.py +++ /dev/null @@ -1,93 +0,0 @@ -#!/usr/bin/env python -u""" -read_SLR_C50.py -Written by Yara Mohajerani and Tyler Sutterley (04/2022) - -Reads monthly degree 5 zonal spherical harmonic data files from SLR - -Dataset distributed by CSR - ftp://ftp.csr.utexas.edu/pub/slr/degree_5/ - CSR_Monthly_5x5_Gravity_Harmonics.txt -Dataset distributed by GSFC - https://earth.gsfc.nasa.gov/geo/data/slr - gsfc_slr_5x5c61s61.txt - -CALLING SEQUENCE: - SLR_C50 = read_SLR_C50(SLR_file) - -INPUTS: - SLR_file: - GSFC: gsfc_slr_5x5c61s61.txt - CSR: CSR_Monthly_5x5_Gravity_Harmonics.txt - -OUTPUTS: - data: SLR degree 5 order 0 cosine stokes coefficients (C50) - error: SLR degree 5 order 0 cosine stokes coefficient error (eC50) - month: GRACE/GRACE-FO month of measurement (April 2002 = 004) - time: date of SLR measurement - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - C50_MEAN: mean C50 to add to LARES C50 anomalies - DATE: mid-point of monthly solution for calculating 28-day arc averages - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - read_SLR_harmonics.py: low-degree spherical harmonic coefficients from SLR - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - include utf-8 encoding in reads to be windows compliant - Updated 12/2021: use function for converting from 7-day arcs - Updated 11/2021: reader for new weekly 5x5+6,1 fields from NASA GSFC - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: simplified program similar to other SLR readers - define int/float precision to prevent deprecation warning - Updated 04/2021: using utilities from time module - Updated 08/2020: flake8 compatible regular expression strings - Updated 07/2020: added function docstrings - Written 11/2019 -""" -import warnings -import gravity_toolkit.SLR - -# PURPOSE: read Degree 5 zonal data from Satellite Laser Ranging (SLR) -def read_SLR_C50(*args, **kwargs): - """ - Reads C50 spherical harmonic coefficients from SLR measurements - - Parameters - ---------- - SLR_file: str - Satellite Laser Ranging file - C50_MEAN: float, default 0.0 - Mean C50 to add to LARES C50 anomalies - DATE: float or NoneType, default None - Mid-point of monthly solution for calculating 28-day arc averages - HEADER: bool, default True - File contains header text to be skipped - - Returns - ------- - data: float - SLR degree 5 order 0 cosine stokes coefficients - error: float - SLR degree 5 order 0 cosine stokes coefficient error - month: int - GRACE/GRACE-FO month of measurement - time: float - date of SLR measurement - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.SLR instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.SLR.C50(*args,**kwargs) diff --git a/gravity_toolkit/read_SLR_CS2.py b/gravity_toolkit/read_SLR_CS2.py deleted file mode 100644 index 77a3e56a..00000000 --- a/gravity_toolkit/read_SLR_CS2.py +++ /dev/null @@ -1,117 +0,0 @@ -#!/usr/bin/env python -u""" -read_SLR_CS2.py -Written by Hugo Lecomte and Tyler Sutterley (04/2022) - -Reads monthly degree 2,m (figure axis and azimuthal dependence) - spherical harmonic data files from satellite laser ranging (SLR) - -Dataset distributed by CSR - http://download.csr.utexas.edu/pub/slr/degree_2/ - C21_S21_RL06.txt or C22_S22_RL06.txt -Dataset distributed by GFZ - ftp://isdcftp.gfz-potsdam.de/grace/GravIS/GFZ/Level-2B/aux_data/ - GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat -Dataset distributed by GSFC - https://earth.gsfc.nasa.gov/geo/data/slr - -CALLING SEQUENCE: - SLR_2m = read_SLR_CS2(SLR_file) - -INPUTS: - SLR_file: - CSR 2,1: C21_S21_RL06.txt - CSR 2,2: C22_S22_RL06.txt - GFZ: GRAVIS-2B_GFZOP_GRACE+SLR_LOW_DEGREES_0002.dat - GSFC: GSFC_C21_S21.txt - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - DATE: mid-point of monthly solution for calculating 28-day arc averages - -OUTPUTS: - C2m: SLR degree 2 order m cosine stokes coefficients - S2m: SLR degree 2 order m sine stokes coefficients - eC2m: SLR degree 2 order m cosine stokes coefficient error - eS2m: SLR degree 2 order m sine stokes coefficient error - month: GRACE/GRACE-FO month of measurement (Apr. 2002 = 004) - time: date of SLR measurement - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - read_SLR_harmonics.py: low-degree spherical harmonic coefficients from SLR - -REFERENCES: - Cheng et al., " Variations of the Earth's figure axis from satellite - laser ranging and GRACE", Journal of Geophysical Research, - 116, B01409, (2011). https://doi.org/10.1029/2010JB000850 - Dahle et al., "The GFZ GRACE RL06 Monthly Gravity Field Time Series: - Processing Details, and Quality Assessment", Remote Sensing, - 11(18), 2116, (2019). https://doi.org/10.3390/rs11182116 - Dahle and Murboeck, "Post-processed GRACE/GRACE-FO Geopotential - GSM Coefficients GFZ RL06 (Level-2B Product)." - V. 0002. GFZ Data Services, (2019). - https://doi.org/10.5880/GFZ.GRAVIS_06_L2B - Chen el al., "Assessment of degree-2 order-1 gravitational changes - from GRACE and GRACE Follow-on, Earth rotation, satellite laser - ranging, and models", Journal of Geodesy, 95(38), (2021). - https://doi.org/10.1007/s00190-021-01492-x - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - include utf-8 encoding in reads to be windows compliant - Updated 11/2021: reader for new weekly 5x5+6,1 fields from NASA GSFC - Updated 09/2021: use functions for converting to and from GRACE months - Updated 08/2021: output empty spherical harmonic errors for GSFC - Updated 06/2021: added GSFC 7-day SLR figure axis solutions - Updated 05/2021: added GFZ GravIS GRACE/SLR low degree solutions - Updated 04/2021: use adjust_months function to fix special months cases - Written 11/2020 -""" -import warnings -import gravity_toolkit.SLR - -# PURPOSE: read Degree 2,m data from Satellite Laser Ranging (SLR) -def read_SLR_CS2(*args, **kwargs): - """ - Reads CS2,m spherical harmonic coefficients from SLR measurements - - Parameters - ---------- - SLR_file: str - Satellite Laser Ranging file - ORDER: int, default 1 - Spherical harmonic order to extract from low-degree fields - DATE: float or NoneType, default None - Mid-point of monthly solution for calculating 28-day arc averages - HEADER: bool, default True - File contains header text to be skipped - - Returns - ------- - C2m: float - SLR degree 2 order m cosine stokes coefficients - S2m: float - SLR degree 2 order m sine stokes coefficients - eC2m: float - SLR degree 2 order m cosine stokes coefficient error - eS2m: float - SLR degree 2 order m sine stokes coefficient error - month: int - GRACE/GRACE-FO month of measurement - time: float - date of SLR measurement - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.SLR instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.SLR.CS2(*args,**kwargs) diff --git a/gravity_toolkit/read_SLR_geocenter.py b/gravity_toolkit/read_SLR_geocenter.py deleted file mode 100644 index 3cbad2d1..00000000 --- a/gravity_toolkit/read_SLR_geocenter.py +++ /dev/null @@ -1,181 +0,0 @@ -#!/usr/bin/env python -u""" -read_SLR_geocenter.py -Written by Tyler Sutterley (04/2022) - -Reads monthly geocenter files from satellite laser ranging provided by CSR - http://download.csr.utexas.edu/pub/slr/geocenter/ - RL04: GCN_RL04.txt - RL05: GCN_RL05.txt -New CF-CM geocenter dataset to reflect the true degree-1 mass variations - http://download.csr.utexas.edu/pub/slr/geocenter/geocenter/README_L1_L2 - http://download.csr.utexas.edu/pub/slr/geocenter/GCN_L1_L2_30d_CF-CM.txt -New geocenter solutions from Minkang Cheng - http://download.csr.utexas.edu/outgoing/cheng/gct2est.220_5s - -CALLING SEQUENCE: - geocenter = read_SLR_geocenter(geocenter_file) - -INPUTS: - geocenter_file: degree 1 file - -OPTIONS: - RADIUS: Earth's radius for calculating spherical harmonics from SLR data - HEADER: rows of data to skip when importing data - COLUMNS: column names of ascii file - time: date in decimal-years - X: X-component of geocenter variation - Y: Y-component of geocenter variation - Z: Z-component of geocenter variation - X_sigma: X-component uncertainty - Y_sigma: Y-component uncertainty - Z_sigma: Z-component uncertainty - -OUTPUTS: - C10: cosine d1/o0 spherical harmonic coefficients - C11: cosine d1/o1 spherical harmonic coefficients - S11: sine d1/o1 spherical harmonic coefficients - eC10: cosine d1/o0 spherical harmonic coefficient error - eC11: cosine d1/o1 spherical harmonic coefficient error - eS11: sine d1/o1 spherical harmonic coefficient error - month: GRACE/GRACE-FO month - time: date of each month in year-decimal - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - geocenter.py: converts between spherical harmonics and geocenter variations - time.py: utilities for calculating time operations - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 11/2021: function deprecated. merged with gravity_toolkit.geocenter - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 04/2021: use file not found exceptions - Updated 02/2021: use adjust_months function to fix special months cases - Updated 12/2020: added option COLUMNS to generalize the ascii data format - replaced numpy loadtxt with generic read using regular expressions - using utilities from time module for operations - Updated 08/2020: flake8 compatible regular expression strings - Updated 07/2020: added function docstrings - Updated 08/2019: add catch to verify input geocenter file exists - Updated 06/2019: added option RADIUS for setting the Earth's radius - Updated 08/2018: using full release string (RL05 instead of 5) - Updated 04/2017: parallels updates to geocenter function INVERSE option - use enumerate to iterate over dates. added option HEADER for headers - Updated 06/2016: using __future__ print function - Updated 05/2016: use geocenter files from 6-hour AOD1b glo Ylms calculated - in aod1b_geocenter.py - Updated 09/2015: added second function for AOD corrected geocenter values - Updated 04/2015: using Julian dates to calculate GRACE/GRACE-FO month - Written 08/2013 -""" -from __future__ import print_function - -import warnings -import gravity_toolkit.geocenter - -# PURPOSE: read geocenter data from Satellite Laser Ranging (SLR) -def read_SLR_geocenter(geocenter_file, RADIUS=None, HEADER=0, - COLUMNS=['time', 'X', 'Y', 'Z', 'X_sigma', 'Y_sigma', 'Z_sigma']): - """ - Reads monthly geocenter files from satellite laser ranging - - Parameters - ---------- - geocenter_file: str - Satellite Laser Ranging file - RADIUS: float or NoneType, default None - Earth's radius for calculating spherical harmonics from SLR data - HEADER: int, default 0 - Rows of data to skip when importing data - COLUMNS: list - Column names of ascii file - - - ``'time'``: date in decimal-years - - ``'X'``: X-component of geocenter variation - - ``'Y'``: Y-component of geocenter variation - - ``'Z'``: Z-component of geocenter variation - - ``'X_sigma'``: X-component uncertainty - - ``'Y_sigma'``: Y-component uncertainty - - ``'Z_sigma'``: Z-component uncertainty - - Returns - ------- - C10: float - cosine d1/o0 spherical harmonic coefficients - C11: float - cosine d1/o1 spherical harmonic coefficients - S11: float - sine d1/o1 spherical harmonic coefficients - eC10: float - cosine d1/o0 spherical harmonic coefficient error - eC11: float - cosine d1/o1 spherical harmonic coefficient error - eS11: float - sine d1/o1 spherical harmonic coefficient error - month: int - GRACE/GRACE-FO month - time: float - date of each month in year-decimal - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.harmonics instead", - DeprecationWarning) - # call renamed version to not break workflows - DEG1 = gravity_toolkit.geocenter(radius=RADIUS).from_SLR(geocenter_file, - AOD=False, header=HEADER, columns=COLUMNS) - # return the SLR-derived geocenter solutions - return DEG1.to_dict() - - -# special function for outputting AOD corrected SLR geocenter values -# need to run aod1b_geocenter.py to calculate the monthly geocenter dealiasing -def aod_corrected_SLR_geocenter(geocenter_file, DREL, RADIUS=None, HEADER=0, - COLUMNS=[]): - """ - Reads monthly geocenter files from satellite laser ranging corrected - for non-tidal ocean and atmospheric variation - - Parameters - ---------- - geocenter_file: Satellite Laser Ranging file - DREL: GRACE/GRACE-FO/Swarm data release - - RADIUS: Earth's radius for calculating spherical harmonics from SLR data - HEADER: rows of data to skip when importing data - COLUMNS: column names of ascii file - time: date in decimal-years - X: X-component of geocenter variation - Y: Y-component of geocenter variation - Z: Z-component of geocenter variation - X_sigma: X-component uncertainty - Y_sigma: Y-component uncertainty - Z_sigma: Z-component uncertainty - - Returns - ------- - C10: cosine d1/o0 spherical harmonic coefficients - C11: cosine d1/o1 spherical harmonic coefficients - S11: sine d1/o1 spherical harmonic coefficients - eC10: cosine d1/o0 spherical harmonic coefficient error - eC11: cosine d1/o1 spherical harmonic coefficient error - eS11: sine d1/o1 spherical harmonic coefficient error - month: GRACE/GRACE-FO month - time: date of each month in year-decimal - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.geocenter instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - DEG1 = gravity_toolkit.geocenter(radius=RADIUS).from_SLR(geocenter_file, - AOD=True, release=DREL, header=HEADER, columns=COLUMNS) - # return the SLR-derived geocenter solutions - return DEG1.to_dict() diff --git a/gravity_toolkit/read_gravis_geocenter.py b/gravity_toolkit/read_gravis_geocenter.py deleted file mode 100644 index ddbe92d7..00000000 --- a/gravity_toolkit/read_gravis_geocenter.py +++ /dev/null @@ -1,90 +0,0 @@ -#!/usr/bin/env python -u""" -read_gravis_geocenter.py -Written by Tyler Sutterley (04/2022) - -Reads monthly geocenter spherical harmonic data files from - GFZ GravIS calculated using GRACE/GRACE-FO measurements - and Ocean Models of degree 1 - -Dataset distributed by GFZ - ftp://isdcftp.gfz-potsdam.de/grace/GravIS/GFZ/Level-2B/aux_data/ - GRAVIS-2B_GFZOP_GEOCENTER_0002.dat - -CALLING SEQUENCE: - geocenter = read_gravis_geocenter(geocenter_file) - -INPUTS: - geocenter_file: degree 1 file - -OUTPUTS: - C10: cosine d1/o0 spherical harmonic coefficients - C11: cosine d1/o1 spherical harmonic coefficients - S11: sine d1/o1 spherical harmonic coefficients - month: GRACE/GRACE-FO month - time: date of each month in year-decimal - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - -REFERENCES: - Dahle and Murboeck, "Post-processed GRACE/GRACE-FO Geopotential - GSM Coefficients GFZ RL06 (Level-2B Product)." - V. 0002. GFZ Data Services, (2019). - https://doi.org/10.5880/GFZ.GRAVIS_06_L2B - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 11/2021: function deprecated. merged with gravity_toolkit.geocenter - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: define int/float precision to prevent deprecation warning - Written 05/2021 -""" -import warnings -import gravity_toolkit.geocenter - -# PURPOSE: read geocenter data from GFZ GravIS SLR/GRACE solutions -def read_gravis_geocenter(geocenter_file, HEADER=True): - """ - Reads monthly geocenter spherical harmonic data files from - GFZ GravIS calculated using GRACE/GRACE-FO measurements - and Ocean Models of degree 1 - - Parameters - ---------- - geocenter_file: str - degree 1 file - HEADER: bool, default True - File contains header text to be skipped - - Returns - ------- - C10: float - cosine d1/o0 spherical harmonic coefficients - C11: float - cosine d1/o1 spherical harmonic coefficients - S11: float - sine d1/o1 spherical harmonic coefficients - month: int - GRACE/GRACE-FO month - time: float - date of each month in year-decimal - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.geocenter instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - DEG1 = gravity_toolkit.geocenter().from_gravis(geocenter_file,header=HEADER) - # return the GFZ GravIS geocenter solutions - return DEG1.to_dict() diff --git a/gravity_toolkit/read_swenson_geocenter.py b/gravity_toolkit/read_swenson_geocenter.py deleted file mode 100755 index 199b88e1..00000000 --- a/gravity_toolkit/read_swenson_geocenter.py +++ /dev/null @@ -1,98 +0,0 @@ -#!/usr/bin/env python -u""" -read_swenson_geocenter.py -Written by Tyler Sutterley (04/2022) - -Reads monthly geocenter coefficients from GRACE measurements and - Ocean Models of Degree 1 provided by Sean Swenson in mm w.e. - -Swenson, S., D. Chambers, and J. Wahr, "Estimating geocenter variations - from a combination of GRACE and ocean model output", J. Geophys. Res., - 113(B08410), 2008. doi:10.1029/2007JB005338 - -CALLING SEQUENCE: - geocenter = read_swenson_geocenter(geocenter_file) - -INPUTS: - geocenter_file: degree 1 file - -OUTPUTS: - C10: cosine d1/o0 spherical harmonic coefficients - C11: cosine d1/o1 spherical harmonic coefficients - S11: sine d1/o1 spherical harmonic coefficients - month: GRACE/GRACE-FO month - time: date of each month in year-decimal - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 11/2021: function deprecated. merged with gravity_toolkit.geocenter - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 04/2021: use file not found exceptions - Updated 02/2021: use adjust_months function to fix special months cases - Updated 12/2020: using utilities from time module - Updated 08/2019: add catch to verify input geocenter file exists - Updated 10/2018: verify integers for python3 compatibility - Updated 03/2017: added catch for HEADER flag - Updated 06/2016: added option HEADER for files that do not have header text - Updated 10/2015: checks if months are included as the last column - and will import. else will calculate GRACE month from the date - Updated regex operator to include pure integers as well - Updated 04/2015: using Julian date to determine GRACE month - Updated 03/2015: minor update to read and regular expression - Updated 10/2014: rewrote with general code updates - using regular expressions to extract data - using better algorithm to find grace month - Written 11/2012 -""" -import warnings -import gravity_toolkit.geocenter - -# PURPOSE: read geocenter data from Sean Swenson -def read_swenson_geocenter(geocenter_file, HEADER=True): - """ - Reads monthly geocenter files computed by Sean Swenson using - GRACE/GRACE-FO measurements and Ocean Models of degree 1 - - Parameters - ---------- - geocenter_file: str - degree 1 file - HEADER: bool, default True - file contains header text to be skipped - - Returns - ------- - C10: float - cosine d1/o0 spherical harmonic coefficients - C11: float - cosine d1/o1 spherical harmonic coefficients - S11: float - sine d1/o1 spherical harmonic coefficients - month: int - GRACE/GRACE-FO month - time: float - date of each month in year-decimal - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.geocenter instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - DEG1 = gravity_toolkit.geocenter().from_swenson(geocenter_file, - header=HEADER) - # return the geocenter solutions from Sean Swenson - return DEG1.to_dict() diff --git a/gravity_toolkit/read_tellus_geocenter.py b/gravity_toolkit/read_tellus_geocenter.py deleted file mode 100644 index c91634f0..00000000 --- a/gravity_toolkit/read_tellus_geocenter.py +++ /dev/null @@ -1,125 +0,0 @@ -#!/usr/bin/env python -u""" -read_tellus_geocenter.py -Written by Tyler Sutterley (04/2022) - -Reads monthly geocenter spherical harmonic data files from GRACE Tellus - Technical Notes (TN-13) calculated using GRACE/GRACE-FO measurements and - Ocean Models of Degree 1 - -Datasets distributed by NASA PO.DAAC -https://podaac-tools.jpl.nasa.gov/drive/files/allData/tellus/L2/degree_1 - -Swenson, S., D. Chambers, and J. Wahr, "Estimating geocenter variations - from a combination of GRACE and ocean model output", J. Geophys. Res., - 113(B08410), 2008. doi:10.1029/2007JB005338 - -Sun, Y., R. Riva, and P. Ditmar, "Observed changes in the Earth's dynamic - oblateness from GRACE data and geophysical models", J. Geodesy., - 90(1), 81-89, 2016. doi:10.1007/s00190-015-0852-y - -Sun, Y., R. Riva, and P. Ditmar, "Optimizing estimates of annual - variations and trends in geocenter motion and J2 from a combination - of GRACE data and geophysical models", J. Geophys. Res. Solid Earth, - 121, 2016. doi:10.1002/2016JB013073 - -CALLING SEQUENCE: - geocenter = read_tellus_geocenter(geocenter_file) - -INPUTS: - geocenter_file: degree 1 file - -OUTPUTS: - C10: cosine d1/o0 spherical harmonic coefficients - C11: cosine d1/o1 spherical harmonic coefficients - S11: sine d1/o1 spherical harmonic coefficients - eC10: cosine d1/o0 spherical harmonic coefficient error - eC11: cosine d1/o1 spherical harmonic coefficient error - eS11: sine d1/o1 spherical harmonic coefficient error - month: GRACE/GRACE-FO month - time: date of each month in year-decimal - -OPTIONS: - HEADER: file contains header text to be skipped (default: True) - JPL: use JPL TN-13 geocenter files with self-attraction and loading - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python - https://numpy.org - https://numpy.org/doc/stable/user/numpy-for-matlab-users.html - dateutil: powerful extensions to datetime - https://dateutil.readthedocs.io/en/stable/ - -PROGRAM DEPENDENCIES: - time.py: utilities for calculating time operations - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 11/2021: function deprecated. merged with gravity_toolkit.geocenter - Updated 09/2021: use functions for converting to and from GRACE months - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 04/2021: use file not found exceptions - Updated 02/2021: use adjust_months function to fix special months cases - Updated 12/2020: using utilities from time module - Updated 08/2020: flake8 compatible regular expression strings - Updated 07/2020: added function docstrings - Updated 08/2019: add catch to verify input geocenter file exists - Updated 07/2019: month adjustments for new TN-13 geocenter files - calculate GRACE/GRACE-FO month based on mean time for JPL TN-13 data files - Updated 06/2019: can use the new JPL TN-13 geocenter files from Tellus - Updated 10/2018: using future division for python3 Compatibility - Updated 06/2016: added option HEADER for files that do not have header text - Updated 04/2015: added time output with convert_calendar_decimal - Updated 03/2015: minor update to read and regular expression - Updated 10/2014: rewrote with general code updates. - using regular expressions to extract data - Updated 03/2013: changed outputs to be C10, C11, S11 instead of C1, S1 -""" -from __future__ import print_function, division - -import warnings -import gravity_toolkit.geocenter - -# PURPOSE: read geocenter data from PO.DAAC -def read_tellus_geocenter(geocenter_file, HEADER=True, JPL=False): - """ - Reads monthly geocenter files computed by JPL Tellus using - GRACE/GRACE-FO measurements and Ocean Models of degree 1 - - Parameters - ---------- - geocenter_file: str - degree 1 file - HEADER: bool, default True - file contains header text to be skipped - JPL: bool, default True - Use JPL TN-13 geocenter files with self-attraction and loading - - Returns - ------- - C10: float - cosine d1/o0 spherical harmonic coefficients - C11: float - cosine d1/o1 spherical harmonic coefficients - S11: float - sine d1/o1 spherical harmonic coefficients - eC10: float - cosine d1/o0 spherical harmonic coefficient error - eC11: float - cosine d1/o1 spherical harmonic coefficient error - eS11: float - sine d1/o1 spherical harmonic coefficient error - month: int - GRACE/GRACE-FO month - time: float - date of each month in year-decimal - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.geocenter instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - DEG1 = gravity_toolkit.geocenter().from_tellus(geocenter_file, - header=HEADER, JPL=JPL) - # return the JPL GRACE Tellus geocenter solutions - return DEG1.to_dict() diff --git a/gravity_toolkit/savitzky_golay.py b/gravity_toolkit/savitzky_golay.py deleted file mode 100644 index 2213f4bf..00000000 --- a/gravity_toolkit/savitzky_golay.py +++ /dev/null @@ -1,109 +0,0 @@ -#!/usr/bin/env python -u""" -savitzky_golay.py -Written by Tyler Sutterley (04/2022) -Adapted from Numerical Recipes, Third Edition - -Smooth and optionally differentiate data of non-uniform sampling - with a Savitzky-Golay filter - -Can preserves the original shape and features of the signal better - than moving averages techniques - -CALLING SEQUENCE: - y_out = savitzky_golay(t_in, y_in, WINDOW=13, ORDER=2) - -INPUTS: - t_in: input time array - y_in: input data array - -OPTIONS: - WINDOW: length of the window (such as 13 for annual). - Must be an odd integer number. - ORDER: order of the polynomial used in the filtering. - Must be less than (window_size - 1) - DERIV: the order of the derivative to compute - default = 0 means only smoothing - RATE: scaling factor for output data and error - DATA_ERR: estimated data error of known and equal value - -OUTPUTS: - data: smoothed time-series (or n-th derivative) - error: estimated error at time points - time: time points for window - -NOTES: - The Savitzky-Golay is a type of low-pass filter, particularly - suited for smoothing noisy data. The main idea behind this - approach is to make for each point a least-square fit with a - polynomial of high order over an odd-sized window centered at - the point. - - A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of - Data by Simplified Least Squares Procedures. Analytical - Chemistry, 1964, 36 (8), pp 1627-1639. - Numerical Recipes 3rd Edition: The Art of Scientific Computing - W.H. Press, S.A. Teukolsky, W. T. Vetterling, B.P. Flannery - Cambridge University Press - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 07/2020: added function docstrings - Updated 08/2019: importing factorial from scipy.special.factorial - Updated 11/2018: using future division for python3 compatibility - Updated 08/2015: changed sys.exit to raise ValueError - Written 06/2014 -""" -import warnings -import gravity_toolkit.time_series - -def savitzky_golay(*args, **kwargs): - """ - Smooth and optionally differentiate data with a Savitzky-Golay - filter [Savitzky1964]_ [Press2007]_ - - Parameters - ---------- - t_in: float - time array - y_in: float - data magnitude array - WINDOW: int or NoneType, default None - Length of the window - - Must be an odd integer - ORDER: int, defualt 2 - Order of the polynomial used in the filtering - - Must be less than (window_size - 1) - DERIV: int, default 0 - Order of the derivative to compute - RATE: float, default 1 - Scaling factor for output data and error - DATA_ERR: float, default 0 - Estimated data error of known and equal value - - Returns - ------- - data: float - Smoothed signal (or n-th derivative) - error: float - Estimated error at time points - time: float - Time points for window - - References - ---------- - .. [Savitzky1964] A. Savitzky, M. J. E. Golay, "Smoothing and - Differentiation of Data by Simplified Least Squares Procedures". - *Analytical Chemistry*, 36(8), 1627--1639, (1964). - .. [Press2007] *Numerical Recipes 3rd Edition: The Art of Scientific - Computing*, W.H. Press, S.A. Teukolsky, W. T. Vetterling, - B.P. Flannery. Cambridge University Press, (2007). - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.time_series instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.time_series.savitzky_golay(*args,**kwargs) diff --git a/gravity_toolkit/tsamplitude.py b/gravity_toolkit/tsamplitude.py deleted file mode 100755 index a304fbd3..00000000 --- a/gravity_toolkit/tsamplitude.py +++ /dev/null @@ -1,56 +0,0 @@ -#!/usr/bin/env python -u""" -tsamplitude.py -Written by Tyler Sutterley (04/2022) - -Calculate the amplitude and phase of a harmonic function from calculated - sine and cosine of a series of measurements - -CALLING SEQUENCE: - ampl,ph = tsamplitude(bsin, bcos) - -INPUTS: - bsin: amplitude of the calculated sine values - bcos: amplitude of the calculated cosine values - -OUTPUTS: - ampl: amplitude from the harmonic functions - ph: phase from the harmonic functions in degrees - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 07/2020: added function docstrings - Updated 10/2019: output both amplitude and phase - Updated 05/2013: converted to python - Written 07/2012: -""" -import warnings -import gravity_toolkit.time_series - -def tsamplitude(*args): - """ - Calculate the amplitude and phase of a harmonic function - - Parameters - ---------- - bsin: float - amplitude of the calculated sine values - bcos: float - amplitude of the calculated cosine values - - Returns - ------- - ampl: float - amplitude from the harmonic functions - ph: float - phase from the harmonic functions in degrees - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.time_series instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.time_series.amplitude(*args) diff --git a/gravity_toolkit/tsregress.py b/gravity_toolkit/tsregress.py deleted file mode 100755 index 0b5b0d0d..00000000 --- a/gravity_toolkit/tsregress.py +++ /dev/null @@ -1,191 +0,0 @@ -#!/usr/bin/env python -u""" -tsregress.py -Written by Tyler Sutterley (04/2022) - -Fits a synthetic signal to data over a time period by ordinary or weighted - least-squares - -Fit significance derivations are based on Burnham and Anderson (2002) - Model Selection and Multimodel Inference - -CALLING SEQUENCE: - tsbeta = tsregress(t_in, d_in, ORDER=1, CYCLES=[0.5,1.0], CONF=0.95) - -INPUTS: - t_in: input time array - d_in: input data array - -OUTPUTS: - beta: regressed coefficients array - error: regression fit error for each coefficient for an input deviation - STDEV: standard deviation of output error - CONF: confidence interval of output error - std_err: standard error for each coefficient - R2: coefficient of determination (r**2). - Proportion of variability accounted by the model - R2Adj: adjusted r**2. adjusts the r**2 for the number of terms in the model - MSE: mean square error - WSSE: Weighted sum of squares error - NRMSE: normalized root mean square error - AIC: Akaike information criterion (Second-Order, AICc) - BIC: Bayesian information criterion (Schwarz criterion) - model: modeled timeseries - simple: modeled timeseries without oscillating components - residual: model residual - DOF: degrees of freedom - N: number of terms used in fit - cov_mat: covariance matrix - -OPTIONS: - DATA_ERR: data precision - single value if equal - array if unequal for weighted least squares - WEIGHT: Set if measurement errors for use in weighted least squares - RELATIVE: relative period - ORDER: maximum polynomial order in fit (0=constant, 1=linear, 2=quadratic) - CYCLES: list of cyclical terms (0.5=semi-annual, 1=annual) - STDEV: standard deviation of output error - CONF: confidence interval of output error - AICc: use second order AIC - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - scipy: Scientific Tools for Python (https://docs.scipy.org/doc/) - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 07/2020: added function docstrings - Updated 10/2019: changing Y/N flags to True/False - Updated 12/2018: put transpose of design matrix within FIT_TYPE if statement - Updated 08/2018: import packages before function definition - Updated 10/2017: output a seasonal model (will be 0 if no oscillating terms) - Updated 09/2017: using rcond=-1 in numpy least-squares algorithms - Updated 03/2017: added a catch for zero error in weighted least-squares - Updated 08/2015: changed sys.exit to raise ValueError - Updated 09/2014: made AICc option for second order AIC - previously was default with no option for standard AIC - Updated 07/2014: output the covariance matrix Hinv - import scipy.stats and scipy.special - Updated 06/2014: changed message to sys.exit - new output for number of terms - Updated 04/2014: added parameter RELATIVE for the relative time - Updated 02/2014: minor update to if statements. output simple regression - Updated 10/2013: - Added calculation for AICc (corrected for small sample size) - Added output DOF (degrees of freedom, nu) - Updated 09/2013: updated weighted least-squares and added AIC, BIC and - LOGLIK options for parameter evaluation - Fixed case with known standard deviations - Changed Weight flag to Y/N - Minor change: P_cons, P_lin and P_quad changed to P_x0, P_x1 and P_x2 - Updated 07/2013: added output for the modelled time-series - Updated 05/2013: converted to Python - Updated 02/2013: added in case for equal data error - and made associated edits.. added option WEIGHT - Updated 10/2012: added in additional FIT_TYPES that do not have a trend - added option for the Schwarz Criterion for model selection - Updated 08/2012: added option for changing the confidence interval or - adding a standard deviation of the error - changed 'std' to data_err for measurement errors - Updated 06/2012: added options for MSE and NRMSE. adjusted rsq_adj - Updated 06/2012: changed design matrix creation to 'FIT_TYPE' - added r_square to calcuating goodness of fit compared to others - Updated 03/2012: combined polynomial and harmonic regression functions - Updated 01/2012: added std weighting for a error weighted least-squares - Written 10/2011 -""" -import warnings -import gravity_toolkit.time_series - -def tsregress(*args, **kwargs): - """ - Fits a synthetic signal to data over a time period by - ordinary or weighted least-squares - - Parameters - ---------- - t_in: float - input time array - d_in: float - input data array - ORDER: int, default 1 - maximum polynomial order in fit - - * ``0``: constant - * ``1``: linear - * ``2``: quadratic - CYCLES: list, default [0.5, 1.0] - list of cyclical terms - DATA_ERR: float or list - Data precision - - - single value if equal - - array if unequal for weighted least squares - WEIGHT: bool, default False - Use weighted least squares with measurement errors - RELATIVE: float or List, default Ellipsis - Epoch for calculating relative dates - - - float: use exact value as epoch - - list: use mean from indices of available times - - ``Ellipsis``: use mean of all available times - STDEV: float, default 0 - Standard deviation of output error - CONF: float, default 0 - Confidence interval of output error - AICc: bool, default False - Use second order AIC for small sample sizes [Burnham2002]_ - - Returns - ------- - beta: float - regressed coefficients array - error: float - regression fit error for each coefficient for an input deviation - - - ``STDEV``: standard deviation of output error - - ``CONF``: confidence interval of output error - std_err: float - standard error for each coefficient - R2: float - coefficient of determination (r\ :sup:`2`) - R2Adj: float - r\ :sup:`2` adjusted for the number of terms in the model - MSE: float - mean square error - WSSE: float - Weighted sum of squares error - NRMSE: float - normalized root mean square error - AIC: float - Akaike information criterion - BIC: float - Bayesian information criterion (Schwarz criterion) - model: float - modeled timeseries - simple: float - modeled timeseries without oscillating components - residual: float - model residual - DOF: int - degrees of freedom - N: int - number of terms used in fit - cov_mat: float - covariance matrix - - Reference - --------- - .. [Burnham2002] K. P. Burnham and D. R. Anderson, - *Model Selection and Multimodel Inference*, - 2nd Edition, 488 pp., (2002). - `doi: 10.1007/b97636 `_ - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.time_series instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.time_series.regress(*args, **kwargs) diff --git a/gravity_toolkit/tssmooth.py b/gravity_toolkit/tssmooth.py deleted file mode 100755 index 4a917a5a..00000000 --- a/gravity_toolkit/tssmooth.py +++ /dev/null @@ -1,146 +0,0 @@ -#!/usr/bin/env python -u""" -tssmooth.py -Written by Tyler Sutterley (04/2022) - -Computes a moving average of a time-series using three possible routines: - 1) centered moving average - 2) 13-month Loess filter (default) - 3) 13-month Loess filter weighted and outputs for all dates - -Note: due to the missing months in the GRACE/GRACE-FO time series, - a standard moving average will have problems if the - missing months are not interpolated. - -CALLING SEQUENCE: - smth = tssmooth(t_in, d_in, HFWTH=6) - -INPUTS: - t_in: input time array - d_in: input data array - -OUTPUTS: - time: time after removing start and end half-windows - data: smoothed time-series - season: seasonal component calculated by the Loess filter - annual: annual component calculated by the Loess filter - semiann: semi-annual component calculated by the Loess filter - trend: instantaneous trend calculated by the Loess filter - error: estimated error of the instantaneous trend - noise: noise component after removing the Loess trend and seasonal components - reduce: original time series after removing start and end half-windows - -OPTIONS: - MOVING: calculates centered moving average using mean of window - mean of: (January up to December) and (February up to January) - WEIGHT: smoothing algorithm that backward models dates before - half-width and forward models dates after half-width - 0: use unweighted Loess filter - 1: use linear weights with Loess filter - 2: use gaussian weights with Loess filter - HFWTH: half-width of the moving average (default = 6 for 13-month Loess) - DATA_ERR: input error for known and equal errors (single value) - STDEV: standard deviation of output error - CONF: confidence interval of output error (default is for 95%) - -PYTHON DEPENDENCIES: - numpy: Scientific Computing Tools For Python (https://numpy.org) - scipy: Scientific Tools for Python (https://docs.scipy.org/doc/) - -UPDATE HISTORY: - Updated 04/2022: updated docstrings to numpy documentation format - Updated 05/2021: define int/float precision to prevent deprecation warning - Updated 07/2020: added function docstrings - Updated 10/2019: changing Y/N flags to True/False. output amplitudes - Updated 09/2019: calculate and output annual and semi-annual phase - Updated 08/2018: use implicit import of scipy stats and special - Updated 09/2017: using rcond=-1 in numpy least-squares algorithms - Updated 09/2014: added output for the reduced time-series - Updated 06/2014: added parameter DATA_ERR for known and equal errors - Updated 03/2014: created a new smoothing algorithm - similar to Loess-type least-squares algorithm but has - backward models dates before HFWTH and forward models dates after - if all dates are available will be centrally weighted - need to update to include error and should find a reference - as i came up with this algorithm at 4:30am - Updated 02/2014: minor update to if statements - Updated 09/2013: switched MOVING flag to Y/N - Minor change: P_cons, and P_lin changed to P_x0, and P_x1 - Updated 06/2013: added error for instantaneous trend - Updated 04/2013: converted to python and added more outputs - Updated 02/2013: using a centered moving average - added seasonal option to compute the smooth seasonal variation - calculated by the loess filter program. - added option to compute the noise component after removing the - smoothed trend and the seasonal component - Updated 03/2012: added Loess smoothing following Velicogna (2009) - Written 12/2011 -""" -import warnings -import gravity_toolkit.time_series - -def tssmooth(*args, **kwargs): - """ - Computes the moving average of a time-series - - 1) centered moving average - 2) 13-month Loess filter [Velicogna2009]_ - 3) 13-month Loess filter weighted and outputs for all dates - - Parameters - ---------- - t_in: float - input time array - d_in: float - input data array - HFWTH: int - half-width of the moving average - MOVING: bool, default False - calculates centered moving average using mean of window - WEIGHT: smoothing algorithm that backward models dates before - half-width and forward models dates after half-width - - - ``0``: use unweighted Loess filter - - ``1``: use linear weights with Loess filter - - ``2``: use gaussian weights with Loess filter - DATA_ERR: float or list - input error for known and equal errors - STDEV: float, default 0 - Standard deviation of output error - CONF: float, default 0 - Confidence interval of output error - - Returns - ------- - time: float - time after removing start and end half-windows - data: float - smoothed time-series - season: float - seasonal component calculated by the Loess filter - annual: float - annual component calculated by the Loess filter - semiann: float - semi-annual component calculated by the Loess filter - trend: float - instantaneous trend calculated by the Loess filter - error: float - estimated error of the instantaneous trend - noise: float - noise component after removing the Loess trend and seasonal components - reduce: float - original time series after removing start and end half-windows - - References - ---------- - .. [Velicogna2009] I. Velicogna, "Increasing rates of ice mass loss - from the Greenland and Antarctic ice sheets revealed by GRACE", - *Geophysical Research Letters*, 36(L19503), - `doi: 10.1029/2009GL040222 `_ - """ - warnings.filterwarnings("module") - warnings.warn("Deprecated. Please use gravity_toolkit.time_series instead", - DeprecationWarning) - warnings.filterwarnings("ignore") - # call renamed version to not break workflows - return gravity_toolkit.time_series.smooth(*args, **kwargs) diff --git a/version.txt b/version.txt index 9084fa2f..26aaba0e 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -1.1.0 +1.2.0