refactor: use np.einsum for spherical harmonic summations

.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

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

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

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}},

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,

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",

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",