Notebook

PySHBundle Tutorials
Tutorial 1: Introduction to Spherical Harmonics using PySHBundle

by: Abhishek Mhamane, IIT Kanpur
Date: July 1, 2023

In [1]:

import pandas as pd
import os
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pyshbundle

Initializing PySHbundle v0.3.0

In [3]:

from pyshbundle.shutils import plm
from pyshbundle.viz_utils import ylm, ylm_plot

What are Spherical Harmonics?¶

Solid and Surface spherical harmonics

Harmonics - solutions of Laplace equation
Spherical - they are in spherical coordinates ($\theta$, $\lambda$, r)
Solid - cause they are 3D functions
Surface - when we remove the radial part (or consider a unit sphere)

Representing Gravitational potential as spherical harmonic expansion¶

$$ V(r, \theta, \lambda) = \frac{GM}{R} \sum_{l=0}^{\infty} \left(\frac{R}{r} \right)^{l+1} \sum_{m=0}^{l} P_{lm}(\cos{\theta}) (C_{lm}\cos{m\lambda} + S_{lm}\sin{m\lambda}) $$

Where,

l is the spherical harmonic degree
m is the spherical harmonic order

Associated Legendre Polynomials¶

In [6]:

help(plm)

Help on function plm in module pyshbundle.shutils:

plm(l: <built-in function array>, m: int, thetaRAD, nargin, nargout)
    plm Fully normalized associated Legendre functions for a selected order M
    
    Args:
        l (numpy.array): Degree, but not necessarily monotonic.
               For l < m a vector of zeros will be returned.
        m (int): order. If absent, m = 0 is assumed.
        thetaRAD (numpy.array): co-latitude [rad] (vector)
        nargin (int): number of input argument
        nargout (int): number of output argument
    Returns:
        p : array
            plm fully normalized Legendre functions
        dp : array
            first derivative of plm
        ddp : array
            second derivative of plm
    
    Author:
        Vivek Kumar Yadav, Interdisciplinary Center for Water Research (ICWaR), Indian Institute of Science (IISc)

In [7]:

#
l = np.zeros((1,1))
l[0] = 7

l2 = np.zeros((1,1))
l2[0] = 7
n = [0, 7, 14, 20]
colatitude = np.deg2rad(np.arange(0, 180.5, 0.5))
lon = np.deg2rad(np.arange(-180, 180.5, 0.5))

In [8]:

print(l)

[[7.]]

In [10]:

# computing polynomial 
plm_7_0 = plm(l=l, m=0, thetaRAD=colatitude, nargin=1, nargout=1)
plm_7_4 = plm(l=l2, m=4, thetaRAD=colatitude, nargin=1, nargout=1)
plm_7_7 = plm(l=l2, m=7, thetaRAD=colatitude, nargin=1, nargout=1)

In [11]:

plm_7_0.shape

Out[11]:

(361, 1, 1)

In [12]:

# Plotting the plm 
fig, ax = plt.subplots(figsize=(8, 5))
plt.plot(np.rad2deg(colatitude), plm_7_0[:, 0, 0], label="PLM(7, 0)")
plt.plot(np.rad2deg(colatitude), plm_7_4[:, 0, 0], label="PLM(7, 4)")
plt.plot(np.rad2deg(colatitude), plm_7_7[:, 0, 0], label="PLM(7, 7)")

plt.hlines(y=0, xmin=0, xmax=180, linestyles='--', colors='k')
plt.xlim([0, 180])
plt.legend()

plt.title("Visualizing the PLM function -> number of zeros = degree - order")
plt.xlabel("Co-latitude")
plt.ylabel("plm(l, m)")

plt.show()

Check Orthogonality of Spherical Harmonic Coefficients¶

Refering to Stack Overflow post - https://stackoverflow.com/questions/39537794/orthogonality-issue-in-scipys-legendre-polynomials

Great GitHub Repo by @markmbaum for 2-D Orthogonal Spherical Harmonic. Find the docs for orthopoly - https://markmbaum.github.io/orthopoly/

Visualizing the orthogonality of Legendre Polynomial¶

In [14]:

# the plm implementation using a co-latitude vector (geodetic use-case)
l1 = np.zeros((1,1))
l1[0] = 2

l2 = np.zeros((1,1))
l2[0] = 4


colatitude = np.deg2rad(np.arange(0, 180.5, 0.5))

plm_2_0 = plm(l=l1, m=0, thetaRAD=colatitude, nargin=1, nargout=1)
plm_4_0 = plm(l=l2, m=0, thetaRAD=colatitude, nargin=1, nargout=1)

product_2_4 = plm_2_0*plm_4_0

# Plotting the plm 
fig, ax = plt.subplots(figsize=(9, 5))
plt.plot(np.rad2deg(colatitude), plm_2_0[:, 0, 0], label="PLM(2, 0)")
plt.plot(np.rad2deg(colatitude), plm_4_0[:, 0, 0], label="PLM(4, 0)")

plt.plot(np.rad2deg(colatitude), product_2_4[:, 0, 0], label="PLM(2, 0)*PLM(4, 0)")

