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
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 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 FITSImageData)
No.    Name         Type      Cards   Dimensions   Format     Units
  0  TEMPERATURE  PrimaryHDU      22   (2048, 2048)   float64   K

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                                                  
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 [7]:
prj_fits["temperature"].data
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["temperature"].data
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)
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 1.3 and will be removed in a future version. Use argument "overwrite" instead. [astropy.utils.decorators]

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

In [10]:
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 [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 = yt.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 = 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 [14]:
fid = yt.FITSImageData.from_file("sloshing.fits")
fid.info()
Filename: sloshing.fits
No.    Name         Type      Cards   Dimensions   Format     Units
  0  TEMPERATURE  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 [15]:
prj_fits2 = yt.FITSProjection(ds, "z", ["density"])
prj_fits3 = yt.FITSImageData.from_images([prj_fits, prj_fits2])
prj_fits3.info()
Filename: (No file associated with this FITSImageData)
No.    Name         Type      Cards   Dimensions   Format     Units
  0  TEMPERATURE  PrimaryHDU      22   (2048, 2048)   float64   R
  1  DENSITY     ImageHDU        23   (2048, 2048)   float64   g/cm**2

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

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

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

In [17]:
dens_fits.info()
Filename: (No file associated with this FITSImageData)
No.    Name         Type      Cards   Dimensions   Format     Units
  0  DENSITY     ImageHDU        23   (2048, 2048)   float64   g/cm**2

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

In [18]:
prj_fits3.info()
Filename: (No file associated with this FITSImageData)
No.    Name         Type      Cards   Dimensions   Format     Units
  0  TEMPERATURE  PrimaryHDU      22   (2048, 2048)   float64   R

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 [19]:
prj_fits["temperature"].header
Out[19]:
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  =             1.953125                                                  
CDELT2  =             1.953125                                                  
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 [20]:
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 [21]:
prj_fits["temperature"].header
Out[21]:
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.0013563368055556                                                  
CDELT2  =   0.0013563368055556                                                  
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 [22]:
prj_fits3.create_sky_wcs(sky_center, sky_scale, ctype=["RA---TAN","DEC--TAN"], replace_old_wcs=False)

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

In [23]:
prj_fits3["temperature"].header
Out[23]:
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.0013563368055556                                                  
CDELT2  =   0.0013563368055556                                                  
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    '                                                            
WCSAXESA=                    2                                                  
CRPIX1A =               1024.5                                                  
CRPIX2A =               1024.5                                                  
CDELT1A =  -0.0013563368055556                                                  
CDELT2A =   0.0013563368055556                                                  
CUNIT1A = 'deg     '                                                            
CUNIT2A = 'deg     '                                                            
CTYPE1A = 'RA---TAN'                                                            
CTYPE2A = 'DEC--TAN'                                                            
CRVAL1A =                 30.0                                                  
CRVAL2A =                 45.0                                                  
LONPOLEA=                180.0                                                  
LATPOLEA=                 45.0                                                  
WCSNAMEA= 'celestial'                                                           
RADESYSA= 'ICRS    '                                                            

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 [24]:
fid_frb.update_header("all", "time", 0.1) # Update all the fields
fid_frb.update_header("temperature", "scale", "Rankine") # Update just one field
In [25]:
print (fid_frb["density"].header["time"])
print (fid_frb["temperature"].header["scale"])
0.1
Rankine

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