fix geoid gravity functional

Author: strawpants

src/shtns_backend/shtns.py

@@ -368,7 +368,9 @@ def lonlat_grid(self,nmax=None,lon=None,lat=None,gtype="gauss"):
                 dslonlat=xr.Dataset(coords=dict(lon=("nlonlat",lon),lat=("nlonlat",lat)))
             else:
                 #we need to check whether the proposed grid is compatible with SHTns
-
+                #rough estimate of nmax needed
+                dres=min(np.median(np.diff(lon)),np.median(np.diff(lat)))
+                nmax=int(180/dres)
                 sh = shtns.sht(nmax, nmax, norm=shtns.sht_orthonormal+shtns.SHT_NO_CS_PHASE)
                 sh.set_grid(len(lat),len(lon),shtns_type)
                 dslonlat=xr.Dataset(coords=dict(lon=lon,lat=lat))

src/shxarray/earth/ellipsoid.py

@@ -0,0 +1,29 @@
+
+# This file is part of the shxarray software which is licensed
+# under the Apache License version 2.0 (see the LICENSE file in the main repository)
+# Copyright Roelof Rietbroek (r.rietbroek@utwente.nl), 2025
+#
+
+import xarray as xr
+import shxarray
+
+def get_GRS80_Stokes():
+    """
+    Return the GRS80 ellipsoid as Stokes coefficients
+    """
+    a = 6378137.0
+    f = 1.0 / 298.257222101
+    nmax=8
+    dsgrs80=xr.DataArray.sh.zeros(nmax)
+    dsgrs80.attrs['title']='GRS80'
+    dsgrs80.attrs['a_earth']=a
+    dsgrs80.attrs['f_earth']=f
+    #fill out derived values
+    dsgrs80.loc[dict(n=0,m=0)]=1
+    dsgrs80.loc[dict(n=2,m=0)]=-0.48416685489612e-03 
+    dsgrs80.loc[dict(n=4,m=0)]=0.79030407333333e-06 
+    dsgrs80.loc[dict(n=6,m=0)]=-0.16872510013651e-08 
+    dsgrs80.loc[dict(n=8,m=0)]=-0.34609833692685e-11 
+
+    return dsgrs80
+    

src/shxarray/kernels/gravfunctionals.py

@@ -33,7 +33,7 @@ class Stokes2Geoid(IsoKernelBase):
     transform=("stokes","geoid")
     def __init__(self,nmax):
         super().__init__()
-        self._dsiso=a_earth*xr.DataArray(np.ones([121]),dims=['n'],coords=dict(n=np.arange(nmax+1)))
+        self._dsiso=a_earth*xr.DataArray(np.ones([nmax+1]),dims=['n'],coords=dict(n=np.arange(nmax+1)))
 
 class Load2Geoid(IsoKernelBase):
     """Provides an isotropic kernel representing the transformation of a surface load (in m) to geoid height in meter"""