plt.fill_between(np.rad2deg(colatitude),product_2_4[:, 0, 0],color='green',alpha=0.5)

plt.hlines(y=0, xmin=0, xmax=180, linestyles='--', colors='k')
plt.xlim([0, 180])
plt.legend()

plt.title("Orthogonality of Legendre Polynomials - different degree with same order")
plt.xlabel("Co-latitude")
plt.ylabel("plm(l, m)")

plt.show()

It can be observed that the product of PLM(2,0) and PLM4(4,0) yields a output symmetric about X-axis and also about co-latitude = 90 deg. Thus area under the plot turns out to be zero -> both are orthogonal functions

In [15]:

# the plm implementation using a co-latitude vector (geodetic use-case)

l3 = np.zeros((1,1))
l3[0] = 3

colatitude = np.deg2rad(np.arange(0, 180.5, 0.5))

plm_3_0 = plm(l=l3, m=0, thetaRAD=colatitude, nargin=1, nargout=1)

product_3_3 = plm_3_0*plm_3_0

# Plotting the plm 
fig, ax = plt.subplots(figsize=(9, 5))
plt.plot(np.rad2deg(colatitude), plm_2_0[:, 0, 0], label="PLM(3, 0)")

plt.plot(np.rad2deg(colatitude), product_3_3[:, 0, 0], label="PLM(3, 0)*PLM(3, 0)", color='green')

plt.fill_between(np.rad2deg(colatitude),product_3_3[:, 0, 0],color='green',alpha=0.5)

plt.hlines(y=0, xmin=0, xmax=180, linestyles='--', colors='k')
plt.xlim([0, 180])
plt.legend()

plt.title("Orthogonality of Legendre Polynomials (same order and same degree)")
plt.xlabel("Co-latitude")
plt.ylabel("plm(l, m)")

plt.show()

In this case the integral cannot be zero, its some finite positive quantity.

In [16]:

# the plm implementation using a co-latitude vector (geodetic use-case)
l1 = np.zeros((1,1))
l1[0] = 3

l2 = np.zeros((1,1))
l2[0] = 3


colatitude = np.deg2rad(np.arange(0, 180.5, 0.5))

plm_3_1 = plm(l=l1, m=1, thetaRAD=colatitude, nargin=1, nargout=1)
plm_3_0 = plm(l=l2, m=2, thetaRAD=colatitude, nargin=1, nargout=1)

product_33_10 = plm_3_1*plm_3_0

# Plotting the plm 
fig, ax = plt.subplots(figsize=(9, 5))
plt.plot(np.rad2deg(colatitude), plm_3_1[:, 0, 0], label="PLM(3, 1)")
plt.plot(np.rad2deg(colatitude), plm_3_0[:, 0, 0], label="PLM(3, 2)")

plt.plot(np.rad2deg(colatitude), product_33_10[:, 0, 0], label="PLM(3, 1)*PLM(3, 2)")

plt.fill_between(np.rad2deg(colatitude),product_33_10[:, 0, 0],color='green',alpha=0.5)

plt.hlines(y=0, xmin=0, xmax=180, linestyles='--', colors='k')
plt.xlim([0, 180])
plt.legend()

plt.title("Orthogonality of Legendre Polynomials - Same degree with different order")
plt.xlabel("Co-latitude")
plt.ylabel("plm(l, m)")

plt.show()

In [17]:

# the plm implementation using a co-latitude vector (geodetic use-case)
l1 = np.zeros((1,1))
l1[0] = 3

l2 = np.zeros((1,1))
l2[0] = 4


colatitude = np.deg2rad(np.arange(0, 180.5, 0.5))

plm_3_1 = plm(l=l1, m=1, thetaRAD=colatitude, nargin=1, nargout=1)
plm_3_0 = plm(l=l2, m=2, thetaRAD=colatitude, nargin=1, nargout=1)

product_33_10 = plm_3_1*plm_3_0

# Plotting the plm 
fig, ax = plt.subplots(figsize=(9, 5))
plt.plot(np.rad2deg(colatitude), plm_3_1[:, 0, 0], label="PLM(3, 1)")
plt.plot(np.rad2deg(colatitude), plm_3_0[:, 0, 0], label="PLM(4, 2)")

plt.plot(np.rad2deg(colatitude), product_33_10[:, 0, 0], label="PLM(3, 1)*PLM(3, 2)")

plt.fill_between(np.rad2deg(colatitude),product_33_10[:, 0, 0],color='green',alpha=0.5)

plt.hlines(y=0, xmin=0, xmax=180, linestyles='--', colors='k')
plt.xlim([0, 180])
plt.legend()

plt.title("Orthogonality of Legendre Polynomials - Different degree with different order")
plt.xlabel("Co-latitude")
plt.ylabel("plm(l, m)")

plt.show()

The shaded regions if symmetric (such that sum=0) indicates that both the polynomials are orthogonal. Else both the functions are not symmetric.

Surface Spherical Harmonics¶

similar to the plm and plmplot in SHBundle

