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.
%matplotlib inline import matplotlib.pyplot as plt import numpy as np import pyshtools as pysh
pysh.utils.figstyle(rel_width=0.75) %config InlineBackend.figure_format = 'retina' # if you are not using a retina display, comment this line
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:
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:
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)
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:
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.
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()
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:
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()