X-ray FITS Images

%matplotlib inline
import numpy as np

import yt

This notebook shows how to use yt to make plots and examine FITS X-ray images and events files.

Sloshing, Shocks, and Bubbles in Abell 2052

This example uses data provided by Scott Randall, presented originally in Blanton, E.L., Randall, S.W., Clarke, T.E., et al. 2011, ApJ, 737, 99. They consist of two files, a “flux map” in counts/s/pixel between 0.3 and 2 keV, and a spectroscopic temperature map in keV.

ds = yt.load(
    "xray_fits/A2052_merged_0.3-2_match-core_tmap_bgecorr.fits",
    auxiliary_files=["xray_fits/A2052_core_tmap_b1_m2000_.fits"],
)

Since the flux and projected temperature images are in two different files, we had to use one of them (in this case the “flux” file) as a master file, and pass in the “temperature” file with the auxiliary_files keyword to load.

Next, let’s derive some new fields for the number of counts, the “pseudo-pressure”, and the “pseudo-entropy”:

def _counts(field, data):
    exposure_time = data.get_field_parameter("exposure_time")
    return data["fits", "flux"] * data["fits", "pixel"] * exposure_time


ds.add_field(
    ("gas", "counts"),
    function=_counts,
    sampling_type="cell",
    units="counts",
    take_log=False,
)


def _pp(field, data):
    return np.sqrt(data["gas", "counts"]) * data["fits", "projected_temperature"]


ds.add_field(
    ("gas", "pseudo_pressure"),
    function=_pp,
    sampling_type="cell",
    units="sqrt(counts)*keV",
    take_log=False,
)


def _pe(field, data):
    return data["fits", "projected_temperature"] * data["gas", "counts"] ** (-1.0 / 3.0)


ds.add_field(
    ("gas", "pseudo_entropy"),
    function=_pe,
    sampling_type="cell",
    units="keV*(counts)**(-1/3)",
    take_log=False,
)

