2.3. Input and Output Operations and Extension Manipulation - MEF
AstroData
is not intended to be Multi-Extension FITS (MEF) centric. The core
is independent of the file format. At Gemini, our data model uses MEF.
Therefore we have implemented a FITS handler that maps a MEF to the
internal AstroData
representation. A different handler can be implemented
for a different file format.
In this chapter, we present examples that will help the reader understand how
to access the information stored in a MEF with the AstroData
object and
understand that mapping.
Try it yourself
Download the data package (Try it yourself) if you wish to follow along and run the examples. Then
$ cd <path>/ad_usermanual/playground
$ python
2.3.1. Imports
Before doing anything, you need to import AstroData
and the Gemini instrument
configuration package gemini_instruments
.
>>> import astrodata
>>> import gemini_instruments
2.3.2. Open and access existing dataset
2.3.2.1. Read in the dataset
The file on disk is loaded into the AstroData
class associated with the
instrument the data is from. This association is done automatically based on
header content.
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> type(ad)
<class 'gemini_instruments.gmos.adclass.AstroDataGmos'>
From now on, ad
knows it is GMOS data. It knows how to access its headers
and when using the Recipe System (recipe_system
), it will trigger the
selection of the GMOS primitives and recipes.
The original path and filename are stored in the object. If you were to write
the AstroData
object to disk without specifying anything, those path and
filename would be used.
>>> ad.path
'../playdata/N20170609S0154.fits'
>>> ad.filename
'N20170609S0154.fits'
2.3.2.2. Accessing the content of a MEF file
Accessing pixel data, headers, and tables will be covered in detail in the following chapters. Here we just introduce the basic content interface.
For details on the AstroData
structure, please refer to the
previous chapter.
AstroData
uses NDData
as the core of its structure. Each FITS extension
becomes a NDAstroData
object, subclassed from NDData
, and is added to
a list.
2.3.2.2.1. Pixel data
To access pixel data, the list index and the .data
attribute are used. That
returns a numpy.ndarray
. The list of NDAstroData
is zero-indexed.
Extension number 1 in a MEF is index 0 in an |AstroData| object.
>>> ad = astrodata.open('../playdata/N20170609S0154_varAdded.fits')
>>> data = ad[0].data
>>> type(data)
<class 'numpy.ndarray'>
>>> data.shape
(2112, 256)
Remember that in a ndarray
the y-axis is the first number.
The variance and data quality planes, the VAR and DQ planes in Gemini MEF
files, are represented by the .variance
and .mask
attributes,
respectively. They are not their own “extension”, they don’t have their
own index in the list, unlike in a MEF. They are attached to the pixel data,
packaged together by the NDAstroData
object. They are represented as
numpy.ndarray
just like the pixel data
>>> var = ad[0].variance
>>> dq = ad[0].mask
2.3.2.2.2. Tables
Tables in the MEF file will also be loaded into the AstroData
object. If a table
is associated with a specific science extension through the EXTVER header keyword, that
table will be packaged within the same AstroData extension as the pixel data.
The AstroData
“extension” is the NDAstroData
object plus any table or other pixel
array. If the table is not associated with a specific extension and applies
globally, it will be added to the AstroData object as a global addition. No
indexing will be required to access it. In the example below, one OBJCAT
is
associated with each extension, while the REFCAT
has a global scope
>>> ad.info()
Filename: ../playdata/N20170609S0154_varAdded.fits
Tags: ACQUISITION GEMINI GMOS IMAGE NORTH OVERSCAN_SUBTRACTED OVERSCAN_TRIMMED
PREPARED SIDEREAL
Pixels Extensions
Index Content Type Dimensions Format
[ 0] science NDAstroData (2112, 256) float32
.variance ndarray (2112, 256) float32
.mask ndarray (2112, 256) uint16
.OBJCAT Table (6, 43) n/a
.OBJMASK ndarray (2112, 256) uint8
[ 1] science NDAstroData (2112, 256) float32
.variance ndarray (2112, 256) float32
.mask ndarray (2112, 256) uint16
.OBJCAT Table (8, 43) n/a
.OBJMASK ndarray (2112, 256) uint8
[ 2] science NDAstroData (2112, 256) float32
.variance ndarray (2112, 256) float32
.mask ndarray (2112, 256) uint16
.OBJCAT Table (7, 43) n/a
.OBJMASK ndarray (2112, 256) uint8
[ 3] science NDAstroData (2112, 256) float32
.variance ndarray (2112, 256) float32
.mask ndarray (2112, 256) uint16
.OBJCAT Table (5, 43) n/a
.OBJMASK ndarray (2112, 256) uint8
Other Extensions
Type Dimensions
.REFCAT Table (245, 16)
The tables are stored internally as astropy.table.Table
objects.
>>> ad[0].OBJCAT
<Table length=6>
NUMBER X_IMAGE Y_IMAGE ... REF_MAG_ERR PROFILE_FWHM PROFILE_EE50
int32 float32 float32 ... float32 float32 float32
------ ------- ------- ... ----------- ------------ ------------
1 283.461 55.4393 ... 0.16895 -999.0 -999.0
...
>>> type(ad[0].OBJCAT)
<class 'astropy.table.table.Table'>
>>> refcat = ad.REFCAT
>>> type(refcat)
<class 'astropy.table.table.Table'>
2.3.2.2.3. Headers
Headers are stored in the NDAstroData
.meta
attribute as astropy.io.fits.Header
objects,
which is a form of Python ordered dictionaries. Headers associated with extensions
are stored with the corresponding NDAstroData
object. The MEF Primary Header
Unit (PHU) is stored “globally” in the AstroData
object. Note that when slicing an AstroData
object,
for example copying over just the first extension, the PHU will follow. The
slice of an AstroData
object is an AstroData
object.
Headers can be accessed directly, or for some predefined concepts, the use of
Descriptors is preferred. See the chapters on headers for details.
Using Descriptors:
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> ad.filter_name()
'open1-6&g_G0301'
>>> ad.filter_name(pretty=True)
'g'
Using direct header access:
>>> ad.phu['FILTER1']
'open1-6'
>>> ad.phu['FILTER2']
'g_G0301'
Accessing the extension headers:
>>> ad.hdr['CCDSEC']
['[1:512,1:4224]', '[513:1024,1:4224]', '[1025:1536,1:4224]', '[1537:2048,1:4224]']
>>> ad[0].hdr['CCDSEC']
'[1:512,1:4224]'
With descriptors:
>>> ad.array_section(pretty=True)
['[1:512,1:4224]', '[513:1024,1:4224]', '[1025:1536,1:4224]', '[1537:2048,1:4224]']
2.3.3. Modify Existing MEF Files
Before you start modify the structure of an AstroData
object, you should be
familiar with it. Please make sure that you have read the previous chapter
on the structure of the AstroData object.
2.3.3.1. Appending an extension
In this section, we take an extension from one AstroData
object and append it
to another.
Here is an example appending a whole AstroData extension, with pixel data, variance, mask and tables.
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> advar = astrodata.open('../playdata/N20170609S0154_varAdded.fits')
>>> ad.info()
Filename: ../playdata/N20170609S0154.fits
Tags: ACQUISITION GEMINI GMOS IMAGE NORTH RAW SIDEREAL UNPREPARED
Pixels Extensions
Index Content Type Dimensions Format
[ 0] science NDAstroData (2112, 288) uint16
[ 1] science NDAstroData (2112, 288) uint16
[ 2] science NDAstroData (2112, 288) uint16
[ 3] science NDAstroData (2112, 288) uint16
>>> ad.append(advar[3])
>>> ad.info()
Filename: ../playdata/N20170609S0154.fits
Tags: ACQUISITION GEMINI GMOS IMAGE NORTH RAW SIDEREAL UNPREPARED
Pixels Extensions
Index Content Type Dimensions Format
[ 0] science NDAstroData (2112, 288) uint16
[ 1] science NDAstroData (2112, 288) uint16
[ 2] science NDAstroData (2112, 288) uint16
[ 3] science NDAstroData (2112, 288) uint16
[ 4] science NDAstroData (2112, 256) float32
.variance ndarray (2112, 256) float32
.mask ndarray (2112, 256) int16
.OBJCAT Table (5, 43) n/a
.OBJMASK ndarray (2112, 256) uint8
>>> ad[4].hdr['EXTVER']
4
>>> advar[3].hdr['EXTVER']
4
As you can see above, the fourth extension of advar
, along with everything
it contains was appended at the end of the first AstroData
object. However,
note that, because the EXTVER of the extension in advar
was 4, there are
now two extensions in ad
with this EXTVER. This is not a problem because
EXTVER is not used by AstroData
(it uses the index instead) and it is handled
only when the file is written to disk.
In this next example, we are appending only the pixel data, leaving behind the other associated data. One can attach the headers too, like we do here.
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> advar = astrodata.open('../playdata/N20170609S0154_varAdded.fits')
>>> ad.append(advar[3].data, header=advar[3].hdr)
>>> ad.info()
Filename: ../playdata/N20170609S0154.fits
Tags: ACQUISITION GEMINI GMOS IMAGE NORTH RAW SIDEREAL UNPREPARED
Pixels Extensions
Index Content Type Dimensions Format
[ 0] science NDAstroData (2112, 288) uint16
[ 1] science NDAstroData (2112, 288) uint16
[ 2] science NDAstroData (2112, 288) uint16
[ 3] science NDAstroData (2112, 288) uint16
[ 4] science NDAstroData (2112, 256) float32
Notice how a new extension was created but variance
, mask
, the OBJCAT
table and OBJMASK image were not copied over. Only the science pixel data was
copied over.
Please note, there is no implementation for the “insertion” of an extension.
2.3.3.2. Removing an extension or part of one
Removing an extension or a part of an extension is straightforward. The
Python command del()
is used on the item to remove. Below are a few
examples, but first let us load a file
>>> ad = astrodata.open('../playdata/N20170609S0154_varAdded.fits')
>>> ad.info()
As you go through these examples, check the new structure with ad.info()
after every removal to see how the structure has changed.
Deleting a whole AstroData
extension, the fourth one
>>> del ad[3]
Deleting only the variance array from the second extension
>>> ad[1].variance = None
Deleting a table associated with the first extension
>>> del ad[0].OBJCAT
Deleting a global table, not attached to a specific extension
>>> del ad.REFCAT
2.3.4. Writing back to disk
The AstroData
layer takes care of converting
the AstroData
object back to a MEF file on disk. When writing to disk,
one should be aware of the path and filename information associated
with the AstroData
object.
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> ad.path
'../playdata/N20170609S0154.fits'
>>> ad.filename
'N20170609S0154.fits'
2.3.4.1. Writing to a new file
There are various ways to define the destination for the new FITS file. The most common and natural way is
>>> ad.write('new154.fits')
>>> ad.write('new154.fits', overwrite=True)
This will write a FITS file named ‘new154.fits’ in the current directory.
With overwrite=True
, it will overwrite the file if it already exists.
A path can be prepended to the filename if the current directory is not
the destination.
Note that ad.filename
and ad.path
have not changed, we have just
written to the new file, the AstroData
object is in no way associated
with that new file.
>>> ad.path
'../playdata/N20170609S0154.fits'
>>> ad.filename
'N20170609S0154.fits'
If you want to create that association, the ad.filename
and ad.path
needs to be modified first. For example:
>>> ad.filename = 'new154.fits'
>>> ad.write(overwrite=True)
>>> ad.path
'../playdata/new154.fits'
>>> ad.filename
'new154.fits'
Changing ad.filename
also changes the filename in the ad.path
. The
sequence above will write ‘new154.fits’ not in the current directory but
rather to the directory that is specified in ad.path
.
WARNING: ad.write()
has an argument named filename
. Setting filename
in the call to ad.write()
, as in ad.write(filename='new154.fits')
will NOT
modify ad.filename
or ad.path
. The two “filenames”, one a method argument
the other a class attribute have no association to each other.
2.3.4.2. Updating an existing file on disk
Updating an existing file on disk requires explicitly allowing overwrite.
If you have not written ‘new154.fits’ to disk yet (from previous section)
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> ad.write('new154.fits', overwrite=True)
Now let’s open ‘new154.fits’, and write to it
>>> adnew = astrodata.open('new154.fits')
>>> adnew.write(overwrite=True)
2.3.4.3. A note on FITS header keywords
When writing an AstroData
object to disk as a FITS file, it is necessary to add or
update header keywords to represent some of the internally-stored information. Any
extensions that did not originally belong to this AstroData
will be assigned new
EXTVER keywords to avoid conflicts with existing extensions, and the internal WCS is
converted to the appropriate FITS keywords. Note that in some cases it may not be
possible for standard FITS keywords to accurately represent the true WCS. In such
cases, the FITS keywords are written as an approximation to the true WCS, together
with an additional keyword
FITS-WCS= 'APPROXIMATE' / FITS WCS is approximate
to indicate this. The accurate WCS is written as an additional FITS extension with
EXTNAME='WCS'
that AstroData will recognize when the file is read back in. The
WCS
extension will not be written to disk if there is an accurate FITS
representation of the WCS (e.g., for a simple image).
2.3.5. Create New MEF Files
A new MEF file can be created from an existing, maybe modified, file or it can be created from scratch. We discuss both cases here.
2.3.5.1. Create New Copy of MEF Files
To create a new copy of a MEF file, modified or not, the user has already been given most of the tools in the sections above. Yet, let’s throw a couple examples for completeness.
2.3.5.1.1. Basic example
As seen above, a MEF file can be opened with astrodata
, the AstroData
object can be modified (or not), and then written back to disk under a
new name.
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
... optional modifications here ...
>>> ad.write('newcopy.fits')
2.3.5.1.2. Needing true copies in memory
Sometimes it is a true copy in memory that is needed. This is not specific
to MEF. In Python, doing something like adnew = ad
does not create a
new copy of the AstrodData object; it just gives it a new name. If you
modify adnew
you will be modifying ad
too. They point to the same
block of memory.
To create a true independent copy, the deepcopy
utility needs to be used.
>>> from copy import deepcopy
>>> ad = astrodata.open('../playdata/N20170609S0154.fits')
>>> adcopy = deepcopy(ad)
Be careful using deepcopy
, your memory could balloon really fast. Use it
only when truly needed.
2.3.5.2. Create New MEF Files from Scratch
Before one creates a new MEF file on disk, one has to create the AstroData
object that will be eventually written to disk. The AstroData
object
created also needs to know that it will have to be written using the MEF
format. This is fortunately handled fairly transparently by astrodata
.
The key to associating the FITS data to the AstroData
object is simply to
create the AstroData
object from astropy.io.fits
header objects. Those
will be recognized by astrodata
as FITS and the constructor for FITS will be
used. The user does not need to do anything else special. Here is how it is
done.
2.3.5.2.1. Create a MEF with basic header and data array set to zeros
>>> import numpy as np
>>> from astropy.io import fits
>>> phu = fits.PrimaryHDU()
>>> pixel_data = np.zeros((100,100))
>>> hdu = fits.ImageHDU()
>>> hdu.data = pixel_data
>>> ad = astrodata.create(phu)
>>> ad.append(hdu, name='SCI')
or another way to do the last two blocs:
>>> hdu = fits.ImageHDU(data=pixel_data, name='SCI')
>>> ad = astrodata.create(phu, [hdu])
Then it is just a matter of calling ad.write('somename.fits')
on that
new Astrodata
object.
2.3.5.2.2. Associate a pixel array with a science pixel array
Only main science (“SCI”) pixel arrays are added as slices to an astrodata object. It not uncommon to have pixels information associated with those main science pixels, for example an object mask where marked pixels in the mask are directly associated with sources in the science array.
Such pixel arrays are added to specific slice of the astrodata object they are associated with.
Building on the astrodata object we created in the previous subsection, one would add a pixel array to the first slice of the astrodata object as follows:
>>> extra_data = np.ones((100, 100))
>>> ad[0].EXTRADATA = extra_data
When the file is written to disk as a MEF, an extension will be created with
EXTNAME = EXTRADATA
and an EXTVER
that matches the slice’s EXTVER
,
in this case is would be 1
.
2.3.5.2.3. Represent a table as a FITS binary table in an AstroData
object
One first needs to create a table, either an astropy.table.Table
or a BinTableHDU
. See the Astropy documentation
on tables and this manual’s section dedicated to tables for
more information.
In the first example, we assume that my_astropy_table
is
a Table
ready to be attached to an AstroData
object. (Warning: we have not created my_astropy_table
therefore the
example below will not run, though this is how it would be done.)
>>> phu = fits.PrimaryHDU()
>>> ad = astrodata.create(phu)
>>> astrodata.add_header_to_table(my_astropy_table)
>>> ad.append(my_astropy_table, name='SMAUG')
In the second example, we start with a FITS BinTableHDU
and attach it to a new AstroData
object. (Again, we have not created
my_fits_table
so the example will not run.)
>>> phu = fits.PrimaryHDU()
>>> ad = astrodata.create(phu)
>>> ad.append(my_fits_table, name='DROGON')
As before, once the AstroData
object is constructed, the ad.write()
method can be used to write it to disk as a MEF file.