From af3859c20f7530277e1b2ead8a210c3704231ab7 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Tue, 18 Aug 2026 12:13:33 -0700 Subject: [PATCH 1/3] refactor: add ocean and land geocenter functions to class refactor: change assertions to value errors test: add geocenter tests --- doc/source/_assets/gravity-refs.bib | 14 ++ gravity_toolkit/geocenter.py | 143 +++++++++++++----- gravity_toolkit/harmonic_summation.py | 9 +- gravity_toolkit/read_SLR_harmonics.py | 8 +- gravity_toolkit/scripts/calc_degree_one.py | 60 ++------ .../scripts/monte_carlo_degree_one.py | 58 ++----- gravity_toolkit/time.py | 4 +- gravity_toolkit/time_series/lomb_scargle.py | 6 +- test/test_geocenter.py | 62 ++++++++ 9 files changed, 223 insertions(+), 141 deletions(-) create mode 100644 test/test_geocenter.py diff --git a/doc/source/_assets/gravity-refs.bib b/doc/source/_assets/gravity-refs.bib index 5c6714d..28c30d9 100644 --- a/doc/source/_assets/gravity-refs.bib +++ b/doc/source/_assets/gravity-refs.bib @@ -93,6 +93,20 @@ @article{Chao:1987fq pages = {569--596}, publisher = {Oxford University Press (OUP)}, } +@article{Chen:1999ki, +author = {Chen, J L and Wilson, C R and Eanes, R J and Nerem, R S}, +title = {{Geophysical interpretation of observed geocenter variations}}, +journal = {Journal of Geophysical Research: Solid Earth}, +year = {1999}, +month = feb, +volume = {104}, +number = {B2}, +issn = {0148-0227}, +url = {https://doi.org/10.1029/1998jb900019}, +doi = {10.1029/1998JB900019}, +pages = {2683--2690}, +publisher = {American Geophysical Union (AGU)}, +} @article{Cheng:2011hh, author = {Cheng, M and Ries, J C and Tapley, B D}, title = {{Variations of the Earth's figure axis from satellite laser ranging and GRACE}}, diff --git a/gravity_toolkit/geocenter.py b/gravity_toolkit/geocenter.py index 9b119b9..4d00d33 100644 --- a/gravity_toolkit/geocenter.py +++ b/gravity_toolkit/geocenter.py @@ -16,6 +16,8 @@ UPDATE HISTORY: Updated 08/2026: add file logger for reading and writing files + add ocean and land seasonal geocenter functions + turn some functions into classmethods to create geocenter objects Updated 07/2026: add dunder (magic) methods for mathematical operations add HTML representation of geocenter class Updated 06/2024: use wrapper to importlib for optional dependencies @@ -95,26 +97,22 @@ def __init__(self, **kwargs): # WGS84 ellipsoid parameters a_axis = 6378137.0 # [m] semimajor axis of the ellipsoid flat = 1.0 / 298.257223563 # flattening of the ellipsoid - # Mean Earth's Radius in mm having the same volume as WGS84 ellipsoid - kwargs.setdefault( - 'radius', 1000.0 * a_axis * (1.0 - flat) ** (1.0 / 3.0) - ) - # cartesian coordinates - kwargs.setdefault('X', None) - kwargs.setdefault('Y', None) - kwargs.setdefault('Z', None) + # Earth's mean radius having the same volume as WGS84 ellipsoid + radius = 1000.0 * a_axis * (1.0 - flat) ** (1.0 / 3.0) # set default class attributes - self.C10 = None - self.C11 = None - self.S11 = None - self.X = copy.copy(kwargs['X']) - self.Y = copy.copy(kwargs['Y']) - self.Z = copy.copy(kwargs['Z']) - self.time = None - self.month = None - self.filename = None + self.C10 = kwargs.get('C10', None) + self.C11 = kwargs.get('C11', None) + self.S11 = kwargs.get('S11', None) + # cartesian coordinates + self.X = kwargs.get('X', None) + self.Y = kwargs.get('Y', None) + self.Z = kwargs.get('Z', None) + # time and month variables + self.time = kwargs.get('time', None) + self.month = kwargs.get('month', None) + self.filename = kwargs.get('filename', None) # Average Radius of the Earth [mm] - self.radius = copy.copy(kwargs['radius']) + self.radius = kwargs.get('radius', radius) # create logger for class self.logger = logging.getLogger(__name__) # iterator @@ -195,8 +193,10 @@ def from_AOD1B(self, release, year, month, product='glo'): # first column: ISO-formatted date and time cal_date = time.strptime(line_contents[0], r'%Y-%m-%dT%H:%M:%S') # verify that dates are within year and month - assert cal_date.tm_year == year - assert cal_date.tm_mon == month + if cal_date.tm_year != year: + raise ValueError(f'Years mismatch: {cal_date.tm_year} {year}') + if cal_date.tm_mon != month: + raise ValueError(f'Months mismatch: {cal_date.tm_mon} {month}') # second-fourth columns: X, Y and Z geocenter variations temp.X[i], temp.Y[i], temp.Z[i] = np.array( line_contents[1:], dtype='f' @@ -824,6 +824,72 @@ def from_netCDF4(self, geocenter_file, group=None, **kwargs): # return the geocenter harmonics return self.from_dict(DEG1) + @classmethod + def land_seasonal(cls, time): + """ + Model the seasonal component of geocenter motion induced by land water + variations as estimated by :cite:t:`Chen:1999ki` + + Parameters + ---------- + time: np.ndarray + decimal years for geocenter motion + """ + # Annual amplitudes of geocenter components (mm) + AA = dict(X=1.28, Y=0.52, Z=3.30) + # Annual phase of geocenter components (radians) + AP = dict(X=np.radians(44), Y=np.radians(182), Z=np.radians(43)) + # Semi-Annual amplitudes of geocenter components (mm) + SAA = dict(X=0.15, Y=0.56, Z=0.50) + # Semi-Annual phase of geocenter components (radians) + SAP = dict(X=np.radians(331), Y=np.radians(312), Z=np.radians(75)) + # dictionary with geocenter components + XYZ = {} + for key in ['X', 'Y', 'Z']: + # annual and semi-annual components of the geocenter motion + ANN = AA[key] * np.sin(2.0 * np.pi * time + AP[key]) + SEMI = SAA[key] * np.sin(4.0 * np.pi * time + SAP[key]) + XYZ[key] = ANN + SEMI + # create geocenter object from the components + temp = cls(time=time, **XYZ).from_cartesian() + # remove the mean geocenter motion from the seasonal component + temp.mean(apply=True) + # return the seasonal geocenter motion + return temp + + @classmethod + def ocean_seasonal(cls, time): + """ + Model the seasonal component of geocenter motion induced by ocean + variations as estimated by :cite:t:`Chen:1999ki` + + Parameters + ---------- + time: np.ndarray + decimal years for geocenter motion + """ + # Annual amplitudes of geocenter components (mm) + AA = dict(X=0.96, Y=0.97, Z=0.49) + # Annual phase of geocenter components (radians) + AP = dict(X=np.radians(73), Y=np.radians(52), Z=np.radians(3)) + # Semi-Annual amplitudes of geocenter components (mm) + SAA = dict(X=0.86, Y=0.73, Z=0.25) + # Semi-Annual phase of geocenter components (radians) + SAP = dict(X=np.radians(187), Y=np.radians(173), Z=np.radians(232)) + # dictionary with geocenter components + XYZ = {} + for key in ['X', 'Y', 'Z']: + # annual and semi-annual components of the geocenter motion + ANN = AA[key] * np.sin(2.0 * np.pi * time + AP[key]) + SEMI = SAA[key] * np.sin(4.0 * np.pi * time + SAP[key]) + XYZ[key] = ANN + SEMI + # create geocenter object from the components + temp = cls(time=time, **XYZ).from_cartesian() + # remove the mean geocenter motion from the seasonal component + temp.mean(apply=True) + # return the seasonal geocenter motion + return temp + def copy(self, **kwargs): """ Copy a ``geocenter`` object to a new ``geocenter`` object @@ -899,7 +965,8 @@ def from_dict(self, temp, **kwargs): pass return self - def from_harmonics(self, temp, **kwargs): + @classmethod + def from_harmonics(cls, temp, **kwargs): """ Convert a ``harmonics`` object to a ``geocenter`` object @@ -914,26 +981,29 @@ def from_harmonics(self, temp, **kwargs): temp.update_dimensions() # set default keyword arguments kwargs.setdefault('fields', ['time', 'month', 'filename']) - # try to assign variables to self + # try to get attributes from the harmonics object + attributes = {} for key in kwargs['fields']: try: val = getattr(temp, key) - setattr(self, key, np.copy(val)) except AttributeError: pass - # get spherical harmonic objects + else: + attributes[key] = np.copy(val) + # get spherical harmonics if temp.ndim == 2: - self.C10 = np.copy(temp.clm[1, 0]) - self.C11 = np.copy(temp.clm[1, 1]) - self.S11 = np.copy(temp.slm[1, 1]) + C10 = np.copy(temp.clm[1, 0]) + C11 = np.copy(temp.clm[1, 1]) + S11 = np.copy(temp.slm[1, 1]) elif temp.ndim == 3: - self.C10 = np.copy(temp.clm[1, 0, :]) - self.C11 = np.copy(temp.clm[1, 1, :]) - self.S11 = np.copy(temp.slm[1, 1, :]) + C10 = np.copy(temp.clm[1, 0, :]) + C11 = np.copy(temp.clm[1, 1, :]) + S11 = np.copy(temp.slm[1, 1, :]) # return the geocenter object - return self + return cls(C10=C10, C11=C11, S11=S11, **attributes) - def from_matrix(self, clm, slm): + @classmethod + def from_matrix(cls, clm, slm, **kwargs): """ Converts spherical harmonic matrices to a ``geocenter`` object @@ -948,10 +1018,11 @@ def from_matrix(self, clm, slm): clm = np.atleast_3d(clm) slm = np.atleast_3d(slm) # output geocenter object - self.C10 = np.copy(clm[1, 0, :]) - self.C11 = np.copy(clm[1, 1, :]) - self.S11 = np.copy(slm[1, 1, :]) - return self + C10 = np.copy(clm[1, 0, :]) + C11 = np.copy(clm[1, 1, :]) + S11 = np.copy(slm[1, 1, :]) + # return the geocenter object + return cls(C10=C10, C11=C11, S11=S11, **kwargs) def to_dict(self, **kwargs): """ diff --git a/gravity_toolkit/harmonic_summation.py b/gravity_toolkit/harmonic_summation.py index 47c92b8..366d134 100755 --- a/gravity_toolkit/harmonic_summation.py +++ b/gravity_toolkit/harmonic_summation.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ harmonic_summation.py -Written by Tyler Sutterley (07/2026) +Written by Tyler Sutterley (08/2026) Returns the spatial field for a series of spherical harmonics @@ -30,6 +30,7 @@ units.py: class for converting spherical harmonic data to specific units UPDATE HISTORY: + Updated 08/2026: change assertions to value errors Updated 07/2026: use np.einsum for spherical harmonic summations use np.radians to convert from degrees to radians Updated 04/2023: allow love numbers to be None for custom units case @@ -159,8 +160,10 @@ def harmonic_transform( MMAX = np.copy(LMAX) # verify that longitudes cover the complete sphere - assert np.isclose(np.min(lon), 0.0) - assert np.isclose(np.max(lon), 360.0) + if not np.isclose(np.min(lon), 0.0): + raise ValueError(f'Longitude array must start at 0 degrees') + if not np.isclose(np.max(lon), 360.0): + raise ValueError(f'Longitude array must end at 360 degrees') # number of longitudinal points phimax = len(np.squeeze(lon)) # colatitude in radians diff --git a/gravity_toolkit/read_SLR_harmonics.py b/gravity_toolkit/read_SLR_harmonics.py index dc575eb..c1edf5b 100644 --- a/gravity_toolkit/read_SLR_harmonics.py +++ b/gravity_toolkit/read_SLR_harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ read_SLR_harmonics.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Reads in low-degree spherical harmonic coefficients calculated from Satellite Laser Ranging (SLR) measurements @@ -50,6 +50,7 @@ time.py: utilities for calculating time operations UPDATE HISTORY: + Updated 08/2026: change assertions to value errors Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: improve typing for variables in docstrings Updated 11/2022: use f-strings for formatting verbose or ascii output @@ -208,7 +209,10 @@ def read_CSR_monthly_6x1(SLR_file, SCALE=1e-10, HEADER=True): line_contents = file_contents[count].split() # verify arc number from iteration and file IARC = int(line_contents[0]) - assert IARC == (d + 1) + if IARC != (d + 1): + raise ValueError( + f'Arc number mismatch: expected {d + 1}, got {IARC}' + ) # modified Julian date of the middle of the month Ylms['MJD'][d] = np.mean(np.array(line_contents[5:7], dtype=np.float64)) # date of the mid-point of the arc given in years diff --git a/gravity_toolkit/scripts/calc_degree_one.py b/gravity_toolkit/scripts/calc_degree_one.py index df98628..99d2d50 100755 --- a/gravity_toolkit/scripts/calc_degree_one.py +++ b/gravity_toolkit/scripts/calc_degree_one.py @@ -4,7 +4,7 @@ Written by Tyler Sutterley (08/2026) Calculates degree 1 variations using GRACE coefficients of degree 2 and greater, - and ocean bottom pressure variations from ECCO and OMCT/MPIOM + and modeled ocean bottom pressure variations Relation between Geocenter Motion in mm and Normalized Geoid Coefficients X = sqrt(3)*a*C11 @@ -177,6 +177,7 @@ UPDATE HISTORY: Updated 08/2026: use default file logger for valid and failed program runs + moved seasonal land water degree-1 model to geocenter class Updated 07/2026: use np.einsum for spherical harmonic summations can remove sets of harmonic files from the GRACE/GRACE-FO data use np.radians to convert from degrees to radians @@ -377,43 +378,6 @@ def load_AOD(base_dir, PROC, DREL, DSET, START, END, MISSING, LMAX): return gravtk.harmonics().from_dict(grace_Ylms) -# PURPOSE: model the seasonal component of an initial degree 1 model -# using preliminary estimates of annual and semi-annual variations from LWM -# as calculated in Chen et al. (1999), doi:10.1029/1998JB900019 -# NOTE: this is to get an accurate assessment of the land water mass for the -# eustatic component (not for the ocean component from GRACE) -def model_seasonal_geocenter(grace_date): - # Annual amplitudes of (Soil Moisture + Snow) geocenter components (mm) - AAx = 1.28 - AAy = 0.52 - AAz = 3.30 - # Annual phase of (Soil Moisture + Snow) geocenter components (degrees) - APx = 44.0 - APy = 182.0 - APz = 43.0 - # Semi-Annual amplitudes of (Soil Moisture + Snow) geocenter components - SAAx = 0.15 - SAAy = 0.56 - SAAz = 0.50 - # Semi-Annual phase of (Soil Moisture + Snow) geocenter components - SAPx = 331.0 - SAPy = 312.0 - SAPz = 75.0 - # calculate each geocenter component from the amplitude and phase - # converting the phase from degrees to radians - X = AAx * np.sin( - 2.0 * np.pi * grace_date + np.radians(APx) - ) + SAAx * np.sin(4.0 * np.pi * grace_date + np.radians(SAPx)) - Y = AAy * np.sin( - 2.0 * np.pi * grace_date + np.radians(APy) - ) + SAAy * np.sin(4.0 * np.pi * grace_date + np.radians(SAPy)) - Z = AAz * np.sin( - 2.0 * np.pi * grace_date + np.radians(APz) - ) + SAAz * np.sin(4.0 * np.pi * grace_date + np.radians(SAPz)) - DEG1 = gravtk.geocenter(X=X - X.mean(), Y=Y - Y.mean(), Z=Z - Z.mean()) - return DEG1.from_cartesian() - - # PURPOSE: calculate a geocenter time-series def calc_degree_one( base_dir, @@ -637,7 +601,7 @@ def calc_degree_one( GAD_Ylms.mean(apply=True) GAC_Ylms.mean(apply=True) # convert GAC to geocenter object - GAC = gravtk.geocenter().from_harmonics(GAC_Ylms) + GAC = gravtk.geocenter.from_harmonics(GAC_Ylms) # filter GRACE/GRACE-FO coefficients if DESTRIPE: # destriping GRACE GSM and GAD coefficients @@ -666,7 +630,7 @@ def calc_degree_one( GIA_Ylms = GIA_Ylms_rate.drift(GSM_Ylms.time, epoch=2003.3) GIA_Ylms.month[:] = np.copy(GSM_Ylms.month) # save geocenter coefficients of monthly GIA variability - gia = gravtk.geocenter().from_harmonics(GIA_Ylms) + gia = gravtk.geocenter.from_harmonics(GIA_Ylms) # GRACE GAD degree 1 GAD = gravtk.geocenter() @@ -704,7 +668,7 @@ def calc_degree_one( # truncate to degree and order LMAX/MMAX ATM_Ylms = ATM_Ylms.truncate(lmax=LMAX, mmax=MMAX) # save geocenter coefficients of the atmospheric jump corrections - atm = gravtk.geocenter().from_harmonics(ATM_Ylms) + atm = gravtk.geocenter.from_harmonics(ATM_Ylms) # read bottom pressure model if applicable if MODEL not in ('OMCT', 'MPIOM'): @@ -775,7 +739,7 @@ def calc_degree_one( # redistributing the mass over the ocean if specified remove_Ylms.add(Ylms) # save geocenter coefficients of the auxiliary corrections - remove = gravtk.geocenter().from_harmonics(remove_Ylms) + remove = gravtk.geocenter.from_harmonics(remove_Ylms) # Calculating cos/sin of phi arrays # output [m,phi] @@ -841,7 +805,7 @@ def calc_degree_one( # get seasonal variations of an initial geocenter correction # for use in the land water mass calculation - seasonal_geocenter = model_seasonal_geocenter(tdec) + seasonal_geocenter = gravtk.geocenter.land_seasonal(tdec) # iterate solutions: if not single iteration n_iter = 0 @@ -1661,14 +1625,12 @@ def print_global(fid, PROC, DREL, MODEL, AOD, GIA, SLR, S21, month): fid.write(' {0:22}: {1}\n'.format('creator_institution', inst)) # date range and date created calendar_year, calendar_month = gravtk.time.grace_to_calendar(month) - start_time = '{0:4.0f}-{1:02.0f}'.format( - calendar_year[0], calendar_month[0] - ) + start_time = f'{calendar_year[0]:4.0f}-{calendar_month[0]:02.0f}' fid.write(' {0:22}: {1}\n'.format('time_coverage_start', start_time)) - end_time = '{0:4.0f}-{1:02.0f}'.format( - calendar_year[-1], calendar_month[-1] - ) + end_time = f'{calendar_year[-1]:4.0f}-{calendar_month[-1]:02.0f}' fid.write(' {0:22}: {1}\n'.format('time_coverage_end', end_time)) + duration = f'{month[-1] - month[0]:02.0f} months' + fid.write(' {0:22}: {1}\n'.format('time_coverage_duration', duration)) today = time.strftime('%Y-%m-%d', time.localtime()) fid.write(' {0:22}: {1}\n'.format('date_created', today)) fid.write('\n') diff --git a/gravity_toolkit/scripts/monte_carlo_degree_one.py b/gravity_toolkit/scripts/monte_carlo_degree_one.py index 7ea8aac..7ef8be0 100644 --- a/gravity_toolkit/scripts/monte_carlo_degree_one.py +++ b/gravity_toolkit/scripts/monte_carlo_degree_one.py @@ -4,7 +4,7 @@ Written by Tyler Sutterley (08/2026) Calculates degree 1 errors using GRACE coefficients of degree 2 and greater, - and ocean bottom pressure variations from OMCT/MPIOM in a Monte Carlo scheme + and modeled ocean bottom pressure variations in a Monte Carlo scheme Relation between Geocenter Motion in mm and Normalized Geoid Coefficients X = sqrt(3)*a*C11 @@ -158,6 +158,7 @@ UPDATE HISTORY: Updated 08/2026: use default file logger for valid and failed program runs + moved seasonal land water degree-1 model to geocenter class Updated 07/2026: use np.einsum for spherical harmonic summations can remove sets of harmonic files from the GRACE/GRACE-FO data use np.radians to convert from degrees to radians @@ -248,43 +249,6 @@ def info(args): logger.info(f'process id: {os.getpid():d}') -# PURPOSE: model the seasonal component of an initial degree 1 model -# using preliminary estimates of annual and semi-annual variations from LWM -# as calculated in Chen et al. (1999), doi:10.1029/1998JB900019 -# NOTE: this is to get an accurate assessment of the land water mass for the -# eustatic component (not for the ocean component from GRACE) -def model_seasonal_geocenter(grace_date): - # Annual amplitudes of (Soil Moisture + Snow) geocenter components (mm) - AAx = 1.28 - AAy = 0.52 - AAz = 3.30 - # Annual phase of (Soil Moisture + Snow) geocenter components (degrees) - APx = 44.0 - APy = 182.0 - APz = 43.0 - # Semi-Annual amplitudes of (Soil Moisture + Snow) geocenter components - SAAx = 0.15 - SAAy = 0.56 - SAAz = 0.50 - # Semi-Annual phase of (Soil Moisture + Snow) geocenter components - SAPx = 331.0 - SAPy = 312.0 - SAPz = 75.0 - # calculate each geocenter component from the amplitude and phase - # converting the phase from degrees to radians - X = AAx * np.sin( - 2.0 * np.pi * grace_date + np.radians(APx) - ) + SAAx * np.sin(4.0 * np.pi * grace_date + np.radians(SAPx)) - Y = AAy * np.sin( - 2.0 * np.pi * grace_date + np.radians(APy) - ) + SAAy * np.sin(4.0 * np.pi * grace_date + np.radians(SAPy)) - Z = AAz * np.sin( - 2.0 * np.pi * grace_date + np.radians(APz) - ) + SAAz * np.sin(4.0 * np.pi * grace_date + np.radians(SAPz)) - DEG1 = gravtk.geocenter(X=X - X.mean(), Y=Y - Y.mean(), Z=Z - Z.mean()) - return DEG1.from_cartesian() - - # PURPOSE: calculate the satellite error for a geocenter time-series def monte_carlo_degree_one( base_dir, @@ -535,7 +499,7 @@ def monte_carlo_degree_one( GIA_Ylms = GIA_Ylms_rate.drift(GSM_Ylms.time, epoch=2003.3) GIA_Ylms.month[:] = np.copy(GSM_Ylms.month) # save geocenter coefficients of monthly GIA variability - gia = gravtk.geocenter().from_harmonics(GIA_Ylms) + gia = gravtk.geocenter.from_harmonics(GIA_Ylms) # read atmospheric jump corrections from Fagiolini et al. (2015) ATM_Ylms = GSM_Ylms.zeros_like() @@ -550,7 +514,7 @@ def monte_carlo_degree_one( # truncate to degree and order LMAX/MMAX ATM_Ylms = ATM_Ylms.truncate(lmax=LMAX, mmax=MMAX) # save geocenter coefficients of the atmospheric jump corrections - atm = gravtk.geocenter().from_harmonics(ATM_Ylms) + atm = gravtk.geocenter.from_harmonics(ATM_Ylms) # input spherical harmonic datafiles to be removed from the GRACE data # Remove sets of Ylms from the GRACE data before returning @@ -603,7 +567,7 @@ def monte_carlo_degree_one( # redistributing the mass over the ocean if specified remove_Ylms.add(Ylms) # save geocenter coefficients of the auxiliary corrections - remove = gravtk.geocenter().from_harmonics(remove_Ylms) + remove = gravtk.geocenter.from_harmonics(remove_Ylms) # input spherical harmonic datafiles to be used in monte carlo error_Ylms = [] @@ -752,7 +716,7 @@ def monte_carlo_degree_one( # get seasonal variations of an initial geocenter correction # for use in the land water mass calculation - seasonal_geocenter = model_seasonal_geocenter(tdec) + seasonal_geocenter = gravtk.geocenter.land_seasonal(tdec) # degree 1 iterations for each monte carlo run iteration = gravtk.geocenter() @@ -1370,14 +1334,12 @@ def print_global(fid, PROC, DREL, MODEL, GIA, SLR, S21, month): fid.write(' {0:22}: {1}\n'.format('creator_institution', inst)) # date range and date created calendar_year, calendar_month = gravtk.time.grace_to_calendar(month) - start_time = '{0:4.0f}-{1:02.0f}'.format( - calendar_year[0], calendar_month[0] - ) + start_time = f'{calendar_year[0]:4.0f}-{calendar_month[0]:02.0f}' fid.write(' {0:22}: {1}\n'.format('time_coverage_start', start_time)) - end_time = '{0:4.0f}-{1:02.0f}'.format( - calendar_year[-1], calendar_month[-1] - ) + end_time = f'{calendar_year[-1]:4.0f}-{calendar_month[-1]:02.0f}' fid.write(' {0:22}: {1}\n'.format('time_coverage_end', end_time)) + duration = f'{month[-1] - month[0]:02.0f} months' + fid.write(' {0:22}: {1}\n'.format('time_coverage_duration', duration)) today = time.strftime('%Y-%m-%d', time.localtime()) fid.write(' {0:22}: {1}\n'.format('date_created', today)) fid.write('\n') diff --git a/gravity_toolkit/time.py b/gravity_toolkit/time.py index 6eef200..55d2f4e 100644 --- a/gravity_toolkit/time.py +++ b/gravity_toolkit/time.py @@ -15,6 +15,7 @@ UPDATE HISTORY: Updated 08/2026: output numpy.datetime64 objects from file parsers added astype option to calendar_days function (can now be int) + change assertions to value errors Updated 07/2026: added function to convert times to numpy datetimes Updated 06/2024: remove timescale class and leap second calculations Updated 06/2023: add timescale class for converting between time scales @@ -226,7 +227,8 @@ def to_datetime(timedelta: np.ndarray, attributes: str, unit: str = 's'): ``numpy.datetime64`` array """ # verify that units are numpy datetime compatible - assert unit in ['D', 'h', 'm', 's', 'ms', 'us', 'ns'] + if unit not in ['D', 'h', 'm', 's', 'ms', 'us', 'ns']: + raise ValueError(f'Invalid unit: {unit}') # get the epoch and units from the attributes epoch, to_secs = parse_date_string(attributes) # convert epoch to datetime variable diff --git a/gravity_toolkit/time_series/lomb_scargle.py b/gravity_toolkit/time_series/lomb_scargle.py index 03ccc3a..3372c82 100755 --- a/gravity_toolkit/time_series/lomb_scargle.py +++ b/gravity_toolkit/time_series/lomb_scargle.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ lomb_scargle.py -Written by Tyler Sutterley (06/2024) +Written by Tyler Sutterley (08/2026) Wrapper function for computing Lomb-Scargle periodograms @@ -49,6 +49,7 @@ https://doi.org/10.1086/164037 UPDATE HISTORY: + Updated 08/2026: change assertions to value errors Updated 06/2024: add function docstrings Updated 01/2015: added centroid output Written 08/2013 @@ -126,7 +127,8 @@ def lomb_scargle(t_in, d_in, **kwargs): # number of angular frequencies N = int(kwargs['N']) - assert len(OMEGA) == 2, 'Angular frequency range must have 2 values' + if len(OMEGA) != 2: + raise ValueError('Angular frequency range must have 2 values') # array of angular frequencies angular_freq = np.linspace(OMEGA[0], OMEGA[1], N) diff --git a/test/test_geocenter.py b/test/test_geocenter.py new file mode 100644 index 0000000..594b134 --- /dev/null +++ b/test/test_geocenter.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python +""" +test_geocenter.py (08/2026) +""" +import pytest +import numpy as np +import gravity_toolkit as gravtk + +# PURPOSE: model the seasonal component of an initial degree 1 model +# using preliminary estimates of annual and semi-annual variations from LWM +# as calculated in Chen et al. (1999), doi:10.1029/1998JB900019 +# NOTE: this is to get an accurate assessment of the land water mass for the +# eustatic component (not for the ocean component from GRACE) +def test_seasonal_geocenter(): + # create a range of test dates + grace_date = np.arange(2002.25, 2020.25, 1.0/12.0) + # Annual amplitudes of (Soil Moisture + Snow) geocenter components (mm) + AAx = 1.28 + AAy = 0.52 + AAz = 3.30 + # Annual phase of (Soil Moisture + Snow) geocenter components (degrees) + APx = 44.0 + APy = 182.0 + APz = 43.0 + # Semi-Annual amplitudes of (Soil Moisture + Snow) geocenter components + SAAx = 0.15 + SAAy = 0.56 + SAAz = 0.50 + # Semi-Annual phase of (Soil Moisture + Snow) geocenter components + SAPx = 331.0 + SAPy = 312.0 + SAPz = 75.0 + # calculate each geocenter component from the amplitude and phase + # converting the phase from degrees to radians + X = AAx * np.sin( + 2.0 * np.pi * grace_date + np.radians(APx) + ) + SAAx * np.sin(4.0 * np.pi * grace_date + np.radians(SAPx)) + Y = AAy * np.sin( + 2.0 * np.pi * grace_date + np.radians(APy) + ) + SAAy * np.sin(4.0 * np.pi * grace_date + np.radians(SAPy)) + Z = AAz * np.sin( + 2.0 * np.pi * grace_date + np.radians(APz) + ) + SAAz * np.sin(4.0 * np.pi * grace_date + np.radians(SAPz)) + valid = gravtk.geocenter(X=X - X.mean(), Y=Y - Y.mean(), Z=Z - Z.mean()) + valid.from_cartesian() + # calculate using direct function + DEG1 = gravtk.geocenter.land_seasonal(grace_date) + # compare geocenter and degree one components + for key in ['X', 'Y', 'Z', 'C10', 'C11', 'S11']: + assert np.allclose(valid[key], DEG1[key]) + +# PURPOSE: test the from_harmonics class method +def test_from_harmonics(): + rng = np.random.default_rng() + Ylms = gravtk.harmonics(lmax=1).zeros() + Ylms.clm[1, 0] = rng.random() + Ylms.clm[1, 1] = rng.random() + Ylms.slm[1, 1] = rng.random() + geocenter = gravtk.geocenter.from_harmonics(Ylms) + assert geocenter.C10 == Ylms.clm[1, 0] + assert geocenter.C11 == Ylms.clm[1, 1] + assert geocenter.S11 == Ylms.slm[1, 1] From 87465218c03a0cbb9c9945ee99bc45c37f42230d Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Tue, 18 Aug 2026 13:04:39 -0700 Subject: [PATCH 2/3] fix: copilot finds --- gravity_toolkit/geocenter.py | 20 +++++++++++++++- test/test_geocenter.py | 46 ++++++++++++++++++++++++++++++++++-- 2 files changed, 63 insertions(+), 3 deletions(-) diff --git a/gravity_toolkit/geocenter.py b/gravity_toolkit/geocenter.py index 4d00d33..cf1bb48 100644 --- a/gravity_toolkit/geocenter.py +++ b/gravity_toolkit/geocenter.py @@ -977,6 +977,7 @@ def from_harmonics(cls, temp, **kwargs): fields: list default keys in ``harmonics`` object """ + # reshape to matrices if flattened # assign degree and order fields temp.update_dimensions() # set default keyword arguments @@ -991,11 +992,22 @@ def from_harmonics(cls, temp, **kwargs): else: attributes[key] = np.copy(val) # get spherical harmonics - if temp.ndim == 2: + if temp.flattened: + # find the indices of the degree 1 spherical harmonics + # l0 should be 1 and l1 should be lmax + 1 + (l0,) = np.flatnonzero((temp.l == 1) & (temp.m == 0)) + (l1,) = np.flatnonzero((temp.l == 1) & (temp.m == 1)) + # extract from flattened arrays + C10 = temp.clm[l0] + C11 = temp.clm[l1] + S11 = temp.slm[l1] + elif temp.ndim == 2: + # extract from matrix C10 = np.copy(temp.clm[1, 0]) C11 = np.copy(temp.clm[1, 1]) S11 = np.copy(temp.slm[1, 1]) elif temp.ndim == 3: + # extract from matrix with time dimension C10 = np.copy(temp.clm[1, 0, :]) C11 = np.copy(temp.clm[1, 1, :]) S11 = np.copy(temp.slm[1, 1, :]) @@ -1242,12 +1254,18 @@ def mean(self, apply=False, indices=Ellipsis): temp.C10 = np.mean(self.C10[indices]) temp.C11 = np.mean(self.C11[indices]) temp.S11 = np.mean(self.S11[indices]) + temp.to_cartesian() # calculating the time-variable gravity field by removing # the static component of the gravitational field if apply: + # remove the mean degree-one components self.C10 -= temp.C10 self.C11 -= temp.C11 self.S11 -= temp.S11 + # remove the mean geocenter motion + self.X -= temp.X + self.Y -= temp.Y + self.Z -= temp.Z # calculate mean of temporal variables for key in ['time', 'month']: try: diff --git a/test/test_geocenter.py b/test/test_geocenter.py index 594b134..3faae9c 100644 --- a/test/test_geocenter.py +++ b/test/test_geocenter.py @@ -2,18 +2,20 @@ """ test_geocenter.py (08/2026) """ + import pytest import numpy as np import gravity_toolkit as gravtk + # PURPOSE: model the seasonal component of an initial degree 1 model # using preliminary estimates of annual and semi-annual variations from LWM # as calculated in Chen et al. (1999), doi:10.1029/1998JB900019 # NOTE: this is to get an accurate assessment of the land water mass for the # eustatic component (not for the ocean component from GRACE) -def test_seasonal_geocenter(): +def test_land_seasonal(): # create a range of test dates - grace_date = np.arange(2002.25, 2020.25, 1.0/12.0) + grace_date = np.arange(2002.25, 2020.25, 1.0 / 12.0) # Annual amplitudes of (Soil Moisture + Snow) geocenter components (mm) AAx = 1.28 AAy = 0.52 @@ -49,6 +51,46 @@ def test_seasonal_geocenter(): for key in ['X', 'Y', 'Z', 'C10', 'C11', 'S11']: assert np.allclose(valid[key], DEG1[key]) + +def test_ocean_seasonal(): + # create a range of test dates + grace_date = np.arange(2002.25, 2020.25, 1.0 / 12.0) + # Annual amplitudes of ocean (TOPEX) geocenter components (mm) + AAx = 0.96 + AAy = 0.97 + AAz = 0.49 + # Annual phase of ocean (TOPEX) geocenter components (degrees) + APx = 73.0 + APy = 52.0 + APz = 3.0 + # Semi-Annual amplitudes of ocean (TOPEX) geocenter components + SAAx = 0.86 + SAAy = 0.73 + SAAz = 0.25 + # Semi-Annual phase of ocean (TOPEX) geocenter components + SAPx = 187.0 + SAPy = 173.0 + SAPz = 232.0 + # calculate each geocenter component from the amplitude and phase + # converting the phase from degrees to radians + X = AAx * np.sin( + 2.0 * np.pi * grace_date + np.radians(APx) + ) + SAAx * np.sin(4.0 * np.pi * grace_date + np.radians(SAPx)) + Y = AAy * np.sin( + 2.0 * np.pi * grace_date + np.radians(APy) + ) + SAAy * np.sin(4.0 * np.pi * grace_date + np.radians(SAPy)) + Z = AAz * np.sin( + 2.0 * np.pi * grace_date + np.radians(APz) + ) + SAAz * np.sin(4.0 * np.pi * grace_date + np.radians(SAPz)) + valid = gravtk.geocenter(X=X - X.mean(), Y=Y - Y.mean(), Z=Z - Z.mean()) + valid.from_cartesian() + # calculate using direct function + DEG1 = gravtk.geocenter.ocean_seasonal(grace_date) + # compare geocenter and degree one components + for key in ['X', 'Y', 'Z', 'C10', 'C11', 'S11']: + assert np.allclose(valid[key], DEG1[key]) + + # PURPOSE: test the from_harmonics class method def test_from_harmonics(): rng = np.random.default_rng() From b812e3c6c15a1161001d0492108ba16c2c088319 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Tue, 18 Aug 2026 13:06:28 -0700 Subject: [PATCH 3/3] Update geocenter.py --- gravity_toolkit/geocenter.py | 1 + 1 file changed, 1 insertion(+) diff --git a/gravity_toolkit/geocenter.py b/gravity_toolkit/geocenter.py index cf1bb48..9992d5e 100644 --- a/gravity_toolkit/geocenter.py +++ b/gravity_toolkit/geocenter.py @@ -1262,6 +1262,7 @@ def mean(self, apply=False, indices=Ellipsis): self.C10 -= temp.C10 self.C11 -= temp.C11 self.S11 -= temp.S11 + if apply and np.any(self.X) and np.any(self.Y) and np.any(self.Z): # remove the mean geocenter motion self.X -= temp.X self.Y -= temp.Y