More about localized spectral analysis on the sphere

A localised spectrum estimate is the spectrum of a global function multiplied by a localization window. In a 2D Cartesian coordinate system, this multiplication becomes a convolution in the Fourier domain. As a result of the convolution, a localized spectral estimate for a given wavelength depends upon a range of coefficients of both the global function and localization window. One says that the window couples the global Fourier coefficients to the localized coefficients.

Windowing on the sphere leads to similar effects: Multiplication of a global function by a window leads to coupling of the spherical harmonic coefficients of the globally defined function to the localized spectrum. In the case of an isotropic and stationary stochastic function, a coupling matrix describes the relation between the global power and the expectation of the localized power, and at high degrees, this matrix operation also resembles a convolution.

In [1]:
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pyshtools as pysh
In [2]:
%config InlineBackend.figure_format = 'retina'  # if you are not using a retina display, comment this line

Generate a random model using the SHCoeffs class interface

We first generate a random model using the SHCoeffs class interface. This class simplifies coefficient generation and provides a simple access to most SHTOOLS functions. We start by defining a global power spectrum that follows a power-law-like decay:

In [3]:
lmax = 100
a = 30  # scale length
degrees = np.arange(lmax+1, dtype=np.float)
power = 1. / (1. + (degrees / a) ** 2) ** 1.5

Next, we generate random coefficients from this input power spectrum, plot the power spectrum of the realization, and expand the coeffificients on a grid:

In [4]:
coeffs_global = pysh.SHCoeffs.from_random(power, seed=12345)

power_global = coeffs_global.spectrum()
fig, ax = coeffs_global.plot_spectrum(unit='per_dlogl', show=False)    # show=False is used to avoid a warning when plotting in inline mode

grid_global = coeffs_global.expand(grid='DH2')
fig2, ax2 = grid_global.plot(colorbar='right', show=False)

Generate a box window function

We next generate a window function that picks a few local regions from the globally defined model. In this example, the window contains sharp boundaries, so the spherical harmonic bandwidth of the function is infinite. In general, this could severely bias the localized power spectrum away from its global value, especially when the power spectrum follows a power law. To combat this negative characteristic, we will show in a different tutorial how to construct windows with a specified spherical harmonic bandwith that are optimally concentrated within the provided region.

This example makes use of the class SHGrid, which is the counterpart to SHCoeffs. Let's start by making an arbitrary mask:

In [5]:
latgrid, longrid = np.meshgrid(grid_global.lats(), grid_global.lons(), indexing='ij')

window = (-40 < latgrid) & (latgrid < -30) & (10 < longrid) & (longrid < 30)
window += (0 < latgrid) & (latgrid < 30) & (60 < longrid) & (longrid < 80)
window += (-70 < latgrid) & (latgrid < -40) & (130 < longrid) & (longrid < 150)
window += (20 < latgrid) & (latgrid < 40) & (125 < longrid) & (longrid < 145)
window += (10 < latgrid) & (latgrid < 30) & (220 < longrid) & (longrid < 250)

Next, generate an SHGrid instance from the input array, plot the grid, expand it in spherical harmonics, and calculate and plot the power spectrum. Note that the first element of the grid corresponds to 0 degrees east and 90 degrees north.

In [6]:
grid_window = pysh.SHGrid.from_array(window.astype(np.float64))
fig, ax = grid_window.plot(colorbar='right', show=False)

coeffs_window = grid_window.expand()
fig2, ax2 = coeffs_window.plot_spectrum(unit='per_dlogl', show=False)

power_window = coeffs_window.spectrum()

Multiply the random model by the window

Multiplication of the global model with the window function localizes the data, and its expansion in spherical harmonics gives a local power spectrum estimate. The interaction of the window function with the random model distorts the spectrum. In particular, the ouput spectrum at degree l is influenced by the input spherical harmonic degrees from l - Lwin to l + Lwin, where Lwin in the spherical harmonic bandwidth of the windowing function.

Let's multiply the data by the window, expand the result in spherical harmonics, and plot the resulting localized power spectrum:

In [7]:
grid_local = grid_global * grid_window
grid_local.plot(colorbar='right', show=False)

coeffs_local = grid_local.expand()
coeffs_local.plot_spectrum(unit='per_dlogl', show=False)

power_local = coeffs_local.spectrum()