$$ $$

In [18]:

ylmc, ylms = ylm(l=2, m=1)

In [19]:

ylmc.shape

Out[19]:

(37, 1, 73)

In [20]:

np.linspace(0,np.pi,37)

Out[20]:

array([0.        , 0.08726646, 0.17453293, 0.26179939, 0.34906585,
       0.43633231, 0.52359878, 0.61086524, 0.6981317 , 0.78539816,
       0.87266463, 0.95993109, 1.04719755, 1.13446401, 1.22173048,
       1.30899694, 1.3962634 , 1.48352986, 1.57079633, 1.65806279,
       1.74532925, 1.83259571, 1.91986218, 2.00712864, 2.0943951 ,
       2.18166156, 2.26892803, 2.35619449, 2.44346095, 2.53072742,
       2.61799388, 2.70526034, 2.7925268 , 2.87979327, 2.96705973,
       3.05432619, 3.14159265])

In [21]:

# basic plotting of SH coeff on plane
l = 0
m = 0

ylmc_00, ylms_00 = ylm(l=l, m=m)

fig = plt.figure(figsize=(15, 7.5))
ax = plt.axes(projection = ccrs.PlateCarree())

lons = np.linspace(-180, 180, 73)
lats = np.linspace(-90, 90, 37)

x, y = np.meshgrid(lons, lats)

if m >=0 :

    img_extent = (-180, 180, -90, 90)

    # plot the data
    im = ax.imshow(ylmc_00[:, 0, :], origin='upper', extent=img_extent, transform=ccrs.PlateCarree(),)
else:
    plt.contourf(x, y, ylms_00[:, 0, :], cmap='RdYlBu_r')


# setting gridlines
gl = ax.gridlines(crs = ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, color='gray', alpha=0.9, linestyle='--')
# remove top x label
gl.top_labels = False
# change x label styles - font size ad colour
gl.xlabel_style = {'size':12,}
# left and right labels
gl.left_labels = True
gl.right_labels = False
# coastlines
ax.coastlines()

plt.colorbar(im, orientation='vertical', shrink=0.85, pad=0.02,label=f"[...]")

plt.title(f"Visualization of Spherical Harmonics - degree: {l} order: {m}")
plt.show()

Instead we can rely on the 'ylm_plot' utility in the visualization module

In [22]:

# basic plotting of SH coeff on plane - Zonal
ylmc_40, ylms_40 = ylm(l=4, m=0)

ylm_plot(l=4, m=0)

In [23]:

# basic plotting of SH coeff on plane - Teseral
l = 5
m = 5
ylmc_55, ylms_55 = ylm(l=l, m=m)

ylm_plot(l, m)

In [24]:

# basic plotting of SH coeff on plane - Sectoral
l = 11
m = 8
ylmc_11_8, ylms_11_8 = ylm(l=l, m=m)

ylm_plot(l, m)

In [25]:

# Summing up (Superposition) SH Coeff C00 + C10 + S0-1 + C01
# basic plotting of SH coeff on plane - Sectoral

# for explaination purpose
C00 = 1.0
C10 = 10.0
C55 = 15.0
C11_8 = 5.0

ylmc = C00*ylmc_00 + C10*ylmc_40 + C55*ylmc_55 + C11_8*ylmc_11_8

In [26]:

# Summing SH Coeff C00 + C10 + S0-1 + C01
# basic plotting of SH coeff on plane - Sectoral
l = 2
m = 2

ylmc = C00*ylmc_00 + C10*ylmc_40 + C55*ylmc_55 + C11_8*ylmc_11_8

fig = plt.figure(figsize=(15, 7.5))
ax = plt.axes(projection = ccrs.Orthographic(central_longitude=74.0, central_latitude=0.0, globe=None))

lons = np.linspace(-180, 180, 73)
lats = np.linspace(-90, 90, 37)

# max_level = 6
# min_level = -4.5
# step_level = 0.1
# levels = np.arange(min_level, max_level + step_level, step_level)
x, y = np.meshgrid(lons, lats)

if m >=0 :
    # plt.contourf(x, y, ylmc[:, 0, :], cmap='RdYlBu_r', )
    img_extent = (-180, 180, -90, 90)

    # plot the data
    im = ax.imshow(ylmc[:, 0, :], origin='upper', extent=img_extent, transform=ccrs.PlateCarree(),)
else:
    plt.contourf(x, y, ylms[:, 0, :], cmap='RdYlBu_r')


# setting gridlines
gl = ax.gridlines(crs = ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, color='gray', alpha=0.9, linestyle='--')
# remove top x label
gl.top_labels = False
# change x label styles - font size ad colour
gl.xlabel_style = {'size':12,}
# left and right labels
gl.left_labels = True
gl.right_labels = False
# coastlines
ax.coastlines()

plt.colorbar(im, orientation='vertical', shrink=0.84, pad=0.02,label=f"[...]")

plt.title(f"Surface Spherical Harmonics - Summation - C00 + C40 + C55 + C11_8")
plt.show()

Thank you for reading