Introduction to pyshtools (Part 1): Grids and Coeffs

pyshtools is a Python module that provides functions for working with spherical harmonics. The base functions are time-tested routines from the Fortran 95 SHTOOLS package. pyshtools provides easy access to these routines by use of Python-wrapper functions and a few simplified classes for spherical harmonic coefficients, grids, and localization windows.

To get started, import the standard matplotlib library for graphics, numpy for mathematical extensions to Python, and pyshtools:

In [1]:
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pyshtools

The pyshtools namespace is composed of three classes (SHCoeffs, SHGrid, and SHWindow) a submodule shtools that contains all of the Python-wrapped Fortran functions, and a series of submodules that contain subsets of these functions (constant, legendre, expand, shio, spectralanalysis, localizedspectralanalysis, rotate, gravmag, and other). The vast majority of pyshtools functionality is implemented in the classes SHCoeffs and SHGrid, which we will demonstrate here.

Let's start by first creating a power spectrum that follows a power law with exponent -2, up to and including degree 100. To avoid a division by zero, we will set the degree 0 term to zero:

In [2]:
degrees = np.arange(101, dtype=float)
degrees[0] = np.inf
power = degrees**(-2)

Next, we will create a random realization of spherical harmonic coefficients whose expectation of the power spectrum is given by the spectrum we just created:

In [3]:
clm = pyshtools.SHCoeffs.from_random(power)

This creates a new class instance of SHCoeffs that contains several attributes and method functions that will be explored below. By default, pyshtools assumes that the coefficients are real and that they are normalized using the '4pi' convention that exlcudes the Condon-Shortley phase factor. This is the standard normalization in Geodesy and many fields of geophysics and spectral analysis. Nevertheless, another normalization can be specified explicitly by providing the optional parameter normalization, which can be either '4pi', 'ortho', or 'schmidt'. The Condon-Shortley phase can be included by setting csphase=-1, and if you wanted complex coefficients, you could set kind='complex'. In general you don't need to worry about these parameters, because most output is provided in a normalization independent manner.

This is just one way to create a set of spherical harmonic coefficients. The other constructor methods are from_file() to read the coefficients from a file using the SHRead routine, from_zeros() if you just want all the coefficients to be set to zero, and from_array() if you already have a numpy array of the coefficients.

Next, let's calculate the power spectrum (i.e., the power per spherical harmonic degree) and plot it. pyshtools provides a built in plotting function to do this, and as we will see below, the power spectrum can also be returned as a numpy array using the get_powerperdegree() method.

In [4]:
fig, ax = clm.plot_powerperdegree()

The coefficients can be converted easily to a different normalization using the return_coeffs() method. Here, we will convert the coeffients to the orthonormalized convention using the Condon-Shortley phase, which is common in many fields of physics and seismology. Also, let's pretend that we only need the first 50 degrees of the function:

In [5]:
clm_ortho = clm.return_coeffs(normalization='ortho', csphase=-1, lmax=50)

If you ever forget how your data are normalized, you can use the built in info() method to remind you:

In [6]:
kind = 'real'
normalization = 'ortho'
csphase = -1
lmax = 50

We will next calculate the power spectrum of our two functions and plot them along with our theoretical input spectrum:

In [7]:
fix, ax = plt.subplots(1,1)
ax.plot(clm.get_degrees(), power, '-k', clm.get_degrees(), clm.get_powerperdegree(), '-r', clm_ortho.get_degrees(), clm_ortho.get_powerperdegree(), '-b')
ax.set(yscale='log', xlabel='degree', ylabel='power')
 <matplotlib.text.Text at 0x10b2626a0>,
 <matplotlib.text.Text at 0x10b10ceb8>]

As you can see, the power spectrum of the random realization follows closely the input spectrum. Furthermore, given that the power spectrum is independent of the employed normalization, the power spectrum using the '4pi' and 'ortho' normalizations are identical. You will also note that we used the method get_degrees() to return a numpy list of the degrees from 0 up to the maximum value.

Next, let's expand the data onto a grid and plot it. We first use the expand() method, which returns a new instance of the class SHGrid. Then we plot the data using the built in method plot_rawdata():

In [8]:
clm_grid = clm.expand()

When expanding data onto a grid, pyshtools provides three options. The default is to use grid='DH', which is an equally-sampled grid in latitude and longitude that conforms to Driscoll and Healy's (1994) sampling theorem. If you would like a grid that is oversampled in longitude by a factor or two, such that the number of longitude bands is twice that as the number of latitude bands, use grid='DH2' instead.

