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()