# Writing FITS Images¶

Notebook

yt has capabilities for writing 2D and 3D uniformly gridded data generated from datasets to FITS files. This is via the FITSImageData class. We'll test these capabilities out on an Athena dataset.

In [1]:
import yt
from yt.utilities.fits_image import FITSImageData, FITSProjection

In [2]:
ds = yt.load("MHDSloshing/virgo_low_res.0054.vtk", parameters={"length_unit":(1.0,"Mpc"),
"mass_unit":(1.0e14,"Msun"),
"time_unit":(1.0,"Myr")})

WARNING:yt:Supplying unit conversions from the parameters dict is deprecated, and will be removed in a future release. Use units_override instead.
WARNING:yt:Overriding code units. This is an experimental and potentially dangerous option that may yield inconsistent results, and must be used very carefully, and only if you know what you want from it.


## Creating FITS images from Slices and Projections¶

There are several ways to make a FITSImageData instance. The most intuitive ways are to use the FITSSlice, FITSProjection, FITSOffAxisSlice, and FITSOffAxisProjection classes to write slices and projections directly to FITS. To demonstrate a useful example of creating a FITS file, let's first make a ProjectionPlot:

In [3]:
prj = yt.ProjectionPlot(ds, "z", ["temperature"], weight_field="density", width=(500.,"kpc"))
prj.show()


Suppose that we wanted to write this projection to a FITS file for analysis and visualization in other programs, such as ds9. We can do that using FITSProjection:

In [4]:
prj_fits = FITSProjection(ds, "z", ["temperature"], weight_field="density")

/usr/lib64/python3.4/site-packages/IPython/kernel/__init__.py:13: ShimWarning: The IPython.kernel package has been deprecated. You should import from ipykernel or jupyter_client instead.
"You should import from ipykernel or jupyter_client instead.", ShimWarning)


which took the same parameters as ProjectionPlot except the width, because FITSProjection and FITSSlice always make slices and projections of the width of the domain size, at the finest resolution available in the simulation, in a unit determined to be appropriate for the physical size of the dataset.

We can call a number of the AstroPy HDUList class's methods from a FITSImageData object. For example, info shows us the contents of the virtual FITS file:

In [5]:
prj_fits.info()

Filename: (No file associated with this HDUList)
No.    Name         Type      Cards   Dimensions   Format
0    TEMPERATURE  PrimaryHDU      21   (2048, 2048)   float64


We can also look at the header for a particular field:

In [6]:
prj_fits["temperature"].header

Out[6]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -64 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                 2048
NAXIS2  =                 2048
EXTEND  =                    T
EXTNAME = 'TEMPERATURE'        / extension name
BTYPE   = 'temperature'
BUNIT   = 'K       '
WCSAXES =                    2
CRPIX1  =               1024.5
CRPIX2  =               1024.5
CDELT1  =             1.953125
CDELT2  =             1.953125
CUNIT1  = 'kpc     '
CUNIT2  = 'kpc     '
CTYPE1  = 'LINEAR  '
CTYPE2  = 'LINEAR  '
CRVAL1  =                  0.0
CRVAL2  =                  0.0
LATPOLE =                 90.0                                                  

where we can see that the temperature units are in Kelvin and the cell widths are in kiloparsecs. If we want the raw image data with units, we can call get_data:

In [7]:
prj_fits.get_data("temperature")

Out[7]:
YTArray([[ 29791393.9040929 ,  29791393.9040929 ,  29791393.9040929 , ...,
31976315.9281826 ,  31976315.9281826 ,  31976315.9281826 ],
[ 29791393.9040929 ,  29791393.9040929 ,  29791393.9040929 , ...,
31976315.9281826 ,  31976315.9281826 ,  31976315.9281826 ],
[ 29791393.9040929 ,  29791393.9040929 ,  29791393.9040929 , ...,
31976315.9281826 ,  31976315.9281826 ,  31976315.9281826 ],
...,
[ 30272987.3013772 ,  30272987.3013772 ,  30272987.3013772 , ...,
30351651.93357775,  30351651.93357775,  30351651.93357775],
[ 30272987.3013772 ,  30272987.3013772 ,  30272987.3013772 , ...,
30351651.93357775,  30351651.93357775,  30351651.93357775],
[ 30272987.3013772 ,  30272987.3013772 ,  30272987.3013772 , ...,
30351651.93357775,  30351651.93357775,  30351651.93357775]]) K

We can use the set_unit method to change the units of a particular field:

In [8]:
prj_fits.set_unit("temperature","R")
prj_fits.get_data("temperature")

