Loading Spherical Data

Support in yt for non-Cartesian geometries is partial and is still being extended, but here is an example of how to load spherical data from a regularly-spaced grid. For irregularly spaced grids, a similar setup can be used, either by using supplying the cell_widths parameter to load_uniform_grid (see Stretched Grid Data for more details), or by using the load_hexahedral_mesh method.

Note that in yt, “spherical” means that it is ordered \(r\), \(\theta\), \(\phi\), where \(\theta\) is the declination from the azimuth (running from \(0\) to \(\pi\) and \(\phi\) is the angle around the zenith (running from \(0\) to \(2\pi\)).

We first start out by loading yt.

import numpy as np

import yt

Now, we create a few derived fields. The first three are just straight translations of the Cartesian coordinates, so that we can see where we are located in the data, and understand what we’re seeing. The final one is just a fun field that is some combination of the three coordinates, and will vary in all dimensions.

@yt.derived_field(name="sphx", units="cm", take_log=False, sampling_type="cell")
def sphx(field, data):
    return np.cos(data["phi"]) * np.sin(data["theta"]) * data["r"]


@yt.derived_field(name="sphy", units="cm", take_log=False, sampling_type="cell")
def sphy(field, data):
    return np.sin(data["phi"]) * np.sin(data["theta"]) * data["r"]


@yt.derived_field(name="sphz", units="cm", take_log=False, sampling_type="cell")
def sphz(field, data):
    return np.cos(data["theta"]) * data["r"]


@yt.derived_field(name="funfield", units="cm", take_log=False, sampling_type="cell")
def funfield(field, data):
    return (np.sin(data["phi"]) ** 2 + np.cos(data["theta"]) ** 2) * (
        1.0 * data["r"].uq + data["r"]
    )

Loading Data

Now we can actually load our data. We use the load_uniform_grid function here. Normally, the first argument would be a dictionary of field data, where the keys were the field names and the values the field data arrays. Here, we’re just going to look at derived fields, so we supply an empty one.

The next few arguments are the number of dimensions, the bounds, and we then specify the geometry as spherical.

ds = yt.load_uniform_grid(
    {},
    [128, 128, 128],
    bbox=np.array([[0.0, 1.0], [0.0, np.pi], [0.0, 2 * np.pi]]),
    geometry="spherical",
)

Looking at Data

Now we can take slices. The first thing we will try is making a slice of data along the “phi” axis, here \(\pi/2\), which will be along the y axis in the positive direction. We use the .slice attribute, which creates a slice, and then we convert this into a plot window. Note that here 2 is used to indicate the third axis (0-indexed) which for spherical data is \(\phi\).

This is the manual way of creating a plot – below, we’ll use the standard, automatic ways. Note that the coordinates run from \(-r\) to \(r\) along the \(z\) axis and from \(0\) to \(r\) along the \(R\) axis. We use the capital \(R\) to indicate that it’s the \(R\) along the \(x-y\) plane.

s = ds.slice(2, np.pi / 2)
p = s.to_pw("funfield", origin="native")
p.set_zlim("all", 0.0, 4.0)
p.show()

We can also slice along \(r\). For now, this creates a regular grid with incorrect units for phi and theta. We are currently exploring two other options – a simple aitoff projection, and fixing it to use the correct units as-is.

s = yt.SlicePlot(ds, "r", "funfield")
s.set_zlim("all", 0.0, 4.0)
s.show()

We can also slice at constant \(\theta\). But, this is a weird thing! We’re slicing at a constant declination from the azimuth. What this means is that when thought of in a Cartesian domain, this slice is actually a cone. The axes have been labeled appropriately, to indicate that these are not exactly the \(x\) and \(y\) axes, but instead differ by a factor of \(\sin(\theta))\).

s = yt.SlicePlot(ds, "theta", "funfield")
s.set_zlim("all", 0.0, 4.0)
s.show()

We’ve seen lots of the funfield plots, but we can also look at the Cartesian axes. This next plot plots the Cartesian \(x\), \(y\) and \(z\) values on a \(\theta\) slice. Because we’re not supplying an argument to the center parameter, yt will place it at the center of the \(\theta\) axis, which will be at \(\pi/2\), where it will be aligned with the \(x-y\) plane. The slight change in sphz results from the cells themselves migrating, and plotting the center of those cells.

s = yt.SlicePlot(ds, "theta", ["sphx", "sphy", "sphz"])
s.show()

We can do the same with the \(\phi\) axis.

s = yt.SlicePlot(ds, "phi", ["sphx", "sphy", "sphz"])
s.show()