Download this Jupyter notebook from github

---

Create spherical harmonic coefficients from an ocean mask

Author: R. Rietbroek Nov 2024 (r.rietbroek@utwente.nl)

The ocean function is defined to be one over the ocean and 0 over land areas. This mask function can also be converted to a set of spherical hamronic coefficients. These may for example be used in applying the Sea level equation in the spherical harmonic domain.

1. Obtain a (fine ocean grid)

As a starting point we can take a masked grid which was generated using using an example recipe from geoslurp.

import requests
import os
import xarray as xr
import shxarray
import matplotlib.pyplot as mpl

url="https://github.com/strawpants/geoshapes/raw/refs/heads/master/raster/ocean/ne_10m_oceangrid_0.125.nc"
#temporary storage
outfile="/tmp/"+os.path.basename(url)
if not os.path.exists(outfile):
    r = requests.get(url)
    print(f"Downloading ocean grid {outfile}")
    with open(outfile,'wb') as fid:
        fid.write(r.content)
else:
    print(f"{outfile}, already downloaded")


dsoce=xr.open_dataset(outfile).fillna(0.0)
display(dsoce)

Downloading ocean grid /tmp/ne_10m_oceangrid_0.125.nc

<xarray.Dataset> Size: 33MB
Dimensions:      (latitude: 1440, longitude: 2880)
Coordinates:
  * longitude    (longitude) float64 23kB -179.9 -179.8 -179.7 ... 179.8 179.9
  * latitude     (latitude) float64 12kB -89.94 -89.81 -89.69 ... 89.81 89.94
    spatial_ref  int64 8B ...
Data variables:
    oceanfunc    (latitude, longitude) float64 33MB 0.0 0.0 0.0 ... 1.0 1.0 1.0
Attributes:
    Conventions:  CF-1.9
    title:        Gridded Ocean mask generated from natural earth 10m dataset
    institution:  Roelof Rietbroek <r.rietbroek@utwente.nl>, Faculty of Geoin...
    source:       geoslurp
    history:      2024-11-26 14:53:20.439413 geoslurp
    references:   https://geoslurp.wobbly.earth/en/latest/notebooks/OceanFunc...
    comment:      Auto generated

xarray.Dataset

• Dimensions:
  • latitude: 1440
  • longitude: 2880
• Coordinates: (3)
  • longitude
    (longitude)
    float64
    -179.9 -179.8 ... 179.8 179.9

    units :

    degree_east

    standard_name :

    longitude

    array([-179.9375, -179.8125, -179.6875, ...,  179.6875,  179.8125,  179.9375])

  • latitude
    (latitude)
    float64
    -89.94 -89.81 ... 89.81 89.94

    units :

    degree_north

    standard_name :

    latitude

    array([-89.9375, -89.8125, -89.6875, ...,  89.6875,  89.8125,  89.9375])

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    [1 values with dtype=int64]

• Data variables: (1)
  • oceanfunc
    (latitude, longitude)
    float64
    0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [1., 1., 1., ..., 1., 1., 1.],
           [1., 1., 1., ..., 1., 1., 1.],
           [1., 1., 1., ..., 1., 1., 1.]])

• Indexes: (2)
  • longitude
    PandasIndex

    PandasIndex(Index([-179.9375, -179.8125, -179.6875, -179.5625, -179.4375, -179.3125,
           -179.1875, -179.0625, -178.9375, -178.8125,
           ...
            178.8125,  178.9375,  179.0625,  179.1875,  179.3125,  179.4375,
            179.5625,  179.6875,  179.8125,  179.9375],
          dtype='float64', name='longitude', length=2880))

  • latitude
    PandasIndex

    PandasIndex(Index([-89.9375, -89.8125, -89.6875, -89.5625, -89.4375, -89.3125, -89.1875,
           -89.0625, -88.9375, -88.8125,
           ...
            88.8125,  88.9375,  89.0625,  89.1875,  89.3125,  89.4375,  89.5625,
            89.6875,  89.8125,  89.9375],
          dtype='float64', name='latitude', length=1440))