/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1949: RuntimeWarning: divide by zero encountered in power
  out_arr = func(

Here, we’re deriving a “counts” field from the “flux” field by passing it a field_parameter for the exposure time of the time and multiplying by the pixel scale. Second, we use the fact that the surface brightness is strongly dependent on density (\(S_X \propto \rho^2\)) to use the counts in each pixel as a “stand-in”. Next, we’ll grab the exposure time from the primary FITS header of the flux file and create a YTQuantity from it, to be used as a field_parameter:

exposure_time = ds.quan(ds.primary_header["exposure"], "s")

Now, we can make the SlicePlot object of the fields we want, passing in the exposure_time as a field_parameter. We’ll also set the width of the image to 250 pixels.

slc = yt.SlicePlot(
    ds,
    "z",
    [
        ("fits", "flux"),
        ("fits", "projected_temperature"),
        ("gas", "pseudo_pressure"),
        ("gas", "pseudo_entropy"),
    ],
    origin="native",
    field_parameters={"exposure_time": exposure_time},
)
slc.set_log(("fits", "flux"), True)
slc.set_zlim(("fits", "flux"), 1e-5)
slc.set_log(("gas", "pseudo_pressure"), False)
slc.set_log(("gas", "pseudo_entropy"), False)
slc.set_width(250.0)
slc.show()

/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1949: RuntimeWarning: invalid value encountered in power
  out_arr = func(
/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1949: RuntimeWarning: invalid value encountered in multiply
  out_arr = func(
/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1824: RuntimeWarning: invalid value encountered in sqrt
  out_arr = func(np.asarray(inp), out=out_func, **kwargs)

To add the celestial coordinates to the image, we can use PlotWindowWCS, if you have a recent version of AstroPy (>= 1.3) installed:

from yt.frontends.fits.misc import PlotWindowWCS

wcs_slc = PlotWindowWCS(slc)
wcs_slc.show()

We can make use of yt’s facilities for profile plotting as well.

v, c = ds.find_max(("fits", "flux"))  # Find the maximum flux and its center
my_sphere = ds.sphere(c, (100.0, "code_length"))  # Radius of 150 pixels
my_sphere.set_field_parameter("exposure_time", exposure_time)

Such as a radial profile plot:

radial_profile = yt.ProfilePlot(
    my_sphere,
    "radius",
    ["counts", "pseudo_pressure", "pseudo_entropy"],
    n_bins=30,
    weight_field="ones",
)
radial_profile.set_log("counts", True)
radial_profile.set_log("pseudo_pressure", True)
radial_profile.set_log("pseudo_entropy", True)
radial_profile.set_xlim(3, 100.0)
radial_profile.show()

Or a phase plot:

phase_plot = yt.PhasePlot(
    my_sphere, "pseudo_pressure", "pseudo_entropy", ["counts"], weight_field=None
)
phase_plot.show()

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

from yt.frontends.fits.misc import ds9_region

reg_file = [
    "# Region file format: DS9 version 4.1\n",
    "global color=green dashlist=8 3 width=3 include=1 source=1\n",
    "FK5\n",
    'circle(15:16:44.817,+7:01:19.62,34.6256")',
]
f = open("circle.reg", "w")
f.writelines(reg_file)
f.close()
circle_reg = ds9_region(
    ds, "circle.reg", field_parameters={"exposure_time": exposure_time}
)

This region may now be used to compute derived quantities:

print(
    circle_reg.quantities.weighted_average_quantity("projected_temperature", "counts")
)

2.1562820167812964 keV

Or used in projections:

prj = yt.ProjectionPlot(
    ds,
    "z",
    [
        ("fits", "flux"),
        ("fits", "projected_temperature"),
        ("gas", "pseudo_pressure"),
        ("gas", "pseudo_entropy"),
    ],
    origin="native",
    field_parameters={"exposure_time": exposure_time},
    data_source=circle_reg,
    method="sum",
)
prj.set_log(("fits", "flux"), True)
prj.set_zlim(("fits", "flux"), 1e-5)
prj.set_log(("gas", "pseudo_pressure"), False)
prj.set_log(("gas", "pseudo_entropy"), False)
prj.set_width(250.0)
prj.show()

/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1949: RuntimeWarning: divide by zero encountered in power
  out_arr = func(
/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1949: RuntimeWarning: invalid value encountered in multiply
  out_arr = func(
/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1824: RuntimeWarning: invalid value encountered in sqrt
  out_arr = func(np.asarray(inp), out=out_func, **kwargs)
/tmp/yt/.venv/lib/python3.11/site-packages/unyt/array.py:1949: RuntimeWarning: invalid value encountered in power
  out_arr = func(

The Bullet Cluster

This example uses an events table file from a ~100 ks exposure of the “Bullet Cluster” from the Chandra Data Archive. In this case, the individual photon events are treated as particle fields in yt. However, you can make images of the object in different energy bands using the setup_counts_fields function.

from yt.frontends.fits.api import setup_counts_fields

load will handle the events file as FITS image files, and will set up a grid using the WCS information in the file. Optionally, the events may be reblocked to a new resolution. by setting the "reblock" parameter in the parameters dictionary in load. "reblock" must be a power of 2.

ds2 = yt.load("xray_fits/acisf05356N003_evt2.fits.gz", parameters={"reblock": 2})

setup_counts_fields will take a list of energy bounds (emin, emax) in keV and create a new field from each where the photons in that energy range will be deposited onto the image grid.

ebounds = [(0.1, 2.0), (2.0, 5.0)]
setup_counts_fields(ds2, ebounds)

The “x”, “y”, “energy”, and “time” fields in the events table are loaded as particle fields. Each one has a name given by “event_” plus the name of the field:

dd = ds2.all_data()
print(dd["io", "event_x"])
print(dd["io", "event_y"])

[4104.83642578 4726.79150391 4440.02783203 ... 3000.33911133 5874.03808594
 5138.02148438] code_length
[3971.82617188 4977.97851562 4273.42919922 ... 4864.12988281 3419.90551758
 3655.55444336] code_length

Now, we’ll make a plot of the two counts fields we made, and pan and zoom to the bullet:

slc = yt.SlicePlot(
    ds2, "z", [("gas", "counts_0.1-2.0"), ("gas", "counts_2.0-5.0")], origin="native"
)
slc.pan((100.0, 100.0))
slc.set_width(500.0)
slc.show()

The counts fields can take the field parameter "sigma" and use AstroPy’s convolution routines to smooth the data with a Gaussian:

slc = yt.SlicePlot(
    ds2,
    "z",
    [("gas", "counts_0.1-2.0"), ("gas", "counts_2.0-5.0")],
    origin="native",
    field_parameters={"sigma": 2.0},
)  # This value is in pixel scale
slc.pan((100.0, 100.0))
slc.set_width(500.0)
slc.set_zlim(("gas", "counts_0.1-2.0"), 0.01, 100.0)
slc.set_zlim(("gas", "counts_2.0-5.0"), 0.01, 50.0)
slc.show()