The third grid is constructed explicitly for use with Gauss-Legendre quadrature integration techniques. This grid contains about half as many latitude bands as an equivalent DH grid, but the latitudes are unequally sampled at the zeros of the Legendre function of degree lmax. The following commands show how to expand the spherical harmonic coefficents onto a GLQ grid, plot it, and output lists that contain the latitudes and longitudes (in degrees) for each row and column of the grid:

In [9]:
clm_glq_grid= clm.expand(grid='GLQ')

lats = clm_glq_grid.get_lats()
lons = clm_glq_grid.get_lons()

Once again, if you ever forget how your grid was constructed, the info() method provides you with everything you need to know:

In [10]:
kind = 'real'
grid = 'GLQ'
nlat = 101
nlon = 201
lmax = 100
In [11]:
kind = 'real'
grid = 'DH'
sampling = 1
nlat = 202
nlon = 202
lmax = 100

Sometimes you need to set individual spherical harmonic coefficients to a specified value. For example, for planetary topography, the degree 2 order 0 term is usually large as a result of the planet's rotation. Let's set this term equal to zero and replot it using an over-sampled DH grid. When using the set_coeffs() method, you set the cosine terms using positive values for the orders, and the sine terms using negative orders:

In [12]:
value = 0.

clm.set_coeffs(value, l, m)
clm_grid_dh2 = clm.expand(grid='DH2')
kind = 'real'
grid = 'DH'
sampling = 2
nlat = 202
nlon = 404
lmax = 100

Of course, you could also have input arrays of equal length for l, m, and value. If you want to extract the grid as a numpy array, this can be done using the get_grid() method. If instead you want to extract the coefficients as a numpy array, this can be done using the get_coeffs() method. Here, we will just return the first few degrees of the spherical harmonic coefficients and verify that the coefficient was indeed set to zero:

In [13]:
coeffs = clm.get_coeffs(lmax=4)
print(coeffs, end='\n\n')
print('c20 = {:f}'.format(coeffs[0,2,0]))
[[[ -0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
  [ -7.03293569e-01   8.88501971e-02   0.00000000e+00   0.00000000e+00
  [  0.00000000e+00  -2.92033184e-01  -5.79078284e-01   0.00000000e+00
  [ -1.76405372e-02  -2.45524713e-01  -4.34332126e-02   4.22255366e-02
  [  1.01620014e-02  -2.85416844e-02   2.92043770e-02   1.99774057e-02

 [[  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
  [  0.00000000e+00  -1.74070704e-01   0.00000000e+00   0.00000000e+00
  [  0.00000000e+00  -1.62044479e-01  -7.14799452e-02   0.00000000e+00
  [  0.00000000e+00  -1.18153463e-01  -4.51717452e-04   1.55294960e-01
  [  0.00000000e+00   2.43668012e-02   5.47790343e-02  -1.09160268e-01

c20 = 0.000000

It is also easy to rotate either the physical body or the underlying coordinate frame of the function expressed in spherical harmonics. This is accomplished using three Euler angles, alpha, beta and gamma. There are many different notations for these angles, so please read the documentation before doing this blindly! Here, we will rotate the point at the north pole to 60 degrees north latitude and 180 degrees east longitude:

In [14]:
clat = 60.
clon = 180.

alpha = 0.
beta = -(90.-clat)
gamma = -clon

clm_rotated = clm.rotate(alpha, beta, gamma, degrees=True)
grid_rotated = clm_rotated.expand()

By default, this routine expects the angles to be in degrees, which we specified redundantly by the optional parameter. If the angles were in radians, you would instead set this to False.

Finally, pyshtools provides many physical constants that are useful when working with planetary data. As an example, here is the GM of Earth, along with information on its source:

In [15]:
print(pyshtools.constant.gm_earth, end='\n\n')


Gravitational constant times the mass of the Earth.


398.6004415 10^12 m3 s-2


N. K. Pavlis, S. A. Holmes, S. C. Kenyon, and J. K. Factor (2012). The
development and evaluation of the Earth Gravitational Model 2008 (EGM2008). J.
Geophys. Res. 117, B04406, doi:10.1029/2011JB008916.

All of the operations in this introduction could have also been performed using the raw Python-wrapped Fortran 95 routines. In general, when using these routines, one needs to pay attention to the specific properties of the grids and normalizations of the coefficients. The SHCoeffs and SHGrid classes simpify the access to these routines given that all the metadata is stored directly in the class attributes. Nevertheless, in some cases the wrapped SHTOOLS routines might be preferable to use as they might be more computationally efficient, or perhaps use less memory. Not all pyshtools routines are accessible directly from the class interfaces, and in those cases, it is often necessary to input the raw coefficients and gridded data, which can be obtained from the methods get_coeffs() and get_grid().