Analyzing FITS Radio Cubes

%matplotlib inline
import yt

This notebook demonstrates some of the capabilities of yt on some FITS “position-position-spectrum” cubes of radio data.

Note that it depends on some external dependencies, including astropy and regions.

M33 VLA Image

The dataset "m33_hi.fits" has NaNs in it, so we’ll mask them out by setting nan_mask = 0:

ds = yt.load("radio_fits/m33_hi.fits", nan_mask=0.0)

First, we’ll take a slice of the data along the z-axis, which is the velocity axis of the FITS cube:

slc = yt.SlicePlot(ds, "z", ("fits", "intensity"), origin="native")
slc.show()

The x and y axes are in units of the image pixel. When making plots of FITS data, to see the image coordinates as they are in the file, it is helpful to set the keyword origin = "native". If you want to see the celestial coordinates along the axes, you can import the PlotWindowWCS class and feed it the SlicePlot. For this to work, a version of AstroPy >= 1.3 needs to be installed.

from yt.frontends.fits.misc import PlotWindowWCS

PlotWindowWCS(slc)

Generally, it is best to get the plot in the shape you want it before feeding it to PlotWindowWCS. Once it looks the way you want, you can save it just like a normal PlotWindow plot:

slc.save()

['m33_hi.fits_Slice_z_intensity.png']

We can also take slices of this dataset at a few different values along the “z” axis (corresponding to the velocity), so let’s try a few. To pick specific velocity values for slices, we will need to use the dataset’s spec2pixel method to determine which pixels to slice on:

import yt.units as u

new_center = ds.domain_center
new_center[2] = ds.spec2pixel(-250000.0 * u.m / u.s)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[6], line 4
      1 import yt.units as u
      3 new_center = ds.domain_center
----> 4 new_center[2] = ds.spec2pixel(-250000.0 * u.m / u.s)

File /tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1760, in unyt_array.__setitem__(self, item, value)
   1758     if value.units != self.units and value.units != NULL_UNIT:
   1759         value = value.to(self.units)
-> 1760 super().__setitem__(item, value)

ValueError: assignment destination is read-only

Now we can use this new center to create a new slice:

slc = yt.SlicePlot(ds, "z", ("fits", "intensity"), center=new_center, origin="native")
slc.show()

We can do this a few more times for different values of the velocity:

new_center[2] = ds.spec2pixel(-100000.0 * u.m / u.s)
slc = yt.SlicePlot(ds, "z", ("fits", "intensity"), center=new_center, origin="native")
slc.show()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 1
----> 1 new_center[2] = ds.spec2pixel(-100000.0 * u.m / u.s)
      2 slc = yt.SlicePlot(ds, "z", ("fits", "intensity"), center=new_center, origin="native")
      3 slc.show()

File /tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1760, in unyt_array.__setitem__(self, item, value)
   1758     if value.units != self.units and value.units != NULL_UNIT:
   1759         value = value.to(self.units)
-> 1760 super().__setitem__(item, value)

ValueError: assignment destination is read-only

new_center[2] = ds.spec2pixel(-150000.0 * u.m / u.s)
slc = yt.SlicePlot(ds, "z", ("fits", "intensity"), center=new_center, origin="native")
slc.show()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[9], line 1
----> 1 new_center[2] = ds.spec2pixel(-150000.0 * u.m / u.s)
      2 slc = yt.SlicePlot(ds, "z", ("fits", "intensity"), center=new_center, origin="native")
      3 slc.show()

File /tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1760, in unyt_array.__setitem__(self, item, value)
   1758     if value.units != self.units and value.units != NULL_UNIT:
   1759         value = value.to(self.units)
-> 1760 super().__setitem__(item, value)

ValueError: assignment destination is read-only

These slices demonstrate the intensity of the radio emission at different line-of-sight velocities.

We can also make a projection of all the emission along the line of sight:

prj = yt.ProjectionPlot(ds, "z", ("fits", "intensity"), origin="native")
prj.show()

We can also look at the slices perpendicular to the other axes, which will show us the structure along the velocity axis:

slc = yt.SlicePlot(ds, "x", ("fits", "intensity"), origin="native", window_size=(8, 8))
slc.show()

slc = yt.SlicePlot(ds, "y", ("fits", "intensity"), origin="native", window_size=(8, 8))
slc.show()

In these cases, we needed to explicitly declare a square window_size to get a figure that looks good.

\(^{13}\)CO GRS Data

This next example uses one of the cubes from the Boston University Galactic Ring Survey.

ds = yt.load("radio_fits/grs-50-cube.fits", nan_mask=0.0)

We can use the quantities methods to determine derived quantities of the dataset. For example, we could find the maximum and minimum temperature:

dd = ds.all_data()  # A region containing the entire dataset
extrema = dd.quantities.extrema(("fits", "temperature"))
print(extrema)

[-41.05050659  79.63768005] K

We can compute the average temperature along the “velocity” axis for all positions by making a ProjectionPlot:

prj = yt.ProjectionPlot(
    ds, "z", ("fits", "temperature"), origin="native", weight_field=("index", "ones")
)  # "ones" weights each cell by 1
prj.set_zlim(("fits", "temperature"), zmin=(1e-3, "K"))
prj.set_log(("fits", "temperature"), True)
prj.show()

We can also make a histogram of the temperature field of this region:

pplot = yt.ProfilePlot(
    dd, ("fits", "temperature"), [("index", "ones")], weight_field=None, n_bins=128
)
pplot.show()

We can see from this histogram and our calculation of the dataset’s extrema that there is a lot of noise. Suppose we wanted to make a projection, but instead make it only of the cells which had a positive temperature value. We can do this by doing a “field cut” on the data:

fc = dd.cut_region(['obj["fits", "temperature"] > 0'])

Now let’s check the extents of this region:

print(fc.quantities.extrema(("fits", "temperature")))

[7.45058060e-09 7.96376801e+01] K

Looks like we were successful in filtering out the negative temperatures. To compute the average temperature of this new region:

fc.quantities.weighted_average_quantity(("fits", "temperature"), ("index", "ones"))

unyt_quantity(0.16913964, 'K')

Now, let’s make a projection of the dataset, using the field cut fc as a data_source:

prj = yt.ProjectionPlot(
    ds,
    "z",
    [("fits", "temperature")],
    data_source=fc,
    origin="native",
    weight_field=("index", "ones"),
)  # "ones" weights each cell by 1
prj.set_log(("fits", "temperature"), True)
prj.show()

Finally, we can also take an existing ds9 region and use it to create a “cut region” as well, using ds9_region (the regions package needs to be installed for this):

from yt.frontends.fits.misc import ds9_region

For this example we’ll create a ds9 region from scratch and load it up:

region = 'galactic;box(+49:26:35.150,-0:30:04.410,1926.1927",1483.3701",0.0)'
box_reg = ds9_region(ds, region)

This region may now be used to compute derived quantities:

print(box_reg.quantities.extrema(("fits", "temperature")))

[-1.12755251  3.1104238 ] K

Or in projections:

prj = yt.ProjectionPlot(
    ds,
    "z",
    ("fits", "temperature"),
    origin="native",
    data_source=box_reg,
    weight_field=("index", "ones"),
)  # "ones" weights each cell by 1
prj.set_zlim(("fits", "temperature"), 1.0e-2, 1.5)
prj.set_log(("fits", "temperature"), True)
prj.show()