• Attributes: (7)

  Conventions :

  CF-1.9

  title :

  Gridded Ocean mask generated from natural earth 10m dataset

  institution :

  Roelof Rietbroek <r.rietbroek@utwente.nl>, Faculty of Geoinformation Science and Earth Observation, University of Twente, Netherlands

  source :

  geoslurp

  history :

  2024-11-26 14:53:20.439413 geoslurp

  references :

  https://geoslurp.wobbly.earth/en/latest/notebooks/OceanFunction.html

  comment :

  Auto generated

dsoce.oceanfunc.plot()

<matplotlib.collections.QuadMesh at 0x76e6edb34110>

2. Convert the ocean grid to spherical harmonic coefficients

nmax=300 # note: a large nmax can take quite a while
# default uses a naive integration approach
dsocesh=dsoce.oceanfunc.sh.analysis(nmax)
display(dsocesh)
print(f"C00 coefficient (~71% of Earth surface is ocean) {dsocesh.sel(n=0).data[0]}")

<xarray.DataArray (nm: 90601)> Size: 725kB
array([ 7.11538971e-01, -6.02066955e-02, -1.22835400e-01, ...,
        3.12961008e-05, -2.97188334e-05, -2.17756304e-04])
Coordinates:
  * nm       (nm) object 725kB MultiIndex
  * n        (nm) int64 725kB 0 1 1 1 2 2 2 2 ... 300 300 300 300 300 300 300
  * m        (nm) int64 725kB 0 -1 0 1 -2 -1 0 1 ... 294 295 296 297 298 299 300

xarray.DataArray

• nm: 90601

• 0.7115 -0.06021 -0.1228 -0.1088 ... 3.13e-05 -2.972e-05 -0.0002178

  array([ 7.11538971e-01, -6.02066955e-02, -1.22835400e-01, ...,
          3.12961008e-05, -2.97188334e-05, -2.17756304e-04])

• Coordinates: (3)
  • nm
    (nm)
    object
    MultiIndex

    array([(0, 0), (1, -1), (1, 0), ..., (300, 298), (300, 299), (300, 300)],
          dtype=object)

  • n
    (nm)
    int64
    0 1 1 1 2 2 ... 300 300 300 300 300

    array([  0,   1,   1, ..., 300, 300, 300])

  • m
    (nm)
    int64
    0 -1 0 1 -2 ... 296 297 298 299 300

    array([  0,  -1,   0, ..., 298, 299, 300])

• Indexes: (1)
  • nm
    n
    m
    PandasMultiIndex

    PandasIndex(MultiIndex([(  0,   0),
                (  1,  -1),
                (  1,   0),
                (  1,   1),
                (  2,  -2),
                (  2,  -1),
                (  2,   0),
                (  2,   1),
                (  2,   2),
                (  3,  -3),
                ...
                (300, 291),
                (300, 292),
                (300, 293),
                (300, 294),
                (300, 295),
                (300, 296),
                (300, 297),
                (300, 298),
                (300, 299),
                (300, 300)],
               name='nm', length=90601))

• Attributes: (0)

C00 coefficient (72% of Earth surface is ocean) 0.7115389707690188

3. Convert the resulting spherical harmonic coefficients back to a grid and visualize the result

dsoce_back=dsocesh.sh.synthesis(lon=dsoce.longitude,lat=dsoce.latitude)

dsoce_back.plot(vmin=0,vmax=1)
ca=mpl.gca()
ca.set_title(f"Ocean function after grid -> sh -> grid, nmax={nmax}, (min,max)=({dsoce_back.min().data.item():.2f},{dsoce_back.max().data.item():.2f})")

Text(0.5, 1.0, 'Ocean function after grid -> sh -> grid, nmax=300, (min,max)=(-0.42,1.30)')

Note that differences up to 0.3 are visible due to Gibb’s phenomena

4. Write the results together with the input grid to the original netcdf file

outfilesh=f"/tmp/ne_10m_oceansh_n{nmax}.nc"
dsocesh.name="oceansh"
dsocesh.reset_index("nm").to_netcdf(outfilesh)