Out[8]:
YTArray([[ 53624509.02736722,  53624509.02736722,  53624509.02736722, ...,
57557368.67072867,  57557368.67072867,  57557368.67072867],
[ 53624509.02736722,  53624509.02736722,  53624509.02736722, ...,
57557368.67072867,  57557368.67072867,  57557368.67072867],
[ 53624509.02736722,  53624509.02736722,  53624509.02736722, ...,
57557368.67072867,  57557368.67072867,  57557368.67072867],
...,
[ 54491377.14247895,  54491377.14247895,  54491377.14247895, ...,
54632973.48043995,  54632973.48043995,  54632973.48043995],
[ 54491377.14247895,  54491377.14247895,  54491377.14247895, ...,
54632973.48043995,  54632973.48043995,  54632973.48043995],
[ 54491377.14247895,  54491377.14247895,  54491377.14247895, ...,
54632973.48043995,  54632973.48043995,  54632973.48043995]]) R

The image can be written to disk using the writeto method:

In [9]:
prj_fits.writeto("sloshing.fits", clobber=True)


Since yt can read FITS image files, it can be loaded up just like any other dataset:

In [10]:
ds2 = yt.load("sloshing.fits")

WARNING:yt:Cannot find time


and we can make a SlicePlot of the 2D image, which shows the same data as the previous image:

In [11]:
slc2 = yt.SlicePlot(ds2, "z", ["temperature"], width=(500.,"kpc"))
slc2.set_log("temperature", True)
slc2.show()


## Using FITSImageData directly¶

If you want more fine-grained control over what goes into the FITS file, you can call FITSImageData directly, with various kinds of inputs. For example, you could use a FixedResolutionBuffer, and specify you want the units in parsecs instead:

In [12]:
slc3 = ds.slice(0, 0.0)
frb = slc3.to_frb((500.,"kpc"), 800)
fid_frb = FITSImageData(frb, fields=["density","temperature"], units="pc")


A 3D FITS cube can also be created from a covering grid:

In [13]:
cvg = ds.covering_grid(ds.index.max_level, [-0.5,-0.5,-0.5], [64, 64, 64], fields=["density","temperature"])
fid_cvg = FITSImageData(cvg, fields=["density","temperature"], units="Mpc")


## Other FITSImageData Methods¶

A FITSImageData instance can be generated from one previously written to disk using the from_file classmethod:

In [14]:
fid = FITSImageData.from_file("sloshing.fits")
fid.info()

Filename: (No file associated with this HDUList)
No.    Name         Type      Cards   Dimensions   Format
0    TEMPERATURE  PrimaryHDU      21   (2048, 2048)   float64


Multiple FITSImageData can be combined to create a new one, provided that the coordinate information is the same:

In [15]:
prj_fits2 = FITSProjection(ds, "z", ["density"])
prj_fits3 = FITSImageData.from_images([prj_fits, prj_fits2])
prj_fits3.info()

Filename: (No file associated with this HDUList)
No.    Name         Type      Cards   Dimensions   Format
0    TEMPERATURE  PrimaryHDU      21   (2048, 2048)   float64
1    DENSITY     ImageHDU        22   (2048, 2048)   float64


Alternatively, individual fields can be popped as well:

In [16]:
dens_fits = prj_fits3.pop("density")

In [17]:
dens_fits.info()

Filename: (No file associated with this HDUList)
No.    Name         Type      Cards   Dimensions   Format
0    DENSITY     PrimaryHDU      21   (2048, 2048)   float64

In [18]:
prj_fits3.info()

Filename: (No file associated with this HDUList)
No.    Name         Type      Cards   Dimensions   Format
0    TEMPERATURE  PrimaryHDU      21   (2048, 2048)   float64


So far, the FITS images we have shown have linear spatial coordinates. One may want to take a projection of an object and make a crude mock observation out of it, with celestial coordinates. For this, we can use the create_sky_wcs method. Specify a center (RA, Dec) coordinate in degrees, as well as a linear scale in terms of angle per distance:

In [19]:
sky_center = [30.,45.] # in degrees
sky_scale = (2.5, "arcsec/kpc") # could also use a YTQuantity
prj_fits.create_sky_wcs(sky_center, sky_scale, ctype=["RA---TAN","DEC--TAN"])


By the default, a tangent RA/Dec projection is used, but one could also use another projection using the ctype keyword. We can now look at the header and see it has the appropriate WCS:

In [20]:
prj_fits["temperature"].header

Out[20]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -64 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                 2048
NAXIS2  =                 2048
EXTEND  =                    T
EXTNAME = 'TEMPERATURE'        / extension name
BTYPE   = 'temperature'
BUNIT   = 'R       '
WCSAXES =                    2
CRPIX1  =               1024.5
CRPIX2  =               1024.5
CDELT1  =    -0.00135633680556
CDELT2  =     0.00135633680556
CUNIT1  = 'deg     '
CUNIT2  = 'deg     '
CTYPE1  = 'RA---TAN'
CTYPE2  = 'DEC--TAN'
CRVAL1  =                 30.0
CRVAL2  =                 45.0
LATPOLE =                 45.0
LONPOLE =                180.0                                                  

Finally, we can add header keywords to a single field or for all fields in the FITS image using update_header:

In [21]:
fid_frb.update_header("all", "time", 0.1) # Update all the fields
fid_frb.update_header("temperature", "scale", "Rankine") # Update just one field

In [22]:
print (fid_frb["density"].header["time"])

0.1