*pyshtools* is a Python package that provides functions for working with spherical harmonics. These routines can be used to perform spherical harmonic transforms and reconstructions, rotations of data expressed in spherical harmonics, and multitaper spectral analyses on the sphere. The base functions are fast routines written in Fortran 95 from the SHTOOLS package, and *pyshtools* provides easy access to these routines by use of Python-wrapper functions and a compact class interface.

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

In [1]:

```
from __future__ import print_function # only necessary if using Python 2.x
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pyshtools
```

The *pyshtools* package contains several classes, including `SHCoeffs`

, `SHGrid`

, and `SHWindow`

, a subpackage `shtools`

that contains the Python-wrapped Fortran functions, and a series of subpackages that contain both native python routines and refernces to routines in `shtools`

(`constant`

, `legendre`

, `expand`

, `shio`

, `spectralanalysis`

, `rotate`

, `gravmag`

, and `util`

). The classes `SHCoeffs`

, `SHGrid`

, and `SHWindow`

bundle the vast majority of *pyshtools* functionality, and in this notebook we will demonstrate how to use the classes related to grids and coefficients.

*pyshtools* can also define *matplotlib* configuration variables for creating publication quality graphics. These can be set using the function `figstyle`

and by specifing the relative width of the image with respect to a journal page. If necessary, one can also specify optional parameters for the resolution of the screen, the maximum width of the image, and the default image aspect ratio.

In [2]:

```
pyshtools.utils.figstyle(rel_width=0.75)
%config InlineBackend.figure_format = 'retina' # if you are not using a retina display, comment this line!
```

Let's start by first creating a power spectrum (i.e., the total power per degree) 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 [3]:

```
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 expected power is given by the spectrum we just created:

In [4]:

```
clm = pyshtools.SHCoeffs.from_random(power, seed=12345)
```

This creates a new class instance of `SHCoeffs`

that contains several attributes and methods. For reproducibility, the optional parameter `seed`

specifies the seed of the *numpy* random number generator. By default, *pyshtools* assumes that the coefficients are real and that they are normalized using the `'4pi'`

convention exlcuding the Condon-Shortley phase factor. This is the standard normalization in Geodesy and many fields of geophysics and spectral analysis. Other normalizations can be specified explicitly by specifying the optional parameter `normalization`

, which can be `'4pi'`

, `'ortho'`

, `'schmidt'`

, or `'unnorm'`

. The Condon-Shortley phase can be included by setting the optional parameter `csphase`

to `-1`

, and if you wanted complex coefficients, you could set `kind='complex'`

.

`from_random()`

is just one way to create a set of spherical harmonic coefficients. The other constructor methods are `from_file()`

to read the coefficients from an `shtools`

or `npy`

formatted file, `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 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 `spectrum()`

method.

In [5]:

```
fig, ax = clm.plot_spectrum(show=False) # show=False is used to avoid a warning when plotting in inline mode
```

By default, this method plots the power in units of power per spherical harmonic degree (`per_l`

), but you can also choose from the average power per coefficient as a function of spherical harmonic degree (`per_lm`

) and the power per log bandwith (`per_dlogl`

). You can visualize the power asociated with each spherical-harmonic coefficient using the `plot_spectrum2d()`

method. In this example, we use the optional parameter `vrange`

to scale the colormap from 0.1 to 1.e-7 times the maximum value.

In [6]:

```
fig, ax = clm.plot_spectrum2d(vrange=(1.e-7,0.1), show=False)
```

The spherical harmonic coefficients of a given `SHCoeffs`

class instance can be converted easily to a different normalization using the `convert()`

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 just return the first 50 degrees of the function:

In [7]:

```
clm_ortho = clm.convert(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 [8]:

```
clm_ortho.info()
```

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

In [9]:

```
fig, ax = clm.plot_spectrum(legend='Orthonormalized', show=False)
clm_ortho.plot_spectrum(ax=ax, linestyle='dashed', legend='4$\pi$ normalized')
ax.plot(clm.degrees(), power, '-k')
limits = ax.set_xlim(0, 100)
```

The method `degrees()`

was here used to return a *numpy* list of the degrees from 0 up to the maximum value of the coefficients. As you can see, the power spectrum of the random realization follows closely the input expectation spectrum. Furthermore, the power spectra are seen to be independent of the chosen normalization.

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`

, and then plot the data using the built in method `plot()`

.

In [10]:

```
grid = clm.expand()
fig, ax = grid.plot(show=False)
```

*pyshtools* makes use of three different grid types. 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 polynomial of the maximum degree of the coefficients. 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 [11]:

```
grid_glq = clm.expand(grid='GLQ')
grid_glq.plot(show=False)
lats = grid_glq.lats()
lons = grid_glq.lons()
```

If you ever forget how your grid was constructed, the `info()`

method provides you with everything you need to know:

In [12]:

```
grid_glq.info()
```

In [13]:

```
grid.info()
```

Sometimes you need to set individual spherical harmonic coefficients to a specified value. For example, for planetary topography, the degree 2 order 0 term can be atypically 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. Though we will modify only a single coefficient here, you could also provide input arrays for `l`

, `m`

, and `value`

. When using the `set_coeffs()`

method, you set the cosine terms using positive values for the orders, and the sine terms using negative orders. We expand the modified coefficients on a grid, and then plot the data with a vertical color bar:

In [14]:

```
clm.set_coeffs(values=0., ls=2, ms=0)
grid_dh2 = clm.expand(grid='DH2')
fig, ax = grid_dh2.plot(colorbar=True, cb_orientation='vertical', cb_label='My data', show=False)
grid_dh2.info()
```