python: example grid plotting script for Meraxes

read_grids.py

#!/usr/bin/env python

"""Simple test script showing how to read / plot slices from Meraxes grids."""

import numpy as np
import matplotlib.pyplot as plt
from ssimpl import meraxes

# path to meraxes file and snapshot to read
fname = "/path/to/meraxes.hdf5"
snapshot = 51

# list the available grids
print("Available grids in {:}:".format(fname))
print(meraxes.list_grids(fname, snapshot))

# read in the delta_T grid
delta_T = meraxes.read_grid(fname, snapshot, "delta_T")

# Create a reusable slice (this is just a way to avoid having to re-write the
# slice everywhere if you change it).
s = np.s_[:, :, delta_T.shape[0]/2]  # centre of box in x,y plane

# Set up the figure
fig, ax = plt.subplots(1, 1)
ax.set_xlabel("x")
ax.set_ylabel("y")

# plot the slice
# N.B. the `.T` which takes the transpose and the `origin="lower"` which
# ensures [0,0] is at the bottom left.
ax.imshow(delta_T[s].T, origin="lower", cmap=plt.cm.gist_earth)

# save the figure
plt.tight_layout()
plt.savefig("test_grid.png")