In [1]:
%matplotlib inline
from __future__ import print_function # only necessary if using Python 2.x
import matplotlib.pyplot as plt
import numpy as np
In this example, we will use the topography of Earth that is provided in the example directory,
In [2]:
infile = '../ExampleDataFiles/topo.dat.gz'
topo = np.loadtxt(infile)
# plot with matplotlib
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.imshow(topo, extent=(-180, 180, -90, 90), cmap='viridis')
ax.set(xlabel='longitude', ylabel='latitude');
Using this data set, we will first demonstrate how to calculate the power spectrum, and second, how to perform simple filtering operations.
The power spectrum of a function describes how the variance of the function is distributed as a function of spherical harmonic degree. Concentration of power (or energy) in spherical harmonics with a particular wavelength can give hints about the origin of the feature. First, let's import the required functions, and expand the grid into spherical harmonics:
In [3]:
from pyshtools.expand import SHExpandDH;
from pyshtools.spectralanalysis import spectrum
coeffs = SHExpandDH(topo, sampling=2)
nl = coeffs.shape[1]
ls = np.arange(nl)
The "power spectrum" can be calculated using different conventions. The default in pyshtools is to calculate the total power of all angular orders as a function of spherical harmonic degee, which corresponds to the option unit='per_l'. The power per degree l is somewhat equivalent to the power at a given wavenumber magnitude |k| in 2D Fourier analyses.
In [4]:
power_per_l = spectrum(coeffs)
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.plot(ls, power_per_l)
ax.set_yscale('log')
ax.set_xscale('log')
ax.grid()
The average power per coefficient as a function of spherical harmonic degree can be calculated by setting the parameter unit equal to 'per_lm'. This is simply the power per degree divided by (2l+1), and is analogous to the power per wavenumber kx and ky in 2D Fourier analyses.
In [5]:
power_per_lm = spectrum(coeffs, unit='per_lm')
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.plot(ls, power_per_lm)
ax.set_yscale('log')
ax.set_xscale('log')
ax.grid()
Finally, the contribution to the total power from all angular orders over an infinitessimal logarithmic degree band dlog_a(l) can be calculated by setting unit='per_dlogl'. In this case, the contrubution in the band
dlog_a(l) is spectrum(l, 'per_dlogl')*dlog_a(l), where a is the base, and where spectrum(l, 'per_dlogl) is equal to spectrum(l, 'per_l')*l*log(a). This power spectrum is useful to analyse dominant scales.
In [6]:
power_per_dlogl = spectrum(coeffs, unit='per_dlogl', base=2.)
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.plot(ls, power_per_dlogl)
ax.set_yscale('log', basey=2)
ax.set_xscale('log', basex=2)
ax.grid()
A global dataset can be filtered isostropically by multiplying the spherical harmonic coefficients by a degree-dependent function. We demonstrate this by setting the coefficients greater or equal to degree 8 equal to zero.
In [7]:
from pyshtools.expand import MakeGridDH;
coeffs_filtered = coeffs.copy()
lmax = 8
coeffs_filtered[:, lmax:, :] = 0.
topo_filtered = MakeGridDH(coeffs_filtered, sampling=2)
Next, we bandpass filter the data to retain only spherical harmonic coefficients between 8 <= l < 20.
In [8]:
coeffs_filtered2 = coeffs.copy()
lmin, lmax = 8, 20
coeffs_filtered2[:, :lmin, :] = 0.
coeffs_filtered2[:, lmax:, :] = 0.
topo_filtered2 = MakeGridDH(coeffs_filtered2, sampling=2)
Finally, let's plot the two filtered data sets:
In [9]:
fig, (row1, row2) = plt.subplots(2, 1, figsize=(10, 10))
row1.imshow(topo_filtered, extent=(-180, 180, -90, 90), cmap='viridis')
row1.set(xlabel='longitude', ylabel='latitude', title='l=0 - 7');
row2.imshow(topo_filtered2, extent=(-180, 180, -90, 90), cmap='viridis')
row2.set(xlabel='longitude', ylabel='latitude', title='l=8 - 19');