Download this Jupyter notebook from github

Convert Geopandas GeoSeries to spherical harmonics#

Author: R. Rietbroek Jan 2025 (r.rietbroek@utwente.nl)*

Often, geographical geometry objects are needed as masks (1 inside,0 outside) expressed in the spherical harmonic domain. This notebook shows how a geopandas Geoseries, containing polygon masks or points can be converted to spherical harmonic coefficients

1. Preparations#

[1]:

#Optionally enable autoreloading for development purposes. Note that this does not automagically reload the binary extensions
%load_ext autoreload
%autoreload 2

[2]:

import xarray as xr
import os
import shxarray
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

import numpy as np
import geopandas as gpd
from shapely import Point

[16]:

#set the maximum desired degree of the output spherical harmonics
nmax=120

/home/roelof/cld_UTwente/Soft/shxarray-git/src/shxarray/core/xr_accessor.py:18: AccessorRegistrationWarning: registration of accessor <class 'shxarray.core.xr_accessor.SHDaAccessor'> under name 'sh' for type <class 'xarray.core.dataarray.DataArray'> is overriding a preexisting attribute with the same name.
  @xr.register_dataarray_accessor("sh")
/home/roelof/cld_UTwente/Soft/shxarray-git/src/shxarray/core/xr_accessor.py:280: AccessorRegistrationWarning: registration of accessor <class 'shxarray.core.xr_accessor.SHDsAccessor'> under name 'sh' for type <class 'xarray.core.dataset.Dataset'> is overriding a preexisting attribute with the same name.
  def synthesis(self,lon=None, lat=None,engine="shlib",**kwargs):

2. Get a set of Antartic grounded drainage basins (as polygons) to use as continental loads#

Note that in this case we load the geometries from a geoslurp database, but the example can be easily modified to load shapes from e.g. geopackage files.

[4]:

# In this example we're loading the drainage basin outlines froma geoslurp database
# The original datasources can be found at https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-greenland-drainage-systems
from geoslurp.manager import GeoslurpManager
import geoslurp.tools.pandas
import pandas as pd
gsman=GeoslurpManager(dbalias="tunnelmarge")
tname="cryo.antarc_ddiv_icesat_grnd"

# One can populate the database like this:
#dsant=gsman.dataset(tname)
#dsant.pull()
#dsant.register()

#load the polygons from the databse in a geopandas dataframe NOTE: convert to polar stereographic projection because the polgyon contains the south Pole
#And this projection is needed when converting to spherical harmonics
dfant=pd.DataFrame.gslrp.load(gsman.conn,f"SELECT basinid,ST_transform(geom::geometry,3031) as geometry FROM {tname}")
# quick plot of the drainage basin outlines
ax=dfant.plot(edgecolor="tab:blue",facecolor="grey")
ax.set_title("Antarctic grounded drainage basins")

[4]:

Text(0.5, 1.0, 'Antarctic grounded drainage basins')

../_images/notebooks_Geometry2sphericalHarmonics_6_1.png

4 Check if the intermediate drainage basin results made sense#

A spatially discontinuous function cannnot be one-to-one mapped to spectrally truncated function, so some Gibb’s phenomena is to be expected

[12]:

# Sum all spherical harmonic coefficient sets for all basins and grid
import cartopy.crs as ccrs
dres=0.25
basinid=18
# basinid=3
daantgrd=dantsh.sel(basinid=basinid).sh.synthesis(lon=np.arange(-180+dres/2,180,dres),lat=np.arange(-90+dres/2,-60,dres))
daantgrd.name="Mask value"
# daantgrd=daantgrd[daantgrd.lat < -50]
ax = plt.axes(projection=ccrs.SouthPolarStereo())
daantgrd.plot(ax=ax,vmax=1,vmin=0,transform=ccrs.PlateCarree())
ax.add_geometries(dfant[dfant.basinid==basinid].geometry,edgecolor='tab:red',crs=ccrs.SouthPolarStereo(),facecolor="none",lw=2)
plt.title(f"Antarctic drainage basin {basinid} from sh nmax={nmax}")

[12]:

Text(0.5, 1.0, 'Antarctic drainage basin 18 from sh nmax=120')

../_images/notebooks_Geometry2sphericalHarmonics_10_1.png

6 Create a visualization of all the disk loads together#

[15]:

dadsks=dadskssh.sum("nlonlat").sh.synthesis()

proj=ccrs.NearsidePerspective(central_longitude=0.0, central_latitude=50, satellite_height=35785831)
fig, ax = plt.subplots(subplot_kw={'projection': proj},figsize=(12,8.5))

ax.coastlines()
dadsks.plot(ax=ax,transform=ccrs.PlateCarree())
gspnts.plot(ax=ax,transform=ccrs.PlateCarree())
ax.set_title("100 randomly positioned spherical harmonic disk loads")
plt.show()

../_images/notebooks_Geometry2sphericalHarmonics_15_0.png