diff --git a/.github/workflows/python-request.yml b/.github/workflows/python-request.yml index fd0c89d6..cecb5289 100644 --- a/.github/workflows/python-request.yml +++ b/.github/workflows/python-request.yml @@ -32,7 +32,7 @@ jobs: with: lfs: true - name: Set up pixi environment - uses: prefix-dev/setup-pixi@v0.9.1 + uses: prefix-dev/setup-pixi@v0.9.6 - name: Lint with flake8 run: | # stop the build if there are Python syntax errors or undefined names diff --git a/access/itsg_graz_grace_sync.py b/access/itsg_graz_grace_sync.py index c765a4cd..fcf0d1be 100755 --- a/access/itsg_graz_grace_sync.py +++ b/access/itsg_graz_grace_sync.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" itsg_graz_grace_sync.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (07/2026) Syncs GRACE/GRACE-FO and auxiliary data from the ITSG GRAZ server CALLING SEQUENCE: @@ -39,6 +39,7 @@ utilities.py: download and management utilities for syncing files UPDATE HISTORY: + Updated 07/2026: ITSG GRACE server moved from outgoing to pub Updated 05/2023: use pathlib to define and operate on paths Updated 12/2022: single implicit import of gravity toolkit Updated 11/2022: use f-strings for formatting verbose or ascii output @@ -82,7 +83,7 @@ def itsg_graz_grace_sync(DIRECTORY, RELEASE=None, LMAX=None, TIMEOUT=0, logging.basicConfig(level=logging.INFO) # ITSG GRAZ server - HOST = ['http://ftp.tugraz.at','outgoing','ITSG','GRACE'] + HOST = ['http://ftp.tugraz.at','pub','ITSG','GRACE'] # open connection with ITSG GRAZ server at remote directory release_directory = f'ITSG-{RELEASE}' # regular expression operators for ITSG data and models diff --git a/dealiasing/dealiasing_global_uplift.py b/dealiasing/dealiasing_global_uplift.py index fd6fa239..a139c3d8 100644 --- a/dealiasing/dealiasing_global_uplift.py +++ b/dealiasing/dealiasing_global_uplift.py @@ -239,7 +239,7 @@ def dealiasing_global_uplift(base_dir, attributes['time']['standard_name'] = 'time' # Computing plms for converting to spatial domain - theta = (90.0 - grid.lat)*np.pi/180.0 + theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # for each tar file diff --git a/doc/source/_assets/gravity-refs.bib b/doc/source/_assets/gravity-refs.bib index cbbf4cce..156723bf 100644 --- a/doc/source/_assets/gravity-refs.bib +++ b/doc/source/_assets/gravity-refs.bib @@ -186,6 +186,20 @@ @book{Dershowitz:2007cc publisher = {Cambridge University Press}, edition = {3}, } +@article{Driscoll:1994bp, +author = {Driscoll, J R and Healy, D M}, +title = {{Computing Fourier Transforms and Convolutions on the 2-Sphere}}, +journal = {Advances in Applied Mathematics}, +year = {1994}, +month = jun, +volume = {15}, +number = {2}, +issn = {0196-8858}, +url = {https://doi.org/10.1006/aama.1994.1008}, +doi = {10.1006/aama.1994.1008}, +pages = {202--250}, +publisher = {Elsevier BV}, +} @article{Dziewonski:1981bz, author = {Dziewonski, A M and Anderson, D L}, title = {{Preliminary reference Earth model}}, @@ -250,6 +264,20 @@ @misc{Gegout:2010gc url = {https://doi.org/10.13140/RG.2.1.1866.7045}, doi = {10.13140/RG.2.1.1866.7045}, } +@article{Gruber:2016hn, +author = {Gruber, C and Abrykosov, O}, +title = {{On computation and use of Fourier coefficients for associated Legendre functions}}, +journal = {Journal of Geodesy}, +year = {2016}, +month = jun, +volume = {90}, +number = {6}, +issn = {0949-7714}, +url = {https://doi.org/10.1007/s00190-016-0891-z}, +doi = {10.1007/s00190-016-0891-z}, +pages = {525--535}, +publisher = {Springer Science and Business Media LLC}, +} @inbook{Han:1989kj, author = {Han, D and Wahr, J}, title = {{Post-Glacial Rebound Analysis for a Rotating Earth}}, @@ -300,6 +328,16 @@ @book{HofmannWellenhof:2006hy doi = {10.1007/978-3-211-33545-1}, publisher = {Springer Vienna}, } +@article{Hofsommer:1960wg, +author = {Hofsommer, D J and Potters, M L}, +title = {{Table of Fourier coefficients of associated Legendre functions}}, +journal = {Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen: Series A: Mathematical Sciences}, +year = {1960}, +volume = {63}, +number = {5}, +issn = {0023-3358}, +pages = {460--480}, +} @article{Holmes:2002ff, author = {Holmes, S A and Featherstone, W E}, title = {{A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions}}, diff --git a/doc/source/notebooks/GRACE-Geostrophic-Maps.ipynb b/doc/source/notebooks/GRACE-Geostrophic-Maps.ipynb index a147a48d..ca4ce4c2 100644 --- a/doc/source/notebooks/GRACE-Geostrophic-Maps.ipynb +++ b/doc/source/notebooks/GRACE-Geostrophic-Maps.ipynb @@ -310,6 +310,7 @@ "widgets.select_corrections()\n", "widgets.select_output()\n", "# display widgets for setting GRACE/GRACE-FO corrections parameters\n", + "widgets.gaussian.value = 600.0\n", "ipywidgets.VBox([\n", " widgets.GIA_file,\n", " widgets.GIA,\n", @@ -339,7 +340,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -363,10 +364,12 @@ "# grid latitude and longitude\n", "grid.lon = np.copy(landsea.lon)\n", "grid.lat = np.copy(landsea.lat)\n", + "# mask equatorial regions due to hydrostrophic inaccuracies\n", + "valid, = np.nonzero((np.abs(grid.lat) > 10))\n", "\n", "# Computing plms for converting to spatial domain\n", - "theta = (90.0 - grid.lat)*np.pi/180.0\n", - "PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta))\n", + "theta = np.radians(90.0 - grid.lat)\n", + "PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta[valid]))\n", "RAD = widgets.gaussian.value\n", "\n", "# read load love numbers file\n", @@ -443,8 +446,6 @@ "# output geostrophic current grid\n", "grid.data = np.zeros((nlat, nlon, 2,nt))\n", "grid.mask = np.ones((nlat, nlon, 2,nt), dtype=bool)\n", - "# mask equatorial regions due to hydrostrophic inaccuracies\n", - "valid, = np.nonzero((np.abs(grid.lat) > 10))\n", "grid.mask[valid,:,:,:] = False\n", "# set land values from land-sea mask to invalid\n", "indy,indx = np.nonzero(np.logical_not(landsea.mask))\n", @@ -575,7 +576,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.8.10 64-bit", + "display_name": "py13", "language": "python", "name": "python3" }, @@ -589,12 +590,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" - }, - "vscode": { - "interpreter": { - "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" - } + "version": "3.13.0" } }, "nbformat": 4, diff --git a/doc/source/notebooks/GRACE-Spatial-Error.ipynb b/doc/source/notebooks/GRACE-Spatial-Error.ipynb index 1b765a04..b2cf733e 100644 --- a/doc/source/notebooks/GRACE-Spatial-Error.ipynb +++ b/doc/source/notebooks/GRACE-Spatial-Error.ipynb @@ -288,13 +288,13 @@ " nlat = len(grid.lat)\n", "\n", "# Computing plms for converting to spatial domain\n", - "theta = (90.0 - grid.lat)*np.pi/180.0\n", + "theta = np.radians(90.0 - grid.lat)\n", "PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta))\n", "# square of legendre polynomials truncated to order MMAX\n", "mm = np.arange(0,MMAX+1)\n", "PLM2 = PLM[:,mm,:]**2\n", "# Calculating cos(m*phi)^2 and sin(m*phi)^2\n", - "phi = grid.lon[np.newaxis,:]*np.pi/180.0\n", + "phi = np.radians(grid.lon[np.newaxis,:])\n", "ccos = np.cos(np.dot(mm[:,np.newaxis],phi))**2\n", "ssin = np.sin(np.dot(mm[:,np.newaxis],phi))**2\n", " \n", diff --git a/doc/source/notebooks/GRACE-Spatial-Maps.ipynb b/doc/source/notebooks/GRACE-Spatial-Maps.ipynb index 2b1b908f..c8802c90 100644 --- a/doc/source/notebooks/GRACE-Spatial-Maps.ipynb +++ b/doc/source/notebooks/GRACE-Spatial-Maps.ipynb @@ -405,7 +405,7 @@ " nlat = len(grid.lat)\n", "\n", "# Computing plms for converting to spatial domain\n", - "theta = (90.0-grid.lat)*np.pi/180.0\n", + "theta = np.radians(90.0-grid.lat)\n", "PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta))\n", "\n", "# read load love numbers file\n", diff --git a/geocenter/calc_degree_one.py b/geocenter/calc_degree_one.py index bf0f5386..beeb9c22 100755 --- a/geocenter/calc_degree_one.py +++ b/geocenter/calc_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" calc_degree_one.py -Written by Tyler Sutterley (01/2025) +Written by Tyler Sutterley (07/2026) Calculates degree 1 variations using GRACE coefficients of degree 2 and greater, and ocean bottom pressure variations from ECCO and OMCT/MPIOM @@ -100,6 +100,15 @@ --ocean-file X: Index file for ocean model harmonics --mean-file X: GRACE/GRACE-FO mean file to remove from the harmonic data --mean-format X: Input data format for GRACE/GRACE-FO mean file + --remove-file X: Monthly files to be removed from the GRACE/GRACE-FO data + --remove-format X: Input data format for files to be removed + ascii + netCDF4 + HDF5 + index-ascii + index-netCDF4 + index-HDF5 + --redistribute-removed: redistribute removed mass fields over the ocean --iterative: Iterate degree one solutions -s X, --solver X: Least squares solver for degree one solutions inv: matrix inversion @@ -167,6 +176,9 @@ https://doi.org/10.1029/2007JB005338 UPDATE HISTORY: + 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 Updated 01/2025: fixed deprecated tick label resizing Updated 06/2024: use wrapper to importlib for optional dependencies Updated 10/2023: generalize mission variable to be GRACE/GRACE-FO @@ -340,12 +352,12 @@ def model_seasonal_geocenter(grace_date): 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 + APx*np.pi/180.0) + \ - SAAx*np.sin(4.0*np.pi*grace_date + SAPx*np.pi/180.0) - Y = AAy*np.sin(2.0*np.pi*grace_date + APy*np.pi/180.0) + \ - SAAy*np.sin(4.0*np.pi*grace_date + SAPy*np.pi/180.0) - Z = AAz*np.sin(2.0*np.pi*grace_date + APz*np.pi/180.0) + \ - SAAz*np.sin(4.0*np.pi*grace_date + SAPz*np.pi/180.0) + 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() @@ -371,6 +383,9 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, DATAFORM=None, MEAN_FILE=None, MEANFORM=None, + REMOVE_FILES=None, + REMOVE_FORMAT=None, + REDISTRIBUTE_REMOVED=False, MODEL_INDEX=None, ITERATIVE=False, SOLVER=None, @@ -484,13 +499,13 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # spatial parameters in radians - dphi = dlon*np.pi/180.0 - dth = dlat*np.pi/180.0 + dphi = np.radians(dlon) + dth = np.radians(dlat) # longitude and colatitude in radians - phi = landsea.lon[np.newaxis,:]*np.pi/180.0 - th = (90.0 - np.squeeze(landsea.lat))*np.pi/180.0 + phi = np.radians(np.squeeze(landsea.lon)) + th = np.radians(90.0 - np.squeeze(landsea.lat)) # create land function - land_function = np.zeros((nlon, nlat),dtype=np.float64) + land_function = np.zeros((nlon, nlat), dtype=np.float64) # extract land function from file # combine land and island levels for land function indx,indy = np.nonzero((landsea.data.T >= 1) & (landsea.data.T <= 3)) @@ -575,17 +590,16 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, GAD = gravtk.geocenter() GAD.time = np.copy(GAD_Ylms.time) GAD.month = np.copy(GAD_Ylms.month) - GAD.C10 = np.zeros((n_files)) - GAD.C11 = np.zeros((n_files)) - GAD.S11 = np.zeros((n_files)) - for t in range(0,n_files): - # converting GAD degree 1 harmonics to mass - # NOTE: following Swenson (2008): do not use the kl Load Love number - # to convert the GAD coefficients into coefficients of mass as - # the GAC and GAD products are computed with a Load Love number of 0 - GAD.C10[t] = rho_e*rad_e*np.squeeze(GAD_Ylms.clm[1,0,t])*(2.0 + 1.0)/3.0 - GAD.C11[t] = rho_e*rad_e*np.squeeze(GAD_Ylms.clm[1,1,t])*(2.0 + 1.0)/3.0 - GAD.S11[t] = rho_e*rad_e*np.squeeze(GAD_Ylms.slm[1,1,t])*(2.0 + 1.0)/3.0 + GAD.C10 = np.empty((n_files)) + GAD.C11 = np.empty((n_files)) + GAD.S11 = np.empty((n_files)) + # converting GAD degree 1 harmonics to mass + # NOTE: following Swenson (2008): do not use the kl Load Love number + # to convert the GAD coefficients into coefficients of mass as + # the GAC and GAD products are computed with a Load Love number of 0 + GAD.C10[:] = rho_e*rad_e*np.squeeze(GAD_Ylms.clm[1,0,:])*(2.0 + 1.0)/3.0 + GAD.C11[:] = rho_e*rad_e*np.squeeze(GAD_Ylms.clm[1,1,:])*(2.0 + 1.0)/3.0 + GAD.S11[:] = rho_e*rad_e*np.squeeze(GAD_Ylms.slm[1,1,:])*(2.0 + 1.0)/3.0 # removing the mean of the GAD OBP coefficients GAD.mean(apply=True) @@ -622,14 +636,66 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, # converting ecco degree 1 harmonics to coefficients of mass OBP = gravtk.geocenter.from_harmonics(OBP_Ylms).scale(dfactor[1]) + # input spherical harmonic datafiles to be removed from the GRACE data + # Remove sets of Ylms from the GRACE data before returning + remove_Ylms = GSM_Ylms.zeros_like() + remove_Ylms.time[:] = np.copy(GSM_Ylms.time) + remove_Ylms.month[:] = np.copy(GSM_Ylms.month) + if REMOVE_FILES: + # extend list if a single format was entered for all files + if len(REMOVE_FORMAT) < len(REMOVE_FILES): + REMOVE_FORMAT = REMOVE_FORMAT*len(REMOVE_FILES) + # for each file to be removed + for REMOVE_FILE,REMOVEFORM in zip(REMOVE_FILES,REMOVE_FORMAT): + if REMOVEFORM in ('ascii','netCDF4','HDF5'): + # ascii (.txt) + # netCDF4 (.nc) + # HDF5 (.H5) + Ylms = gravtk.harmonics().from_file(REMOVE_FILE, + format=REMOVEFORM) + attributes['lineage'].append(Ylms.filename) + elif REMOVEFORM in ('index-ascii','index-netCDF4','index-HDF5'): + # read from index file + _,removeform = REMOVEFORM.split('-') + # index containing files in data format + Ylms = gravtk.harmonics().from_index(REMOVE_FILE, + format=removeform) + attributes['lineage'].extend([f.name for f in Ylms.filename]) + # reduce to GRACE/GRACE-FO months and truncate to degree and order + Ylms = Ylms.subset(GSM_Ylms.month).truncate(lmax=LMAX, mmax=MMAX) + # remove the temporal mean of the coefficients + Ylms.mean(apply=True) + # distribute removed Ylms uniformly over the ocean + if REDISTRIBUTE_REMOVED: + # calculate ratio between total removed mass and + # a uniformly distributed cm of water over the ocean + ratio = Ylms.clm[0,0,:]/ocean_Ylms.clm[0,0] + # for each spherical harmonic + for m in range(0,MMAX+1):# MMAX+1 to include MMAX + for l in range(m,LMAX+1):# LMAX+1 to include LMAX + # remove the ratio*ocean Ylms from Ylms + # note: x -= y is equivalent to x = x - y + Ylms.clm[l,m,:] -= ratio*ocean_Ylms.clm[l,m] + Ylms.slm[l,m,:] -= ratio*ocean_Ylms.slm[l,m] + # filter removed coefficients + if DESTRIPE: + Ylms = Ylms.destripe() + # add data for month t and INDEX_FILE to the total + # remove_clm and remove_slm matrices + # 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) + # Calculating cos/sin of phi arrays # output [m,phi] - m = GSM_Ylms.m[:, np.newaxis] + m = GSM_Ylms.m # Integration factors (solid angle) int_fact = np.sin(th)*dphi*dth - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # 4-pi normalization + norm = 1.0/(4.0*np.pi) + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", m, phi)) # Legendre polynomials for degree 1 P10 = np.squeeze(PLM[1,0,:]) @@ -637,34 +703,31 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, # PLM for spherical harmonic degrees 2+ up to LMAX # converted into mass and smoothed if specified plmout = np.zeros((LMAX+1, MMAX+1, nlat)) - for l in range(1,LMAX+1): - m = np.arange(0,np.min([l,MMAX])+1) - # convert to smoothed coefficients of mass - # Convolving plms with degree dependent factor and smoothing - plmout[l,m,:] = PLM[l,m,:]*dfactor[l]*wt[l] + # convert to smoothed coefficients of mass + # Convolving plms with degree dependent factor and smoothing + plmout[:] = np.einsum("l,l,lmh->lmh", dfactor, wt, PLM[:LMAX+1,:MMAX+1,:]) # Initializing 3x3 I-Parameter matrix - IMAT = np.zeros((3,3)) - # Calculating I-Parameter matrix by integrating over latitudes + # (see equations 12 and 13 of Swenson et al., 2008) + IMAT = np.zeros((3, 3)) # I-Parameter matrix accounts for the fact that the GRACE data only # includes spherical harmonic degrees greater than or equal to 2 - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # C10: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PS11)/(4.0*np.pi) - # C11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PS11)/(4.0*np.pi) - # S11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PS11)/(4.0*np.pi) + # C10, C11, S11 + PC10 = np.einsum("h...,p...->ph...", P10, m_phi[0,:].real) + PC11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].real) + PS11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].imag) + # C10: C10, C11, S11 + IMAT[0,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC10) + IMAT[1,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC11) + IMAT[2,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PS11) + # C11: C10, C11, S11 + IMAT[0,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC10) + IMAT[1,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC11) + IMAT[2,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PS11) + # S11: C10, C11, S11 + IMAT[0,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC10) + IMAT[1,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC11) + IMAT[2,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PS11) # get seasonal variations of an initial geocenter correction # for use in the land water mass calculation @@ -691,39 +754,32 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, DMAT = np.zeros((3,n_files)) # degree 1 iterations iteration = gravtk.geocenter() - iteration.C10 = np.zeros((n_files,max_iter)) - iteration.C11 = np.zeros((n_files,max_iter)) - iteration.S11 = np.zeros((n_files,max_iter)) + iteration.C10 = np.zeros((n_files, max_iter)) + iteration.C11 = np.zeros((n_files, max_iter)) + iteration.S11 = np.zeros((n_files, max_iter)) # calculate non-iterated terms for each file (G-matrix parameters) for t in range(n_files): # calculate geocenter component of ocean mass with GRACE - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] # Summing product of plms and c/slms over all SH degrees >= 2 - # Removing monthly GIA signal and atmospheric correction + # Removing monthly GIA signal, atmospheric correction + # and the auxiliary coefficients Ylms = GSM_Ylms.index(t) Ylms.subtract(GIA_Ylms.index(t)) Ylms.subtract(ATM_Ylms.index(t)) - for i in range(0, nlat): - l = np.arange(2,LMAX+1) - pcos[:,i] = np.sum(plmout[l,:,i]*Ylms.clm[l,:], axis=0) - psin[:,i] = np.sum(plmout[l,:,i]*Ylms.slm[l,:], axis=0) + Ylms.subtract(remove_Ylms.index(t)) + # subset GRACE to degrees 2+ for calculating ocean mass + l2 = slice(2, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l2, :, :], Ylms.ilm[l2, :]) # Multiplying by c/s(phi#m) to get surface density in cmwe (lon,lat) # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - rmass = np.dot(np.transpose(ccos),pcos) + np.dot(np.transpose(ssin),psin) + rmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # calculate G matrix parameters through a summation of each latitude - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # summation of integration factors, Legendre polynomials, - # (convolution of order and harmonics) and the ocean mass at t - G.C10[t] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.C11[t] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.S11[t] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) + # summation of integration factors, Legendre polynomials, + # (convolution of order and harmonics) and the ocean mass at t + G.C10[t] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, rmass) + G.C11[t] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, rmass) + G.S11[t] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, rmass) # calculate degree one solution for each iteration (or single if not) while (eps > eps_max) and (n_iter < max_iter): @@ -743,26 +799,22 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, GSM_Ylms.clm[1,1,t] = iteration.C11[t,n_iter-1] GSM_Ylms.slm[1,1,t] = iteration.S11[t,n_iter-1] - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] # Summing product of plms and c/slms over all SH degrees - # Removing monthly GIA signal and atmospheric correction + # Removing monthly GIA signal, atmospheric correction + # and the auxiliary coefficients Ylms = GSM_Ylms.index(t) Ylms.subtract(GIA_Ylms.index(t)) Ylms.subtract(ATM_Ylms.index(t)) - for i in range(0, nlat): - # for land water: use an initial seasonal geocenter estimate - # from Chen et al. (1999) then the iterative if specified - l = np.arange(1,LMAX+1) - pcos[:,i] = np.sum(plmout[l,:,i]*Ylms.clm[l,:], axis=0) - psin[:,i] = np.sum(plmout[l,:,i]*Ylms.slm[l,:], axis=0) + Ylms.subtract(remove_Ylms.index(t)) + # for land water: use an initial seasonal geocenter estimate + # from Chen et al. (1999) then the iterative if specified + l1 = slice(1, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l1, :, :], Ylms.ilm[l1, :]) # Multiplying by c/s(phi#m) to get surface density in cm w.e. (lonxlat) - # this will be a spatial field similar to outputs from stokes_combine.py # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - lmass = np.dot(np.transpose(ccos),pcos) + np.dot(np.transpose(ssin),psin) + lmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # use sea level fingerprints or eustatic from GRACE land components if FINGERPRINT: @@ -780,8 +832,7 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, sea_level = gravtk.sea_level_equation(land_Ylms.clm, land_Ylms.slm, landsea.lon, landsea.lat, land_function, LMAX=EXPANSION, LOVE=LOVE, BODY_TIDE_LOVE=0, FLUID_LOVE=0, ITERATIONS=3, - POLAR=True, PLM=PLM, ASTYPE=np.float64, SCALE=1e-32, - FILL_VALUE=0) + POLAR=True, PLM=PLM, FILL_VALUE=0) # 3) convert sea level fingerprints into spherical harmonics slf_Ylms = gravtk.gen_stokes(sea_level, landsea.lon, landsea.lat, UNITS=1, LMIN=0, LMAX=1, PLM=PLM[:2,:2,:], LOVE=LOVE) @@ -838,10 +889,14 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, elif SOLVER in ('gelsd', 'gelsy', 'gelss'): DMAT[:,t], res, rnk, s = scipy.linalg.lstsq(IMAT, (CMAT-GMAT), lapack_driver=SOLVER) - # save geocenter for iteration and time t after restoring GIA+ATM - iteration.C10[t,n_iter] = DMAT[0,t]/dfactor[1]+gia.C10[t]+atm.C10[t] - iteration.C11[t,n_iter] = DMAT[1,t]/dfactor[1]+gia.C11[t]+atm.C11[t] - iteration.S11[t,n_iter] = DMAT[2,t]/dfactor[1]+gia.S11[t]+atm.S11[t] + # save geocenter for iteration and time t after restoring fields + iteration.C10[t,n_iter] = DMAT[0,t]/dfactor[1] + \ + gia.C10[t] + atm.C10[t] + remove.C10[t] + iteration.C11[t,n_iter] = DMAT[1,t]/dfactor[1] + \ + gia.C11[t] + atm.C11[t] + remove.C11[t] + iteration.S11[t,n_iter] = DMAT[2,t]/dfactor[1] + \ + gia.S11[t] + atm.S11[t] + remove.S11[t] + # remove mean of each solution for iteration iteration.C10[:,n_iter] -= iteration.C10[:,n_iter].mean() iteration.C11[:,n_iter] -= iteration.C11[:,n_iter].mean() @@ -859,12 +914,13 @@ def calc_degree_one(base_dir, PROC, DREL, MODEL, LMAX, RAD, # Convert inverted solutions into fully normalized spherical harmonics # restore geocenter variation from glacial isostatic adjustment (GIA) # restore atmospheric jump corrections from Fagiolini (2015) if applicable + # restore auxiliary removed spherical harmonics if applicable # for each of the geocenter solutions (C10, C11, S11) # for the iterative case this will be the final iteration DEG1 = gravtk.geocenter() - DEG1.C10 = DMAT[0,:]/dfactor[1] + gia.C10[:] + atm.C10[:] - DEG1.C11 = DMAT[1,:]/dfactor[1] + gia.C11[:] + atm.C11[:] - DEG1.S11 = DMAT[2,:]/dfactor[1] + gia.S11[:] + atm.S11[:] + DEG1.C10 = DMAT[0,:]/dfactor[1] + gia.C10[:] + atm.C10[:] + remove.C10[:] + DEG1.C11 = DMAT[1,:]/dfactor[1] + gia.C11[:] + atm.C11[:] + remove.C11[:] + DEG1.S11 = DMAT[2,:]/dfactor[1] + gia.S11[:] + atm.S11[:] + remove.S11[:] # remove mean of geocenter for each component DEG1.mean(apply=True) # calculate geocenter variations with dealiasing restored @@ -1528,6 +1584,19 @@ def arguments(): parser.add_argument('--mean-format', type=str, default='netCDF4', choices=['ascii','netCDF4','HDF5','gfc'], help='Input data format for GRACE/GRACE-FO mean file') + # monthly files to be removed from the GRACE/GRACE-FO data + parser.add_argument('--remove-file', + type=pathlib.Path, nargs='+', + help='Monthly files to be removed from the GRACE/GRACE-FO data') + choices = [] + choices.extend(['ascii','netCDF4','HDF5']) + choices.extend(['index-ascii','index-netCDF4','index-HDF5']) + parser.add_argument('--remove-format', + type=str, nargs='+', choices=choices, + help='Input data format for files to be removed') + parser.add_argument('--redistribute-removed', + default=False, action='store_true', + help='Redistribute removed mass fields over the ocean') # run with iterative scheme parser.add_argument('--iterative', default=False, action='store_true', @@ -1617,6 +1686,9 @@ def main(): MODEL_INDEX=args.ocean_file, MEAN_FILE=args.mean_file, MEANFORM=args.mean_format, + REMOVE_FILES=args.remove_file, + REMOVE_FORMAT=args.remove_format, + REDISTRIBUTE_REMOVED=args.redistribute_removed, ITERATIVE=args.iterative, SOLVER=args.solver, FINGERPRINT=args.fingerprint, diff --git a/geocenter/delta_degree_one.py b/geocenter/delta_degree_one.py index a621affb..ef6c2824 100644 --- a/geocenter/delta_degree_one.py +++ b/geocenter/delta_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" delta_degree_one.py -Written by Tyler Sutterley (10/2023) +Written by Tyler Sutterley (07/2026) Calculates degree 1 errors using GRACE coefficients of degree 2 and greater, and ocean bottom pressure variations from OMCT/MPIOM @@ -147,6 +147,8 @@ https://doi.org/10.1029/2005GL025305 UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 10/2023: generalize mission variable to be GRACE/GRACE-FO Updated 09/2023: simplify I-matrix and G-matrix calculations Updated 05/2023: use pathlib to define and operate on paths @@ -317,11 +319,11 @@ def delta_degree_one(base_dir, PROC, DREL, LMAX, RAD, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # spatial parameters in radians - dphi = dlon*np.pi/180.0 - dth = dlat*np.pi/180.0 + dphi = np.radians(dlon) + dth = np.radians(dlat) # longitude and colatitude in radians - phi = landsea.lon[np.newaxis,:]*np.pi/180.0 - th = (90.0 - np.squeeze(landsea.lat))*np.pi/180.0 + phi = np.radians(landsea.lon[np.newaxis,:]) + th = np.radians(90.0 - np.squeeze(landsea.lat)) # create land function land_function = np.zeros((nlon, nlat),dtype=np.float64) # extract land function from file @@ -442,47 +444,45 @@ def delta_degree_one(base_dir, PROC, DREL, LMAX, RAD, # Calculating cos/sin of phi arrays # output [m,phi] - m = GSM_Ylms.m[:, np.newaxis] + m = GSM_Ylms.m # Integration factors (solid angle) int_fact = np.sin(th)*dphi*dth - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # 4-pi normalization + norm = 1.0/(4.0*np.pi) + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", m, phi)) # Legendre polynomials for degree 1 P10 = np.squeeze(PLM[1,0,:]) P11 = np.squeeze(PLM[1,1,:]) - # PLM for spherical harmonic degrees 2+ + # PLM for spherical harmonic degrees 2+ up to LMAX # converted into mass and smoothed if specified plmout = np.zeros((LMAX+1, MMAX+1, nlat)) - for l in range(1,LMAX+1): - m = np.arange(0,np.min([l,MMAX])+1) - # convert to smoothed coefficients of mass - # Convolving plms with degree dependent factor and smoothing - plmout[l,m,:] = PLM[l,m,:]*dfactor[l]*wt[l] + # convert to smoothed coefficients of mass + # Convolving plms with degree dependent factor and smoothing + plmout[:] = np.einsum("l,l,lmh->lmh", dfactor, wt, PLM[:LMAX+1,:MMAX+1,:]) # Initializing 3x3 I-Parameter matrix - IMAT = np.zeros((3,3)) - # Calculating I-Parameter matrix by integrating over latitudes + # (see equations 12 and 13 of Swenson et al., 2008) + IMAT = np.zeros((3, 3)) # I-Parameter matrix accounts for the fact that the GRACE data only # includes spherical harmonic degrees greater than or equal to 2 - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # C10: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PS11)/(4.0*np.pi) - # C11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PS11)/(4.0*np.pi) - # S11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PS11)/(4.0*np.pi) + # C10, C11, S11 + PC10 = np.einsum("h...,p...->ph...", P10, m_phi[0,:].real) + PC11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].real) + PS11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].imag) + # C10: C10, C11, S11 + IMAT[0,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC10) + IMAT[1,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC11) + IMAT[2,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PS11) + # C11: C10, C11, S11 + IMAT[0,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC10) + IMAT[1,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC11) + IMAT[2,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PS11) + # S11: C10, C11, S11 + IMAT[0,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC10) + IMAT[1,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC11) + IMAT[2,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PS11) # iterate solutions: if not single iteration n_iter = 0 @@ -499,9 +499,6 @@ def delta_degree_one(base_dir, PROC, DREL, LMAX, RAD, # Allocate for G matrix parameters # G matrix calculates the GRACE ocean mass variations G = gravtk.geocenter() - G.C10 = 0.0 - G.C11 = 0.0 - G.S11 = 0.0 # degree 1 iterations iteration = gravtk.geocenter() iteration.C10 = np.zeros((max_iter)) @@ -509,29 +506,19 @@ def delta_degree_one(base_dir, PROC, DREL, LMAX, RAD, iteration.S11 = np.zeros((max_iter)) # calculate non-iterated terms (G-matrix parameters) # calculate geocenter component of ocean mass with GRACE - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] - # Summing product of plms and c/slms over all SH degrees >= 2 - for i in range(0, nlat): - l = np.arange(2,LMAX+1) - pcos[:,i] = np.sum(((plmout[l,:,i]*delta_Ylms.clm[l,:])**2)/nsmth, axis=0) - psin[:,i] = np.sum(((plmout[l,:,i]*delta_Ylms.slm[l,:])**2)/nsmth, axis=0) + # subset GRACE to degrees 2+ for calculating ocean mass + l2 = slice(2, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l2, :, :], delta_Ylms.ilm[l2, :]) # Multiplying by c/s(phi#m) to get surface density in cmwe (lon,lat) # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - rmass = np.sqrt(np.dot(np.transpose(ccos**2),pcos) + np.dot(np.transpose(ssin**2),psin)) + rmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # calculate G matrix parameters through a summation of each latitude - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # summation of integration factors, Legendre polynomials, - # (convolution of order and harmonics) and the ocean mass at t - G.C10 += np.sum(int_fact[i]*PC10*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.C11 += np.sum(int_fact[i]*PC11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.S11 += np.sum(int_fact[i]*PS11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) + # summation of integration factors, Legendre polynomials, + # (convolution of order and harmonics) and the ocean mass at t + G.C10 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, rmass) + G.C11 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, rmass) + G.S11 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, rmass) # calculate degree one solution for each iteration (or single if not) while (eps > eps_max) and (n_iter < max_iter): @@ -549,22 +536,14 @@ def delta_degree_one(base_dir, PROC, DREL, LMAX, RAD, delta_Ylms.clm[1,1] = iteration.C11[n_iter-1] delta_Ylms.slm[1,1] = iteration.S11[n_iter-1] - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] # Summing product of plms and c/slms over all SH degrees - for i in range(0, nlat): - # for land water: use an initial seasonal geocenter estimate - # from Chen et al. (1999) - l = np.arange(1,LMAX+1) - pcos[:,i] = np.sum(((plmout[l,:,i]*delta_Ylms.clm[l,:])**2)/nsmth, axis=0) - psin[:,i] = np.sum(((plmout[l,:,i]*delta_Ylms.slm[l,:])**2)/nsmth, axis=0) + l1 = slice(1, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l1, :, :], delta_Ylms.ilm[l1, :]) # Multiplying by c/s(phi#m) to get surface density in cm w.e. (lonxlat) - # this will be a spatial field similar to outputs from stokes_combine.py # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - lmass = np.sqrt(np.dot(np.transpose(ccos**2),pcos) + np.dot(np.transpose(ssin**2),psin)) + lmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # use sea level fingerprints or eustatic from GRACE land components if FINGERPRINT: diff --git a/geocenter/geocenter_monte_carlo.py b/geocenter/geocenter_monte_carlo.py index 925627a1..62d578d6 100644 --- a/geocenter/geocenter_monte_carlo.py +++ b/geocenter/geocenter_monte_carlo.py @@ -128,7 +128,7 @@ def geocenter_monte_carlo(grace_dir,PROC,DREL,START_MON,END_MON,MISSING): artist = matplotlib.offsetbox.AnchoredText(axes_labels[key], pad=0., prop=dict(size=16,weight='bold'), frameon=False, loc=2) ax[j].add_artist(artist) - lbl = f'$\sigma$ = {RMS:0.2f} mm\nmax = {max_var:0.2f} mm' + lbl = r'$\sigma$' + f' = {RMS:0.2f} mm\nmax = {max_var:0.2f} mm' artist = matplotlib.offsetbox.AnchoredText(lbl, pad=0., prop=dict(size=12), frameon=False, loc=3) ax[j].add_artist(artist) diff --git a/geocenter/geocenter_spatial_maps.py b/geocenter/geocenter_spatial_maps.py index 1d06bcda..84640abc 100644 --- a/geocenter/geocenter_spatial_maps.py +++ b/geocenter/geocenter_spatial_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" geocenter_spatial_maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (07/2026) Reads in GRACE/GRACE-FO geocenter coefficients and exports trends in the monthly spatial fields in millimeters water equivalent @@ -90,6 +90,7 @@ utilities.py: download and management utilities for files UPDATE HISTORY: + Updated 07/2026: use np.radians to convert from degrees to radians Updated 05/2023: split S2 tidal aliasing terms into GRACE and GRACE-FO eras use pathlib to define and operate on paths Updated 01/2023: refactored time series analysis functions @@ -265,7 +266,7 @@ def geocenter_spatial_maps(base_dir, PROC, DREL, nlat = len(grid.lat) # Computing plms for converting to spatial domain - theta = (90.0-grid.lat)*np.pi/180.0 + theta = np.radians(90.0-grid.lat) PLM, dPLM = gravtk.plm_holmes(2, np.cos(theta)) # fit coefficients fits = ['x1','x2','SS','SC','AS','AC','S2SGRC','S2CGRC','S2SGFO','S2CGFO'] @@ -342,8 +343,7 @@ def geocenter_spatial_maps(base_dir, PROC, DREL, # output phase to file for key in ['SP','AP','S2PGRC','S2PGFO']: # convert phase from -180:180 to 0:360 - ix,iy = np.nonzero(var[key] < 0) - var[key][ix,iy] += 360.0 + var[key] = np.where(var[key] < 0, var[key] + 360.0, var[key]) # copy variables to output grid grid.data = np.copy(var[key]) grid.time = np.copy(MEAN.time) diff --git a/geocenter/kernel_degree_one.py b/geocenter/kernel_degree_one.py index 30d4d194..4f59f1eb 100644 --- a/geocenter/kernel_degree_one.py +++ b/geocenter/kernel_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" kernel_degree_one.py -Written by Tyler Sutterley (06/2024) +Written by Tyler Sutterley (07/2026) Calculates the sensitivity of geocenter calculations for each degree and order Can be used to estimate the geocenter uncertainties for sets of harmonics @@ -107,6 +107,8 @@ https://doi.org/10.1029/2007JB005338 UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 06/2024: use wrapper to importlib for optional dependencies Updated 09/2023: simplify I-matrix and G-matrix calculations Updated 05/2023: use pathlib to define and operate on paths @@ -201,11 +203,11 @@ def kernel_degree_one(base_dir, LMAX, RAD, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # spatial parameters in radians - dphi = dlon*np.pi/180.0 - dth = dlat*np.pi/180.0 + dphi = np.radians(dlon) + dth = np.radians(dlat) # longitude and colatitude in radians - phi = landsea.lon[np.newaxis,:]*np.pi/180.0 - th = (90.0 - np.squeeze(landsea.lat))*np.pi/180.0 + phi = np.radians(landsea.lon[np.newaxis,:]) + th = np.radians(90.0 - np.squeeze(landsea.lat)) # create land function land_function = np.zeros((nlon, nlat),dtype=np.float64) # extract land function from file @@ -236,39 +238,39 @@ def kernel_degree_one(base_dir, LMAX, RAD, # Calculating cos/sin of phi arrays # output [m,phi] - m = np.arange(0,MMAX+1)[:, np.newaxis] + m = np.arange(0,MMAX+1) # Integration factors (solid angle) int_fact = np.sin(th)*dphi*dth - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # 4-pi normalization + norm = 1.0/(4.0*np.pi) + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", m, phi)) # Legendre polynomials for degree 1 P10 = np.squeeze(PLM[1,0,:]) P11 = np.squeeze(PLM[1,1,:]) # Initializing 3x3 I-Parameter matrix - IMAT = np.zeros((3,3)) - # Calculating I-Parameter matrix by integrating over latitudes + # (see equations 12 and 13 of Swenson et al., 2008) + IMAT = np.zeros((3, 3)) # I-Parameter matrix accounts for the fact that the GRACE data only # includes spherical harmonic degrees greater than or equal to 2 - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # C10: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PS11)/(4.0*np.pi) - # C11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PS11)/(4.0*np.pi) - # S11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PS11)/(4.0*np.pi) + # C10, C11, S11 + PC10 = np.einsum("h...,p...->ph...", P10, m_phi[0,:].real) + PC11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].real) + PS11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].imag) + # C10: C10, C11, S11 + IMAT[0,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC10) + IMAT[1,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC11) + IMAT[2,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PS11) + # C11: C10, C11, S11 + IMAT[0,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC10) + IMAT[1,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC11) + IMAT[2,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PS11) + # S11: C10, C11, S11 + IMAT[0,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC10) + IMAT[1,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC11) + IMAT[2,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PS11) # output flag for using sea level fingerprints slf_str = '_SLF' if FINGERPRINT else '' @@ -313,9 +315,9 @@ def kernel_degree_one(base_dir, LMAX, RAD, fit_factor[ii] = wt[l]*(2.0*l + 1.0)/(1.0 + LOVE.kl[l]) # cosine and sine factors if (csharm == 'clm'): - mphi[ii,:] = ccos[m,:] + mphi[ii,:] = m_phi[m,:].real elif (csharm == 'slm'): - mphi[ii,:] = ssin[m,:] + mphi[ii,:] = m_phi[m,:].imag # add 1 to counter ii += 1 @@ -323,10 +325,10 @@ def kernel_degree_one(base_dir, LMAX, RAD, output['covariance'] = np.zeros((n_harm, n_harm, 3)) for i in range(n_harm): # setting kern_i equal to 1 for d/o - kern_i = np.zeros((n_harm,nlat)) + kern_i = np.zeros((n_harm, nlat)) kern_i[i,:] = 1.0*fit_factor[i] # calculate land mass maps - lmass = np.dot(mphi.T,plm*kern_i) + lmass = np.dot(mphi.T, plm*kern_i) # Calculating data matrices # GRACE Eustatic degree 1 from land variations eustatic = gravtk.geocenter() @@ -346,7 +348,7 @@ def kernel_degree_one(base_dir, LMAX, RAD, sea_level = gravtk.sea_level_equation(land_Ylms.clm, land_Ylms.slm, landsea.lon, landsea.lat, land_function, LMAX=EXPANSION, LOVE=LOVE_K1, BODY_TIDE_LOVE=0, FLUID_LOVE=0, ITERATIONS=3, - POLAR=True, PLM=PLM, ASTYPE=np.float64, SCALE=1e-32, FILL_VALUE=0) + POLAR=True, PLM=PLM, FILL_VALUE=0) # 3) convert sea level fingerprints into spherical harmonics slf_Ylms = gravtk.gen_stokes(sea_level, landsea.lon, landsea.lat, UNITS=1, LMIN=0, LMAX=1, PLM=PLM[:2,:2,:], LOVE=LOVE) @@ -385,24 +387,16 @@ def kernel_degree_one(base_dir, LMAX, RAD, if (j >= 3): kern_j[j,:] = 1.0*fit_factor[j] # calculate ocean mass maps - rmass = np.dot(mphi.T,plm*kern_j) + rmass = np.dot(mphi.T, plm*kern_j) # Allocate for G matrix parameters # G matrix calculates the GRACE ocean mass variations G = gravtk.geocenter() - G.C10 = 0.0 - G.C11 = 0.0 - G.S11 = 0.0 # calculate G matrix parameters through a summation of each latitude - for n in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # summation of integration factors, Legendre polynomials, - # (convolution of order and harmonics) and the ocean mass - G.C10 += np.sum(int_fact[n]*PC10*ocean_function[:,n]*rmass[:,n])/(4.0*np.pi) - G.C11 += np.sum(int_fact[n]*PC11*ocean_function[:,n]*rmass[:,n])/(4.0*np.pi) - G.S11 += np.sum(int_fact[n]*PS11*ocean_function[:,n]*rmass[:,n])/(4.0*np.pi) + # summation of integration factors, Legendre polynomials, + # (convolution of order and harmonics) and the ocean mass at t + G.C10 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, rmass) + G.C11 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, rmass) + G.S11 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, rmass) # G Matrix for time t GMAT = np.array([G.C10, G.C11, G.S11]) diff --git a/geocenter/model_degree_one.py b/geocenter/model_degree_one.py index ec788f15..e6aeb516 100755 --- a/geocenter/model_degree_one.py +++ b/geocenter/model_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" model_degree_one.py -Written by Tyler Sutterley (01/2025) +Written by Tyler Sutterley (07/2026) Calculates degree 1 variations using synthetic coefficients of degree 2 and greater for testing the reliability of the algorithm @@ -120,6 +120,8 @@ https://doi.org/10.1029/2007JB005338 UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 01/2025: fixed deprecated tick label resizing Updated 09/2023: simplify I-matrix and G-matrix calculations Updated 05/2023: use pathlib to define and operate on paths @@ -208,12 +210,12 @@ def model_seasonal_geocenter(grace_date): 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 + APx*np.pi/180.0) + \ - SAAx*np.sin(4.0*np.pi*grace_date + SAPx*np.pi/180.0) - Y = AAy*np.sin(2.0*np.pi*grace_date + APy*np.pi/180.0) + \ - SAAy*np.sin(4.0*np.pi*grace_date + SAPy*np.pi/180.0) - Z = AAz*np.sin(2.0*np.pi*grace_date + APz*np.pi/180.0) + \ - SAAz*np.sin(4.0*np.pi*grace_date + SAPz*np.pi/180.0) + 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() @@ -271,11 +273,11 @@ def model_degree_one(input_file, LMAX, RAD, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # spatial parameters in radians - dphi = dlon*np.pi/180.0 - dth = dlat*np.pi/180.0 + dphi = np.radians(dlon) + dth = np.radians(dlat) # longitude and colatitude in radians - phi = landsea.lon[np.newaxis,:]*np.pi/180.0 - th = (90.0 - np.squeeze(landsea.lat))*np.pi/180.0 + phi = np.radians(landsea.lon[np.newaxis,:]) + th = np.radians(90.0 - np.squeeze(landsea.lat)) # create land function land_function = np.zeros((nlon, nlat),dtype=np.float64) # extract land function from file @@ -326,12 +328,13 @@ def model_degree_one(input_file, LMAX, RAD, # Calculating cos/sin of phi arrays # output [m,phi] - m = data_Ylms.m[:, np.newaxis] + m = data_Ylms.m # Integration factors (solid angle) int_fact = np.sin(th)*dphi*dth - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # 4-pi normalization + norm = 1.0/(4.0*np.pi) + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", m, phi)) # Legendre polynomials for degree 1 P10 = np.squeeze(PLM[1,0,:]) @@ -339,34 +342,31 @@ def model_degree_one(input_file, LMAX, RAD, # PLM for spherical harmonic degrees 2+ up to LMAX # converted into mass and smoothed if specified plmout = np.zeros((LMAX+1, MMAX+1, nlat)) - for l in range(1,LMAX+1): - m = np.arange(0,np.min([l,MMAX])+1) - # convert to smoothed coefficients of mass - # Convolving plms with degree dependent factor and smoothing - plmout[l,m,:] = PLM[l,m,:]*dfactor[l]*wt[l] + # convert to smoothed coefficients of mass + # Convolving plms with degree dependent factor and smoothing + plmout[:] = np.einsum("l,l,lmh->lmh", dfactor, wt, PLM[:LMAX+1,:MMAX+1,:]) # Initializing 3x3 I-Parameter matrix - IMAT = np.zeros((3,3)) - # Calculating I-Parameter matrix by integrating over latitudes - # I-Parameter matrix accounts for the fact that the harmonic data only + # (see equations 12 and 13 of Swenson et al., 2008) + IMAT = np.zeros((3, 3)) + # I-Parameter matrix accounts for the fact that the GRACE data only # includes spherical harmonic degrees greater than or equal to 2 - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # C10: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PS11)/(4.0*np.pi) - # C11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PS11)/(4.0*np.pi) - # S11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PS11)/(4.0*np.pi) + # C10, C11, S11 + PC10 = np.einsum("h...,p...->ph...", P10, m_phi[0,:].real) + PC11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].real) + PS11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].imag) + # C10: C10, C11, S11 + IMAT[0,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC10) + IMAT[1,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC11) + IMAT[2,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PS11) + # C11: C10, C11, S11 + IMAT[0,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC10) + IMAT[1,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC11) + IMAT[2,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PS11) + # S11: C10, C11, S11 + IMAT[0,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC10) + IMAT[1,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC11) + IMAT[2,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PS11) # get seasonal variations of an initial geocenter correction # for use in the land water mass calculation @@ -400,29 +400,21 @@ def model_degree_one(input_file, LMAX, RAD, # calculate non-iterated terms for each file (G-matrix parameters) for t in range(n_files): # calculate geocenter component of ocean mass - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] # Summing product of plms and c/slms over all SH degrees >= 2 - for i in range(0, nlat): - l = np.arange(2,LMAX+1) - pcos[:,i] = np.sum(plmout[l,:,i]*data_Ylms.clm[l,:,t], axis=0) - psin[:,i] = np.sum(plmout[l,:,i]*data_Ylms.slm[l,:,t], axis=0) + Ylms = data_Ylms.index(t) + # subset GRACE to degrees 2+ for calculating ocean mass + l2 = slice(2, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l2, :, :], Ylms.ilm[l2, :]) # Multiplying by c/s(phi#m) to get surface density in cmwe (lon,lat) # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - rmass = np.dot(np.transpose(ccos),pcos) + np.dot(np.transpose(ssin),psin) + rmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # calculate G matrix parameters through a summation of each latitude - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # summation of integration factors, Legendre polynomials, - # (convolution of order and harmonics) and the ocean mass at t - G.C10[t] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.C11[t] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.S11[t] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) + # summation of integration factors, Legendre polynomials, + # (convolution of order and harmonics) and the ocean mass at t + G.C10[t] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, rmass) + G.C11[t] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, rmass) + G.S11[t] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, rmass) # calculate degree one solution for each iteration (or single if not) while (eps > eps_max) and (n_iter < max_iter): @@ -442,20 +434,15 @@ def model_degree_one(input_file, LMAX, RAD, data_Ylms.clm[1,1,t] = iteration.C11[t,n_iter-1] data_Ylms.slm[1,1,t] = iteration.S11[t,n_iter-1] - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] # Summing product of plms and c/slms over all SH degrees - for i in range(0, nlat): - l = np.arange(1,LMAX+1) - pcos[:,i] = np.sum(plmout[l,:,i]*data_Ylms.clm[l,:,t], axis=0) - psin[:,i] = np.sum(plmout[l,:,i]*data_Ylms.slm[l,:,t], axis=0) + Ylms = data_Ylms.index(t) + l1 = slice(1, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l1, :, :], Ylms.ilm[l1, :]) # Multiplying by c/s(phi#m) to get surface density in cm w.e. (lonxlat) - # this will be a spatial field similar to outputs from stokes_combine.py # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - lmass = np.dot(np.transpose(ccos),pcos) + np.dot(np.transpose(ssin),psin) + lmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # use sea level fingerprints or eustatic from land components if FINGERPRINT: diff --git a/geocenter/monte_carlo_degree_one.py b/geocenter/monte_carlo_degree_one.py index 7ecb2b6f..d822a76b 100644 --- a/geocenter/monte_carlo_degree_one.py +++ b/geocenter/monte_carlo_degree_one.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" monte_carlo_degree_one.py -Written by Tyler Sutterley (01/2025) +Written by Tyler Sutterley (07/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 @@ -157,6 +157,9 @@ https://doi.org/10.1029/2005GL025305 UPDATE HISTORY: + 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 Updated 01/2025: fixed deprecated tick label resizing Updated 06/2024: use wrapper to importlib for optional dependencies Updated 10/2023: generalize mission variable to be GRACE/GRACE-FO @@ -263,12 +266,12 @@ def model_seasonal_geocenter(grace_date): 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 + APx*np.pi/180.0) + \ - SAAx*np.sin(4.0*np.pi*grace_date + SAPx*np.pi/180.0) - Y = AAy*np.sin(2.0*np.pi*grace_date + APy*np.pi/180.0) + \ - SAAy*np.sin(4.0*np.pi*grace_date + SAPy*np.pi/180.0) - Z = AAz*np.sin(2.0*np.pi*grace_date + APz*np.pi/180.0) + \ - SAAz*np.sin(4.0*np.pi*grace_date + SAPz*np.pi/180.0) + 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() @@ -295,6 +298,9 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, DATAFORM=None, MEAN_FILE=None, MEANFORM=None, + REMOVE_FILES=None, + REMOVE_FORMAT=None, + REDISTRIBUTE_REMOVED=False, ERROR_FILES=[], SOLVER=None, FINGERPRINT=False, @@ -404,11 +410,11 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # spatial parameters in radians - dphi = dlon*np.pi/180.0 - dth = dlat*np.pi/180.0 + dphi = np.radians(dlon) + dth = np.radians(dlat) # longitude and colatitude in radians - phi = landsea.lon[np.newaxis,:]*np.pi/180.0 - th = (90.0 - np.squeeze(landsea.lat))*np.pi/180.0 + phi = np.radians(landsea.lon[np.newaxis,:]) + th = np.radians(90.0 - np.squeeze(landsea.lat)) # create land function land_function = np.zeros((nlon, nlat),dtype=np.float64) # extract land function from file @@ -505,6 +511,57 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, # save geocenter coefficients of the atmospheric jump corrections 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 + remove_Ylms = GSM_Ylms.zeros_like() + remove_Ylms.time[:] = np.copy(GSM_Ylms.time) + remove_Ylms.month[:] = np.copy(GSM_Ylms.month) + if REMOVE_FILES: + # extend list if a single format was entered for all files + if len(REMOVE_FORMAT) < len(REMOVE_FILES): + REMOVE_FORMAT = REMOVE_FORMAT*len(REMOVE_FILES) + # for each file to be removed + for REMOVE_FILE,REMOVEFORM in zip(REMOVE_FILES,REMOVE_FORMAT): + if REMOVEFORM in ('ascii','netCDF4','HDF5'): + # ascii (.txt) + # netCDF4 (.nc) + # HDF5 (.H5) + Ylms = gravtk.harmonics().from_file(REMOVE_FILE, + format=REMOVEFORM) + attributes['lineage'].append(Ylms.filename) + elif REMOVEFORM in ('index-ascii','index-netCDF4','index-HDF5'): + # read from index file + _,removeform = REMOVEFORM.split('-') + # index containing files in data format + Ylms = gravtk.harmonics().from_index(REMOVE_FILE, + format=removeform) + attributes['lineage'].extend([f.name for f in Ylms.filename]) + # reduce to GRACE/GRACE-FO months and truncate to degree and order + Ylms = Ylms.subset(GSM_Ylms.month).truncate(lmax=LMAX, mmax=MMAX) + # remove the temporal mean of the coefficients + Ylms.mean(apply=True) + # distribute removed Ylms uniformly over the ocean + if REDISTRIBUTE_REMOVED: + # calculate ratio between total removed mass and + # a uniformly distributed cm of water over the ocean + ratio = Ylms.clm[0,0,:]/ocean_Ylms.clm[0,0] + # for each spherical harmonic + for m in range(0,MMAX+1):# MMAX+1 to include MMAX + for l in range(m,LMAX+1):# LMAX+1 to include LMAX + # remove the ratio*ocean Ylms from Ylms + # note: x -= y is equivalent to x = x - y + Ylms.clm[l,m,:] -= ratio*ocean_Ylms.clm[l,m] + Ylms.slm[l,m,:] -= ratio*ocean_Ylms.slm[l,m] + # filter removed coefficients + if DESTRIPE: + Ylms = Ylms.destripe() + # add data for month t and INDEX_FILE to the total + # remove_clm and remove_slm matrices + # 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) + # input spherical harmonic datafiles to be used in monte carlo error_Ylms = [] # for each file to be removed @@ -580,47 +637,45 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, # Calculating cos/sin of phi arrays # output [m,phi] - m = GSM_Ylms.m[:, np.newaxis] + m = GSM_Ylms.m # Integration factors (solid angle) int_fact = np.sin(th)*dphi*dth - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # 4-pi normalization + norm = 1.0/(4.0*np.pi) + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", m, phi)) # Legendre polynomials for degree 1 P10 = np.squeeze(PLM[1,0,:]) P11 = np.squeeze(PLM[1,1,:]) - # PLM for spherical harmonic degrees 2+ + # PLM for spherical harmonic degrees 2+ up to LMAX # converted into mass and smoothed if specified plmout = np.zeros((LMAX+1, MMAX+1, nlat)) - for l in range(1,LMAX+1): - m = np.arange(0,np.min([l,MMAX])+1) - # convert to smoothed coefficients of mass - # Convolving plms with degree dependent factor and smoothing - plmout[l,m,:] = PLM[l,m,:]*dfactor[l]*wt[l] + # convert to smoothed coefficients of mass + # Convolving plms with degree dependent factor and smoothing + plmout[:] = np.einsum("l,l,lmh->lmh", dfactor, wt, PLM[:LMAX+1,:MMAX+1,:]) # Initializing 3x3 I-Parameter matrix - IMAT = np.zeros((3,3)) - # Calculating I-Parameter matrix by integrating over latitudes + # (see equations 12 and 13 of Swenson et al., 2008) + IMAT = np.zeros((3, 3)) # I-Parameter matrix accounts for the fact that the GRACE data only # includes spherical harmonic degrees greater than or equal to 2 - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # C10: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,0] += np.sum(int_fact[i]*PC10*ocean_function[:,i]*PS11)/(4.0*np.pi) - # C11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,1] += np.sum(int_fact[i]*PC11*ocean_function[:,i]*PS11)/(4.0*np.pi) - # S11: C10, C11, S11 (see equations 12 and 13 of Swenson et al., 2008) - IMAT[0,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC10)/(4.0*np.pi) - IMAT[1,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PC11)/(4.0*np.pi) - IMAT[2,2] += np.sum(int_fact[i]*PS11*ocean_function[:,i]*PS11)/(4.0*np.pi) + # C10, C11, S11 + PC10 = np.einsum("h...,p...->ph...", P10, m_phi[0,:].real) + PC11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].real) + PS11 = np.einsum("h...,p...->ph...", P11, m_phi[1,:].imag) + # C10: C10, C11, S11 + IMAT[0,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC10) + IMAT[1,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PC11) + IMAT[2,0] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, PS11) + # C11: C10, C11, S11 + IMAT[0,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC10) + IMAT[1,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PC11) + IMAT[2,1] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, PS11) + # S11: C10, C11, S11 + IMAT[0,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC10) + IMAT[1,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PC11) + IMAT[2,2] = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, PS11) # get seasonal variations of an initial geocenter correction # for use in the land water mass calculation @@ -644,67 +699,56 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, Ylms.clm += (1.0-2.0*np.random.rand(LMAX+1,MMAX+1))*eYlms.clm Ylms.slm += (1.0-2.0*np.random.rand(LMAX+1,MMAX+1))*eYlms.slm - # Removing monthly GIA signal and atmospheric correction + # Removing monthly GIA signal, atmospheric correction + # and the auxiliary coefficients GRACE_Ylms = GSM_Ylms.index(t) GRACE_Ylms.subtract(GIA_Ylms.index(t)) GRACE_Ylms.subtract(ATM_Ylms.index(t)) + GRACE_Ylms.subtract(remove_Ylms.index(t)) + # combining GRACE/GRACE-FO with uncertainty for monte carlo run + GRACE_Ylms.add(Ylms) # G matrix calculates the GRACE ocean mass variations G = gravtk.geocenter() G.C10 = 0.0 G.C11 = 0.0 G.S11 = 0.0 - # calculate non-iterated terms (G-matrix parameters) - # calculate geocenter component of ocean mass with GRACE - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] - # Summing product of plms and c/slms over all SH degrees >= 2 - for i in range(0, nlat): - l = np.arange(2,LMAX+1) - pcos[:,i] = np.sum(plmout[l,:,i]*(GRACE_Ylms.clm[l,:]+Ylms.clm[l,:]), axis=0) - psin[:,i] = np.sum(plmout[l,:,i]*(GRACE_Ylms.slm[l,:]+Ylms.slm[l,:]), axis=0) + # subset GRACE to degrees 2+ for calculating ocean mass + l2 = slice(2, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l2, :, :], GRACE_Ylms.ilm[l2, :]) # Multiplying by c/s(phi#m) to get surface density in cmwe (lon,lat) # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - rmass = np.dot(np.transpose(ccos),pcos) + np.dot(np.transpose(ssin),psin) + rmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # calculate G matrix parameters through a summation of each latitude - for i in range(0,nlat): - # C10, C11, S11 - PC10 = P10[i]*ccos[0,:] - PC11 = P11[i]*ccos[1,:] - PS11 = P11[i]*ssin[1,:] - # summation of integration factors, Legendre polynomials, - # (convolution of order and harmonics) and the ocean mass at t - G.C10 += np.sum(int_fact[i]*PC10*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.C11 += np.sum(int_fact[i]*PC11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) - G.S11 += np.sum(int_fact[i]*PS11*ocean_function[:,i]*rmass[:,i])/(4.0*np.pi) + # summation of integration factors, Legendre polynomials, + # (convolution of order and harmonics) and the ocean mass at t + G.C10 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC10, ocean_function, rmass) + G.C11 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PC11, ocean_function, rmass) + G.S11 = norm*np.einsum("h...,ph...,ph...,ph...->...", int_fact, PS11, ocean_function, rmass) # seasonal component of geocenter variation for land water GSM_Ylms.clm[1,0,t] = seasonal_geocenter.C10[t] GSM_Ylms.clm[1,1,t] = seasonal_geocenter.C11[t] GSM_Ylms.slm[1,1,t] = seasonal_geocenter.S11[t] - # Removing monthly GIA signal and atmospheric correction + # Removing monthly GIA signal, atmospheric correction + # and the auxiliary coefficients GRACE_Ylms = GSM_Ylms.index(t) GRACE_Ylms.subtract(GIA_Ylms.index(t)) GRACE_Ylms.subtract(ATM_Ylms.index(t)) + GRACE_Ylms.subtract(remove_Ylms.index(t)) + # combining GRACE/GRACE-FO with uncertainty for monte carlo run + GRACE_Ylms.add(Ylms) - # allocate for product of grace and legendre polynomials - pcos = np.zeros((MMAX+1, nlat))#-[m,lat] - psin = np.zeros((MMAX+1, nlat))#-[m,lat] - # Summing product of plms and c/slms over all SH degrees - for i in range(0, nlat): - # for land water: use an initial seasonal geocenter estimate - # from Chen et al. (1999) - l = np.arange(1,LMAX+1) - pcos[:,i] = np.sum(plmout[l,:,i]*(GRACE_Ylms.clm[l,:]+Ylms.clm[l,:]), axis=0) - psin[:,i] = np.sum(plmout[l,:,i]*(GRACE_Ylms.slm[l,:]+Ylms.slm[l,:]), axis=0) + # for land water: use an initial seasonal geocenter estimate + # from Chen et al. (1999) then the iterative if specified + l1 = slice(1, LMAX+1) + pconv = np.einsum("lmh...,lm...->mh...", plmout[l1, :, :], GRACE_Ylms.ilm[l1, :]) # Multiplying by c/s(phi#m) to get surface density in cm w.e. (lonxlat) - # this will be a spatial field similar to outputs from stokes_combine.py # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta # The summation over spherical harmonic order is in this multiplication - lmass = np.dot(np.transpose(ccos),pcos) + np.dot(np.transpose(ssin),psin) + lmass = np.einsum("mp...,mh...->ph...", m_phi, pconv) # use sea level fingerprints or eustatic from GRACE land components if FINGERPRINT: @@ -761,10 +805,13 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, elif SOLVER in ('gelsd', 'gelsy', 'gelss'): DMAT, res, rnk, s = scipy.linalg.lstsq(IMAT, (CMAT-GMAT), lapack_driver=SOLVER) - # save geocenter for iteration and time t after restoring GIA+ATM - iteration.C10[t,n_iter] = DMAT[0]+gia.C10[t]+atm.C10[t] - iteration.C11[t,n_iter] = DMAT[1]+gia.C11[t]+atm.C11[t] - iteration.S11[t,n_iter] = DMAT[2]+gia.S11[t]+atm.S11[t] + # save geocenter for iteration and time t after restoring fields + iteration.C10[t,n_iter] = DMAT[0]/dfactor[1] + \ + gia.C10[t] + atm.C10[t] + remove.C10[t] + iteration.C11[t,n_iter] = DMAT[1]/dfactor[1] + \ + gia.C11[t] + atm.C11[t] + remove.C11[t] + iteration.S11[t,n_iter] = DMAT[2]/dfactor[1] + \ + gia.S11[t] + atm.S11[t] + remove.S11[t] # remove mean of each solution for iteration iteration.C10[:,n_iter] -= iteration.C10[:,n_iter].mean() iteration.C11[:,n_iter] -= iteration.C11[:,n_iter].mean() @@ -866,6 +913,8 @@ def monte_carlo_degree_one(base_dir, PROC, DREL, LMAX, RAD, nc[key].setncattr(att_name, att_val) # define global attributes + for att_name, att_val in attributes.items(): + fileID.setncattr(att_name, att_val) fileID.date_created = time.strftime('%Y-%m-%d',time.localtime()) # close the output file fileID.close() @@ -1311,6 +1360,19 @@ def arguments(): parser.add_argument('--mean-format', type=str, default='netCDF4', choices=['ascii','netCDF4','HDF5','gfc'], help='Input data format for GRACE/GRACE-FO mean file') + # monthly files to be removed from the GRACE/GRACE-FO data + parser.add_argument('--remove-file', + type=pathlib.Path, nargs='+', + help='Monthly files to be removed from the GRACE/GRACE-FO data') + choices = [] + choices.extend(['ascii','netCDF4','HDF5']) + choices.extend(['index-ascii','index-netCDF4','index-HDF5']) + parser.add_argument('--remove-format', + type=str, nargs='+', choices=choices, + help='Input data format for files to be removed') + parser.add_argument('--redistribute-removed', + default=False, action='store_true', + help='Redistribute removed mass fields over the ocean') # additional error files to be used in the monte carlo run parser.add_argument('--error-file', type=pathlib.Path, @@ -1396,6 +1458,9 @@ def main(): DATAFORM=args.format, MEAN_FILE=args.mean_file, MEANFORM=args.mean_format, + REMOVE_FILES=args.remove_file, + REMOVE_FORMAT=args.remove_format, + REDISTRIBUTE_REMOVED=args.redistribute_removed, ERROR_FILES=args.error_file, SOLVER=args.solver, FINGERPRINT=args.fingerprint, diff --git a/gravity_toolkit/SLR/C20.py b/gravity_toolkit/SLR/C20.py index 992018f6..48369361 100644 --- a/gravity_toolkit/SLR/C20.py +++ b/gravity_toolkit/SLR/C20.py @@ -330,7 +330,7 @@ def C20(SLR_file, AOD=True, HEADER=True): YY,MM,DD,hh,mm,ss = gravity_toolkit.time.convert_julian( MJD+2400000.5, format='tuple') # converting from month, day, year into decimal year - dinput['time'][t] = gravity_toolkit.time.convert_calendar_decimal( + dinput['time'][t], = gravity_toolkit.time.convert_calendar_decimal( YY, MM, day=DD, hour=hh) # Spherical Harmonic data for line dinput['data'][t] = np.float64(line_contents[2]) @@ -390,7 +390,7 @@ def C20(SLR_file, AOD=True, HEADER=True): YY,MM,DD,hh,mm,ss = gravity_toolkit.time.convert_julian( MJD+2400000.5, format='tuple') # converting from month, day, year into decimal year - date_conv[t] = gravity_toolkit.time.convert_calendar_decimal( + date_conv[t], = gravity_toolkit.time.convert_calendar_decimal( YY, MM, day=DD, hour=hh) # Spherical Harmonic data for line C20_input[t] = np.float64(line_contents[2]) diff --git a/gravity_toolkit/SLR/C30.py b/gravity_toolkit/SLR/C30.py index 84c2dda4..3cd9d3b7 100644 --- a/gravity_toolkit/SLR/C30.py +++ b/gravity_toolkit/SLR/C30.py @@ -161,7 +161,7 @@ def C30(SLR_file, C30_MEAN=9.5717395773300e-07, HEADER=True): YY,MM,DD,hh,mm,ss = gravity_toolkit.time.convert_julian( MJD+2400000.5, format='tuple') # converting from month, day, year into decimal year - dinput['time'][t] = gravity_toolkit.time.convert_calendar_decimal( + dinput['time'][t], = gravity_toolkit.time.convert_calendar_decimal( YY, MM, day=DD, hour=hh) # Spherical Harmonic data for line dinput['data'][t] = np.float64(line_contents[5]) diff --git a/gravity_toolkit/SLR/C50.py b/gravity_toolkit/SLR/C50.py index af93dada..53dff904 100644 --- a/gravity_toolkit/SLR/C50.py +++ b/gravity_toolkit/SLR/C50.py @@ -143,7 +143,7 @@ def C50(SLR_file, C50_MEAN=0.0, DATE=None, HEADER=True): YY,MM,DD,hh,mm,ss = gravity_toolkit.time.convert_julian( MJD+2400000.5, format='tuple') # converting from month, day, year into decimal year - dinput['time'][t] = gravity_toolkit.time.convert_calendar_decimal( + dinput['time'][t], = gravity_toolkit.time.convert_calendar_decimal( YY, MM, day=DD, hour=hh) # Spherical Harmonic data for line dinput['data'][t] = np.float64(line_contents[10]) diff --git a/gravity_toolkit/clenshaw_summation.py b/gravity_toolkit/clenshaw_summation.py index c7bbeaff..4315d322 100644 --- a/gravity_toolkit/clenshaw_summation.py +++ b/gravity_toolkit/clenshaw_summation.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" clenshaw_summation.py -Written by Tyler Sutterley (04/2023) +Written by Tyler Sutterley (07/2026) Calculates the spatial field for a series of spherical harmonics for a sequence of ungridded points @@ -49,6 +49,8 @@ Bollettino di Geodesia e Scienze (1982) UPDATE HISTORY: + 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 Updated 03/2023: improve typing for variables in docstrings Updated 02/2023: set custom units as top option in if/else statements @@ -123,8 +125,8 @@ def clenshaw_summation(clm, slm, lon, lat, raise ValueError('Incompatible vector dimensions (lon, lat)') # calculate colatitude and longitude in radians - th = (90.0 - lat)*np.pi/180.0 - phi = np.squeeze(lon*np.pi/180.0) + th = np.radians(90.0 - lat) + phi = np.squeeze(np.radians(lon)) # calculate cos and sin of colatitudes t = np.cos(th) u = np.sin(th) @@ -134,7 +136,7 @@ def clenshaw_summation(clm, slm, lon, lat, # Gaussian Smoothing if (RAD != 0): - wl = 2.0*np.pi*gauss_weights(RAD,LMAX) + wl = 2.0*np.pi*gauss_weights(RAD, LMAX) else: # else = 1 wl = np.ones((LMAX+1)) @@ -154,43 +156,40 @@ def clenshaw_summation(clm, slm, lon, lat, else: raise ValueError(f'Unknown units {UNITS}') + # complex spherical harmonics + ylm = clm - 1j * slm + # smooth degree dependent factors + f = dfactor*wl # calculate arrays for clenshaw summations over colatitudes - s_m_c = np.zeros((npts,LMAX*2+2)) + cs_m = np.zeros((npts, LMAX+1), dtype=np.clongdouble) for m in range(LMAX, -1, -1): # convolve harmonics with unit factors and smoothing - s_m_c[:,2*m:2*m+2] = clenshaw_s_m(t, dfactor*wl, m, clm, slm, - LMAX, ASTYPE=ASTYPE, SCALE=SCALE) + cs_m[:, m] = _clenshaw(t, f, m, ylm, LMAX, SCALE=SCALE) # calculate cos(phi) cos_phi_2 = 2.0*np.cos(phi) - # matrix of cos/sin m*phi summation - cos_m_phi = np.zeros((npts,LMAX+2),dtype=ASTYPE) - sin_m_phi = np.zeros((npts,LMAX+2),dtype=ASTYPE) + # matrix of cos/sin m*phi summation (Euler's form) + m_phi = np.zeros((npts, LMAX+2), dtype=np.clongdouble) # initialize matrix with values at lmax+1 and lmax - cos_m_phi[:,LMAX+1] = np.cos(ASTYPE(LMAX + 1)*phi) - sin_m_phi[:,LMAX+1] = np.sin(ASTYPE(LMAX + 1)*phi) - cos_m_phi[:,LMAX] = np.cos(ASTYPE(LMAX)*phi) - sin_m_phi[:,LMAX] = np.sin(ASTYPE(LMAX)*phi) + m_phi[:,LMAX+1] = np.exp(1j * (LMAX + 1) * phi) + m_phi[:,LMAX] = np.exp(1j * LMAX*phi) # calculate summation for order LMAX - s_m = s_m_c[:,2*LMAX]*cos_m_phi[:,LMAX] + s_m_c[:,2*LMAX+1]*sin_m_phi[:,LMAX] + s_m = (cs_m[:,LMAX]*m_phi[:,LMAX]).real # iterate to calculate complete summation for m in range(LMAX-1, 0, -1): - cos_m_phi[:,m] = cos_phi_2*cos_m_phi[:,m+1] - cos_m_phi[:,m+2] - sin_m_phi[:,m] = cos_phi_2*sin_m_phi[:,m+1] - sin_m_phi[:,m+2] # calculate summation for order m - a_m = np.sqrt((2.0*m+3.0)/(2.0*m+2.0)) - s_m = a_m*u*s_m + s_m_c[:,2*m]*cos_m_phi[:,m] + s_m_c[:,2*m+1]*sin_m_phi[:,m] + m_phi[:,m] = cos_phi_2*m_phi[:,m+1] - m_phi[:,m+2] + a_m = np.sqrt((2.0*m + 3.0)/(2.0*m + 2.0)) + # update summation and discard imaginary component + s_m = a_m*u*s_m + (cs_m[:,m]*m_phi[:,m]).real # calculate spatial field - spatial = np.sqrt(3.0)*u*s_m + s_m_c[:,0] + spatial = np.sqrt(3.0)*u*s_m + cs_m[:,0].real # return the calculated spatial field return spatial -# PURPOSE: compute conditioned arrays for Clenshaw summation from the -# fully-normalized associated Legendre's function for an order m -def clenshaw_s_m(t, f, m, clm1, slm1, lmax, - ASTYPE=np.longdouble, - SCALE=1e-280 - ): +# PURPOSE: compute Clenshaw summation of the fully normalized associated +# Legendre's function for constant order m +def _clenshaw(t, f, m, Ylm1, lmax, SCALE=1e-280): """ Compute conditioned arrays for Clenshaw summation from the fully-normalized associated Legendre's function for an order m @@ -203,14 +202,10 @@ def clenshaw_s_m(t, f, m, clm1, slm1, lmax, degree dependent factors m: int spherical harmonic order - clm1: np.ndarray - cosine spherical harmonics - slm1: np.ndarray - sine spherical harmonics + Ylm1: np.ndarray + complex form of spherical harmonics lmax: int maximum spherical harmonic degree - ASTYPE: np.dtype, default np.longdouble - floating point precision for calculating Clenshaw summation SCALE: float, default 1e-280 scaling factor to prevent underflow in Clenshaw summation @@ -221,49 +216,40 @@ def clenshaw_s_m(t, f, m, clm1, slm1, lmax, """ # allocate for output matrix N = len(t) - s_m = np.zeros((N,2),dtype=ASTYPE) + s_m = np.zeros((N), dtype=np.clongdouble) # scaling to prevent overflow - clm = SCALE*clm1.astype(ASTYPE) - slm = SCALE*slm1.astype(ASTYPE) + ylm = SCALE*Ylm1.astype(np.clongdouble) # convert lmax and m to float - lm = ASTYPE(lmax) - mm = ASTYPE(m) + lm = np.float64(lmax) + mm = np.float64(m) if (m == lmax): - s_m[:,0] = f[lmax]*clm[lmax,lmax] - s_m[:,1] = f[lmax]*slm[lmax,lmax] + s_m[:] = f[lmax]*ylm[lmax,lmax] elif (m == (lmax-1)): a_lm = np.sqrt(((2.0*lm-1.0)*(2.0*lm+1.0))/((lm-mm)*(lm+mm)))*t - s_m[:,0] = a_lm*f[lmax]*clm[lmax,lmax-1] + f[lmax-1]*clm[lmax-1,lmax-1] - s_m[:,1] = a_lm*f[lmax]*slm[lmax,lmax-1] + f[lmax-1]*slm[lmax-1,lmax-1] + s_m[:] = a_lm*f[lmax]*ylm[lmax,lmax-1] + f[lmax-1]*ylm[lmax-1,lmax-1] elif ((m <= (lmax-2)) and (m >= 1)): - s_mm_c_pre_2 = f[lmax]*clm[lmax,m] - s_mm_s_pre_2 = f[lmax]*slm[lmax,m] + s_mm_minus_2 = f[lmax]*ylm[lmax,m] a_lm = np.sqrt(((2.0*lm-1.0)*(2.0*lm+1.0))/((lm-mm)*(lm+mm)))*t - s_mm_c_pre_1 = a_lm*s_mm_c_pre_2 + f[lmax-1]*clm[lmax-1,m] - s_mm_s_pre_1 = a_lm*s_mm_s_pre_2 + f[lmax-1]*slm[lmax-1,m] + s_mm_minus_1 = a_lm*s_mm_minus_2 + f[lmax-1]*ylm[lmax-1,m] for l in range(lmax-2, m-1, -1): - ll = ASTYPE(l) + ll = np.float64(l) a_lm=np.sqrt(((2.0*ll+1.0)*(2.0*ll+3.0))/((ll+1.0-mm)*(ll+1.0+mm)))*t b_lm=np.sqrt(((2.*ll+5.)*(ll+mm+1.)*(ll-mm+1.))/((ll+2.-mm)*(ll+2.+mm)*(2.*ll+1.))) - s_mm_c = a_lm * s_mm_c_pre_1 - b_lm * s_mm_c_pre_2 + f[l]*clm[l,m] - s_mm_s = a_lm * s_mm_s_pre_1 - b_lm * s_mm_s_pre_2 + f[l]*slm[l,m] - s_mm_c_pre_2 = np.copy(s_mm_c_pre_1) - s_mm_s_pre_2 = np.copy(s_mm_s_pre_1) - s_mm_c_pre_1 = np.copy(s_mm_c) - s_mm_s_pre_1 = np.copy(s_mm_s) - s_m[:,0] = np.copy(s_mm_c) - s_m[:,1] = np.copy(s_mm_s) + s_mm_l = a_lm * s_mm_minus_1 - b_lm * s_mm_minus_2 + f[l]*ylm[l,m] + s_mm_minus_2 = np.copy(s_mm_minus_1) + s_mm_minus_1 = np.copy(s_mm_l) + s_m[:] = np.copy(s_mm_l) elif (m == 0): - s_mm_c_pre_2 = f[lmax]*clm[lmax,0] + s_mm_minus_2 = f[lmax]*ylm[lmax,0] a_lm = np.sqrt(((2.0*lm-1.0)*(2.0*lm+1.0))/(lm*lm))*t - s_mm_c_pre_1 = a_lm * s_mm_c_pre_2 + f[lmax-1]*clm[lmax-1,0] + s_mm_minus_1 = a_lm * s_mm_minus_2 + f[lmax-1]*ylm[lmax-1,0] for l in range(lmax-2, m-1, -1): - ll = ASTYPE(l) + ll = np.float64(l) a_lm=np.sqrt(((2.0*ll+1.0)*(2.0*ll+3.0))/((ll+1.0)*(ll+1.0)))*t b_lm=np.sqrt(((2.0*ll+5.0)*(ll+1.0)*(ll+1.0))/((ll+2.0)*(ll+2.0)*(2.0*ll+1.0))) - s_mm_c = a_lm * s_mm_c_pre_1 - b_lm * s_mm_c_pre_2 + f[l]*clm[l,0] - s_mm_c_pre_2 = np.copy(s_mm_c_pre_1) - s_mm_c_pre_1 = np.copy(s_mm_c) - s_m[:,0] = np.copy(s_mm_c) - # return s_m rescaled with scalef + s_mm_l = a_lm * s_mm_minus_1 - b_lm * s_mm_minus_2 + f[l]*ylm[l,0] + s_mm_minus_2 = np.copy(s_mm_minus_1) + s_mm_minus_1 = np.copy(s_mm_l) + s_m[:] = np.copy(s_mm_l) + # return rescaled s_m return s_m/SCALE diff --git a/gravity_toolkit/fourier_legendre.py b/gravity_toolkit/fourier_legendre.py index 6cc8b5cd..44f60826 100755 --- a/gravity_toolkit/fourier_legendre.py +++ b/gravity_toolkit/fourier_legendre.py @@ -2,24 +2,25 @@ u""" fourier_legendre.py Original IDL code gen_plms.pro written by Sean Swenson -Adapted by Tyler Sutterley (03/2023) +Adapted by Tyler Sutterley (07/2026) Computes Fourier coefficients of the associated Legendre functions CALLING SEQUENCE: - plm = fourier_legendre(lmax,mmax) + Almk = fourier_legendre(lmax,mmax) INPUTS: lmax: maximum spherical harmonic degree mmax: maximum spherical harmonic order OUTPUTS: - plm: Fourier coefficients + Almk: Fourier coefficients PYTHON DEPENDENCIES: numpy: Scientific Computing Tools For Python (https://numpy.org) UPDATE HISTORY: + Updated 07/2027: add citations to docstrings Updated 03/2023: improve typing for variables in docstrings Updated 10/2022: add polynomial function for calculating gradients Updated 04/2022: updated docstrings to numpy documentation format @@ -34,6 +35,7 @@ def fourier_legendre(lmax, mmax): """ Computes Fourier coefficients of the associated Legendre functions + :cite:p:`Hofsommer:1960wg,Gruber:2016hn` Parameters ---------- @@ -44,12 +46,12 @@ def fourier_legendre(lmax, mmax): Returns ------- - plm: np.ndarray + Almk: np.ndarray Fourier coefficients """ # allocate for output fourier coefficients - plm = np.zeros((lmax+1,lmax+1,lmax+1)) + Almk = np.zeros((lmax+1,lmax+1,lmax+1)) l_even = np.arange(0,lmax+1,2) l_odd = np.arange(1,lmax,2) m_even = np.arange(0,mmax+1,2) @@ -57,100 +59,100 @@ def fourier_legendre(lmax, mmax): # First compute m=0, m=1 terms # Compute m = 0, l = even terms - plm[l_even,0,0] = 1.0 - p1 = (l_even*(l_even+1.0))*plm[l_even,0,0] - plm[l_even,0,2] = p1 / (l_even*(l_even+1.0)-2.0) + Almk[l_even,0,0] = 1.0 + a1 = (l_even*(l_even+1.0))*Almk[l_even,0,0] + Almk[l_even,0,2] = a1 / (l_even*(l_even+1.0)-2.0) for j in range(2,lmax,2):# equivalent to 2:lmax-2 - p1 = 2.0*(l_even*(l_even+1.0)-j**2.0)*plm[l_even,0,j] - p2 = ((j-2.0)*(j-1.0)-l_even*(l_even+1.0))*plm[l_even,0,j-2] + a1 = 2.0*(l_even*(l_even+1.0)-j**2.0)*Almk[l_even,0,j] + a2 = ((j-2.0)*(j-1.0)-l_even*(l_even+1.0))*Almk[l_even,0,j-2] dfactor = (l_even*(l_even+1.0)-(j+2.0)*(j+1.0)) - plm[l_even,0,j+2] = (p1 + p2) / dfactor + Almk[l_even,0,j+2] = (a1 + a2) / dfactor # Special case for j = 0 fourier coefficient - plm[l_even,0,0] = plm[l_even,0,0]/2.0 + Almk[l_even,0,0] = Almk[l_even,0,0]/2.0 # Normalize overall sum to 2 for m == 0 norm = np.zeros((len(l_even))) for j in range(0,lmax+2,2):# equivalent to 0:lmax - ptemp = np.squeeze(plm[l_even[:, np.newaxis],0,m_even]) + ptemp = np.squeeze(Almk[l_even[:, np.newaxis],0,m_even]) dtemp = 1.0/(1.0-j-m_even) + 1.0/(1.0+j-m_even) + \ 1.0/(1.0-j+m_even) + 1.0/(1.0+j+m_even) - norm[l_even//2] = norm[l_even//2] + plm[l_even,0,j] * \ + norm[l_even//2] = norm[l_even//2] + Almk[l_even,0,j] * \ np.dot(ptemp, dtemp)/2.0 - # normalize plms + # normalize Almks norm = np.sqrt(norm/2.0) for l in range(0,lmax+2,2):# equivalent to 0:lmax - plm[l,0,:] = plm[l,0,:]/norm[l//2] + Almk[l,0,:] = Almk[l,0,:]/norm[l//2] # Compute m = 0, l = odd terms - plm[l_odd,0,1] = 1.0 - p1 = (2.0-l_odd*(l_odd+1.0))*plm[l_odd,0,1] - plm[l_odd,0,3] = p1 / (6.0-l_odd*(l_odd+1.0)) + Almk[l_odd,0,1] = 1.0 + a1 = (2.0-l_odd*(l_odd+1.0))*Almk[l_odd,0,1] + Almk[l_odd,0,3] = a1 / (6.0-l_odd*(l_odd+1.0)) for j in range(3,lmax-1,2):# equivalent to 3:lmax-3 - p1 = 2.0*(l_odd*(l_odd+1.0)-j**2.0)*plm[l_odd,0,j] - p2 = ((j-2.0)*(j-1.0)-l_odd*(l_odd+1.0))*plm[l_odd,0,j-2] + a1 = 2.0*(l_odd*(l_odd+1.0)-j**2.0)*Almk[l_odd,0,j] + a2 = ((j-2.0)*(j-1.0)-l_odd*(l_odd+1.0))*Almk[l_odd,0,j-2] dfactor = (l_odd*(l_odd+1.0)-(j+2.0)*(j+1.0)) - plm[l_odd,0,j+2] = (p1 + p2) / dfactor + Almk[l_odd,0,j+2] = (a1 + a2) / dfactor # Normalize overall sum to 2 for m == 0 norm = np.zeros((len(l_odd))) for j in range(1,lmax+1,2):# equivalent to 1:lmax-1 - ptemp = np.squeeze(plm[l_odd[:, np.newaxis],0,m_odd]) + ptemp = np.squeeze(Almk[l_odd[:, np.newaxis],0,m_odd]) dtemp = 1.0/(1.0-j-m_odd) + 1.0/(1.0+j-m_odd) + \ 1.0/(1.0-j+m_odd) + 1.0/(1.0+j+m_odd) - norm[(l_odd-1)//2] = norm[(l_odd-1)//2] + plm[l_odd,0,j] * \ + norm[(l_odd-1)//2] = norm[(l_odd-1)//2] + Almk[l_odd,0,j] * \ np.dot(ptemp, dtemp)/2.0 - # normalize plms + # normalize Almks norm = np.sqrt(norm/2.0) for l in range(1,lmax+1,2):# equivalent to 1:lmax-1 - plm[l,0,:] = plm[l,0,:]/norm[(l-1)//2] + Almk[l,0,:] = Almk[l,0,:]/norm[(l-1)//2] # Compute m = 1, l = even terms - plm[l_even,1,0] = 0.0 - plm[l_even,1,2] = 1.0 + Almk[l_even,1,0] = 0.0 + Almk[l_even,1,2] = 1.0 for j in range(2,lmax,2):# equivalent to 2:lmax-2 - p1 = 2.0*(l_even*(l_even+1)-j**2.0-2.0)*plm[l_even,1,j] - p2 = ((j-2.0)*(j-1.0)-l_even*(l_even+1))*plm[l_even,1,j-2] + a1 = 2.0*(l_even*(l_even+1)-j**2.0-2.0)*Almk[l_even,1,j] + a2 = ((j-2.0)*(j-1.0)-l_even*(l_even+1))*Almk[l_even,1,j-2] dfactor = (l_even*(l_even+1.0)-(j+2.0)*(j+1.0)) - plm[l_even,1,j+2] = (p1 + p2) / dfactor + Almk[l_even,1,j+2] = (a1 + a2) / dfactor # Normalize overall sum to 4 for m == 1 # different norm than that of the cosine series norm = np.zeros((len(l_even))) for j in range(0,lmax+2,2):# equivalent to 0:lmax - ptemp = np.squeeze(plm[l_even[:, np.newaxis],1,m_even]) + ptemp = np.squeeze(Almk[l_even[:, np.newaxis],1,m_even]) dtemp = -1.0/(1.0-j-m_even) + 1.0/(1+j-m_even) + \ 1.0/(1.0-j+m_even) - 1.0/(1+j+m_even) - norm[l_even//2] = norm[l_even//2] + plm[l_even,1,j] * \ + norm[l_even//2] = norm[l_even//2] + Almk[l_even,1,j] * \ np.dot(ptemp, dtemp)/2.0 - # normalize plms + # normalize Almks norm = np.sqrt(norm/4.0) for l in range(0,lmax+2,2):# equivalent to 0:lmax - plm[l,1,:] = plm[l,1,:]/norm[l//2] + Almk[l,1,:] = Almk[l,1,:]/norm[l//2] # Compute m = 1, l = odd terms - plm[l_odd,1,1] = 1.0 - plm[l_odd,1,3] = 3.0*(l_odd*(l_odd+1)-2)*plm[l_odd,1,1]/(l_odd*(l_odd+1)-6) + Almk[l_odd,1,1] = 1.0 + Almk[l_odd,1,3] = 3.0*(l_odd*(l_odd+1)-2)*Almk[l_odd,1,1]/(l_odd*(l_odd+1)-6) for j in range(3,lmax-1,2):# equivalent to 3:lmax-3 - p1 = 2.0*(l_odd*(l_odd+1.0)-j**2.0-2.0)*plm[l_odd,1,j] - p2 = ((j-2.0)*(j-1.0)-l_odd*(l_odd+1.0))*plm[l_odd,1,j-2] + a1 = 2.0*(l_odd*(l_odd+1.0)-j**2.0-2.0)*Almk[l_odd,1,j] + a2 = ((j-2.0)*(j-1.0)-l_odd*(l_odd+1.0))*Almk[l_odd,1,j-2] dfactor = (l_odd*(l_odd+1.0)-(j+2.0)*(j+1.0)) - plm[l_odd,1,j+2] = (p1 + p2) / dfactor + Almk[l_odd,1,j+2] = (a1 + a2) / dfactor # Normalize overall sum to 4 for m == 1 norm = np.zeros((len(l_odd))) for j in range(1,lmax+1,2):# equivalent to 1:lmax-1 - ptemp = np.squeeze(plm[l_odd[:, np.newaxis],1,m_odd]) + ptemp = np.squeeze(Almk[l_odd[:, np.newaxis],1,m_odd]) dtemp = -1.0/(1.0-j-m_odd) + 1.0/(1.0+j-m_odd) + \ 1.0/(1.0-j+m_odd) - 1.0/(1.0+j+m_odd) - norm[(l_odd-1)//2] = norm[(l_odd-1)//2] + plm[l_odd,1,j] * \ + norm[(l_odd-1)//2] = norm[(l_odd-1)//2] + Almk[l_odd,1,j] * \ np.dot(ptemp, dtemp)/2.0 - # normalize plms + # normalize Almks norm = np.sqrt(norm/4.0) for l in range(1,lmax+1,2):# equivalent to 1:lmax-1 - plm[l,1,:] = plm[l,1,:]/norm[(l-1)//2] + Almk[l,1,:] = Almk[l,1,:]/norm[(l-1)//2] # Compute coefficients for m > 0 @@ -158,58 +160,58 @@ def fourier_legendre(lmax, mmax): m = 0 # m = 0, l = even terms for l in range(m,lmax-1):# equivalent to m:lmax-2 - p1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*plm[l,m,m_even] - p2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*plm[l+2,m,m_even] - p3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0)/2.0)*plm[l,m+2,m_even] + a1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*Almk[l,m,m_even] + a2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*Almk[l+2,m,m_even] + a3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0)/2.0)*Almk[l,m+2,m_even] dfactor = np.sqrt((l+m+4.0)*(l+m+3.0)/(2.0*l+5.0)/2.0) - plm[l+2,m+2,m_even] = (p1 - p2 + p3) / dfactor + Almk[l+2,m+2,m_even] = (a1 - a2 + a3) / dfactor # m = 0, l = odd terms for l in range(m+1,lmax-1):# equivalent to m+1:lmax-2 - p1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*plm[l,m,m_odd] - p2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*plm[l+2,m,m_odd] - p3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0)/2.0)*plm[l,m+2,m_odd] + a1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*Almk[l,m,m_odd] + a2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*Almk[l+2,m,m_odd] + a3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0)/2.0)*Almk[l,m+2,m_odd] dfactor = np.sqrt((l+m+4.0)*(l+m+3.0)/(2.0*l+5.0)/2.0) - plm[l+2,m+2,m_odd] = (p1 - p2 + p3) / dfactor + Almk[l+2,m+2,m_odd] = (a1 - a2 + a3) / dfactor # m = even terms for m in range(2,lmax,2):# equivalent to 2:lmax-2 # m = even, > 2, l = even terms for l in range(m,lmax,2):# equivalent to m:lmax-2 - p1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*plm[l,m,m_even] - p2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*plm[l+2,m,m_even] - p3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*plm[l,m+2,m_even] + a1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*Almk[l,m,m_even] + a2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*Almk[l+2,m,m_even] + a3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*Almk[l,m+2,m_even] dfactor = np.sqrt((l+m+4.0)*(l+m+3.0)/(2.0*l+5.0)) - plm[l+2,m+2,m_even] = (p1 - p2 + p3) / dfactor + Almk[l+2,m+2,m_even] = (a1 - a2 + a3) / dfactor # m = even, > 2, l = odd terms for l in range(m+1,lmax-1,2): - p1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*plm[l,m,m_odd] - p2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*plm[l+2,m,m_odd] - p3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*plm[l,m+2,m_odd] + a1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*Almk[l,m,m_odd] + a2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*Almk[l+2,m,m_odd] + a3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*Almk[l,m+2,m_odd] dfactor = np.sqrt((l+m+4.0)*(l+m+3.0)/(2.0*l+5.0)) - plm[l+2,m+2,m_odd] = (p1 - p2 + p3) / dfactor + Almk[l+2,m+2,m_odd] = (a1 - a2 + a3) / dfactor # m = odd terms for m in range(1,lmax-1,2):# equivalent to 1:lmax-3 # m = odd, > 1, l = even terms for l in range(m+1,lmax-1,2):# equivalent to m+1,lmax-2 - p1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*plm[l,m,m_even] - p2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*plm[l+2,m,m_even] - p3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*plm[l,m+2,m_even] + a1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*Almk[l,m,m_even] + a2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*Almk[l+2,m,m_even] + a3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*Almk[l,m+2,m_even] dfactor = np.sqrt((l+m+4.0)*(l+m+3.0)/(2.0*l+5.0)) - plm[l+2,m+2,m_even] = (p1 - p2 + p3) / dfactor + Almk[l+2,m+2,m_even] = (a1 - a2 + a3) / dfactor # m = odd, > 1, l = odd terms for l in range(m,lmax-1,2):# equivalent to m:lmax-2 - p1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*plm[l,m,m_odd] - p2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*plm[l+2,m,m_odd] - p3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*plm[l,m+2,m_odd] + a1 = np.sqrt((l+m+2.0)*(l+m+1.0)/(2.0*l+1.0))*Almk[l,m,m_odd] + a2 = np.sqrt((l-m+1.0)*(l-m+2.0)/(2.0*l+5.0))*Almk[l+2,m,m_odd] + a3 = np.sqrt((l-m)*(l-m-1.0)/(2.0*l+1.0))*Almk[l,m+2,m_odd] dfactor = np.sqrt((l+m+4.0)*(l+m+3.0)/(2.0*l+5.0)) - plm[l+2,m+2,m_odd] = (p1 - p2 + p3) / dfactor + Almk[l+2,m+2,m_odd] = (a1 - a2 + a3) / dfactor # return the fourier coefficients - return plm + return Almk def legendre_gradient(lmax, mmax): """ @@ -225,67 +227,40 @@ def legendre_gradient(lmax, mmax): Returns ------- - vlm: np.ndarray + Vlmk: np.ndarray Fourier coefficients for meridional gradients - wlm: np.ndarray + Wlmk: np.ndarray Fourier coefficients for zonal gradients """ - - plm = fourier_legendre(lmax, mmax) - vlm = np.zeros((lmax+1,lmax+1,lmax+1)) - wlm = np.zeros((lmax+1,lmax+1,lmax+1)) - - # l=0 zero by definition - lind = np.arange(1,lmax+1) - # m=0 special case - # terms with m=0, m=1 have different coefficients - vlm[lind,0,:] = 2.0*np.dot(np.diag(np.sqrt((lind+1)*lind/2.0)), plm[lind,1,:]) - - # m+1 terms - for l in range(2,lmax+1):# from 2 to lmax - m = np.arange(1,l)# from 1 to l-1 - lplus = np.arange(l+2,2*l+1)# from l+2 to 2*l - lminus = np.arange(l-1,0,-1)# from l-1 to 1 - vlm[l,m,:] = np.dot(np.diag(np.sqrt(lplus*lminus/4.0)), plm[l,m+1,:]) - - # m-1 terms, m-1=0 has different coefficients - vlm[lind,1,:] -= np.dot(np.diag(np.sqrt((lind+1)*lind/2.0)), plm[lind,0,:]) - - for l in range(2,lmax+1): - m = np.arange(2,l+1)# from 2 to l - lplus = np.arange(l+2,2*l+1)# from l+2 to 2*l - lminus = np.arange(l-1,0,-1)# from l-1 to 1 - vlm[l,m,:] -= np.dot(np.diag(np.sqrt(lplus*lminus/4.0)), plm[l,m-1,:]) - # normalizations - for l in range(1,lmax+1): - vlm[l,:,:] /= np.sqrt((l+1)*l) - - # m+1 terms - for l in range(2, lmax+1): - m = np.arange(1,l)# from 1 to l-1 - dfactor = (2.0*l+1.0)/(2.0*l-1.0) - lminus2 = np.arange(l-2,-1,-1)# from l-2 to 0 - lminus1 = np.arange(l-1,0,-1)# from l-1 to 1 - wlm[l,m,:] = np.sqrt(dfactor) * \ - np.dot(np.diag(np.sqrt(lminus2*lminus1/4.0)), plm[l,m+1,:]) - - # m-1 terms, m-1=0 has different coefficients - # m=1 term - for l in range(1, lmax+1): - dfactor = (2.0*l+1.0)/(2.0*l-1.0) - wlm[l,1,:] += np.sqrt(dfactor)*np.sqrt(l*(l+1)/2.0)*plm[l-1,0,:] - - for l in range(2,lmax+1): - m = np.arange(2,l+1) - dfactor = (2.0*l+1.0)/(2.0*l-1.0) - lplus2 = np.arange(l+2,2*l+1)# from l+2 to 2*l - lplus1 = np.arange(l+1,2*l)# from l+1 to (2*l-1) - wlm[l,m,:] += np.sqrt(dfactor) * \ - np.dot(np.diag(np.sqrt(lplus2*lplus1)/4.0), plm[l-1,m-1,:]) - # normalizations + # compute the fourier coefficients of the associated legendre functions + Almk = fourier_legendre(lmax, mmax) + # allocate for output fourier coefficients + Vlmk = np.zeros((lmax+1,lmax+1,lmax+1)) + Wlmk = np.zeros((lmax+1,lmax+1,lmax+1)) + # for each spherical harmonic degree for l in range(1, lmax+1): - wlm[l,:,:] /= np.sqrt((l+1)*l) - # normalize vlm - vlm[:,0,:] /= 2.0 - - return (vlm, wlm) + # degree dependent factor + dfactor = np.sqrt((2.0*l + 1.0)/(2.0*l - 1.0)) + # m=0 special case + Vfact = np.sqrt(l*(l + 1.0)/2.0) + Vlmk[l,0,:] = Vfact * Almk[l,1,:] + for m in range(2, l+1):# from 2 to l + Vfact = np.sqrt((l + m) * (l - m + 1.0) / 4.0) + Wfact = dfactor * np.sqrt((l - m) * (l - m + 1) / 4.0) + Vlmk[l,m-1,:] = Vfact * Almk[l,m,:] + Wlmk[l,m-1,:] = -Wfact * Almk[l-1,m,:] + # m = 1 terms + Vfact = np.sqrt(l*(l + 1.0)/2.0) + Wfact = dfactor * np.sqrt(l*(l + 1.0)/2.0) + Vlmk[l,1,:] -= Vfact*Almk[l,0,:] + Wlmk[l,1,:] += dfactor*Wfact*Almk[l-1,0,:] + for m in range(2, l + 1):# from 2 to l + Vfact = np.sqrt((l + m) * (l - m + 1.0) / 4.0) + Wfact = dfactor * np.sqrt((l + m) * (l + m - 1) / 4.0) + Vlmk[l,m,:] -= Vfact * Almk[l,m-1,:] + Wlmk[l,m,:] += Wfact * Almk[l-1,m-1,:] + # normalizations + Vlmk[l,:,:] /= np.sqrt(l * (l + 1.0)) + Wlmk[l,:,:] /= np.sqrt(l * (l + 1.0)) + # return the coefficients + return (Vlmk, Wlmk) diff --git a/gravity_toolkit/gen_disc_load.py b/gravity_toolkit/gen_disc_load.py index 474b1b3f..ca8de548 100644 --- a/gravity_toolkit/gen_disc_load.py +++ b/gravity_toolkit/gen_disc_load.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" gen_disc_load.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (07/2026) Calculates gravitational spherical harmonic coefficients for a uniform disc load CALLING SEQUENCE: @@ -55,6 +55,8 @@ https://doi.org/10.1007/s00190-011-0522-7 UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 06/2023: modified custom units case to not convert to cmwe Updated 03/2023: simplified unit degree factors using units class improve typing for variables in docstrings @@ -129,8 +131,8 @@ def gen_disc_load(data, lon, lat, area, LMAX=60, MMAX=None, UNITS=2, MMAX = np.copy(LMAX) # convert lon and lat to radians - phi = lon*np.pi/180.0# Longitude in radians - th = (90.0 - lat)*np.pi/180.0# Colatitude in radians + phi = np.radians(lon)# Longitude in radians + th = np.radians(90.0 - lat)# Colatitude in radians # Earth Parameters factors = gravity_toolkit.units(lmax=LMAX) @@ -204,29 +206,20 @@ def gen_disc_load(data, lon, lat, area, LMAX=60, MMAX=None, UNITS=2, # data normally is 1 for a uniform 1cm water equivalent layer # but can be a mass point if reconstructing a spherical harmonic field # NOTE: NOT a matrix multiplication as data (and phi) is a single point - dcos = unit_conv*data*np.cos(m*phi) - dsin = unit_conv*data*np.sin(m*phi) + d = unit_conv*data*np.exp(1j*m*phi) # Multiplying by plm_alpha (F_l from Jacob 2012) plm = np.zeros((LMAX+1, MMAX+1)) - # Initializing preliminary spherical harmonic matrices - yclm = np.zeros((LMAX+1, MMAX+1)) - yslm = np.zeros((LMAX+1, MMAX+1)) # Initializing output spherical harmonic matrices Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX) - Ylms.clm = np.zeros((LMAX+1, MMAX+1)) - Ylms.slm = np.zeros((LMAX+1, MMAX+1)) - for m in range(0,MMAX+1):# MMAX+1 to include MMAX - l = np.arange(m,LMAX+1)# LMAX+1 to include LMAX - # rotate disc load to be centered at lat/lon - plm[l,m] = plmout[l,m]*pl_alpha[l] - # multiplying clm by cos(m*phi) and slm by sin(m*phi) - # to get a field of spherical harmonics - yclm[l,m] = plm[l,m]*dcos[m] - yslm[l,m] = plm[l,m]*dsin[m] - # multiplying by factors to convert to geoid coefficients - Ylms.clm[l,m] = dfactor[l]*yclm[l,m] - Ylms.slm[l,m] = dfactor[l]*yslm[l,m] + # rotate disc load to be centered at lat/lon + plm = np.einsum("lm...,l...->lm...", plmout, pl_alpha) + # multiplying clm by cos(m*phi) and slm by sin(m*phi) + # to get a field of spherical harmonics + ylm = np.einsum("lm...,m...->lm...", plm, d) + # Multiplying by factors to convert to fully normalized coefficients + Ylms.clm = np.einsum("l...,lm...->lm...", dfactor, ylm.real) + Ylms.slm = np.einsum("l...,lm...->lm...", dfactor, ylm.imag) # return the output spherical harmonics object return Ylms diff --git a/gravity_toolkit/gen_harmonics.py b/gravity_toolkit/gen_harmonics.py index 8044f695..f4d561cd 100644 --- a/gravity_toolkit/gen_harmonics.py +++ b/gravity_toolkit/gen_harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" gen_harmonics.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (07/2026) Converts data from the spatial domain to spherical harmonic coefficients Does not compute the solid Earth elastic response or convert units @@ -47,6 +47,8 @@ Associated Legendre Functions", Journal of Geodesy (2002) UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 03/2023: improve typing for variables in docstrings Updated 01/2023: refactored associated legendre polynomials Updated 04/2022: updated docstrings to numpy documentation format @@ -154,30 +156,20 @@ def integration(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs): m: int spherical harmonic order to MMAX """ - - # dimensions of the longitude and latitude arrays - nlon = np.int64(len(lon)) - nlat = np.int64(len(lat)) - # grid step - dlon = np.abs(lon[1]-lon[0]) - dlat = np.abs(lat[1]-lat[0]) - # longitude degree spacing in radians - dphi = dlon*np.pi/180.0 - # colatitude degree spacing in radians - dth = dlat*np.pi/180.0 - - # reformatting longitudes to range 0:360 (if previously -180:180) - if np.count_nonzero(lon < 0): - lon[lon < 0] += 360.0 # calculate longitude and colatitude arrays in radians - phi = np.reshape(lon,(1,nlon))*np.pi/180.0# reshape to 1xnlon - th = (90.0 - np.squeeze(lat))*np.pi/180.0# remove singleton dimensions + phi = np.radians(np.squeeze(lon)) + th = np.radians(90.0 - np.squeeze(lat)) + # reformatting longitudes to range 0:360 (if previously -180:180) + phi = np.where(phi < 0, phi + 2.0*np.pi, phi) + # grid step in radians + dphi = np.abs(phi[1] - phi[0]) + dth = np.abs(th[1] - th[0]) - # Calculating cos/sin of phi arrays (output [m,phi]) # LMAX+1 as there are LMAX+1 elements between 0 and LMAX - m = np.arange(MMAX+1)[:, np.newaxis] - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + ll = np.arange(LMAX+1) + mm = np.arange(MMAX+1) + # Calculating cos/sin of phi arrays (output [m,phi]) + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", mm, phi)) # Multiplying sin(th) with differentials of theta and phi # to calculate the integration factor at each latitude @@ -185,41 +177,24 @@ def integration(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs): coeff = 1.0/(4.0*np.pi) # Calculate polynomials using Holmes and Featherstone (2002) relation - plm = np.zeros((LMAX+1, MMAX+1, nlat)) if (np.ndim(PLM) == 0): - plmout,dplm = plm_holmes(LMAX, np.cos(th)) - else: - # use precomputed plms to improve computational speed - # or to use a different recursion relation for polynomials - plmout = PLM - + PLM, dplm = plm_holmes(LMAX, np.cos(th)) # Multiply plms by integration factors [sin(theta)*dtheta*dphi] # truncate plms to maximum spherical harmonic order if MMAX < LMAX - m = np.arange(MMAX+1) - for j in range(0,nlat): - plm[:,m,j] = plmout[:,m,j]*int_fact[j] - - # Initializing preliminary spherical harmonic matrices - yclm = np.zeros((LMAX+1, MMAX+1)) - yslm = np.zeros((LMAX+1, MMAX+1)) + plm = np.einsum("lmh...,h...->lmh...", PLM[:LMAX+1,:MMAX+1,:], int_fact) # Initializing output spherical harmonic matrices Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX) Ylms.clm = np.zeros((LMAX+1, MMAX+1)) Ylms.slm = np.zeros((LMAX+1, MMAX+1)) # Multiplying gridded data with sin/cos of m#phis (output [m,theta]) # This will sum through all phis in the dot product - dcos = np.dot(ccos,data) - dsin = np.dot(ssin,data) - for l in range(0,LMAX+1): - mm = np.min([MMAX,l])# truncate to MMAX if specified (if l > MMAX) - m = np.arange(0,mm+1)# mm+1 elements between 0 and mm - # Summing product of plms and data over all latitudes - yclm[l,m] = np.sum(plm[l,m,:]*dcos[m,:], axis=1) - yslm[l,m] = np.sum(plm[l,m,:]*dsin[m,:], axis=1) - # convert to output normalization (4-pi normalized harmonics) - Ylms.clm[l,m] = coeff*yclm[l,m] - Ylms.slm[l,m] = coeff*yslm[l,m] - + d = np.einsum("mp...,ph...->mh...", m_phi, data) + # Summing product of plms and data over all latitudes + ylm = np.einsum("lmh...,mh...->lm...", plm, d) + # convert to output normalization (4-pi normalized harmonics) + # truncate to MMAX if specified (if l > MMAX) + Ylms.clm = coeff*ylm.real[:LMAX+1, :MMAX+1] + Ylms.slm = coeff*ylm.imag[:LMAX+1, :MMAX+1] # return the output spherical harmonics object return Ylms @@ -257,139 +232,125 @@ def fourier(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs): # dimensions of the longitude and latitude arrays nlon = np.int64(len(lon)) nlat = np.int64(len(lat)) - # remove singleton dimensions and convert to radians - phi = (np.squeeze(lon)*np.pi/180.0) - # Colatitude in radians - theta = ((90.0 - np.squeeze(lat))*np.pi/180.0) + # calculate longitude and colatitude arrays in radians + phi = np.radians(np.squeeze(lon)) + th = np.radians(90.0 - np.squeeze(lat)) + # reformatting longitudes to range 0:360 (if previously -180:180) + phi = np.where(phi < 0, phi + 2.0*np.pi, phi) + # grid step in radians + dphi = np.abs(phi[1] - phi[0]) + dth = np.abs(th[1] - th[0]) # MMAX+1 to include MMAX - mm = np.arange(MMAX+1)[:, np.newaxis] + mm = np.arange(MMAX+1) # Calculate cos and sin coefficients of signal - ccos = np.cos(np.dot(mm,phi[np.newaxis,:])) - ssin = np.sin(np.dot(mm,phi[np.newaxis,:])) - dcos = np.dot(ccos,data) - dsin = np.dot(ssin,data) - - # Normalize fourier coefficients - dcos[0,:] = dcos[0,:]/nlon - dcos[1:MMAX+1,:] = 2.0*dcos[1:MMAX+1,:]/nlon - dsin[0,:] = dsin[0,:]/nlon - dsin[1:MMAX+1,:] = 2.0*dsin[1:MMAX+1,:]/nlon + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", mm, phi)) + d = np.einsum("mp...,ph...->mh...", m_phi, data) + # normalize coefficients + d[0, :] *= 1.0 / nlon + d[1:, :] *= 2.0 / nlon # Calculate cos and sin coefficients of theta component # Because the function is defined on (0,pi) # it can be expanded in just cosine terms. # this routine assumes that 0 and pi are not included - theta_cc = np.zeros((MMAX+1,MMAX+1)) - theta_sc = np.zeros((MMAX+1,MMAX+1)) - m_even = np.arange(0,MMAX+1,2) - m_odd = np.arange(1,MMAX,2) - n_even = len(m_even) - n_odd = len(m_odd) - - if np.isclose([theta[0],theta[nlat-1]],[0.0,np.pi]).all(): - # non-endpoints - nt = np.dot(mm,theta[1:nlat-1][np.newaxis,:]) - theta_cc[m_even,:] = 2.0*np.dot(dcos[m_even,1:nlat-1],np.cos(nt).T) - theta_sc[m_even,:] = 2.0*np.dot(dsin[m_even,1:nlat-1],np.cos(nt).T) - theta_cc[m_odd,:] = 2.0*np.dot(dcos[m_odd,1:nlat-1],np.sin(nt).T) - theta_sc[m_odd,:] = 2.0*np.dot(dsin[m_odd,1:nlat-1],np.sin(nt).T) + f = np.zeros((MMAX+1,MMAX+1), dtype=np.complex128) + m_even = slice(0, MMAX+1, 2) + m_odd = slice(1, MMAX, 2) + if np.isclose([th[0],th[nlat-1]], [0.0,np.pi]).all(): + # global case (includes poles) + # non-endpoints + k_th = np.exp(1j * np.einsum("h...,k...->kh...", th[1:nlat-1], mm)) + f[m_even,:] = 2.0*np.einsum("mh...,kh...->mk", d[m_even,1:nlat-1],k_th.real) + f[m_odd,:] = 2.0*np.einsum("mh...,kh...->mk", d[m_odd,1:nlat-1],k_th.imag) # endpoints - theta_cc[m_even,:] += np.dot((dcos[m_even,0]*np.cos(theta[0]) + - dcos[m_even,nlat-1]*np.cos(theta[nlat-1]))[:,np.newaxis], mm.T) - theta_sc[m_even,:] += np.dot((dsin[m_even,0]*np.cos(theta[0]) + - dsin[m_even,nlat-1]*np.cos(theta[nlat-1]))[:,np.newaxis], mm.T) - theta_cc[m_odd,:] += np.dot((dcos[m_odd,0]*np.sin(theta[0]) + - dcos[m_odd,nlat-1]*np.sin(theta[nlat-1]))[:,np.newaxis], mm.T) - theta_sc[m_odd,:] += np.dot((dsin[m_odd,0]*np.sin(theta[0]) + - dsin[m_odd,nlat-1]*np.sin(theta[nlat-1]))[:,np.newaxis], mm.T) - - elif not np.isclose([theta[0],theta[nlat-1]],[0.0,np.pi]).any(): - nt = np.dot(mm,theta[np.newaxis,:]) - theta_cc[m_even,:] = 2.0*np.dot(dcos[m_even,:],np.cos(nt).T) - theta_sc[m_even,:] = 2.0*np.dot(dsin[m_even,:],np.cos(nt).T) - theta_cc[m_odd,:] = 2.0*np.dot(dcos[m_odd,:],np.sin(nt).T) - theta_sc[m_odd,:] = 2.0*np.dot(dsin[m_odd,:],np.sin(nt).T) + k_th = np.exp(1j * mm* th[0]) + f[m_even,:] += np.einsum("m...,k...->mk", d[m_even,0], k_th) + f[m_odd,:] += np.einsum("m...,k...->mk", d[m_odd,0], k_th) + k_th = np.exp(1j * mm * th[nlat-1]) + f[m_even,:] += np.einsum("m...,k...->mk", d[m_even,nlat-1], k_th) + f[m_odd,:] += np.einsum("m...,k...->mk", d[m_odd,nlat-1], k_th) + elif not np.isclose([th[0],th[nlat-1]], [0.0,np.pi]).any(): + k_th = np.exp(1j * np.einsum("h...,k...->kh...", th, mm)) + f[m_even,:] = 2.0*np.einsum("mh...,kh...->mk", d[m_even,:],k_th.real) + f[m_odd,:] = 2.0*np.einsum("mh...,kh...->mk", d[m_odd,:],k_th.imag) else: raise ValueError('Latitude coordinates incompatible') # Normalize theta fourier coefficients - theta_cc[:,0] = theta_cc[:,0]/(2.0*nlat) - theta_cc[:,1:MMAX+1] = theta_cc[:,1:MMAX+1]/nlat - theta_sc[:,0] = theta_sc[:,0]/(2.0*nlat) - theta_sc[:,1:MMAX+1] = theta_sc[:,1:MMAX+1]/nlat - + f[:,0] *= 1.0/(2.0*nlat) + f[:,1:MMAX+1] *= 1.0/nlat # Correct normalization for the incomplete coverage of the sphere - delphi = np.abs(phi[1]-phi[0]) - deltheta = np.abs(theta[1]-theta[0]) - norm = nlon*delphi/(2.0*np.pi)*nlat*deltheta/np.pi - theta_cc = theta_cc*norm - theta_sc = theta_sc*norm + f[:] *= nlon*dphi/(2.0*np.pi) * nlat*dth/np.pi # Calculate cos and sin coefficients of Legendre functions # Expand m = even terms in a cosine series # Expand m = odd terms in a sine series # Both are stride 2 if (np.ndim(PLM) == 0): - plm = fourier_legendre(LMAX,MMAX) + Almk = fourier_legendre(LMAX, MMAX) else: - # use precomputed plms to improve computational speed - plm = PLM + # use precomputed alms to improve computational speed + Almk = PLM # Initializing output spherical harmonic matrices Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX) Ylms.clm = np.zeros((LMAX+1, MMAX+1)) Ylms.slm = np.zeros((LMAX+1, MMAX+1)) - # Sum theta fourier coefficients - # temp is the integral of cos(n theta) cos(k theta) dcos(theta) - # over the interval 0 to pi - # n and k must have like parities - - # m = even terms - k_even = np.zeros((n_even,n_even)) - for n in range(0,MMAX+2,2): - k_even[:,n//2] = 0.5*(1.0/(1.0-m_even-n) + 1.0/(1.0+m_even-n) + - 1.0/(1.0-m_even+n) + 1.0/(1.0+m_even+n)) - - k_odd = np.zeros((n_odd,n_odd)) - for n in range(1,MMAX+1,2): - k_odd[:,(n-1)//2] = 0.5*(1.0/(1-m_odd-n) + 1.0/(1+m_odd-n) + - 1.0/(1-m_odd+n) + 1.0/(1+m_odd+n)) - # calculate spherical harmonics for m == even terms - l_even = np.arange(0,LMAX+1,2) - l_odd = np.arange(1,LMAX,2) - for m in range(0,MMAX+2,2): - temp = np.dot(plm[l_even,m,m_even[:,np.newaxis]].T,k_even) - Ylms.clm[l_even,m] = np.dot(theta_cc[m,m_even[:,np.newaxis]].T,temp.T) - Ylms.slm[l_even,m] = np.dot(theta_sc[m,m_even[:,np.newaxis]].T,temp.T) - temp = np.dot(plm[l_odd,m,m_odd[:,np.newaxis]].T,k_odd) - Ylms.clm[l_odd,m] = np.dot(theta_cc[m,m_odd[:,np.newaxis]].T,temp.T) - Ylms.slm[l_odd,m] = np.dot(theta_sc[m,m_odd[:,np.newaxis]].T,temp.T) - - # m = odd terms - k_even = np.zeros((n_even,n_even)) - for n in range(0,MMAX+2,2): - k_even[:,n//2] = 0.5*(-1.0/(1-m_even-n) + 1.0/(1.0+m_even-n) + - 1.0/(1.0-m_even+n) - 1.0/(1.0+m_even+n)) - - k_odd = np.zeros((n_odd,n_odd)) - for n in range(1,MMAX+1,2): - k_odd[:,(n-1)//2] = 0.5*(-1.0/(1-m_odd-n) + 1.0/(1.0+m_odd-n) + - 1.0/(1.0-m_odd+n) - 1.0/(1.0+m_odd+n)) + # even l terms (l even, m even, k even) + l_even = slice(0, LMAX+1, 2) + n_even = np.arange(m_even.start, m_even.stop, m_even.step) + k_even = np.zeros((len(n_even), len(n_even))) + for k in range(0,MMAX+2,2): + k_even[:,k//2] = 0.5*(1.0/(1.0-n_even-k) + 1.0/(1.0+n_even-k) + + 1.0/(1.0-n_even+k) + 1.0/(1.0+n_even+k)) + # calculate summation over coefficients + Aeven = np.einsum("lmk...,kn...->lmn...", Almk[l_even,m_even,m_even], k_even) + Yeven = np.einsum("lmn...,mn...->lm...", Aeven, f[m_even,m_even]) + Ylms.clm[l_even,m_even] = Yeven.real + Ylms.slm[l_even,m_even] = Yeven.imag + + # odd l terms (l odd, m even, k odd) + l_odd = slice(1, LMAX, 2) + n_odd = np.arange(m_odd.start, m_odd.stop, m_odd.step) + k_odd = np.zeros((len(n_odd), len(n_odd))) + for k in range(1, MMAX+1, 2): + k_odd[:,(k-1)//2] = 0.5*(1.0/(1.0-n_odd-k) + 1.0/(1.0+n_odd-k) + + 1.0/(1.0-n_odd+k) + 1.0/(1.0+n_odd+k)) + # calculate summation over coefficients + Aodd = np.einsum("lmk...,kn...->lmn...", Almk[l_odd,m_even,m_odd], k_odd) + Yodd = np.einsum("lmn...,mn...->lm...", Aodd, f[m_even,m_odd]) + Ylms.clm[l_odd,m_even] = Yodd.real + Ylms.slm[l_odd,m_even] = Yodd.imag # calculate spherical harmonics for m == odd terms - l_even = np.arange(2,LMAX+1,2)# do not in include l=0 - l_odd = np.arange(1,LMAX,2) - for m in range(1,MMAX+1,2): - temp = np.dot(plm[l_even,m,m_even[:,np.newaxis]].T,k_even) - Ylms.clm[l_even,m] = np.dot(theta_cc[m,m_even[:,np.newaxis]].T,temp.T) - Ylms.slm[l_even,m] = np.dot(theta_sc[m,m_even[:,np.newaxis]].T,temp.T) - temp = np.dot(plm[l_odd,m,m_odd[:,np.newaxis]].T,k_odd) - Ylms.clm[l_odd,m] = np.dot(theta_cc[m,m_odd[:,np.newaxis]].T,temp.T) - Ylms.slm[l_odd,m] = np.dot(theta_sc[m,m_odd[:,np.newaxis]].T,temp.T) + # even l terms (l even, m odd, k even) + l_even = slice(2, LMAX+1, 2)# do not in include l=0 + n_even = np.arange(m_even.start, m_even.stop, m_even.step) + k_even = np.zeros((len(n_even), len(n_even))) + for k in range(0,MMAX+2,2): + k_even[:,k//2] = 0.5*(-1.0/(1.0-n_even-k) + 1.0/(1.0+n_even-k) + + 1.0/(1.0-n_even+k) - 1.0/(1.0+n_even+k)) + Aeven = np.einsum("lmk...,kn...->lmn...", Almk[l_even,m_odd,m_even], k_even) + Yeven = np.einsum("lmn...,mn...->lm...", Aeven, f[m_odd,m_even]) + Ylms.clm[l_even,m_odd] = Yeven.real + Ylms.slm[l_even,m_odd] = Yeven.imag + + # odd l terms (l odd, m odd, k odd) + l_odd = slice(1, LMAX, 2) + n_odd = np.arange(m_odd.start, m_odd.stop, m_odd.step) + k_odd = np.zeros((len(n_odd), len(n_odd))) + for k in range(1,MMAX+1,2): + k_odd[:,(k-1)//2] = 0.5*(-1.0/(1.0-n_odd-k) + 1.0/(1.0+n_odd-k) + + 1.0/(1.0-n_odd+k) - 1.0/(1.0+n_odd+k)) + # calculate summation over coefficients + Aodd = np.einsum("lmk...,kn...->lmn...", Almk[l_odd,m_odd,m_odd], k_odd) + Yodd = np.einsum("lmn...,mn...->lm...", Aodd, f[m_odd,m_odd]) + Ylms.clm[l_odd,m_odd] = Yodd.real + Ylms.slm[l_odd,m_odd] = Yodd.imag # Divide by Plm normalization Ylms.clm[:,0] /= 2.0 diff --git a/gravity_toolkit/gen_point_load.py b/gravity_toolkit/gen_point_load.py index 78ff3de3..563e2807 100644 --- a/gravity_toolkit/gen_point_load.py +++ b/gravity_toolkit/gen_point_load.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" gen_point_load.py -Written by Tyler Sutterley (04/2023) +Written by Tyler Sutterley (07/2026) Calculates gravitational spherical harmonic coefficients for point masses CALLING SEQUENCE: @@ -47,6 +47,8 @@ https://doi.org/10.1029/JB078i011p01760 UPDATE HISTORY: + 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 Updated 03/2023: improve typing for variables in docstrings Updated 02/2023: set custom units as top option in if/else statements @@ -107,8 +109,8 @@ def gen_point_load(data, lon, lat, LMAX=60, MMAX=None, UNITS=1, LOVE=None): # number of input data points npts = len(data.flatten()) # convert output longitude and latitude into radians - phi = np.pi*lon.flatten()/180.0 - theta = np.pi*(90.0 - lat.flatten())/180.0 + phi = np.radians(lon.flatten()) + theta = np.radians(90.0 - lat.flatten()) # extract degree dependent factor for specific units factors = gravity_toolkit.units(lmax=LMAX) @@ -137,7 +139,7 @@ def gen_point_load(data, lon, lat, LMAX=60, MMAX=None, UNITS=1, LOVE=None): # for each degree l for l in range(LMAX+1): m1 = np.min([l,MMAX]) + 1 - SPH = spherical_harmonic_matrix(l, D, phi, theta, dfactor[l]) + SPH = _complex_harmonics(l, D, phi, theta, dfactor[l]) # truncate to spherical harmonic order and save to output Ylms.clm[l,:m1] = SPH.real[:m1] Ylms.slm[l,:m1] = SPH.imag[:m1] @@ -145,7 +147,7 @@ def gen_point_load(data, lon, lat, LMAX=60, MMAX=None, UNITS=1, LOVE=None): return Ylms # calculate spherical harmonics of degree l evaluated at (theta,phi) -def spherical_harmonic_matrix(l, data, phi, theta, coeff): +def _complex_harmonics(l, data, phi, theta, coeff): """ Calculates the spherical harmonics for a particular degree evaluated from data at coordinates @@ -168,15 +170,15 @@ def spherical_harmonic_matrix(l, data, phi, theta, coeff): Ylms: np.ndarray spherical harmonic coefficients in Eulerian form """ - # calculate normalized legendre polynomials (points, order) - Pl = legendre(l, np.cos(theta), NORMALIZE=True).T + # calculate normalized legendre polynomials (order, points) + Pl = legendre(l, np.cos(theta), NORMALIZE=True) # spherical harmonic orders up to degree l m = np.arange(0, l+1) - # calculate Euler's of spherical harmonic order multiplied by azimuth phi - mphi = np.exp(1j*np.dot(np.squeeze(phi)[:,np.newaxis], m[np.newaxis,:])) - # reshape data to order - D = np.kron(np.ones((1, l+1)), data[:,np.newaxis]) - # calculate spherical harmonics and multiply by coefficients and data - Ylms = coeff*D*Pl*mphi - # calculate the sum over all points and return harmonics for degree l - return np.sum(Ylms, axis=0) + # calculate Euler's of order m multiplied by azimuth phi + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", m, phi)) + # reshape data to (order, points) + D = np.kron(np.ones((l+1, 1)), data[np.newaxis, :]) + # calculate spherical harmonics summing over all points + Yl = np.einsum("mp...,mp...,mp...->m...", D, Pl, m_phi) + # return harmonics for degree l multiplied by coefficients + return coeff*Yl diff --git a/gravity_toolkit/gen_spherical_cap.py b/gravity_toolkit/gen_spherical_cap.py index bea477b9..cc95fd78 100755 --- a/gravity_toolkit/gen_spherical_cap.py +++ b/gravity_toolkit/gen_spherical_cap.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" gen_spherical_cap.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (07/2026) Calculates gravitational spherical harmonic coefficients for a spherical cap Creating a spherical cap with generating angle alpha is a 2 step process: @@ -61,6 +61,8 @@ https://doi.org/10.1007/s00190-011-0522-7 UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 06/2023: modified custom units case to not convert to cmwe Updated 03/2023: simplified unit degree factors using units class improve typing for variables in docstrings @@ -147,8 +149,8 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None, MMAX = np.copy(LMAX) # convert lon and lat to radians - phi = lon*np.pi/180.0# Longitude in radians - th = (90.0 - lat)*np.pi/180.0# Colatitude in radians + phi = np.radians(lon)# Longitude in radians + th = np.radians(90.0 - lat)# Colatitude in radians # Earth Parameters factors = gravity_toolkit.units(lmax=LMAX) @@ -160,7 +162,7 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None, if (RAD_CAP != 0): # if given spherical cap radius in degrees # converting to radians - alpha = RAD_CAP*np.pi/180.0 + alpha = np.radians(RAD_CAP) elif (AREA != 0): # if given spherical cap area in cm^2 # radius in centimeters @@ -231,12 +233,7 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None, # this would be the plm for the center of the spherical cap # used to rotate the spherical cap to point lat/lon if PLM is None: - plmout,_ = plm_holmes(LMAX, np.cos(th)) - # truncate precomputed plms to order - plmout = np.squeeze(plmout[:,:MMAX+1,:]) - else: - # truncate precomputed plms to degree and order - plmout = PLM[:LMAX+1,:MMAX+1] + PLM, _ = plm_holmes(LMAX, np.cos(th)) # calculate array of m values ranging from 0 to MMAX (harmonic orders) # MMAX+1 as there are MMAX+1 elements between 0 and MMAX @@ -247,29 +244,20 @@ def gen_spherical_cap(data, lon, lat, LMAX=60, MMAX=None, # data normally is 1 for a uniform 1cm water equivalent layer # but can be a mass point if reconstructing a spherical harmonic field # NOTE: NOT a matrix multiplication as data (and phi) is a single point - dcos = unit_conv*data*np.cos(m*phi) - dsin = unit_conv*data*np.sin(m*phi) + d = unit_conv*data*np.exp(1j*m*phi) # Multiplying by plm_alpha (F_l from Jacob 2012) plm = np.zeros((LMAX+1, MMAX+1)) - # Initializing preliminary spherical harmonic matrices - yclm = np.zeros((LMAX+1, MMAX+1)) - yslm = np.zeros((LMAX+1, MMAX+1)) # Initializing output spherical harmonic matrices Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX) - Ylms.clm = np.zeros((LMAX+1, MMAX+1)) - Ylms.slm = np.zeros((LMAX+1, MMAX+1)) - for m in range(0,MMAX+1):# MMAX+1 to include MMAX - l = np.arange(m,LMAX+1)# LMAX+1 to include LMAX - # rotate spherical cap to be centered at lat/lon - plm[l,m] = plmout[l,m]*pl_alpha[l] - # multiplying clm by cos(m*phi) and slm by sin(m*phi) - # to get a field of spherical harmonics - yclm[l,m] = plm[l,m]*dcos[m] - yslm[l,m] = plm[l,m]*dsin[m] - # multiplying by factors to convert to geoid coefficients - Ylms.clm[l,m] = dfactor[l]*yclm[l,m] - Ylms.slm[l,m] = dfactor[l]*yslm[l,m] + # rotate spherical cap to be centered at lat/lon + plm = np.einsum("lm...,l...->lm...", PLM[:LMAX+1,:MMAX+1], pl_alpha) + # multiplying clm by cos(m*phi) and slm by sin(m*phi) + # to get a field of spherical harmonics + ylm = np.einsum("lm...,m...->lm...", plm, d) + # Multiplying by factors to convert to fully normalized coefficients + Ylms.clm = np.einsum("l...,lm...->lm...", dfactor, ylm.real) + Ylms.slm = np.einsum("l...,lm...->lm...", dfactor, ylm.imag) # return the output spherical harmonics object return Ylms diff --git a/gravity_toolkit/gen_stokes.py b/gravity_toolkit/gen_stokes.py index 6b43f3fa..271c7e64 100755 --- a/gravity_toolkit/gen_stokes.py +++ b/gravity_toolkit/gen_stokes.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" gen_stokes.py -Written by Tyler Sutterley (04/2023) +Written by Tyler Sutterley (07/2026) Converts data from the spatial domain to spherical harmonic coefficients @@ -43,6 +43,8 @@ and filters the GRACE/GRACE-FO coefficients for striping errors UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 06/2025: copy latitude and longitude as float64 for numpy 2.0 stability Updated 04/2023: allow love numbers to be None for custom units case Updated 03/2023: improve typing for variables in docstrings @@ -129,25 +131,15 @@ def gen_stokes(data, lon, lat, LMIN=0, LMAX=60, MMAX=None, UNITS=1, # grid dimensions nlat = np.int64(len(lat)) - # grid step - dlon = np.abs(lon[1]-lon[0]) - dlat = np.abs(lat[1]-lat[0]) - # longitude degree spacing in radians - dphi = dlon*np.pi/180.0 - # colatitude degree spacing in radians - dth = dlat*np.pi/180.0 - - # convert latitude and longitude to float if integers - lon = lon.astype(np.float64) - lat = lat.astype(np.float64) - # reformatting longitudes to range 0:360 (if previously -180:180) - lon = np.squeeze(lon.copy()) - if np.any(lon < 0): - lon[lon < 0] += 360.0 # Longitude in radians - phi = lon[np.newaxis,:]*np.pi/180.0 - # Colatitude in radians - th = (90.0 - np.squeeze(lat.copy()))*np.pi/180.0 + phi = np.radians(np.squeeze(lon.copy())) + # reformatting longitudes to range 0:360 (if previously -180:180) + phi = np.where(phi < 0, phi + 2.0*np.pi, phi) + # colatitude in radians + th = np.radians(90.0 - np.squeeze(lat.copy())) + # grid step in radians + dphi = np.abs(phi[1] - phi[0]) + dth = np.abs(th[1] - th[0]) # reforming data to lonXlat if input latXlon sz = np.shape(data) @@ -181,9 +173,8 @@ def gen_stokes(data, lon, lat, LMIN=0, LMAX=60, MMAX=None, UNITS=1, # Calculating cos/sin of phi arrays # output [m,phi] - m = np.arange(MMAX+1) - ccos = np.cos(np.dot(m[:,np.newaxis],phi)) - ssin = np.sin(np.dot(m[:,np.newaxis],phi)) + mm = np.arange(MMAX+1) + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", mm, phi)) # Calculating fully-normalized Legendre Polynomials # Output is plm[l,m,th] @@ -193,33 +184,20 @@ def gen_stokes(data, lon, lat, LMIN=0, LMAX=60, MMAX=None, UNITS=1, # if plms are not pre-computed: calculate Legendre polynomials PLM, dPLM = plm_holmes(LMAX, np.cos(th)) - # Multiplying by integration factors [sin(theta)*dtheta*dphi] - # truncate legendre polynomials to spherical harmonic order MMAX - for j in range(0,nlat): - plm[:,m,j] = PLM[:,m,j]*int_fact[j] + # truncate legendre polynomials to degree and order + plm = np.einsum("lmh...,h...->lmh...", PLM[:LMAX+1,:MMAX+1,:], int_fact) - # Initializing preliminary spherical harmonic matrices - yclm = np.zeros((LMAX+1, MMAX+1)) - yslm = np.zeros((LMAX+1, MMAX+1)) # Initializing output spherical harmonic matrices Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX) - Ylms.clm = np.zeros((LMAX+1, MMAX+1)) - Ylms.slm = np.zeros((LMAX+1, MMAX+1)) # Multiplying gridded data with sin/cos of m#phis # This will sum through all phis in the dot product # output [m,theta] - dcos = np.dot(ccos,data) - dsin = np.dot(ssin,data) - for l in range(LMIN,LMAX+1):# equivalent to LMIN:LMAX - mm = np.min([MMAX,l])# truncate to MMAX if specified (if l > MMAX) - m = np.arange(0,mm+1)# mm+1 elements between 0 and mm - # Summing product of plms and data over all latitudes - # axis=1 signifies the direction of the summation - yclm[l,m] = np.sum(plm[l,m,:]*dcos[m,:], axis=1) - yslm[l,m] = np.sum(plm[l,m,:]*dsin[m,:], axis=1) - # Multiplying by factors to convert to fully normalized coefficients - Ylms.clm[l,m] = dfactor[l]*yclm[l,m] - Ylms.slm[l,m] = dfactor[l]*yslm[l,m] + d = np.einsum("mp...,ph...->mh...", m_phi, data) + # Summing product of plms and data over all latitudes + ylm = np.einsum("lmh...,mh...->lm...", plm, d) + # Multiplying by factors to convert to fully normalized coefficients + Ylms.clm = np.einsum("l...,lm...->lm...", dfactor, ylm.real) + Ylms.slm = np.einsum("l...,lm...->lm...", dfactor, ylm.imag) # return the output spherical harmonics object return Ylms \ No newline at end of file diff --git a/gravity_toolkit/geocenter.py b/gravity_toolkit/geocenter.py index 247b77ff..b419f7ca 100644 --- a/gravity_toolkit/geocenter.py +++ b/gravity_toolkit/geocenter.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" geocenter.py -Written by Tyler Sutterley (06/2024) +Written by Tyler Sutterley (07/2026) Data class for reading and processing geocenter data PYTHON DEPENDENCIES: @@ -15,6 +15,7 @@ https://github.com/yaml/pyyaml UPDATE HISTORY: + Updated 07/2026: add dunder (magic) methods for mathematical operations Updated 06/2024: use wrapper to importlib for optional dependencies Updated 05/2024: make subscriptable and allow item assignment Updated 09/2023: add group option to netCDF read function @@ -1218,6 +1219,71 @@ def __str__(self): properties.append(f" end_month: {max(self.month)}") return '\n'.join(properties) + def __add__(self, other): + """Add values to a ``geocenter`` object""" + temp = self.copy() + return temp.add(other) + + def __div__(self, other): + """Divide values from a ``geocenter`` object""" + return self.__truediv__(other) + + def __iadd__(self, other): + """In-place add values to a ``geocenter`` object""" + return self.add(other) + + def __idiv__(self, other): + """In-place divide values from a ``geocenter`` object""" + return self.__itruediv__(other) + + def __imul__(self, other): + """In-place multiply values from a ``geocenter`` object""" + if isinstance(other, (int, float, np.ndarray)): + return self.scale(other) + else: + return self.multiply(other) + + def __ipow__(self, other): + """In-place raise values from a ``geocenter`` object to a power""" + return self.power(other) + + def __isub__(self, other): + """In-place subtract values from a ``geocenter`` object""" + return self.subtract(other) + + def __itruediv__(self, other): + """In-place divide values from a ``geocenter`` object""" + if isinstance(other, (int, float, np.ndarray)): + return self.scale(1.0 / other) + else: + return self.divide(other) + + def __mul__(self, other): + """Multiply values from a ``geocenter`` object""" + temp = self.copy() + if isinstance(other, (int, float, np.ndarray)): + return temp.scale(other) + else: + return temp.multiply(other) + + def __pow__(self, other): + """Raise values from a ``geocenter`` object to a power""" + temp = self.copy() + return temp.power(other) + + def __sub__(self, other): + """Subtract values from a ``geocenter`` object""" + temp = self.copy() + return temp.subtract(other) + + def __truediv__(self, other): + """Divide values from a ``geocenter`` object""" + temp = self.copy() + if isinstance(other, (int, float, np.ndarray)): + return temp.scale(1.0 / other) + else: + return temp.divide(other) + def __len__(self): """Number of months """ diff --git a/gravity_toolkit/grace_input_months.py b/gravity_toolkit/grace_input_months.py index 12cf35b7..3bc40197 100644 --- a/gravity_toolkit/grace_input_months.py +++ b/gravity_toolkit/grace_input_months.py @@ -682,7 +682,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, for i,grace_month in enumerate(months): count = np.count_nonzero(C20_input['month'] == grace_month) if (count != 0): - k, = np.nonzero(C20_input['month'] == grace_month) + k, = np.flatnonzero(C20_input['month'] == grace_month) grace_Ylms['clm'][2,0,i] = np.copy(C20_input['data'][k]) grace_Ylms['eclm'][2,0,i] = np.copy(C20_input['error'][k]) @@ -697,7 +697,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, for i,grace_month in enumerate(months): count = np.count_nonzero(C21_input['month'] == grace_month) if (count != 0) and (grace_month > 176): - k, = np.nonzero(C21_input['month'] == grace_month) + k, = np.flatnonzero(C21_input['month'] == grace_month) grace_Ylms['clm'][2,1,i] = np.copy(C21_input['C2m'][k]) grace_Ylms['slm'][2,1,i] = np.copy(C21_input['S2m'][k]) grace_Ylms['eclm'][2,1,i] = np.copy(C21_input['eC2m'][k]) @@ -714,7 +714,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, for i,grace_month in enumerate(months): count = np.count_nonzero(C22_input['month'] == grace_month) if (count != 0) and (grace_month > 176): - k, = np.nonzero(C22_input['month'] == grace_month) + k, = np.flatnonzero(C22_input['month'] == grace_month) grace_Ylms['clm'][2,2,i] = np.copy(C22_input['C2m'][k]) grace_Ylms['slm'][2,2,i] = np.copy(C22_input['S2m'][k]) grace_Ylms['eclm'][2,2,i] = np.copy(C22_input['eC2m'][k]) @@ -731,7 +731,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, for i,grace_month in enumerate(months): count = np.count_nonzero(C30_input['month'] == grace_month) if (count != 0) and (grace_month > 176): - k, = np.nonzero(C30_input['month'] == grace_month) + k, = np.flatnonzero(C30_input['month'] == grace_month) grace_Ylms['clm'][3,0,i] = np.copy(C30_input['data'][k]) grace_Ylms['eclm'][3,0,i] = np.copy(C30_input['error'][k]) @@ -746,7 +746,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, for i,grace_month in enumerate(months): count = np.count_nonzero(C40_input['month'] == grace_month) if (count != 0) and (grace_month > 176): - k, = np.nonzero(C40_input['month'] == grace_month) + k, = np.flatnonzero(C40_input['month'] == grace_month) grace_Ylms['clm'][4,0,i] = np.copy(C40_input['data'][k]) grace_Ylms['eclm'][4,0,i] = np.copy(C40_input['error'][k]) @@ -761,7 +761,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, for i,grace_month in enumerate(months): count = np.count_nonzero(C50_input['month'] == grace_month) if (count != 0) and (grace_month > 176): - k, = np.nonzero(C50_input['month'] == grace_month) + k, = np.flatnonzero(C50_input['month'] == grace_month) grace_Ylms['clm'][5,0,i] = np.copy(C50_input['data'][k]) grace_Ylms['eclm'][5,0,i] = np.copy(C50_input['error'][k]) @@ -790,7 +790,7 @@ def grace_input_months(base_dir, PROC, DREL, DSET, LMAX, start_mon, end_mon, raise IOError(f'No Matching Geocenter Months ({gm})') # for each considered date for i,grace_month in enumerate(months): - k, = np.nonzero(DEG1_input.month == grace_month) + k, = np.flatnonzero(DEG1_input.month == grace_month) count = np.count_nonzero(DEG1_input.month == grace_month) # Degree 1 is missing for particular month if (count == 0) and kwargs['MODEL_DEG1']: diff --git a/gravity_toolkit/harmonic_gradients.py b/gravity_toolkit/harmonic_gradients.py index 8e0704b2..f306c6eb 100644 --- a/gravity_toolkit/harmonic_gradients.py +++ b/gravity_toolkit/harmonic_gradients.py @@ -2,7 +2,7 @@ u""" harmonic_gradients.py Original IDL code calc_grad.pro written by Sean Swenson -Adapted by Tyler Sutterley (03/2023) +Adapted by Tyler Sutterley (07/2026) Calculates the zonal and meridional gradients of a scalar field from a series of spherical harmonics @@ -29,6 +29,8 @@ Legendre functions UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 03/2023: improve typing for variables in docstrings added geostrophic currents program from Wahr et al. (2002) Updated 10/2022: cleaned up program for public release @@ -48,7 +50,7 @@ def harmonic_gradients(clm1, slm1, lon, lat, LMIN=0, LMAX=60, MMAX=None): """ Calculates the gradient of a scalar field from a series of - spherical harmonics + spherical harmonics :cite:p:`Driscoll:1994bp` Parameters ---------- @@ -81,87 +83,50 @@ def harmonic_gradients(clm1, slm1, lon, lat, MMAX = np.copy(LMAX) # Longitude in radians - phi = (np.squeeze(lon)*np.pi/180.0)[np.newaxis,:] + phi = np.radians(np.squeeze(lon)) # Colatitude in radians - th = (90.0 - np.squeeze(lat))*np.pi/180.0 - thmax = len(np.squeeze(lat)) - phimax = len(np.squeeze(lon)) + th = np.radians(90.0 - np.squeeze(lat)) + thmax = len(th) + # spherical harmonic degree and order + ll = np.arange(0,LMAX+1)# lmax+1 to include lmax + mm = np.arange(0,MMAX+1)# mmax+1 to include mmax + # real (cosine) and imaginary (sine) components + Ylm = np.zeros((LMAX+1, MMAX+1), dtype=np.complex128) # Truncating harmonics to degree and order LMAX # removing coefficients below LMIN and above MMAX - mm = np.arange(0,MMAX+1) - clm = np.zeros((LMAX+1, MMAX+1)) - slm = np.zeros((LMAX+1, MMAX+1)) - clm[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm] - slm[LMIN:LMAX+1,mm] = slm1[LMIN:LMAX+1,mm] - # spherical harmonic degree and order - ll = np.arange(0,LMAX+1)[np.newaxis, :]# lmax+1 to include lmax - mm = np.arange(0,MMAX+1)[:, np.newaxis]# mmax+1 to include mmax + Ylm.real[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm].copy() + Ylm.imag[LMIN:LMAX+1,mm] = -slm1[LMIN:LMAX+1,mm].copy() + dlm = np.einsum("l...,lm...->lm", np.sqrt((ll+1.0)*ll), -1j*Ylm) # generate Vlm coefficients (vlm and wlm) - vlm, wlm = legendre_gradient(LMAX, MMAX) - - dlm = np.zeros((LMAX+1,LMAX+1,2)) - # minus sign is because lat and theta change with opposite sign - for l in range(0,LMAX+1): - dlm[l,:,0] = -clm[l,:]*np.sqrt((l+1.0)*l) - dlm[l,:,1] = -slm[l,:]*np.sqrt((l+1.0)*l) - + Vlmk, Wlmk = legendre_gradient(LMAX, MMAX) + # even and odd spherical harmonic orders m_even = np.arange(0,MMAX+2,2) m_odd = np.arange(1,MMAX,2) + # Euler's formula for theta * k and m * phi + k_th = np.exp(1j * np.einsum("h...,k...->kh...", th, ll)) + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", mm, phi)) # Calculate fourier coefficients from legendre coefficients - d_cos = np.zeros((LMAX+1,thmax,2)) - d_sin = np.zeros((LMAX+1,thmax,2)) - cnk = np.cos(np.dot(th[:,np.newaxis],ll)) - snk = np.sin(np.dot(th[:,np.newaxis],ll)) - - wtmp = np.zeros((len(m_even),LMAX+1,2)) - vtmp = np.zeros((len(m_even),LMAX+1,2)) - # m = even terms (vlm,wlm sine series) - for n in range(0,LMAX+1): - wtmp[:,n,0] = np.sum(wlm[:,m_even,n]*dlm[:,m_even,0],axis=0) - wtmp[:,n,1] = np.sum(wlm[:,m_even,n]*dlm[:,m_even,1],axis=0) - vtmp[:,n,0] = np.sum(vlm[:,m_even,n]*dlm[:,m_even,0],axis=0) - vtmp[:,n,1] = np.sum(vlm[:,m_even,n]*dlm[:,m_even,1],axis=0) - - d_cos[m_even,:,0] = np.dot(wtmp[:,:,1],np.transpose(snk)) - d_sin[m_even,:,0] = np.dot(-wtmp[:,:,0],np.transpose(snk)) - d_cos[m_even,:,1] = np.dot(vtmp[:,:,1],np.transpose(snk)) - d_sin[m_even,:,1] = np.dot(-vtmp[:,:,0],np.transpose(snk)) - - # m = odd terms (vlm,wlm cosine series) - wtmp = np.zeros((len(m_odd),LMAX+1,2)) - vtmp = np.zeros((len(m_odd),LMAX+1,2)) - for n in range(0,LMAX+1): - wtmp[:,n,0] = np.sum(wlm[:,m_odd,n]*dlm[:,m_odd,0],axis=0) - wtmp[:,n,1] = np.sum(wlm[:,m_odd,n]*dlm[:,m_odd,1],axis=0) - vtmp[:,n,0] = np.sum(vlm[:,m_odd,n]*dlm[:,m_odd,0],axis=0) - vtmp[:,n,1] = np.sum(vlm[:,m_odd,n]*dlm[:,m_odd,1],axis=0) - - d_cos[m_odd,:,0] = np.dot(wtmp[:,:,1],np.transpose(cnk)) - d_sin[m_odd,:,0] = np.dot(-wtmp[:,:,0],np.transpose(cnk)) - d_cos[m_odd,:,1] = np.dot(vtmp[:,:,1],np.transpose(cnk)) - d_sin[m_odd,:,1] = np.dot(-vtmp[:,:,0],np.transpose(cnk)) - - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(mm,phi)) - ssin = np.sin(np.dot(mm,phi)) - # Final signal recovery from fourier coefficients - gradients = np.zeros((phimax,thmax,2)) - gradients[:,:,0] = np.dot(np.transpose(ccos), d_cos[:,:,0]) + \ - np.dot(np.transpose(ssin), d_sin[:,:,0]) - gradients[:,:,1] = np.dot(np.transpose(ccos), d_cos[:,:,1]) + \ - np.dot(np.transpose(ssin), d_sin[:,:,1]) - # return the gradient fields - return gradients + d = np.zeros((LMAX+1,thmax,2), dtype=np.complex128) + wtmp = np.einsum("lmk...,lm...->mk", Wlmk, dlm) + vtmp = np.einsum("lmk...,lm...->mk", Vlmk, dlm) + d[m_even,:,0] = np.einsum("mk...,kh...->mh", wtmp[m_even,:], k_th.imag) + d[m_even,:,1] = np.einsum("mk...,kh...->mh", vtmp[m_even,:], k_th.imag) + d[m_odd,:,0] = np.einsum("mk...,kh...->mh", wtmp[m_odd,:], k_th.real) + d[m_odd,:,1] = np.einsum("mk...,kh...->mh", vtmp[m_odd,:], k_th.real) + # calculate the zonal and meridional gradients of the scalar field + gradients = np.einsum("mp...,mhd...->phd...", m_phi, d) + # return the gradient fields and drop imaginary component + return gradients.real def geostrophic_currents(clm1, slm1, lon, lat, LMIN=0, LMAX=60, MMAX=None, RAD=0, DENSITY=1.035, LOVE=None, PLM=None): r""" - Converts data from spherical harmonic coefficients to a - spatial fields of ocean geostrophic currents following + Converts data from spherical harmonic coefficients to spatial + fields of approximate ocean geostrophic currents following :cite:p:`Wahr:1998hy,Wahr:2002ie` Parameters @@ -205,10 +170,10 @@ def geostrophic_currents(clm1, slm1, lon, lat, MMAX = np.copy(LMAX) # Longitude in radians - phi = (np.squeeze(lon)*np.pi/180.0)[np.newaxis,:] + phi = np.radians(np.squeeze(lon)) phmax = len(lon) # colatitude in radians - th = (90.0 - np.squeeze(lat))*np.pi/180.0 + th = np.radians(90.0 - np.squeeze(lat)) thmax = len(th) # Gaussian Smoothing @@ -230,7 +195,8 @@ def geostrophic_currents(clm1, slm1, lon, lat, # smooth harmonics and convert to output units clm = np.zeros((LMAX+1, MMAX+1, 2)) slm = np.zeros((LMAX+1, MMAX+1, 2)) - # zonal flow harmonics + # zonal flow harmonics (equation 3) + # differentiating Legendre polynomials with respect to longitude for l in range(1, LMAX): # truncate to degree and order mm = np.arange(0, np.min([l,MMAX])+1) @@ -240,7 +206,8 @@ def geostrophic_currents(clm1, slm1, lon, lat, np.sqrt(((l+1)**2 - mm**2)*(2.0*l + 3.0)/(2.0*l + 1)) clm[l,mm,0] = coeff*wl[l]*(temp1*clm1[l-1,mm] - temp2*clm1[l+1,mm]) slm[l,mm,0] = coeff*wl[l]*(temp1*slm1[l-1,mm] - temp2*slm1[l+1,mm]) - # meridional flow harmonics + # meridional flow harmonics (equation 4) + # differentiating Legendre polynomials with respect to colatitude for l in range(0, LMAX+1): # truncate to degree and order mm = np.arange(0, np.min([l,MMAX])+1) @@ -251,32 +218,20 @@ def geostrophic_currents(clm1, slm1, lon, lat, # Truncating harmonics to degree and order LMAX # removing coefficients below LMIN and above MMAX mm = np.arange(0, MMAX+1) - clm[LMIN:LMAX+1,mm,:] = clm[LMIN:LMAX+1,mm,:] - slm[LMIN:LMAX+1,mm,:] = slm[LMIN:LMAX+1,mm,:] - - # Calculate fourier coefficients from legendre coefficients - d_cos = np.zeros((MMAX+1,thmax,2))# [m,th] - d_sin = np.zeros((MMAX+1,thmax,2))# [m,th] - for k in range(0,thmax): - # summation over all spherical harmonic degrees - temp = 1.0/(np.cos(th[k])*np.sin(th[k])) - # for each direction of flow - for d in range(2): - d_cos[:,k,d] = temp*np.sum(PLM[:,mm,k]*clm[:,mm,d],axis=0) - d_sin[:,k,d] = temp*np.sum(PLM[:,mm,k]*slm[:,mm,d],axis=0) - - # Final signal recovery from fourier coefficients - m = np.arange(0,MMAX+1)[:,np.newaxis] - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # real (cosine) and imaginary (sine) components + Ylm = clm[LMIN:LMAX+1,:MMAX+1,:] - 1j * slm[LMIN:LMAX+1,:MMAX+1,:] + # convolve legendre polynomials and truncate to degree and order + iint = 1.0/(np.cos(th)*np.sin(th)) + plm = np.einsum("h...,lmh...->lmh...", iint, PLM[LMIN:LMAX+1,:MMAX+1,:]) + # summation over all spherical harmonic degrees + pconv = np.einsum("lmh...,lmd...->mhd...", plm, Ylm) + + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", mm, phi)) # output geostrophic current fields - currents = np.zeros((phmax,thmax,2)) - # for each direction of flow - for d in range(2): - # summation of cosine and sine harmonics - currents[:,:,d] = np.dot(np.transpose(ccos),d_cos[:,:,d]) + \ - np.dot(np.transpose(ssin),d_sin[:,:,d]) + currents = np.empty((phmax,thmax,2)) + # summation of cosine and sine harmonics + currents[:] = np.einsum("mp...,mhd...->phd...", m_phi, pconv) - # return the current fields - return currents + # return the current fields and drop imaginary component + return currents.real diff --git a/gravity_toolkit/harmonic_summation.py b/gravity_toolkit/harmonic_summation.py index 3eed02ce..efb071e1 100755 --- a/gravity_toolkit/harmonic_summation.py +++ b/gravity_toolkit/harmonic_summation.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" harmonic_summation.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (07/2026) Returns the spatial field for a series of spherical harmonics @@ -30,6 +30,8 @@ units.py: class for converting spherical harmonic data to specific units UPDATE HISTORY: + 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 Updated 03/2023: allow units inputs to be strings for named types improve typing for variables in docstrings @@ -85,41 +87,32 @@ def harmonic_summation(clm1, slm1, lon, lat, if MMAX is None: MMAX = np.copy(LMAX) - # Longitude in radians - phi = (np.squeeze(lon)*np.pi/180.0)[np.newaxis,:] + # longitude in radians + phi = np.radians(np.squeeze(lon)) # colatitude in radians - th = (90.0 - np.squeeze(lat))*np.pi/180.0 - thmax = len(th) + th = np.radians(90.0 - np.squeeze(lat)) # if plms are not pre-computed: calculate Legendre polynomials if PLM is None: PLM, dPLM = plm_holmes(LMAX, np.cos(th)) + # spherical harmonic order + mm = np.arange(0,MMAX+1)# mmax+1 to include mmax + # real (cosine) and imaginary (sine) components + Ylm = np.zeros((LMAX+1, MMAX+1), dtype=np.complex128) # Truncating harmonics to degree and order LMAX # removing coefficients below LMIN and above MMAX - mm = np.arange(0, MMAX+1) - clm = np.zeros((LMAX+1, MMAX+1)) - slm = np.zeros((LMAX+1, MMAX+1)) - clm[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm] - slm[LMIN:LMAX+1,mm] = slm1[LMIN:LMAX+1,mm] + Ylm.real[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm] + Ylm.imag[LMIN:LMAX+1,mm] = -slm1[LMIN:LMAX+1,mm] # Calculate fourier coefficients from legendre coefficients - d_cos = np.zeros((MMAX+1,thmax))# [m,th] - d_sin = np.zeros((MMAX+1,thmax))# [m,th] - for k in range(0,thmax): - # summation over all spherical harmonic degrees - d_cos[:,k] = np.sum(PLM[:,mm,k]*clm[:,mm],axis=0) - d_sin[:,k] = np.sum(PLM[:,mm,k]*slm[:,mm],axis=0) - - # Final signal recovery from fourier coefficients - m = np.arange(0,MMAX+1)[:,np.newaxis] - # Calculating cos(m*phi) and sin(m*phi) - ccos = np.cos(np.dot(m,phi)) - ssin = np.sin(np.dot(m,phi)) + # summation over all spherical harmonic degrees + pconv = np.einsum("lmh...,lm...->mh...", PLM[:LMAX+1,:MMAX+1,:], Ylm) + # calculating cos(m*phi) and sin(m*phi) using Euler's formula + m_phi = np.exp(1j * np.einsum("m...,p...->mp...", mm, phi)) # summation of cosine and sine harmonics - s = np.dot(np.transpose(ccos),d_cos) + np.dot(np.transpose(ssin),d_sin) - - # return output data - return s + spatial = np.einsum("mp...,mh...->ph...", m_phi, pconv) + # return output data and drop imaginary component + return spatial.real def harmonic_transform(clm1, slm1, lon, lat, LMIN=0, LMAX=60, MMAX=None, PLM=None): @@ -164,25 +157,23 @@ def harmonic_transform(clm1, slm1, lon, lat, # number of longitudinal points phimax = len(np.squeeze(lon)) # colatitude in radians - th = (90.0 - np.squeeze(lat))*np.pi/180.0 + th = np.radians(90.0 - np.squeeze(lat)) thmax = len(th) # if plms are not pre-computed: calculate Legendre polynomials if PLM is None: PLM, dPLM = plm_holmes(LMAX, np.cos(th)) - # combined Ylms and Fourier coefficients (complex) - Ylms = np.zeros((LMAX+1, MMAX+1),dtype=np.complex128) - delta_M = np.zeros((MMAX+1,thmax),dtype=np.complex128)# [m,th] - # Real (cosine) and imaginary (sine) components + # spherical harmonic order + mm = np.arange(0,MMAX+1)# mmax+1 to include mmax + # real (cosine) and imaginary (sine) components + Ylm = np.zeros((LMAX+1, MMAX+1), dtype=np.complex128) # Truncating harmonics to degree and order LMAX # removing coefficients below LMIN and above MMAX - Ylms[LMIN:LMAX+1,:MMAX+1] = clm1[LMIN:LMAX+1,0:MMAX+1] - \ - slm1[LMIN:LMAX+1,0:MMAX+1]*1j + Ylm.real[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm] + Ylm.imag[LMIN:LMAX+1,mm] = -slm1[LMIN:LMAX+1,mm] # calculate Ylms summation for each theta band - for k in range(0,thmax): - # summation over all spherical harmonic degrees - delta_M[:,k] = np.sum(PLM[:,:,k]*Ylms[:,:],axis=0)/2.0 + delta_M = np.einsum("lmh...,lm...->mh...", PLM, Ylm / 2.0) # output spatial field from FFT transformation s = np.zeros((phimax,thmax)) @@ -269,14 +260,11 @@ def stokes_summation(clm1, slm1, lon, lat, else: raise ValueError(f'Unknown units {UNITS}') - # truncate to degree and order - mm = np.arange(0, MMAX+1) + # spherical harmonic order + mm = np.arange(0,MMAX+1)# mmax+1 to include mmax # smooth harmonics and convert to output units - clm = np.zeros((LMAX+1, MMAX+1)) - slm = np.zeros((LMAX+1, MMAX+1)) - for l in range(0, LMAX+1):# LMAX+1 to include LMAX - clm[l,:] = wl[l]*dfactor[l]*clm1[l,mm] - slm[l,:] = wl[l]*dfactor[l]*slm1[l,mm] + clm = np.einsum("l,l,lm->lm", wl, dfactor, clm1[:, mm]) + slm = np.einsum("l,l,lm->lm", wl, dfactor, slm1[:, mm]) # return the spatial field return harmonic_summation(clm, slm, lon, lat, diff --git a/gravity_toolkit/harmonics.py b/gravity_toolkit/harmonics.py index 877c0edf..f90a1ad0 100644 --- a/gravity_toolkit/harmonics.py +++ b/gravity_toolkit/harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" harmonics.py -Written by Tyler Sutterley (06/2024) +Written by Tyler Sutterley (07/2026) Contributions by Hugo Lecomte Spherical harmonic data class for processing GRACE/GRACE-FO Level-2 data @@ -25,6 +25,7 @@ destripe_harmonics.py: filters spherical harmonics for correlated errors UPDATE HISTORY: + Updated 07/2026: add dunder (magic) methods for mathematical operations Updated 06/2024: use wrapper to importlib for optional dependencies Updated 05/2024: make subscriptable and allow item assignment Updated 10/2023: place time and month variables in try/except block @@ -1893,6 +1894,71 @@ def __str__(self): properties.append(f" end_month: {max(self.month)}") return '\n'.join(properties) + def __add__(self, other): + """Add values to a ``harmonics`` object""" + temp = self.copy() + return temp.add(other) + + def __div__(self, other): + """Divide values from a ``harmonics`` object""" + return self.__truediv__(other) + + def __iadd__(self, other): + """In-place add values to a ``harmonics`` object""" + return self.add(other) + + def __idiv__(self, other): + """In-place divide values from a ``harmonics`` object""" + return self.__itruediv__(other) + + def __imul__(self, other): + """In-place multiply values from a ``harmonics`` object""" + if isinstance(other, (int, float, np.ndarray)): + return self.scale(other) + else: + return self.multiply(other) + + def __ipow__(self, other): + """In-place raise values from a ``harmonics`` object to a power""" + return self.power(other) + + def __isub__(self, other): + """In-place subtract values from a ``harmonics`` object""" + return self.subtract(other) + + def __itruediv__(self, other): + """In-place divide values from a ``harmonics`` object""" + if isinstance(other, (int, float, np.ndarray)): + return self.scale(1.0 / other) + else: + return self.divide(other) + + def __mul__(self, other): + """Multiply values from a ``harmonics`` object""" + temp = self.copy() + if isinstance(other, (int, float, np.ndarray)): + return temp.scale(other) + else: + return temp.multiply(other) + + def __pow__(self, other): + """Raise values from a ``harmonics`` object to a power""" + temp = self.copy() + return temp.power(other) + + def __sub__(self, other): + """Subtract values from a ``harmonics`` object""" + temp = self.copy() + return temp.subtract(other) + + def __truediv__(self, other): + """Divide values from a ``harmonics`` object""" + temp = self.copy() + if isinstance(other, (int, float, np.ndarray)): + return temp.scale(1.0 / other) + else: + return temp.divide(other) + def __len__(self): """Number of months """ diff --git a/gravity_toolkit/mascons.py b/gravity_toolkit/mascons.py index 2a8f08c1..98a01c74 100644 --- a/gravity_toolkit/mascons.py +++ b/gravity_toolkit/mascons.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" mascons.py -Written by Tyler Sutterley (03/2023) +Written by Tyler Sutterley (07/2026) Conversion routines for publicly available GRACE/GRACE-FO mascon solutions PYTHON DEPENDENCIES: @@ -12,6 +12,7 @@ mascon2grid.m written by Felix Landerer and David Wiese (JPL) UPDATE HISTORY: + Updated 07/2026: use np.radians to convert from degrees to radians Updated 03/2023: improve typing for variables in docstrings Updated 11/2022: use lowercase keyword arguments Updated 04/2022: updated docstrings to numpy documentation format @@ -63,15 +64,13 @@ def to_gsfc(gdata, lon, lat, lon_center, lat_center, lon_span, lat_span): # number of mascons nmas = len(lon_center) # convert mascon centers to -180:180 - gt180, = np.nonzero(lon_center > 180) - lon_center[gt180] -= 360.0 + lon_center = np.where(lon_center > 180, lon_center - 360.0, lon_center) # remove singleton dimensions lat = np.squeeze(lat) lon = np.squeeze(lon) # for mascons centered on 180: use 0:360 alon = np.copy(lon) - lt0, = np.nonzero(lon < 0) - alon[lt0] += 360.0 + alon = np.where(alon < 0, alon + 360.0, alon) # loop over each mascon bin and average gdata with the cos-lat weights # for that bin @@ -98,12 +97,12 @@ def to_gsfc(gdata, lon, lat, lon_center, lat_center, lon_span, lat_span): else: ilon = lon.copy() # indices for grid points within the mascon - I, = np.nonzero((lat >= lat_bound[0]) & (lat < lat_bound[1])) - J, = np.nonzero((ilon >= lon_bound[0]) & (ilon < lon_bound[1])) + I, = np.flatnonzero((lat >= lat_bound[0]) & (lat < lat_bound[1])) + J, = np.flatnonzero((ilon >= lon_bound[0]) & (ilon < lon_bound[1])) I,J = (I[np.newaxis,:], J[:,np.newaxis]) # calculate average data for mascon bin - mascon_array['data'][k] = np.mean((np.cos(lat[I]*np.pi/180.0) / - np.mean(np.cos(lat[I]*np.pi/180.0)))*gdata[I,J]/len(I)) + mascon_array['data'][k] = np.mean((np.cos(np.radians(lat[I])) / + np.mean(np.cos(np.radians(lat[I]))))*gdata[I,J]/len(I)) mascon_array['lat_center'][k] = lat_center[k] mascon_array['lon_center'][k] = lon_center[k] @@ -154,13 +153,13 @@ def to_jpl(gdata, lon, lat, lon_bound, lat_bound): mascon_array['lat'] = np.zeros((nmas)) for k in range(0,nmas): # indices for grid points within the mascon - I, = np.nonzero((lat >= lat_bound[k,1]) & (lat < lat_bound[k,0])) - J, = np.nonzero((lon >= lon_bound[k,0]) & (lon < lon_bound[k,2])) + I, = np.flatnonzero((lat >= lat_bound[k,1]) & (lat < lat_bound[k,0])) + J, = np.flatnonzero((lon >= lon_bound[k,0]) & (lon < lon_bound[k,2])) nlt = np.count_nonzero((lat >= lat_bound[k,1]) & (lat < lat_bound[k,0])) I,J = (I[np.newaxis,:], J[:,np.newaxis]) # calculate average data for mascon bin - mascon_array['data'][k] = np.mean((np.cos(lat[I]*np.pi/180.0) / - np.mean(np.cos(lat[I]*np.pi/180.0)))*gdata[I,J]/nlt) + mascon_array['data'][k] = np.mean((np.cos(np.radians(lat[I])) / + np.mean(np.cos(np.radians(lat[I]))))*gdata[I,J]/nlt) # calculate coordinates of mascon center mascon_array['lat'][k] = (lat_bound[k,1]+lat_bound[k,0])/2.0 mascon_array['lon'][k] = (lon_bound[k,1]+lon_bound[k,2])/2.0 @@ -219,8 +218,7 @@ def from_gsfc(mscdata, grid_spacing, lon_center, lat_center, lon_span, lat_span, # number of mascons nmas = len(lon_center) # convert mascon centers to -180:180 - gt180, = np.nonzero(lon_center > 180) - lon_center[gt180] -= 360.0 + lon_center = np.where(lon_center > 180, lon_center - 360.0, lon_center) # Define output latitude and longitude grids lon = np.arange(-180.0+grid_spacing/2.0,180.0+grid_spacing/2.0,grid_spacing) @@ -228,8 +226,7 @@ def from_gsfc(mscdata, grid_spacing, lon_center, lat_center, lon_span, lat_span, nlon, nlat = (len(lon),len(lat)) # for mascons centered on 180: use 0:360 alon = np.copy(lon) - lt0, = np.nonzero(lon < 0) - alon[lt0] += 360.0 + alon = np.where(alon < 0, alon + 360.0, alon) # loop over each mascon bin and assign value to grid points inside bin: mdata = np.zeros((nlat, nlon)) @@ -252,8 +249,8 @@ def from_gsfc(mscdata, grid_spacing, lon_center, lat_center, lon_span, lat_span, else: ilon = lon.copy() # indices for grid points within the mascon - I, = np.nonzero((lat >= lat_bound[0]) & (lat < lat_bound[1])) - J, = np.nonzero((ilon >= lon_bound[0]) & (ilon < lon_bound[1])) + I, = np.flatnonzero((lat >= lat_bound[0]) & (lat < lat_bound[1])) + J, = np.flatnonzero((ilon >= lon_bound[0]) & (ilon < lon_bound[1])) I,J = (I[np.newaxis,:], J[:,np.newaxis]) mdata[I,J] = mscdata[k] @@ -310,8 +307,8 @@ def from_jpl(mscdata, grid_spacing, lon_bound, lat_bound, **kwargs): # loop over each mascon bin and assign value to grid points inside bin: mdata = np.zeros((nlat, nlon)) for k in range(0, nmas): - I, = np.nonzero((lat >= lat_bound[k,1]) & (lat < lat_bound[k,0])) - J, = np.nonzero((lon >= lon_bound[k,0]) & (lon < lon_bound[k,2])) + I, = np.flatnonzero((lat >= lat_bound[k,1]) & (lat < lat_bound[k,0])) + J, = np.flatnonzero((lon >= lon_bound[k,0]) & (lon < lon_bound[k,2])) I,J = (I[np.newaxis,:], J[:,np.newaxis]) mdata[I,J] = mscdata[k] diff --git a/gravity_toolkit/sea_level_equation.py b/gravity_toolkit/sea_level_equation.py index 24b53231..1e1e3906 100644 --- a/gravity_toolkit/sea_level_equation.py +++ b/gravity_toolkit/sea_level_equation.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -sea_level_equation.py (06/2025) +sea_level_equation.py (07/2026) Solves the sea level equation with the option of including polar motion feedback Uses a Clenshaw summation to calculate the spherical harmonic summation @@ -37,7 +37,6 @@ ITERATIONS: maximum number of iterations for the solver PLM: input Legendre polynomials FILL_VALUE: value used over land points - ASTYPE: floating point precision for calculating Clenshaw summation SCALE: scaling factor to prevent underflow in Clenshaw summation PYTHON DEPENDENCIES: @@ -91,6 +90,8 @@ https://doi.org/10.1029/JB090iB11p09363 UPDATE HISTORY: + Updated 07/2026: use np.einsum for spherical harmonic summations + use np.radians to convert from degrees to radians Updated 06/2025: added option to set the density of sea water (g/cm^3) Updated 03/2023: improve typing for variables in docstrings Updated 01/2023: refactored associated legendre polynomials @@ -131,8 +132,7 @@ # PURPOSE: Computes Sea Level Fingerprints including polar motion feedback def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, LOVE=None, BODY_TIDE_LOVE=0, FLUID_LOVE=0, DENSITY=1.0, POLAR=True, - ITERATIONS=6, PLM=None, FILL_VALUE=0, ASTYPE=np.longdouble, SCALE=1e-280, - **kwargs): + ITERATIONS=6, PLM=None, FILL_VALUE=0, SCALE=1e-280, **kwargs): r""" Solves the sea level equation with the option of including polar motion feedback :cite:p:`Farrell:1976hm,Kendall:2005ds,Mitrovica:2003cq` @@ -180,8 +180,6 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, Legendre polynomials FILL_VALUE: float, default 0 Invalid value used over land points - ASTYPE: np.dtype, default np.longdouble - Floating point precision for calculating Clenshaw summation SCALE: float, default 1e-280 Scaling factor to prevent underflow in Clenshaw summation @@ -192,17 +190,17 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, """ # dimensions of land function - nphi,nth = np.shape(land_function) + nphi, nth = np.shape(land_function) # calculate colatitude and longitude in radians - th = (90.0 - glat)*np.pi/180.0 - phi = np.squeeze(glon*np.pi/180.0) + th = np.radians(90.0 - glat) + phi = np.radians(np.squeeze(glon)) # calculate ocean function from land function ocean_function = 1.0 - land_function # indices of the ocean function ii,jj = np.nonzero(ocean_function) # extract arrays of kl, hl, and ll Love Numbers - hl,kl,ll = LOVE + hl, kl, ll = LOVE # density of water [g/cm^3] rho_water = np.float64(DENSITY) # Earth Parameters @@ -271,21 +269,25 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, # calculate coefh and coefp for each degree and order # see equation 11 from Tamisiea et al (2010) - coefh = np.zeros((LMAX+1,LMAX+1)) - coefp = np.zeros((LMAX+1,LMAX+1)) + coefh = np.zeros((LMAX+1, LMAX+1)) + coefp = np.zeros((LMAX+1, LMAX+1)) for l in range(LMAX+1): + m = np.arange(0, l+1) + # tilt factor for degree l + gamma_l = (1.0 + kl[l] - hl[l]) # coefh and coefp will be the same for all orders except for degree 2 # and order 1 (if POLAR motion feedback is included) - m = np.arange(0,l+1) - coefh[l,m] = 3.0*rho_water*(1.0 + kl[l] - hl[l])/rho_e/np.float64(2*l+1) - coefp[l,m] = (1.0 + kl[l] - hl[l])/(kl[l] + 1.0) + coefh[l,m] = 3.0*rho_water*gamma_l/rho_e/np.float64(2*l+1) + coefp[l,m] = gamma_l/(kl[l] + 1.0) # if degree 2 and POLAR parameter is set if (l == 2) and POLAR: + # tilt factor for body tides + gamma_2b = (1.0 + k2b - h2b) # calculate coefficient for polar motion feedback and add to coefs # For small perturbations in rotation vector: driving potential # will be dominated by degree two and order one polar wander # effects (quadrantal geometry effects) (Kendall et al., 2005) - coefpmf = (1.0 + k2b - h2b)*(1.0 + kl[l])/(klf - k2b) + coefpmf = gamma_2b*(1.0 + kl[l])/(klf - k2b) # add effects of polar motion feedback to order 1 coefficients coefh[l,1] += 3.0*rho_water*coefpmf/rho_e/np.float64(2*l+1) coefp[l,1] += coefpmf/(kl[l] + 1.0) @@ -295,15 +297,15 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, # calculate Legendre polynomials using Holmes and Featherstone relation PLM, dPLM = plm_holmes(LMAX, np.cos(th)) # calculate sin of colatitudes - gth,gphi = np.meshgrid(th, phi) + gth, gphi = np.meshgrid(th, phi) u = np.sin(gth[ii,jj]) # indices of spherical harmonics for calculating eps - l1,m1 = np.tril_indices(LMAX+1) + l1, m1 = np.tril_indices(LMAX+1) # total mass of the surface mass load [g] from harmonics tmass = 4.0*np.pi*(rad_e**3.0)*rho_e*loadClm[0,0]/3.0 # convert ocean function into a series of spherical harmonics - ocean_Ylms = gen_harmonics(ocean_function,glon,glat,LMAX=LMAX,PLM=PLM) + ocean_Ylms = gen_harmonics(ocean_function, glon, glat, LMAX=LMAX, PLM=PLM) # total area of ocean calculated by integrating the ocean function ocean_area = 4.0*np.pi*ocean_Ylms.clm[0,0] @@ -316,8 +318,12 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, logging.info(f'Total Ocean Area: {ocean_area:0.10g}') logging.info(f'Uniform Ocean Height: {sea_height:0.10g}') + # allocate for output sea level field + sea_level = np.empty((nphi, nth)) + # complex load spherical harmonics + loadYlms = loadClm - 1j*loadSlm # distribute sea height over ocean harmonics - height_Ylms = ocean_Ylms.scale(sea_height) + height_Ylms = ocean_Ylms * sea_height # iterate solutions until convergence or reaching total iterations n_iter = 1 # use maximum eps values from Mitrovica and Peltier (1991) @@ -325,42 +331,38 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, eps = np.inf eps_max = 1e-4 while (eps > eps_max) and (n_iter <= ITERATIONS): - # allocate for sea level field of iteration - sea_level = np.zeros((nphi,nth)) + # zero out the sea level field for this iteration + sea_level[:,:] = 0.0 # calculate combined spherical harmonics for Clenshaw summation - clm1 = coefh*height_Ylms.clm + rad_e*coefp*loadClm - slm1 = coefh*height_Ylms.slm + rad_e*coefp*loadSlm + Ylm1 = coefh*height_Ylms.ilm + rad_e*coefp*loadYlms # calculate clenshaw summations over colatitudes - s_m_c = np.zeros((nth,LMAX*2+2)) + cs_m = np.zeros((nth, LMAX+1), dtype=np.clongdouble) for m in range(LMAX, -1, -1): - s_m_c[:,2*m:2*m+2] = clenshaw_s_m(np.cos(th), m, clm1, slm1, LMAX, - ASTYPE=ASTYPE, SCALE=SCALE) + cs_m[:,m] = _clenshaw(np.cos(th), m, Ylm1, LMAX, SCALE=SCALE) # calculate cos(phi) cos_phi_2 = 2.0*np.cos(phi) # matrix of cos/sin m*phi summation - cos_m_phi = np.zeros((nphi,LMAX+2),dtype=ASTYPE) - sin_m_phi = np.zeros((nphi,LMAX+2),dtype=ASTYPE) + m_phi = np.zeros((nphi, LMAX+2), dtype=np.clongdouble) # initialize matrix with values at lmax+1 and lmax - cos_m_phi[:,LMAX+1] = np.cos(ASTYPE(LMAX + 1)*phi) - sin_m_phi[:,LMAX+1] = np.sin(ASTYPE(LMAX + 1)*phi) - cos_m_phi[:,LMAX] = np.cos(ASTYPE(LMAX)*phi) - sin_m_phi[:,LMAX] = np.sin(ASTYPE(LMAX)*phi) + m_phi[:,LMAX+1] = np.exp(1j * (LMAX + 1) * phi) + m_phi[:,LMAX] = np.exp(1j * LMAX*phi) # calculate summation - gc=np.multiply(s_m_c[np.newaxis,:,2*LMAX],cos_m_phi[:,np.newaxis,LMAX]) - gs=np.multiply(s_m_c[np.newaxis,:,2*LMAX+1],sin_m_phi[:,np.newaxis,LMAX]) - s_m = gc[ii,jj] + gs[ii,jj] + g = np.einsum("h...,p...->ph...", cs_m[:,LMAX], m_phi[:,LMAX]) + # discard imaginary component + s_m = g[ii,jj].real # iterate to calculate complete summation for m in range(LMAX-1, 0, -1): - cos_m_phi[:,m] = cos_phi_2*cos_m_phi[:,m+1] - cos_m_phi[:,m+2] - sin_m_phi[:,m] = cos_phi_2*sin_m_phi[:,m+1] - sin_m_phi[:,m+2] - a_m = np.sqrt((2.0*m+3.0)/(2.0*m+2.0)) - gc=np.multiply(s_m_c[np.newaxis,:,2*m],cos_m_phi[:,np.newaxis,m]) - gs=np.multiply(s_m_c[np.newaxis,:,2*m+1],sin_m_phi[:,np.newaxis,m]) - s_m = a_m*u*s_m + gc[ii,jj] + gs[ii,jj] + # calculate summation for order m + m_phi[:,m] = cos_phi_2*m_phi[:,m+1] - m_phi[:,m+2] + a_m = np.sqrt((2.0*m + 3.0)/(2.0*m + 2.0)) + g = np.einsum("h...,p...->ph...", cs_m[:,m], m_phi[:,m]) + # update summation and discard imaginary component + s_m = a_m*u*s_m + g[ii,jj].real + # add the l=0/m=0 term + gs_m = np.kron(np.ones((nphi, 1)), cs_m.real[:, 0]) # calculate new sea level for iteration - gsmc,gcmp = np.meshgrid(s_m_c[:,0],cos_m_phi[:,0]) - sea_level[ii,jj] = np.sqrt(3.0)*u*s_m + gsmc[ii,jj] + sea_level[ii,jj] = np.sqrt(3.0)*u*s_m + gs_m[ii,jj] # calculate spherical harmonic field for iteration Ylms = gen_harmonics(sea_level, glon, glat, LMAX=LMAX, PLM=PLM) @@ -383,11 +385,10 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, # Equation 48 of Mitrovica and Peltier (1991) # add difference to total sea level field to force mass conservation sea_level += sea_height*ocean_function[:,:] - uniform_Ylms = ocean_Ylms.scale(sea_height) - Ylms.add(uniform_Ylms) + Ylms += ocean_Ylms * sea_height # calculate eps to determine if solution is appropriately converged - mod1 = np.sqrt(height_Ylms.clm**2 + height_Ylms.slm**2) - mod2 = np.sqrt(Ylms.clm**2 + Ylms.slm**2) + mod1 = np.hypot(height_Ylms.clm, height_Ylms.slm) + mod2 = np.hypot(Ylms.clm, Ylms.slm) eps = np.abs(np.sum(mod2[l1,m1] - mod1[l1,m1])/np.sum(mod1[l1,m1])) # save height harmonics for use in the next iteration height_Ylms = Ylms.copy() @@ -397,8 +398,8 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, # calculate final total mass for sanity check omass = 4.0*np.pi*(rad_e**2.0)*rho_water*height_Ylms.clm[0,0] # if verbose output: sanity check of masses - logging.info('Original Total Ocean Mass: {0:0.10g}'.format(-tmass/1e15)) - logging.info('Final Iterated Ocean Mass: {0:0.10g}'.format(omass/1e15)) + logging.info(f'Original Total Ocean Mass: {-tmass/1e15:0.10g}') + logging.info(f'Final Iterated Ocean Mass: {omass/1e15:0.10g}') # set final invalid points to fill value if applicable if (FILL_VALUE != 0): @@ -410,9 +411,7 @@ def sea_level_equation(loadClm, loadSlm, glon, glat, land_function, LMAX=0, # PURPOSE: compute Clenshaw summation of the fully normalized associated # Legendre's function for constant order m -def clenshaw_s_m(t, m, clm1, slm1, lmax, - ASTYPE=np.longdouble, SCALE=1e-280 - ): +def _clenshaw(t, m, Ylm1, lmax, SCALE=1e-280): """ Compute conditioned arrays for Clenshaw summation from the fully-normalized associated Legendre's function for an order m @@ -423,14 +422,10 @@ def clenshaw_s_m(t, m, clm1, slm1, lmax, elements ranging from -1 to 1, typically cos(th) m: int spherical harmonic order - clm1: np.ndarray - cosine spherical harmonics - slm1: np.ndarray - sine spherical harmonics + Ylm1: np.ndarray + complex form of spherical harmonics lmax: int maximum spherical harmonic degree - ASTYPE: np.dtype, default np.longdouble - floating point precision for calculating Clenshaw summation SCALE: float, default 1e-280 scaling factor to prevent underflow in Clenshaw summation @@ -441,49 +436,40 @@ def clenshaw_s_m(t, m, clm1, slm1, lmax, """ # allocate for output matrix N = len(t) - s_m = np.zeros((N,2),dtype=ASTYPE) + s_m = np.zeros((N), dtype=np.clongdouble) # scaling to prevent overflow - clm = SCALE*clm1.astype(ASTYPE) - slm = SCALE*slm1.astype(ASTYPE) + ylm = SCALE*Ylm1.astype(np.clongdouble) # convert lmax and m to float - lm = ASTYPE(lmax) - mm = ASTYPE(m) + lm = np.float64(lmax) + mm = np.float64(m) if (m == lmax): - s_m[:,0] = np.copy(clm[lmax,lmax]) - s_m[:,1] = np.copy(slm[lmax,lmax]) + s_m[:] = np.copy(ylm[lmax,lmax]) elif (m == (lmax-1)): a_lm = np.sqrt(((2.0*lm-1.0)*(2.0*lm+1.0))/((lm-mm)*(lm+mm)))*t - s_m[:,0] = a_lm*clm[lmax,lmax-1] + clm[lmax-1,lmax-1] - s_m[:,1] = a_lm*slm[lmax,lmax-1] + slm[lmax-1,lmax-1] + s_m[:] = a_lm*ylm[lmax,lmax-1] + ylm[lmax-1,lmax-1] elif ((m <= (lmax-2)) and (m >= 1)): - s_mm_c_pre_2 = np.copy(clm[lmax,m]) - s_mm_s_pre_2 = np.copy(slm[lmax,m]) + s_mm_minus_2 = np.copy(ylm[lmax,m]) a_lm = np.sqrt(((2.0*lm-1.0)*(2.0*lm+1.0))/((lm-mm)*(lm+mm)))*t - s_mm_c_pre_1 = a_lm*s_mm_c_pre_2 + clm[lmax-1,m] - s_mm_s_pre_1 = a_lm*s_mm_s_pre_2 + slm[lmax-1,m] + s_mm_minus_1 = a_lm*s_mm_minus_2 + ylm[lmax-1,m] for l in range(lmax-2, m-1, -1): - ll = ASTYPE(l) + ll = np.float64(l) a_lm=np.sqrt(((2.0*ll+1.0)*(2.0*ll+3.0))/((ll+1.0-mm)*(ll+1.0+mm)))*t b_lm=np.sqrt(((2.*ll+5.)*(ll+mm+1.)*(ll-mm+1.))/((ll+2.-mm)*(ll+2.+mm)*(2.*ll+1.))) - s_mm_c = a_lm * s_mm_c_pre_1 - b_lm * s_mm_c_pre_2 + clm[l,m] - s_mm_s = a_lm * s_mm_s_pre_1 - b_lm * s_mm_s_pre_2 + slm[l,m] - s_mm_c_pre_2 = np.copy(s_mm_c_pre_1) - s_mm_s_pre_2 = np.copy(s_mm_s_pre_1) - s_mm_c_pre_1 = np.copy(s_mm_c) - s_mm_s_pre_1 = np.copy(s_mm_s) - s_m[:,0] = np.copy(s_mm_c) - s_m[:,1] = np.copy(s_mm_s) + s_mm_l = a_lm * s_mm_minus_1 - b_lm * s_mm_minus_2 + ylm[l,m] + s_mm_minus_2 = np.copy(s_mm_minus_1) + s_mm_minus_1 = np.copy(s_mm_l) + s_m[:] = np.copy(s_mm_l) elif (m == 0): - s_mm_c_pre_2 = np.copy(clm[lmax,0]) + s_mm_minus_2 = np.copy(ylm[lmax,0]) a_lm = np.sqrt(((2.0*lm-1.0)*(2.0*lm+1.0))/(lm*lm))*t - s_mm_c_pre_1 = a_lm * s_mm_c_pre_2 + clm[lmax-1,0] + s_mm_minus_1 = a_lm * s_mm_minus_2 + ylm[lmax-1,0] for l in range(lmax-2, m-1, -1): - ll = ASTYPE(l) + ll = np.float64(l) a_lm=np.sqrt(((2.0*ll+1.0)*(2.0*ll+3.0))/((ll+1.0)*(ll+1.0)))*t b_lm=np.sqrt(((2.0*ll+5.0)*(ll+1.0)*(ll+1.0))/((ll+2.0)*(ll+2.0)*(2.0*ll+1.0))) - s_mm_c = a_lm * s_mm_c_pre_1 - b_lm * s_mm_c_pre_2 + clm[l,0] - s_mm_c_pre_2 = np.copy(s_mm_c_pre_1) - s_mm_c_pre_1 = np.copy(s_mm_c) - s_m[:,0] = np.copy(s_mm_c) + s_mm_l = a_lm * s_mm_minus_1 - b_lm * s_mm_minus_2 + ylm[l,0] + s_mm_minus_2 = np.copy(s_mm_minus_1) + s_mm_minus_1 = np.copy(s_mm_l) + s_m[:] = np.copy(s_mm_l) # return rescaled s_m return s_m/SCALE diff --git a/gravity_toolkit/spatial.py b/gravity_toolkit/spatial.py index adc80ea5..ded4d4bc 100644 --- a/gravity_toolkit/spatial.py +++ b/gravity_toolkit/spatial.py @@ -20,6 +20,8 @@ time.py: utilities for calculating time operations UPDATE HISTORY: + Updated 07/2026: add dunder (magic) methods for mathematical operations + add option to change the output format for ascii files Updated 10/2024: allow 2D and 3D arrays in output netCDF4 files Updated 06/2024: use wrapper to importlib for optional dependencies Updated 05/2024: make subscriptable and allow item assignment @@ -734,6 +736,8 @@ def to_ascii(self, filename, **kwargs): full path of output ascii file date: bool, default True ``spatial`` objects contain date information + float_format: str, default '12.4f' + format for floating point numbers in output ascii file verbose: bool, default False Output file and variable information """ @@ -744,14 +748,16 @@ def to_ascii(self, filename, **kwargs): logging.info(str(self.filename)) # open the output file fid = self.filename.open(mode='w', encoding='utf8') + file_format = '{0:10.4f} {1:10.4f} ' + float_format = kwargs.get('float_format', '12.4f') if hasattr(self, 'error') and kwargs['date']: - file_format = '{0:10.4f} {1:10.4f} {2:12.4f} {3:12.4f} {4:10.4f}' + file_format += ' '.join([f'{{{i}:' + float_format + '}' for i in (2,3,4)]) elif hasattr(self, 'error'): - file_format = '{0:10.4f} {1:10.4f} {2:12.4f} {3:12.4f}' + file_format += ' '.join([f'{{{i}:' + float_format + '}' for i in (2,3)]) elif kwargs['date']: - file_format = '{0:10.4f} {1:10.4f} {2:12.4f} {4:10.4f}' + file_format += ' '.join([f'{{{i}:' + float_format + '}' for i in (2,4)]) else: - file_format = '{0:10.4f} {1:10.4f} {2:12.4f}' + file_format += ' '.join([f'{{{i}:' + float_format + '}' for i in (2,)]) # write to file for each valid latitude and longitude ii,jj = np.nonzero((self.data != self.fill_value) & (~self.mask)) for i,j in zip(ii,jj): @@ -1707,6 +1713,59 @@ def __str__(self): properties.append(f" end_month: {max(self.month)}") return '\n'.join(properties) + def __add__(self, other): + """Add values to a ``spatial`` object""" + temp = self.copy() + return temp.offset(other) + + def __div__(self, other): + """Divide values from a ``spatial`` object""" + return self.__truediv__(other) + + def __iadd__(self, other): + """In-place add values to a ``spatial`` object""" + return self.offset(other) + + def __idiv__(self, other): + """In-place divide values from a ``spatial`` object""" + return self.__itruediv__(other) + + def __imul__(self, other): + """In-place multiply values from a ``spatial`` object""" + return self.scale(other) + + def __ipow__(self, other): + """In-place raise values from a ``spatial`` object to a power""" + return self.power(other) + + def __isub__(self, other): + """In-place subtract values from a ``spatial`` object""" + return self.offset(-other) + + def __itruediv__(self, other): + """In-place divide values from a ``spatial`` object""" + return self.scale(1.0 / other) + + def __mul__(self, other): + """Multiply values from a ``spatial`` object""" + temp = self.copy() + return temp.scale(other) + + def __pow__(self, other): + """Raise values from a ``spatial`` object to a power""" + temp = self.copy() + return temp.power(other) + + def __sub__(self, other): + """Subtract values from a ``spatial`` object""" + temp = self.copy() + return temp.offset(-other) + + def __truediv__(self, other): + """Divide values from a ``spatial`` object""" + temp = self.copy() + return temp.scale(1.0 / other) + def __len__(self): """Number of months """ diff --git a/gravity_toolkit/time_series/amplitude.py b/gravity_toolkit/time_series/amplitude.py index 7d174cb2..ad4df441 100755 --- a/gravity_toolkit/time_series/amplitude.py +++ b/gravity_toolkit/time_series/amplitude.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" amplitude.py -Written by Tyler Sutterley (01/2023) +Written by Tyler Sutterley (07/2026) Calculate the amplitude and phase of a harmonic function from calculated sine and cosine of a series of measurements @@ -21,6 +21,7 @@ numpy: Scientific Computing Tools For Python (https://numpy.org) UPDATE HISTORY: + Updated 07/2026: use np.radians to convert from degrees to radians Updated 01/2023: refactored time series analysis functions Updated 04/2022: updated docstrings to numpy documentation format Updated 07/2020: added function docstrings @@ -49,5 +50,5 @@ def amplitude(bsin, bcos): phase from the harmonic functions in degrees """ ampl = np.sqrt(bsin**2.0 + bcos**2.0) - ph = 180.0*np.arctan2(bcos, bsin)/np.pi - return (ampl,ph) + ph = np.degrees(np.arctan2(bcos, bsin)) + return (ampl, ph) diff --git a/gravity_toolkit/time_series/piecewise.py b/gravity_toolkit/time_series/piecewise.py index e2349007..66042163 100755 --- a/gravity_toolkit/time_series/piecewise.py +++ b/gravity_toolkit/time_series/piecewise.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" piecewise.py -Written by Tyler Sutterley (04/2023) +Written by Tyler Sutterley (07/2026) Fits a synthetic signal to data over a time period by ordinary or weighted least-squares for breakpoint analysis @@ -61,6 +61,7 @@ scipy: Scientific Tools for Python (https://docs.scipy.org/doc/) UPDATE HISTORY: + Updated 07/2026: use np.hypot to calculate the sum of two squares Updated 04/2023: option to include extra fit terms in the design matrix Updated 01/2023: refactored time series analysis functions Updated 04/2022: updated docstrings to numpy documentation format @@ -293,7 +294,7 @@ def piecewise(t_in, d_in, BREAK_TIME=None, BREAKPOINT=None, # Recalculating beta2 error beta_err = np.copy(temp_err) - beta_err[2] = np.sqrt(temp_err[1]**2 + temp_err[2]**2) + beta_err[2] = np.hypot(temp_err[1], temp_err[2]) # Weighted sum of squares Error WSSE = np.dot(np.transpose(wi*(d_in[0:nmax] - np.dot(DMAT,beta_mat))), wi*(d_in[0:nmax] - np.dot(DMAT,beta_mat)))/np.float64(nu) @@ -314,7 +315,7 @@ def piecewise(t_in, d_in, BREAK_TIME=None, BREAKPOINT=None, temp_err[i] = np.sum((NORMEQ[i,:]*P_err)**2) # Recalculating beta2 error beta_err = np.copy(temp_err) - beta_err[2] = np.sqrt(temp_err[1]**2 + temp_err[2]**2) + beta_err[2] = np.hypot(temp_err[1], temp_err[2]) # Mean square error MSE = np.dot(np.transpose(d_in[0:nmax] - np.dot(DMAT,beta_mat)), (d_in[0:nmax] - np.dot(DMAT,beta_mat)))/np.float64(nu) @@ -359,10 +360,10 @@ def piecewise(t_in, d_in, BREAK_TIME=None, BREAKPOINT=None, # Recalculating standard error for beta2 st_err = np.copy(temp_std) - st_err[2] = np.sqrt(temp_std[1]**2 + temp_std[2]**2) + st_err[2] = np.hypot(temp_std[1], temp_std[2]) # Recalculating beta2 error beta_err = np.copy(temp_err) - beta_err[2] = np.sqrt(temp_err[1]**2 + temp_err[2]**2) + beta_err[2] = np.hypot(temp_err[1], temp_err[2]) return {'beta':beta_out, 'error':beta_err, 'std_err':st_err, 'R2':rsquare, 'R2Adj':rsq_adj, 'MSE':MSE, 'NRMSE':NRMSE, 'AIC':AIC, 'BIC':BIC, diff --git a/gravity_toolkit/time_series/smooth.py b/gravity_toolkit/time_series/smooth.py index 247e934c..d80050e6 100755 --- a/gravity_toolkit/time_series/smooth.py +++ b/gravity_toolkit/time_series/smooth.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" smooth.py -Written by Tyler Sutterley (01/2023) +Written by Tyler Sutterley (07/2026) Computes a moving average of a time-series using three possible routines: 1) centered moving average @@ -48,6 +48,7 @@ scipy: Scientific Tools for Python (https://docs.scipy.org/doc/) UPDATE HISTORY: + Updated 07/2026: use np.hypot to calculate the sum of two squares Updated 01/2023: refactored time series analysis functions Updated 04/2022: updated docstrings to numpy documentation format Updated 05/2021: define int/float precision to prevent deprecation warning @@ -224,13 +225,13 @@ def smooth(t_in, d_in, HFWTH=6, MOVING=False, DATA_ERR=0, WEIGHT=0, # annual component AS,AC = beta_mat[SEAS:SEAS+2] dannual[ran] += wi*np.dot(TMAT[:,SEAS:SEAS+2],[AS,AC]) - annamp[ran] += wi*np.sqrt(AS**2 + AC**2) - annphase[ran] += wi*np.arctan2(AC,AS)*180.0/np.pi + annamp[ran] += wi*np.hypot(AS, AC) + annphase[ran] += wi*np.degrees(np.arctan2(AC, AS)) # semi-annual component SS,SC = beta_mat[SEAS+2:SEAS+4] dsemian[ran] += wi*np.dot(TMAT[:,SEAS+2:SEAS+4],[SS,SC]) - semiamp[ran] += wi*np.sqrt(SS**2 + SC**2) - semiphase[ran] += wi*np.arctan2(SC,SS)*180.0/np.pi + semiamp[ran] += wi*np.hypot(SS, SC) + semiphase[ran] += wi*np.degrees(np.arctan2(SC, SS)) # add weights weight[ran] += wi # divide weighted smoothed time-series by weights @@ -326,13 +327,13 @@ def smooth(t_in, d_in, HFWTH=6, MOVING=False, DATA_ERR=0, WEIGHT=0, # annual component AS,AC = beta_mat[SEAS:SEAS+2] dannual[i] = np.dot(TMAT[HFWTH,SEAS:SEAS+2],[AS,AC]) - annphase[i] = np.arctan2(AC,AS)*180.0/np.pi - annamp[i] = np.sqrt(AS**2 + AC**2) + annphase[i] = np.degrees(np.arctan2(AC, AS)) + annamp[i] = np.hypot(AS, AC) # semi-annual component SS,SC = beta_mat[SEAS+2:SEAS+4] dsemian[i] = np.dot(TMAT[HFWTH,SEAS+2:SEAS+4],[SS,SC]) - semiamp[i] = np.sqrt(SS**2 + SC**2) - semiphase[i] = np.arctan2(SC,SS)*180.0/np.pi + semiamp[i] = np.hypot(SS, SC) + semiphase[i] = np.degrees(np.arctan2(SC, SS)) # noise component dnoise[i] = d_in[i+HFWTH] - dsmth[i] - dseason[i] # reduced time-series diff --git a/gravity_toolkit/tools.py b/gravity_toolkit/tools.py index 36c93f0e..16cc786a 100644 --- a/gravity_toolkit/tools.py +++ b/gravity_toolkit/tools.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" tools.py -Written by Tyler Sutterley (11/2024) +Written by Tyler Sutterley (07/2026) Jupyter notebook, user interface and plotting tools PYTHON DEPENDENCIES: @@ -27,6 +27,7 @@ utilities.py: download and management utilities for files UPDATE HISTORY: + Updated 07/2026: use np.radians and np.degrees for angle conversions Updated 11/2024: fix deprecated widget object copies Updated 04/2024: add widget for setting endpoint for accessing PODAAC data place colormap registration within try/except to check for existing @@ -1123,9 +1124,9 @@ def wrap_longitudes(lon): lon: np.ndarray longitude (degrees east) """ - phi = np.arctan2(np.sin(lon*np.pi/180.0), np.cos(lon*np.pi/180.0)) + phi = np.arctan2(np.sin(np.radians(lon)), np.cos(np.radians(lon))) # convert phi from radians to degrees - return phi*180.0/np.pi + return np.degrees(phi) # PURPOSE: parallels the matplotlib basemap shiftgrid function def shift_grid(lon0, data, lon, CYCLIC=360.0): diff --git a/mapping/plot_AIS_GrIS_maps.py b/mapping/plot_AIS_GrIS_maps.py index a082fb43..381638cf 100644 --- a/mapping/plot_AIS_GrIS_maps.py +++ b/mapping/plot_AIS_GrIS_maps.py @@ -578,9 +578,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave=np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_AIS_grid_3maps.py b/mapping/plot_AIS_grid_3maps.py index 352ea0db..1e12ed3c 100644 --- a/mapping/plot_AIS_grid_3maps.py +++ b/mapping/plot_AIS_grid_3maps.py @@ -473,9 +473,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave=np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_AIS_grid_4maps.py b/mapping/plot_AIS_grid_4maps.py index 550ef497..8aa43cf8 100644 --- a/mapping/plot_AIS_grid_4maps.py +++ b/mapping/plot_AIS_grid_4maps.py @@ -473,11 +473,11 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average - ave=np.sum(area*data[indy,indx])/np.sum(area) + ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average ax1.contour(lon,lat,data,[ave],colors='blue',linestyles='solid', linewidths=1.5,transform=ccrs.PlateCarree()) diff --git a/mapping/plot_AIS_grid_maps.py b/mapping/plot_AIS_grid_maps.py index 0125d3b1..e73e5a90 100644 --- a/mapping/plot_AIS_grid_maps.py +++ b/mapping/plot_AIS_grid_maps.py @@ -461,9 +461,9 @@ def plot_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_AIS_grid_movie.py b/mapping/plot_AIS_grid_movie.py index ce956b41..0b0fe5fa 100644 --- a/mapping/plot_AIS_grid_movie.py +++ b/mapping/plot_AIS_grid_movie.py @@ -583,9 +583,9 @@ def animate_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_AIS_regional_maps.py b/mapping/plot_AIS_regional_maps.py index 4196dcd9..139f5974 100644 --- a/mapping/plot_AIS_regional_maps.py +++ b/mapping/plot_AIS_regional_maps.py @@ -544,9 +544,9 @@ def plot_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_AIS_regional_movie.py b/mapping/plot_AIS_regional_movie.py index ac966ddd..6555577e 100644 --- a/mapping/plot_AIS_regional_movie.py +++ b/mapping/plot_AIS_regional_movie.py @@ -672,9 +672,9 @@ def animate_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_GrIS_grid_3maps.py b/mapping/plot_GrIS_grid_3maps.py index 5e80b57f..9fb522ab 100644 --- a/mapping/plot_GrIS_grid_3maps.py +++ b/mapping/plot_GrIS_grid_3maps.py @@ -490,9 +490,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave=np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_GrIS_grid_5maps.py b/mapping/plot_GrIS_grid_5maps.py index 9d804ab4..382569f9 100644 --- a/mapping/plot_GrIS_grid_5maps.py +++ b/mapping/plot_GrIS_grid_5maps.py @@ -474,9 +474,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave=np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_GrIS_grid_maps.py b/mapping/plot_GrIS_grid_maps.py index d42373bc..ca3de2ac 100644 --- a/mapping/plot_GrIS_grid_maps.py +++ b/mapping/plot_GrIS_grid_maps.py @@ -474,9 +474,9 @@ def plot_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_GrIS_grid_movie.py b/mapping/plot_GrIS_grid_movie.py index b2ebc63e..63fa29e8 100644 --- a/mapping/plot_GrIS_grid_movie.py +++ b/mapping/plot_GrIS_grid_movie.py @@ -605,9 +605,9 @@ def animate_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_QML_grid_3maps.py b/mapping/plot_QML_grid_3maps.py index ae50db05..71eca64c 100644 --- a/mapping/plot_QML_grid_3maps.py +++ b/mapping/plot_QML_grid_3maps.py @@ -479,9 +479,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(data.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave=np.sum(area*data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_global_grid_3maps.py b/mapping/plot_global_grid_3maps.py index d8736213..5594e476 100644 --- a/mapping/plot_global_grid_3maps.py +++ b/mapping/plot_global_grid_3maps.py @@ -346,9 +346,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(dinput.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*dinput.data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_global_grid_4maps.py b/mapping/plot_global_grid_4maps.py index 96fcf4c4..517db035 100644 --- a/mapping/plot_global_grid_4maps.py +++ b/mapping/plot_global_grid_4maps.py @@ -346,9 +346,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(dinput.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*dinput.data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_global_grid_5maps.py b/mapping/plot_global_grid_5maps.py index 5deef916..d63e1bb0 100644 --- a/mapping/plot_global_grid_5maps.py +++ b/mapping/plot_global_grid_5maps.py @@ -352,9 +352,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(dinput.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*dinput.data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_global_grid_9maps.py b/mapping/plot_global_grid_9maps.py index 4406f4ce..2b1dba5e 100644 --- a/mapping/plot_global_grid_9maps.py +++ b/mapping/plot_global_grid_9maps.py @@ -347,9 +347,9 @@ def plot_grid(base_dir, FILENAMES, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(dinput.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*dinput.data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_global_grid_maps.py b/mapping/plot_global_grid_maps.py index 761abd59..da57df62 100644 --- a/mapping/plot_global_grid_maps.py +++ b/mapping/plot_global_grid_maps.py @@ -320,9 +320,9 @@ def plot_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(dinput.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*dinput.data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/mapping/plot_global_grid_movie.py b/mapping/plot_global_grid_movie.py index eff08e2e..40b98eb3 100644 --- a/mapping/plot_global_grid_movie.py +++ b/mapping/plot_global_grid_movie.py @@ -420,9 +420,9 @@ def animate_grid(base_dir, FILENAME, # plot line contour for global average if MEAN_CONTOUR and CONTOURS: # calculate areas of each grid cell - dphi,dth = (dlon*np.pi/180.0,dlat*np.pi/180.0) + dphi,dth = (np.radians(dlon), np.radians(dlat)) indy,indx = np.nonzero(np.logical_not(subset.mask)) - area = (rad_e**2)*dth*dphi*np.cos(lat[indy,indx]*np.pi/180.0) + area = (rad_e**2)*dth*dphi*np.cos(np.radians(lat[indy,indx])) # calculate average ave = np.sum(area*subset.data[indy,indx])/np.sum(area) # plot line contour of global average diff --git a/scripts/calc_sensitivity_kernel.py b/scripts/calc_sensitivity_kernel.py index f3563e2f..840d7e56 100644 --- a/scripts/calc_sensitivity_kernel.py +++ b/scripts/calc_sensitivity_kernel.py @@ -441,7 +441,7 @@ def calc_sensitivity_kernel(LMAX, RAD, n_lat = len(grid.lat) # Computing plms for converting to spatial domain - theta = (90.0-grid.lat)*np.pi/180.0 + theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # for each mascon diff --git a/scripts/combine_HEX_spherical_caps.py b/scripts/combine_HEX_spherical_caps.py index c1d8dc2c..f2874129 100644 --- a/scripts/combine_HEX_spherical_caps.py +++ b/scripts/combine_HEX_spherical_caps.py @@ -305,7 +305,7 @@ def combine_HEX_spherical_caps(PROC, DREL, DSET, LMAX, RAD, # output GRACE regional time-series total_mass = grace_reg[i] total_thick = 1e15*total_mass/area_reg[i] - total_error = np.sqrt(grace_error[i]**2 + statistical_reg[i]**2) + total_error = np.hypot(grace_error[i], statistical_reg[i]) thick_error = 1e15*total_error/area_reg[i] # total area in kilometers^2 area_km = area_reg[i]/1e10 diff --git a/scripts/combine_harmonics.py b/scripts/combine_harmonics.py index 17df8b35..714cff2f 100644 --- a/scripts/combine_harmonics.py +++ b/scripts/combine_harmonics.py @@ -273,7 +273,7 @@ def combine_harmonics(INPUT_FILE, OUTPUT_FILE, attributes['ROOT']['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2' # Computing plms for converting to spatial domain - theta = (90.0 - grid.lat)*np.pi/180.0 + theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # converting harmonics to truncated, smoothed coefficients in output units diff --git a/scripts/combine_sea_level_data.py b/scripts/combine_sea_level_data.py index f1000855..0465c8e1 100755 --- a/scripts/combine_sea_level_data.py +++ b/scripts/combine_sea_level_data.py @@ -109,7 +109,7 @@ def combine_sea_level_data(index_file, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # longitude and colatitude in radians - th = (90.0 - np.squeeze(landsea.lat))*np.pi/180.0 + th = np.radians(90.0 - np.squeeze(landsea.lat)) # Calculating Legendre Polynomials using Holmes and Featherstone relation PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(th)) diff --git a/scripts/convert_harmonics.py b/scripts/convert_harmonics.py index 482ade8c..c6d7eeee 100644 --- a/scripts/convert_harmonics.py +++ b/scripts/convert_harmonics.py @@ -187,7 +187,7 @@ def convert_harmonics(INPUT_FILE, OUTPUT_FILE, attributes['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2' # calculate associated Legendre polynomials - th = (90.0 - input_spatial.lat)*np.pi/180.0 + th = np.radians(90.0 - input_spatial.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(th)) # create list of harmonics objects diff --git a/scripts/grace_spatial_error.py b/scripts/grace_spatial_error.py index f20672b5..149e9270 100755 --- a/scripts/grace_spatial_error.py +++ b/scripts/grace_spatial_error.py @@ -419,11 +419,11 @@ def grace_spatial_error(base_dir, PROC, DREL, DSET, LMAX, RAD, delta.attributes['ROOT'] = attributes # Computing plms for converting to spatial domain - phi = delta.lon[np.newaxis,:]*np.pi/180.0 - theta = (90.0 - delta.lat)*np.pi/180.0 + phi = np.radians(delta.lon[np.newaxis,:]) + theta = np.radians(90.0 - delta.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # square of legendre polynomials truncated to order MMAX - mm = np.arange(0,MMAX+1) + mm = np.arange(0, MMAX+1) PLM2 = PLM[:,mm,:]**2 # Calculating cos(m*phi)^2 and sin(m*phi)^2 diff --git a/scripts/grace_spatial_maps.py b/scripts/grace_spatial_maps.py index ea65903a..d4544897 100755 --- a/scripts/grace_spatial_maps.py +++ b/scripts/grace_spatial_maps.py @@ -426,7 +426,7 @@ def grace_spatial_maps(base_dir, PROC, DREL, DSET, LMAX, RAD, nlat = len(grid.lat) # Computing plms for converting to spatial domain - theta = (90.0-grid.lat)*np.pi/180.0 + theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # output spatial units diff --git a/scripts/piecewise_grace_maps.py b/scripts/piecewise_grace_maps.py index f936dc1f..39df4b75 100755 --- a/scripts/piecewise_grace_maps.py +++ b/scripts/piecewise_grace_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" piecewise_grace_maps.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (07/2026) Reads in GRACE/GRACE-FO spatial files and fits a piecewise regression model at each grid point for breakpoint analysis @@ -60,6 +60,7 @@ utilities.py: download and management utilities for files UPDATE HISTORY: + Updated 07/2026: use np.hypot to calculate the sum of two squares Updated 06/2023: append amplitude and phase titles when creating flags more tidal aliasing periods using values from Ray and Luthcke (2006) Updated 05/2023: split S2 tidal aliasing terms into GRACE and GRACE-FO eras @@ -216,7 +217,7 @@ def piecewise_grace_maps(LMAX, RAD, # find index of breakpoint within GRACE/GRACE-FO months if BREAKPOINT not in grid.month: raise ValueError(f'{BREAKPOINT} not found in GRACE/GRACE-FO months') - breakpoint_index, = np.nonzero(grid.month == BREAKPOINT) + breakpoint_index, = np.flatnonzero(grid.month == BREAKPOINT) # Setting output parameters coef_str = ['x0', 'px1', 'px1'] @@ -366,16 +367,15 @@ def piecewise_grace_maps(LMAX, RAD, out.data[:,:,j], out.data[:,:,j+1] ) # convert phase from -180:180 to 0:360 - ii,jj = np.nonzero(ph.data < 0) - ph.data[ii,jj] += 360.0 + ph.data = np.where(ph.data < 0, ph.data + 360.0, ph.data) # Amplitude Error comp1 = out.error[:,:,j]*out.data[:,:,j]/amp.data comp2 = out.error[:,:,j+1]*out.data[:,:,j+1]/amp.data - amp.error = np.sqrt(comp1**2 + comp2**2) + amp.error = np.hypot(comp1, comp2) # Phase Error (degrees) comp1 = out.error[:,:,j]*out.data[:,:,j+1]/(amp.data**2) comp2 = out.error[:,:,j+1]*out.data[:,:,j]/(amp.data**2) - ph.error = (180.0/np.pi)*np.sqrt(comp1**2 + comp2**2) + ph.error = np.degrees(np.hypot(comp1, comp2)) # output file names for amplitude, phase and errors f2 = (FILE_PREFIX, units, LMAX, order_str, diff --git a/scripts/plot_SLR_azimuthal.py b/scripts/plot_SLR_azimuthal.py index 21d55ac3..55746faf 100644 --- a/scripts/plot_SLR_azimuthal.py +++ b/scripts/plot_SLR_azimuthal.py @@ -103,7 +103,7 @@ def plot_SLR_azimuthal(base_dir, PROC, DREL, START_MON, END_MON, MISSING): 'CSR SLR','GSFC SLR','GFZ GravIS'] plot_zorder = [0,3,1,2] fig_text = {'C21':'a)','S21':'b)'} - plot_ylabel = {'C21':'C$\mathregular{_{21}}$','S21':'S$\mathregular{_{21}}$'} + plot_ylabel = {'C21':r'C$\mathregular{_{21}}$','S21':r'S$\mathregular{_{21}}$'} CSR = dict(time=CSR_CS2['time'],month=CSR_CS2['month']) GSFC = dict(time=GSFC_CS2['time'],month=GSFC_CS2['month']) GFZ = dict(time=GFZ_CS2['time'],month=GFZ_CS2['month']) diff --git a/scripts/regional_spherical_caps.py b/scripts/regional_spherical_caps.py index bf03edde..7a332e8a 100644 --- a/scripts/regional_spherical_caps.py +++ b/scripts/regional_spherical_caps.py @@ -169,7 +169,7 @@ def regional_spherical_caps(LMAX, rad_e = dfactor.rad_e# Average Radius of the Earth [cm] # colatitude of the spherical cap centers - th = (90.0 - lat)*np.pi/180.0 + th = np.radians(90.0 - lat) # Legendre polynomials of the spherical caps PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(th)) @@ -195,7 +195,7 @@ def regional_spherical_caps(LMAX, # area = (volume*density)/(mass/area) ar = 4.0*np.pi*(rad_e**3.0)*rho_e*np.squeeze(Ylms.clm[0,0])/3.0 # calculate equivalent radius (should equal RAD_CAP) - rad = np.sqrt(ar/np.pi)/rad_e*180.0/np.pi + rad = np.degrees(np.sqrt(ar/np.pi) / rad_e) # if verbose output: sanity check of radii args = (cap_number, RAD_CAP[i], rad) logging.info('{0:4d} {1:10.4f} {2:10.4f}'.format(*args)) diff --git a/scripts/regress_grace_maps.py b/scripts/regress_grace_maps.py index 62be313b..1e65a043 100755 --- a/scripts/regress_grace_maps.py +++ b/scripts/regress_grace_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" regress_grace_maps.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (07/2026) Reads in GRACE/GRACE-FO spatial files and fits a regression model at each grid point @@ -60,6 +60,7 @@ utilities.py: download and management utilities for files UPDATE HISTORY: + Updated 07/2026: use np.hypot to calculate the sum of two squares Updated 06/2023: append amplitude and phase titles when creating flags more tidal aliasing periods using values from Ray and Luthcke (2006) Updated 05/2023: split S2 tidal aliasing terms into GRACE and GRACE-FO eras @@ -358,16 +359,15 @@ def regress_grace_maps(LMAX, RAD, out.data[:,:,j], out.data[:,:,j+1] ) # convert phase from -180:180 to 0:360 - ii,jj = np.nonzero(ph.data < 0) - ph.data[ii,jj] += 360.0 + ph.data = np.where(ph.data < 0, ph.data + 360.0, ph.data) # Amplitude Error comp1 = out.error[:,:,j]*out.data[:,:,j]/amp.data comp2 = out.error[:,:,j+1]*out.data[:,:,j+1]/amp.data - amp.error = np.sqrt(comp1**2 + comp2**2) + amp.error = np.hypot(comp1, comp2) # Phase Error (degrees) comp1 = out.error[:,:,j]*out.data[:,:,j+1]/(amp.data**2) comp2 = out.error[:,:,j+1]*out.data[:,:,j]/(amp.data**2) - ph.error = (180.0/np.pi)*np.sqrt(comp1**2 + comp2**2) + ph.error = np.degrees(np.hypot(comp1, comp2)) # output file names for amplitude, phase and errors f2 = (FILE_PREFIX, units, LMAX, order_str, diff --git a/scripts/run_sea_level_equation.py b/scripts/run_sea_level_equation.py index e4a0e27d..fdeb7f46 100644 --- a/scripts/run_sea_level_equation.py +++ b/scripts/run_sea_level_equation.py @@ -178,7 +178,7 @@ def run_sea_level_equation(INPUT_FILE, OUTPUT_FILE, nth,nphi = landsea.shape land_function = np.zeros((nth, nphi), dtype=np.float64) # calculate colatitude in radians - th = (90.0 - landsea.lat)*np.pi/180.0 + th = np.radians(90.0 - landsea.lat) # extract land function from file # combine land and island levels for land function indx,indy = np.nonzero((landsea.data >= 1) & (landsea.data <= 3)) diff --git a/scripts/scale_grace_maps.py b/scripts/scale_grace_maps.py index 3f29117e..62552bc9 100644 --- a/scripts/scale_grace_maps.py +++ b/scripts/scale_grace_maps.py @@ -478,11 +478,11 @@ def scale_grace_maps(base_dir, PROC, DREL, DSET, LMAX, RAD, grid.mask = np.zeros((nlat, nlon, nfiles),dtype=bool) # Computing plms for converting to spatial domain - phi = grid.lon[np.newaxis,:]*np.pi/180.0 - theta = (90.0-grid.lat)*np.pi/180.0 + phi = np.radians(grid.lon[np.newaxis,:]) + theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # square of legendre polynomials truncated to order MMAX - mm = np.arange(0,MMAX+1) + mm = np.arange(0, MMAX+1) PLM2 = PLM[:,mm,:]**2 # dfactor is the degree dependent coefficients diff --git a/scripts/sea_level_error.py b/scripts/sea_level_error.py index 40d83bde..8681073b 100755 --- a/scripts/sea_level_error.py +++ b/scripts/sea_level_error.py @@ -140,7 +140,7 @@ def sea_level_error(PROC, DREL, DSET, LMAX, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # calculate Fully-Normalized Legendre Polynomials - th = (90.0 - landsea.lat)*np.pi/180.0 + th = np.radians(90.0 - landsea.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(th)) # create index file for least_squares_mascons.py diff --git a/scripts/sea_level_regress.py b/scripts/sea_level_regress.py index 33136b55..bac512d7 100644 --- a/scripts/sea_level_regress.py +++ b/scripts/sea_level_regress.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" sea_level_regress.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (07/2026) Reads in sea level grid files and calculates the trends at each grid point following an input regression model @@ -57,6 +57,7 @@ spatial.py: spatial data class for reading, writing and processing data UPDATE HISTORY: + Updated 07/2026: use np.hypot to calculate the sum of two squares Updated 06/2023: append amplitude and phase titles when creating flags more tidal aliasing periods using values from Ray and Luthcke (2006) Updated 05/2023: split S2 tidal aliasing terms into GRACE and GRACE-FO eras @@ -345,16 +346,17 @@ def sea_level_regress(PROC, DREL, DSET, LMAX, out.data[indy,indx,j], out.data[indy,indx,j+1] ) # convert phase from -180:180 to 0:360 - ii,jj = np.nonzero((ph.data < 0) & np.logical_not(ph.mask)) - ph.data[ii,jj] += 360.0 + ph.data = np.where( + (ph.data < 0) & np.logical_not(ph.mask), ph.data + 360.0, ph.data + ) # Amplitude Error comp1=out.error[indy,indx,j]*out.data[indy,indx,j]/amp.data[indy,indx] comp2=out.error[indy,indx,j+1]*out.data[indy,indx,j+1]/amp.data[indy,indx] - amp.error[indy,indx] = np.sqrt(comp1**2 + comp2**2) + amp.error[indy,indx] = np.hypot(comp1, comp2) # Phase Error (degrees) comp1=out.error[indy,indx,j]*out.data[indy,indx,j+1]/(amp.data[indy,indx]**2) comp2=out.error[indy,indx,j+1]*out.data[indy,indx,j]/(amp.data[indy,indx]**2) - ph.error[indy,indx] = (180.0/np.pi)*np.sqrt(comp1**2 + comp2**2) + ph.error[indy,indx] = np.degrees(np.hypot(comp1, comp2)) # output file names for amplitude and phase f2 = (ITERATION, dset_str, gia_str, ocean_str, EXPANSION, diff --git a/scripts/sea_level_stokes.py b/scripts/sea_level_stokes.py index 235e6be7..67545289 100755 --- a/scripts/sea_level_stokes.py +++ b/scripts/sea_level_stokes.py @@ -172,7 +172,7 @@ def sea_level_stokes(PROC, DREL, DSET, LMAX, dlon,dlat = landsea.spacing nlat, nlon = landsea.shape # calculate Fully-Normalized Legendre Polynomials - th = (90.0 - landsea.lat)*np.pi/180.0 + th = np.radians(90.0 - landsea.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(th)) # create index file for calc_mascon.py diff --git a/test/land.fcn.1_deg.gz b/test/land.fcn.1_deg.gz new file mode 100755 index 00000000..fc977c42 Binary files /dev/null and b/test/land.fcn.1_deg.gz differ diff --git a/test/out.slf.green_ice.1_deg.2008.60.gz b/test/out.slf.green_ice.1_deg.2008.60.gz new file mode 100644 index 00000000..0f13a933 Binary files /dev/null and b/test/out.slf.green_ice.1_deg.2008.60.gz differ diff --git a/test/test_download_and_read.py b/test/test_download_and_read.py index eb3a311e..cd83b698 100644 --- a/test/test_download_and_read.py +++ b/test/test_download_and_read.py @@ -3,6 +3,7 @@ test_download_and_read.py (11/2021) Tests the read program to verify that coefficients are being extracted """ +import pytest import pathlib import posixpath import gravity_toolkit as gravtk @@ -42,6 +43,7 @@ def test_gfz_http_download_and_read(): assert (Ylms['clm'][2,0] == -0.484169355584e-03) # PURPOSE: Download a GRACE file from GFZ and check that read program runs +@pytest.mark.skip(reason="Deprecated GFZ FTP server") def test_gfz_ftp_download_and_read(): HOST=['isdcftp.gfz-potsdam.de','grace','Level-2','CSR','RL06', 'GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gz'] @@ -90,7 +92,7 @@ def test_esa_swarm_download_and_read(): # PURPOSE: Download a GRACE ITSG GRAZ file and check that read program runs def test_itsg_graz_download_and_read(): - HOST=['http://ftp.tugraz.at','outgoing','ITSG','GRACE', + HOST=['http://ftp.tugraz.at','pub','ITSG','GRACE', 'ITSG-Grace_operational','monthly','monthly_n60', 'ITSG-Grace_operational_n60_2018-06.gfc'] # download and read as virtual file object diff --git a/test/test_harmonics.py b/test/test_harmonics.py index 62c931b1..34de7e0d 100755 --- a/test/test_harmonics.py +++ b/test/test_harmonics.py @@ -39,7 +39,7 @@ def test_harmonics(): RAD = 250.0 # calculate colatitudes of input distribution - theta = (90.0 - input_distribution.lat)*np.pi/180.0 + theta = np.radians(90.0 - input_distribution.lat) # use fortran thresholds for colatitude bounds theta[theta > np.arccos(-0.9999999)] = np.arccos(-0.9999999) theta[theta < np.arccos(0.9999999)] = np.arccos(0.9999999) @@ -68,13 +68,14 @@ def test_harmonics(): dfactor = gravtk.units(lmax=LMAX).harmonic(*LOVE) wt = 2.0*np.pi*gravtk.gauss_weights(RAD,LMAX) smooth_Ylms.convolve(dfactor.cmwe*wt) + transform_Ylms = smooth_Ylms.copy() # convert harmonics back to spatial domain at same grid spacing test_distribution = gravtk.harmonic_summation(smooth_Ylms.clm, smooth_Ylms.slm, input_distribution.lon, input_distribution.lat, LMAX=LMAX, PLM=PLM).T # convert harmonics using fast-fourier transform method - test_transform = gravtk.harmonic_transform(smooth_Ylms.clm, - smooth_Ylms.slm, input_distribution.lon, input_distribution.lat, + test_transform = gravtk.harmonic_transform(transform_Ylms.clm, + transform_Ylms.slm, input_distribution.lon, input_distribution.lat, LMAX=LMAX, PLM=PLM).T # convert harmonics back to spatial domain using wrapper function test_combine = gravtk.stokes_summation(test_Ylms.clm, diff --git a/test/test_masks.py b/test/test_masks.py index 70c1f923..d04cb873 100644 --- a/test/test_masks.py +++ b/test/test_masks.py @@ -25,10 +25,10 @@ def test_lsmask(LANDMASK): nlat, nlon = landsea.shape # colatitude in radians gridlon, gridlat = np.meshgrid(landsea.lon, landsea.lat) - th = (90.0 - gridlat)*np.pi/180.0 + th = np.radians(90.0 - gridlat) # grid spacing in radians - dphi = np.pi*np.abs(dlon)/180.0 - dth = np.pi*np.abs(dlat)/180.0 + dphi = np.radians(np.abs(dlon)) + dth = np.radians(np.abs(dlat)) # create land function land_function = np.zeros((nlat, nlon), dtype=np.float64) # combine land and island levels for land function diff --git a/test/test_sea_level.py b/test/test_sea_level.py new file mode 100644 index 00000000..ebe093ad --- /dev/null +++ b/test/test_sea_level.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python +u""" +test_sea_level.py (07/2026) +""" +import inspect +import pathlib +import numpy as np +import gravity_toolkit as gravtk + +# path to test files +filename = inspect.getframeinfo(inspect.currentframe()).filename +filepath = pathlib.Path(filename).absolute().parent + +# PURPOSE: test sea level equation programs +def test_sea_level(): + # path to load Love numbers file + love_numbers_file = gravtk.utilities.get_data_path( + ['data','love_numbers']) + # read load Love numbers + LOVE = gravtk.read_love_numbers(love_numbers_file, FORMAT='class') + + # read land function file + LANDMASK = filepath.joinpath('land.fcn.1_deg.gz') + landsea = gravtk.spatial().from_ascii(LANDMASK, + date=False, spacing=[1.0, 1.0], nlat=180, nlon=360, + extent=[0.5,359.5,-89.5,89.5], compression='gzip') + + # spherical harmonic parameters + # maximum spherical harmonic degree + LMAX = 60 + # read harmonics from file + harmonics_file = filepath.joinpath('out.geoid.green_ice.0.5.2008.60.gz') + Ylms = gravtk.harmonics(lmax=LMAX, mmax=LMAX).from_ascii( + harmonics_file, date=False, compression='gzip') + # calculate the legendre functions using Martin Mohlenkamp's relation + th = np.radians(90.0 - landsea.lat) + PLM, dPLM = gravtk.plm_mohlenkamp(LMAX, np.cos(th)) + + # run pseudo-spectral sea level equation solver + sea_level = gravtk.sea_level_equation(Ylms.clm, Ylms.slm, + landsea.lon, landsea.lat, landsea.data.T, LMAX=LMAX, + LOVE=LOVE, BODY_TIDE_LOVE=0, + FLUID_LOVE=0, DENSITY=1.0, POLAR=True, + PLM=PLM, ITERATIONS=2, FILL_VALUE=np.nan).T + + # check that sea level data is equal to file precision + valid_file = filepath.joinpath('out.slf.green_ice.1_deg.2008.60.gz') + validation = gravtk.spatial().from_ascii(valid_file, + date=False, spacing=[1.0, 1.0], nlat=180, nlon=360, + extent=[0.5,359.5,-89.5,89.5], compression='gzip') + # check differences + difference = validation.data - sea_level + valid_difference = difference[np.isfinite(difference)] + assert np.all(np.abs(valid_difference) < 1e-8) + +def test_harmonics(): + # read land function file + LANDMASK = filepath.joinpath('land.fcn.1_deg.gz') + landsea = gravtk.spatial().from_ascii(LANDMASK, + date=False, spacing=[1.0, 1.0], nlat=180, nlon=360, + extent=[0.5,359.5,-89.5,89.5], compression='gzip') + # calculate ocean function from land function + land_function = landsea.data.T + ocean_function = 1.0 - land_function + # maximum spherical harmonic degree + LMAX = 60 + # calculate spherical harmonics using integration and fourier methods + YlmI = gravtk.gen_harmonics(ocean_function, landsea.lon, landsea.lat, + LMAX=LMAX, METHOD="integration") + YlmF = gravtk.gen_harmonics(ocean_function, landsea.lon, landsea.lat, + LMAX=LMAX, METHOD="fourier") + # check that amplitudes of harmonic data are nearly equal + difference = YlmI.amplitude - YlmF.amplitude + harmonic_eps = np.finfo(np.float16).eps + assert np.all(np.abs(difference) < harmonic_eps) diff --git a/utilities/quick_mascon_regress.py b/utilities/quick_mascon_regress.py index 2fcad343..77524d39 100755 --- a/utilities/quick_mascon_regress.py +++ b/utilities/quick_mascon_regress.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" quick_mascon_regress.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (07/2026) Creates a regression summary file for a mascon time series file COMMAND LINE OPTIONS: @@ -22,6 +22,7 @@ time_series.amplitude.py: calculates the amplitude and phase of a harmonic UPDATE HISTORY: + Updated 07/2026: use np.hypot to calculate the sum of two squares Updated 06/2023: can choose different tidal aliasing periods Updated 05/2023: split S2 tidal aliasing terms into GRACE and GRACE-FO eras allow fit to be piecewise using a known breakpoint GRACE/GRACE-FO month @@ -185,11 +186,11 @@ def run_regress(input_file, # Amplitude Errors comp1 = fit[err][j]*fit['beta'][j]/amp['beta'][i] comp2 = fit[err][j+1]*fit['beta'][j+1]/amp['beta'][i] - amp[err][i] = np.sqrt(comp1**2 + comp2**2) + amp[err][i] = np.hypot(comp1, comp2) # Phase Error (degrees) comp1 = fit[err][j]*fit['beta'][j+1]/(amp['beta'][i]**2) comp2 = fit[err][j+1]*fit['beta'][j]/(amp['beta'][i]**2) - ph[err][i] = (180.0/np.pi)*np.sqrt(comp1**2 + comp2**2) + ph[err][i] = np.degrees(np.hypot(comp1, comp2)) # Amplitude, Error and Statistical Significance args = ('Ampl.','Estimate','Std. Error','95% Conf.','Units')