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
/usr/lib64/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
In [2]:
ds = yt.load("MHDSloshing/virgo_low_res.0054.vtk", units_override={"length_unit":(1.0,"Mpc"),
                                                                   "mass_unit":(1.0e14,"Msun"),
                                                                   "time_unit":(1.0,"Myr")})

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 = yt.FITSProjection(ds, "z", ["temperature"], weight_field="density")

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 also set width manually in FITSProjection. For example, set the width to 500 kiloparsec to get FITS file of the same projection plot as discussed above.

In [5]:
prj_fits = yt.FITSProjection(ds, "z", ["temperature"], weight_field="density", width=(500., "kpc"))

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 [6]:
prj_fits.info()
Filename: (No file associated with this FITSImageData)
No.    Name      Ver    Type      Cards   Dimensions   Format     Units
  0  TEMPERATURE    1 PrimaryHDU      22   (2048, 2048)   float64   K

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

In [7]:
prj_fits["temperature"].header
Out[7]:
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  =          0.244140625                                                  
CDELT2  =          0.244140625                                                  
CUNIT1  = 'kpc     '                                                            
CUNIT2  = 'kpc     '                                                            
CTYPE1  = 'LINEAR  '                                                            
CTYPE2  = 'LINEAR  '                                                            
CRVAL1  =                  0.0                                                  
CRVAL2  =                  0.0                                                  
LATPOLE =                 90.0                                                  
WCSNAME = 'yt      '                                                            

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 use the data attribute of this field:

In [8]:
prj_fits["temperature"].data
Out[8]:
YTArray([[34842120.10635341, 34842120.10635341, 34842120.10635341, ...,
          30841671.98158602, 30841671.98158602, 30841671.98158602],
         [34842120.10635341, 34842120.10635341, 34842120.10635341, ...,
          30841671.98158602, 30841671.98158602, 30841671.98158602],
         [34842120.10635341, 34842120.10635341, 34842120.10635341, ...,
          30841671.98158602, 30841671.98158602, 30841671.98158602],
         ...,
         [35518604.60016977, 35518604.60016977, 35518604.60016977, ...,
          31064918.70792199, 31064918.70792199, 31064918.70792199],
         [35518604.60016977, 35518604.60016977, 35518604.60016977, ...,
          31064918.70792199, 31064918.70792199, 31064918.70792199],
         [35518604.60016977, 35518604.60016977, 35518604.60016977, ...,
          31064918.70792199, 31064918.70792199, 31064918.70792199]]) K

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

In [9]:
prj_fits.set_unit("temperature","R")
prj_fits["temperature"].data
Out[9]:
YTArray([[62715816.19143613, 62715816.19143613, 62715816.19143613, ...,
          55515009.56685483, 55515009.56685483, 55515009.56685483],
         [62715816.19143613, 62715816.19143613, 62715816.19143613, ...,
          55515009.56685483, 55515009.56685483, 55515009.56685483],
         [62715816.19143613, 62715816.19143613, 62715816.19143613, ...,
          55515009.56685483, 55515009.56685483, 55515009.56685483],
         ...,
         [63933488.28030558, 63933488.28030558, 63933488.28030558, ...,
          55916853.67425957, 55916853.67425957, 55916853.67425957],
         [63933488.28030558, 63933488.28030558, 63933488.28030558, ...,
          55916853.67425957, 55916853.67425957, 55916853.67425957],
         [63933488.28030558, 63933488.28030558, 63933488.28030558, ...,
          55916853.67425957, 55916853.67425957, 55916853.67425957]]) R

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

In [10]:
prj_fits.writeto("sloshing.fits", clobber=True)
/tmp/yt/yt/utilities/parallel_tools/parallel_analysis_interface.py:319: VisibleDeprecationWarning: The "clobber" keyword argument is deprecated. Use the "overwrite" argument, which has the same effect, instead.
  return func(*args, **kwargs)

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

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

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

In [12]:
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 [13]:
slc3 = ds.slice(0, 0.0)
frb = slc3.to_frb((500.,"kpc"), 800)
fid_frb = yt.FITSImageData(frb, fields=["density","temperature"], units="pc")

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

In [14]:
cvg = ds.covering_grid(ds.index.max_level, [-0.5,-0.5,-0.5], [64, 64, 64], fields=["density","temperature"])
fid_cvg = yt.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 [15]:
fid = yt.FITSImageData.from_file("sloshing.fits")
fid.info()
Filename: sloshing.fits
No.    Name      Ver    Type      Cards   Dimensions   Format     Units
  0  TEMPERATURE    1 PrimaryHDU      22   (2048, 2048)   float64   R

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

In [16]:
prj_fits2 = yt.FITSProjection(ds, "z", ["density"])
prj_fits3 = yt.FITSImageData.from_images([prj_fits, prj_fits2])
prj_fits3.info()
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-16-359b44499d7b> in <module>()
      1 prj_fits2 = yt.FITSProjection(ds, "z", ["density"])
