example_spherical_harmonics

In [1]:

import terratools
from terratools.terra_model import read_netcdf
from terratools.example_data import example_terra_model
import matplotlib.pyplot as plt

import terratools from terratools.terra_model import read_netcdf from terratools.example_data import example_terra_model import matplotlib.pyplot as plt

First we download and read in the example mantle convection model:

In [2]:

# Download and cache model
path = example_terra_model()

# read in the model
model = read_netcdf([path])

# Download and cache model path = example_terra_model() # read in the model model = read_netcdf([path])

We can convert the composition histograms of harzburgite, lherzolite and basalt into a bulk composition scalar field using model.calc_bulk_composition()

In [3]:

model.calc_bulk_composition()

model.calc_bulk_composition()

We can perform spherical harmonic decompositions of scalar field. In terratools we make use of healpy (python wrapper for healpix) to do this. For example we would do:

In [4]:

model.calc_spherical_harmonics("c")

model.calc_spherical_harmonics("c")

to calculate the spherical harmonic coefficients and power spectra of the bulk composition field. We can also include optional arguments nside to change the resolution of the healpix grid (default=2**6), lmax to set the maximum spherical harmonic degree to calculate (default=16), and savemap toggle saving of a healpix map (default=False).

In [5]:

model.calc_spherical_harmonics("t", nside=2**6, lmax=20, savemap=True)

model.calc_spherical_harmonics("t", nside=2**6, lmax=20, savemap=True)

There is a 'getter' function for accessing the spherical harmonic coefficients etc

In [6]:

sph_temp = model.get_spherical_harmonics("t")
sph_comp = model.get_spherical_harmonics("c")

sph_temp = model.get_spherical_harmonics("t") sph_comp = model.get_spherical_harmonics("c")

We can see, for example, that changing lmax when calculating the spherical harmonics for the temperature field changes the length of 'power_per_l' at each radial layer.

In [7]:

print(len(sph_comp[0]["power_per_l"]))
print(len(sph_temp[0]["power_per_l"]))

print(len(sph_comp[0]["power_per_l"])) print(len(sph_temp[0]["power_per_l"]))

17
21

We can plot the power spectra of the for fields as a function of depth using:

In [8]:

fig, ax = model.plot_spectral_heterogeneity("c", lyrmin=0, lyrmax=-1, show=False)
plt.show()

fig, ax = model.plot_spectral_heterogeneity("c", lyrmin=0, lyrmax=-1, show=False) plt.show()

We can also plot depth slices of the spectrally filtered fields using:

In [9]:

model.plot_hp_map("t", index=1, nside=2**5, show=False)

plt.show()

model.plot_hp_map("t", index=1, nside=2**5, show=False) plt.show()

Note that either index or radii must be passed in determine which layer to plot.