feat: add ability to use a custom weighting function (#181)

useful for integrating over gaussian grids

Author: tsutterley

doc/source/_assets/geoid_height.svg

@@ -112,7 +112,9 @@ def filter(self, record):
 pygments_style = 'native'
 bibtex_bibfiles = ['_assets/gravity-refs.bib']
 bibtex_default_style = 'plain'
+plot_formats = ['png']
 plot_html_show_formats = False
+plot_html_show_source_link = False
 numfig = True
 numfig_secnum_depth = 1
 

doc/source/conf.py

@@ -49,6 +49,7 @@
 UPDATE HISTORY:
     Updated 07/2026: use np.einsum for spherical harmonic summations
         use np.radians to convert from degrees to radians
+        added custom weighting function for gridded data
     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
@@ -129,7 +130,16 @@ def gen_harmonics(data, lon, lat, **kwargs):
     return Ylms
 
 
-def integration(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs):
+def integration(
+    data,
+    lon,
+    lat,
+    LMAX=60,
+    MMAX=None,
+    WEIGHT=None,
+    PLM=0,
+    **kwargs,
+):
     """
     Converts data from the spatial domain to spherical harmonic coefficients
 
@@ -145,6 +155,8 @@ def integration(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs):
         Upper bound of Spherical Harmonic Degrees
     MMAX: int or NoneType, default None
         Upper bound of Spherical Harmonic Orders
+    WEIGHT: np.ndarray or NoneType, default None
+        Custom latitudinal weighting function for gridded data
     PLM: float, default 0
         input Legendre polynomials
 
@@ -164,19 +176,28 @@ def integration(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs):
     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])
+    # grid dimensions
+    nlat = np.int64(len(th))
 
     # LMAX+1 as there are LMAX+1 elements between 0 and LMAX
     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
-    int_fact = np.sin(th) * dphi * dth
+    # use an integration factor for gridded data or
+    # calculate from sin(theta)*dtheta*dphi
+    int_fact = np.zeros((nlat))
+    if WEIGHT is not None:
+        # Weighting function for integrating gridded data
+        int_fact[:] = np.broadcast_to(np.atleast_1d(WEIGHT), nlat)
+    else:
+        # Multiplying sin(th) with differentials of theta and phi
+        # to calculate the integration factor at each latitude
+        dphi = np.abs(phi[1] - phi[0])
+        dth = np.abs(th[1] - th[0])
+        int_fact[:] = np.sin(th) * dphi * dth
+    # normalizing coefficients
     coeff = 1.0 / (4.0 * np.pi)
 
     # Calculate polynomials using Holmes and Featherstone (2002) relation
@@ -204,7 +225,15 @@ def integration(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs):
     return Ylms
 
 
-def fourier(data, lon, lat, LMAX=60, MMAX=None, PLM=0, **kwargs):
+def fourier(
+    data,
+    lon,
+    lat,
+    LMAX=60,
+    MMAX=None,
+    PLM=0,
+    **kwargs,
+):
     """
     Computes the spherical harmonic coefficients of a spatial field
 

gravity_toolkit/gen_harmonics.py

@@ -29,6 +29,7 @@
         2: Gigatonnes of mass
         3: kg/m^2
         list: custom degree-dependent unit conversion factor
+    WEIGHT: custom latitudinal weighting function for gridded data
     PLM: input Legendre polynomials
     LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
 
@@ -45,6 +46,7 @@
 UPDATE HISTORY:
     Updated 07/2026: use np.einsum for spherical harmonic summations
         use np.radians to convert from degrees to radians
+        added custom weighting function for gridded data
     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
@@ -82,7 +84,16 @@
 
 
 def gen_stokes(
-    data, lon, lat, LMIN=0, LMAX=60, MMAX=None, UNITS=1, PLM=None, LOVE=None
+    data,
+    lon,
+    lat,
+    LMIN=0,
+    LMAX=60,
+    MMAX=None,
+    UNITS=1,
+    WEIGHT=None,
+    PLM=None,
+    LOVE=None,
 ):
     r"""
     Converts data from the spatial domain to spherical harmonic
@@ -109,6 +120,8 @@ def gen_stokes(
             - ``2``: gigatonnes of mass (Gt)
             - ``3``:  mm water equivalent thickness (mm w.e., kg/m\ :sup:`2`)
             - list: custom degree-dependent unit conversion factor
+    WEIGHT: np.ndarray or NoneType, default None
+        Custom latitudinal weighting function for gridded data
     PLM: np.ndarray or NoneType, default None
         Input Legendre polynomials
     LOVE: tuple or NoneType, default None
@@ -132,8 +145,6 @@ def gen_stokes(
     # upper bound of spherical harmonic orders (default = LMAX)
     MMAX = np.copy(LMAX) if (MMAX is None) else MMAX
 
-    # grid dimensions
-    nlat = np.int64(len(lat))
     # Longitude in radians
     phi = np.radians(np.squeeze(lon.copy()))
     # reformatting longitudes to range 0:360 (if previously -180:180)
@@ -143,6 +154,8 @@ def gen_stokes(
     # grid step in radians
     dphi = np.abs(phi[1] - phi[0])
     dth = np.abs(th[1] - th[0])
+    # grid dimensions
+    nlat = np.int64(len(th))
 
     # reforming data to lonXlat if input latXlon
     sz = np.shape(data)
@@ -173,6 +186,10 @@ def gen_stokes(
         int_fact[:] = np.sin(th) * dphi * dth
     else:
         raise ValueError(f'Unknown units {UNITS}')
+    # use a custom weighting function for gridded data if provided
+    if WEIGHT is not None:
+        # Weighting function for integrating gridded data
+        int_fact[:] = np.broadcast_to(np.atleast_1d(WEIGHT), nlat)
 
     # Calculating cos/sin of phi arrays
     # output [m,phi]