----> 2 prj_fits3 = yt.FITSImageData.from_images([prj_fits, prj_fits2])
      3 prj_fits3.info()

/tmp/yt/yt/visualization/fits_image.py in from_images(cls, image_list)
    452         first = True
    453         for fid in image_list:
--> 454             assert_same_wcs(w, fid.wcs)
    455             if img_shape != fid.shape:
    456                 raise RuntimeError("Images do not have the same shape!")

/tmp/yt/yt/visualization/fits_image.py in assert_same_wcs(wcs1, wcs2)
    604         assert wcs1.wcs.ctype[i] == wcs2.wcs.ctype[i]
    605     assert_allclose(wcs1.wcs.crpix, wcs2.wcs.crpix)
--> 606     assert_allclose(wcs1.wcs.cdelt, wcs2.wcs.cdelt)
    607     assert_allclose(wcs1.wcs.crval, wcs2.wcs.crval)
    608     crota1 = getattr(wcs1.wcs, "crota", None)

~/.local/lib64/python3.6/site-packages/numpy/testing/nose_tools/utils.py in assert_allclose(actual, desired, rtol, atol, equal_nan, err_msg, verbose)
   1394     header = 'Not equal to tolerance rtol=%g, atol=%g' % (rtol, atol)
   1395     assert_array_compare(compare, actual, desired, err_msg=str(err_msg),
-> 1396                          verbose=verbose, header=header, equal_nan=equal_nan)
   1397 
   1398 

~/.local/lib64/python3.6/site-packages/numpy/testing/nose_tools/utils.py in assert_array_compare(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf)
    777                                 verbose=verbose, header=header,
    778                                 names=('x', 'y'), precision=precision)
--> 779             raise AssertionError(msg)
    780     except ValueError:
    781         import traceback

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

(mismatch 100.0%)
 x: array([0.244141, 0.244141])
 y: array([1.953125, 1.953125])

Alternatively, individual fields can be popped as well to produce new instances of FITSImageData:

In [17]:
dens_fits = prj_fits3.pop("density")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-17-ff6a28f7ca73> in <module>()
----> 1 dens_fits = prj_fits3.pop("density")

NameError: name 'prj_fits3' is not defined

So this new instance would only have the "density" field:

In [18]:
dens_fits.info()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-18-b1d61db9f0f1> in <module>()
----> 1 dens_fits.info()

NameError: name 'dens_fits' is not defined

and the old one has the "density" field removed:

In [19]:
prj_fits3.info()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-19-f62c9a80cdd4> in <module>()
----> 1 prj_fits3.info()

NameError: name 'prj_fits3' is not defined

So far, the FITS images we have shown have linear spatial coordinates. We can see this by looking at the header for one of the fields, and examining the CTYPE1 and CTYPE2 keywords:

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.244140625                                                  
CDELT2  =          0.244140625                                                  
CUNIT1  = 'kpc     '                                                            
CUNIT2  = 'kpc     '                                                            
CTYPE1  = 'LINEAR  '                                                            
CTYPE2  = 'LINEAR  '                                                            
CRVAL1  =                  0.0                                                  
CRVAL2  =                  0.0                                                  
LATPOLE =                 90.0                                                  
WCSNAME = 'yt      '                                                            

The WCSNAME keyword is set to "yt" by default.

However, 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 [21]:
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 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 [22]:
prj_fits["temperature"].header
Out[22]:
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.00016954210069444                                                  
CDELT2  =  0.00016954210069444                                                  
CUNIT1  = 'deg     '                                                            
CUNIT2  = 'deg     '                                                            
CTYPE1  = 'RA---TAN'                                                            
CTYPE2  = 'DEC--TAN'                                                            
CRVAL1  =                 30.0                                                  
CRVAL2  =                 45.0                                                  
LATPOLE =                 45.0                                                  
WCSNAME = 'celestial'                                                           
LONPOLE =                180.0                                                  
RADESYS = 'ICRS    '                                                            

and now the WCSNAME has been set to "celestial". If you don't want to override the default WCS but to add another one, then you can make the call to create_sky_wcs and set replace_old_wcs=False:

In [23]:
prj_fits3.create_sky_wcs(sky_center, sky_scale, ctype=["RA---TAN","DEC--TAN"], replace_old_wcs=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-23-69c9adb331e9> in <module>()
----> 1 prj_fits3.create_sky_wcs(sky_center, sky_scale, ctype=["RA---TAN","DEC--TAN"], replace_old_wcs=False)

NameError: name 'prj_fits3' is not defined

We now can see that there are two WCSes in the header, with the celestial WCS keywords having the "A" designation:

In [24]:
prj_fits3["temperature"].header
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-24-40f4cf4633d6> in <module>()
----> 1 prj_fits3["temperature"].header

NameError: name 'prj_fits3' is not defined

Any further WCSes that are added will have "B", "C", etc.

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

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

(FITSImageData.ipynb; FITSImageData_evaluated.ipynb; FITSImageData.py)