Pyfits, numpy, matplotlib and tables are the four most useful tools for analyzing astronomical data. numpy gives Python the ability to do array arithmetic like IDL. Pyfits allows us to get the data from FITS files into these arrays. Matplotlib is a simple-to-use yet powerful plotting package, and tables provides powerful tools for manipulating tabular data.
You should have numpy, pyfits and matplotlib available to you through ureka, while the tables modules are in astropy.
From within Python, you have to import a package to make it available to you; this is rather different from IDL where as long as something is on your IDL search path, it's available without you doing anything else.
import pyfits
pyfits.__version__
'3.2.dev1944'
import numpy as np
np.version.version
'1.6.2'
plot()
[]
The first thing you need to be able to do is get the data from a FITS file and put it into an array. First let's find a FITS file. IPython allows you to use several shell commands from within Python, so we can type 'cd' and 'ls' and IPython executes them. The regular Python interpreter doesn't do this, PyRAF does to a certain extent.
pwd
u'/home/robert/work/jwst'
ls
DQInit/ pyfits/ dqinit.tar Session2_RIJ.ipynb Exposure_Log_FM_Phase5_Latest_29Jul11.csv* test/ fits_generator/ Ureka/ install_ureka_v1.0beta3* Ureka_linux-rhe5_64_dev.tar NRS1_dark88.fits
We have NRS1_dark88.fits. This is a NIRSpec dark exposure made up of 88 groups, each group holds a 10.5 second exposure and the data is read non-destructively, like NICMOS and WFC3IR. We can inspect the contents of this without actually opening it:
pyfits.info('NRS1_dark88.fits')
Filename: NRS1_dark88.fits No. Name Type Cards Dimensions Format 0 PRIMARY PrimaryHDU 200 () uint8 1 SCI ImageHDU 13 (2048, 2048, 88, 1) int16 2 ERR ImageHDU 10 (2048, 2048, 88, 1) float32 3 DQ ImageHDU 10 (2048, 2048, 88, 1) int16
The most flexible way of opening this file and getting at the data and header information is to use the pyfits open command
h = pyfits.open('NRS1_dark88.fits')
h is a handle or pointer to all the information in the FITS file.
help(h)
Help on HDUList in module pyfits.hdu.hdulist object: class HDUList(__builtin__.list, pyfits.verify._Verify) | HDU list class. This is the top-level FITS object. When a FITS | file is opened, a `HDUList` object is returned. | | Method resolution order: | HDUList | __builtin__.list | pyfits.verify._Verify | __builtin__.object | | Methods defined here: | | __delitem__(self, key) | Delete an HDU from the `HDUList`, indexed by number or name. | | __delslice__(self, start, stop) | Delete a slice of HDUs from the `HDUList`, indexed by number only. | | __enter__(self) | # Support the 'with' statement | | __exit__(self, type, value, traceback) | | __getitem__(self, key) | Get an HDU from the `HDUList`, indexed by number or name. | | __getslice__(self, start, end) | | __init__(self, hdus=[], file=None) | Construct a `HDUList` object. | | Parameters | ---------- | hdus : sequence of HDU objects or single HDU, optional | The HDU object(s) to comprise the `HDUList`. Should be | instances of `_BaseHDU`. | | file : file object, optional | The opened physical file associated with the `HDUList`. | | __iter__(self) | | __setitem__(self, key, hdu) | Set an HDU to the `HDUList`, indexed by number or name. | | append(self, hdu) | Append a new HDU to the `HDUList`. | | Parameters | ---------- | hdu : instance of _BaseHDU | HDU to add to the `HDUList`. | | close(self, output_verify='exception', verbose=False, closed=True) | Close the associated FITS file and memmap object, if any. | | Parameters | ---------- | output_verify : str | Output verification option. Must be one of ``"fix"``, | ``"silentfix"``, ``"ignore"``, ``"warn"``, or | ``"exception"``. See :ref:`verify` for more info. | | verbose : bool | When `True`, print out verbose messages. | | closed : bool | When `True`, close the underlying file object. | | fileinfo(self, index) | Returns a dictionary detailing information about the locations | of the indexed HDU within any associated file. The values are | only valid after a read or write of the associated file with | no intervening changes to the `HDUList`. | | Parameters | ---------- | index : int | Index of HDU for which info is to be returned. | | Returns | ------- | fileinfo : dict or None | | The dictionary details information about the locations of | the indexed HDU within an associated file. Returns `None` | when the HDU is not associated with a file. | | Dictionary contents: | | ========== ======================================================== | Key Value | ========== ======================================================== | file File object associated with the HDU | filename Name of associated file object | filemode Mode in which the file was opened (readonly, | copyonwrite, update, append, denywrite, ostream) | resized Flag that when `True` indicates that the data has been | resized since the last read/write so the returned values | may not be valid. | hdrLoc Starting byte location of header in file | datLoc Starting byte location of data block in file | datSpan Data size including padding | ========== ======================================================== | | filename(self) | Return the file name associated with the HDUList object if one exists. | Otherwise returns None. | | Returns | ------- | filename : a string containing the file name associated with the | HDUList object if an association exists. Otherwise returns | None. | | flush(*args, **kwargs) | Force a write of the `HDUList` back to the file (for append and | update modes only). | | Parameters | ---------- | output_verify : str | Output verification option. Must be one of ``"fix"``, | ``"silentfix"``, ``"ignore"``, ``"warn"``, or | ``"exception"``. See :ref:`verify` for more info. | | verbose : bool | When `True`, print verbose messages | | index_of(self, key) | Get the index of an HDU from the `HDUList`. | | Parameters | ---------- | key : int, str or tuple of (string, int) | The key identifying the HDU. If `key` is a tuple, it is of | the form (`key`, `ver`) where `ver` is an ``EXTVER`` value | that must match the HDU being searched for. | | Returns | ------- | index : int | The index of the HDU in the `HDUList`. | | info(self, output=None) | Summarize the info of the HDUs in this `HDUList`. | | Note that this function prints its results to the console---it | does not return a value. | | Parameters | ---------- | output : file, optional | A file-like object to write the output to. If False, does not | output to a file and instead returns a list of tuples representing | the HDU info. Writes to sys.stdout by default. | | insert(self, index, hdu) | Insert an HDU into the `HDUList` at the given `index`. | | Parameters | ---------- | index : int | Index before which to insert the new HDU. | | hdu : _BaseHDU instance | The HDU object to insert | | readall(self) | Read data of all HDUs into memory. | | update_extend(self) | Make sure that if the primary header needs the keyword | ``EXTEND`` that it has it and it is correct. | | writeto(self, fileobj, output_verify='exception', clobber=False, checksum=False) | Write the `HDUList` to a new file. | | Parameters | ---------- | fileobj : file path, file object or file-like object | File to write to. If a file object, must be opened for | append (ab+). | | output_verify : str | Output verification option. Must be one of ``"fix"``, | ``"silentfix"``, ``"ignore"``, ``"warn"``, or | ``"exception"``. See :ref:`verify` for more info. | | clobber : bool | When `True`, overwrite the output file if exists. | | checksum : bool | When `True` adds both ``DATASUM`` and ``CHECKSUM`` cards | to the headers of all HDU's written to the file. | | ---------------------------------------------------------------------- | Class methods defined here: | | fromfile(cls, fileobj, mode='readonly', memmap=False, save_backup=False, **kwargs) from __builtin__.type | Creates an HDUList instance from a file-like object. | | The actual implementation of :func:`fitsopen`, and generally shouldn't | be used directly. Use :func:`pyfits.open` instead (and see its | documentation for details of the parameters accepted by this method). | | fromstring(cls, data, **kwargs) from __builtin__.type | Creates an HDUList instance from a string or other in-memory data | buffer containing an entire FITS file. Similar to | :meth:`HDUList.fromfile`, but does not accept the mode or memmap | arguments, as they are only relevant to reading from a file on disk. | | This is useful for interfacing with other libraries such as CFITSIO, | and may also be useful for streaming applications. | | Parameters | ---------- | data : str, buffer, memoryview, etc. | A string or other memory buffer containing an entire FITS file. It | should be noted that if that memory is read-only (such as a Python | string) the returned :class:`HDUList`'s data portions will also be | read-only. | | kwargs : dict | Optional keyword arguments. See :func:`pyfits.open` for details. | | Returns | ------- | hdul : HDUList | An :class:`HDUList` object representing the in-memory FITS file. | | ---------------------------------------------------------------------- | Data descriptors defined here: | | __dict__ | dictionary for instance variables (if defined) | | __weakref__ | list of weak references to the object (if defined) | | ---------------------------------------------------------------------- | Methods inherited from __builtin__.list: | | __add__(...) | x.__add__(y) <==> x+y | | __contains__(...) | x.__contains__(y) <==> y in x | | __eq__(...) | x.__eq__(y) <==> x==y | | __ge__(...) | x.__ge__(y) <==> x>=y | | __getattribute__(...) | x.__getattribute__('name') <==> x.name | | __gt__(...) | x.__gt__(y) <==> x>y | | __iadd__(...) | x.__iadd__(y) <==> x+=y | | __imul__(...) | x.__imul__(y) <==> x*=y | | __le__(...) | x.__le__(y) <==> x<=y | | __len__(...) | x.__len__() <==> len(x) | | __lt__(...) | x.__lt__(y) <==> x<y | | __mul__(...) | x.__mul__(n) <==> x*n | | __ne__(...) | x.__ne__(y) <==> x!=y | | __repr__(...) | x.__repr__() <==> repr(x) | | __reversed__(...) | L.__reversed__() -- return a reverse iterator over the list | | __rmul__(...) | x.__rmul__(n) <==> n*x | | __setslice__(...) | x.__setslice__(i, j, y) <==> x[i:j]=y | | Use of negative indices is not supported. | | __sizeof__(...) | L.__sizeof__() -- size of L in memory, in bytes | | count(...) | L.count(value) -> integer -- return number of occurrences of value | | extend(...) | L.extend(iterable) -- extend list by appending elements from the iterable | | index(...) | L.index(value, [start, [stop]]) -> integer -- return first index of value. | Raises ValueError if the value is not present. | | pop(...) | L.pop([index]) -> item -- remove and return item at index (default last). | Raises IndexError if list is empty or index is out of range. | | remove(...) | L.remove(value) -- remove first occurrence of value. | Raises ValueError if the value is not present. | | reverse(...) | L.reverse() -- reverse *IN PLACE* | | sort(...) | L.sort(cmp=None, key=None, reverse=False) -- stable sort *IN PLACE*; | cmp(x, y) -> -1, 0, 1 | | ---------------------------------------------------------------------- | Data and other attributes inherited from __builtin__.list: | | __hash__ = None | | __new__ = <built-in method __new__ of type object> | T.__new__(S, ...) -> a new object with type S, a subtype of T | | ---------------------------------------------------------------------- | Methods inherited from pyfits.verify._Verify: | | run_option(self, option='warn', err_text='', fix_text='Fixed.', fix=None, fixable=True) | Execute the verification with selected option. | | verify(self, option='warn') | Verify all values in the instance. | | Parameters | ---------- | option : str | Output verification option. Must be one of ``"fix"``, | ``"silentfix"``, ``"ignore"``, ``"warn"``, or | ``"exception"``. See :ref:`verify` for more info.
type(h)
pyfits.hdu.hdulist.HDUList
pyfits.open returns something called an HDUList, which behaves like a list of header/data units.
h
[<pyfits.hdu.image.PrimaryHDU at 0x3e62e50>, <pyfits.hdu.image.ImageHDU at 0x3e9d150>, <pyfits.hdu.image.ImageHDU at 0x3e9d450>, <pyfits.hdu.image.ImageHDU at 0x3e9d750>]
A list is Python data structure that is an ordered collection of objects, enclosed in square brackets and with its members separated by commas.
This list has four members, a PrimaryHDU object and three ImageHDU objects. These are what we saw when we did pyfits.info(), a single Primary HDU and 3 Image HDUs, the SCI, ERR and DQ HDUs.
We can refer to the individual members of a list using C-like indexing, starting at 0.
h[0]
<pyfits.hdu.image.PrimaryHDU at 0x3e62e50>
h[2]
<pyfits.hdu.image.ImageHDU at 0x3e9d450>
h[4]
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-12-0b2c789b69f7> in <module>() ----> 1 h[4] /home/robert/work/jwst/Ureka/variants/common/lib/python/pyfits/hdu/hdulist.pyc in __getitem__(self, key) 174 175 idx = self.index_of(key) --> 176 return super(HDUList, self).__getitem__(idx) 177 178 def __setitem__(self, key, hdu): IndexError: list index out of range
The index goes from 0 to 3, so if we try to access element #4, Python will give us an error, or 'throw an exception'
In this case, we are throwing an IndexError exception, because we tried to access a list element with an index outside the range of indices of h.
The data we are interested in is in extension #1. We can't get at it using the 'SCI' name like we can with IRAF, we have to use its index.
cube = h[1].data
Now I have the 4-dimensional pixel data from the SCI extension in the 4-dimension numpy array 'cube'. This cube has 1 integration of 88 groups, each with a 2048 by 2048 2-dimensional array
We could have read the data into our array by using the convenience functions as mentioned in the previous lecture:
cube = pyfits.getdata('NRS1_dark88.fits')
which avoids using the 'Open' function of pyfits. Similarly, you can get the header from a FITS file by doing
header = pyfits.getheader('NRS1_dark88.fits')
The pyfits documentation at http://packages.python.org/pyfits/ will give more information.
cube.shape
(1, 88, 2048, 2048)
We can investigate the shape of this array by inspecting the 'shape' attribute of the cube object. All numpy arrays have a shape attribute. The ordering is C-like rather than FORTRAN-like, with the most rapidly-varying index last rather than first. Notice how the ordering in the pyfits.info() function earlier is IRAF/Fortran like, probably from a historical desire to make the header listing somewhat similar to IRAF's FITS header listings.
numpy arrays have many attributes that you can inspect. Some of these are variables, others are functions. You can find out all the atributes of an object by using the Python built-in 'dir' function:
dir(cube)
['T', '__abs__', '__add__', '__and__', '__array__', '__array_finalize__', '__array_interface__', '__array_prepare__', '__array_priority__', '__array_struct__', '__array_wrap__', '__class__', '__contains__', '__copy__', '__deepcopy__', '__delattr__', '__delitem__', '__delslice__', '__div__', '__divmod__', '__doc__', '__eq__', '__float__', '__floordiv__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__getslice__', '__gt__', '__hash__', '__hex__', '__iadd__', '__iand__', '__idiv__', '__ifloordiv__', '__ilshift__', '__imod__', '__imul__', '__index__', '__init__', '__int__', '__invert__', '__ior__', '__ipow__', '__irshift__', '__isub__', '__iter__', '__itruediv__', '__ixor__', '__le__', '__len__', '__long__', '__lshift__', '__lt__', '__mod__', '__mul__', '__ne__', '__neg__', '__new__', '__nonzero__', '__oct__', '__or__', '__pos__', '__pow__', '__radd__', '__rand__', '__rdiv__', '__rdivmod__', '__reduce__', '__reduce_ex__', '__repr__', '__rfloordiv__', '__rlshift__', '__rmod__', '__rmul__', '__ror__', '__rpow__', '__rrshift__', '__rshift__', '__rsub__', '__rtruediv__', '__rxor__', '__setattr__', '__setitem__', '__setslice__', '__setstate__', '__sizeof__', '__str__', '__sub__', '__subclasshook__', '__truediv__', '__xor__', 'all', 'any', 'argmax', 'argmin', 'argsort', 'astype', 'base', 'byteswap', 'choose', 'clip', 'compress', 'conj', 'conjugate', 'copy', 'ctypes', 'cumprod', 'cumsum', 'data', 'diagonal', 'dot', 'dtype', 'dump', 'dumps', 'fill', 'flags', 'flat', 'flatten', 'getfield', 'imag', 'item', 'itemset', 'itemsize', 'max', 'mean', 'min', 'nbytes', 'ndim', 'newbyteorder', 'nonzero', 'prod', 'ptp', 'put', 'ravel', 'real', 'repeat', 'reshape', 'resize', 'round', 'searchsorted', 'setasflat', 'setfield', 'setflags', 'shape', 'size', 'sort', 'squeeze', 'std', 'strides', 'sum', 'swapaxes', 'take', 'tofile', 'tolist', 'tostring', 'trace', 'transpose', 'var', 'view']
Note that it tells you ALL of the attributes of the cube object, including those preceded by one or more underscores, which are conventionally meant to be 'private'.
IPython lets you do this more simply by just hitting the TAB key after typing cube and then a period. You only see the public attributes, not the private ones.
cube.
File "<ipython-input-32-82ef5b289c03>", line 1 cube. ^ SyntaxError: invalid syntax
Let's say we want to plot the ramp for the pixel in column 200 and row 100. We specify what to plot using the following notation:
cube[0,:,99,199]
array([ 16792., 16812., 16776., 16804., 16786., 16781., 16729., 16682., 16744., 16743., 16727., 16793., 16742., 16736., 16743., 16813., 16761., 16785., 16732., 16787., 16752., 16743., 16782., 16689., 16801., 16809., 16763., 16747., 16729., 16779., 16785., 16770., 16846., 16757., 16781., 16736., 16782., 16803., 16858., 16749., 16780., 16761., 16739., 16773., 16771., 16769., 16719., 16765., 16812., 16783., 16793., 16791., 16789., 16821., 16858., 16848., 16794., 16796., 16858., 16825., 16864., 16847., 16885., 16887., 16821., 16782., 16843., 16796., 16840., 16796., 16808., 16837., 16862., 16848., 16828., 16830., 16807., 16875., 16812., 16905., 16824., 16848., 16843., 16828., 16897., 16828., 16831., 16856.], dtype=float32)
Here's what you need to notice:
You can find out more about numpy by going to the documentation at http://docs.scipy.org/doc/numpy/reference/
pylab.plot(cube[0,:,99,199])
[<matplotlib.lines.Line2D at 0x40f35d0>]
We can overplot the ramp for the pixel to the right of this:
pylab.plot(cube[0,:,99,200])
[<matplotlib.lines.Line2D at 0x432ead0>]
As you can see, the ramps have almost identical bumps and wiggles, but are vertically offset by about 600 counts. There is a lot of high-frequency (i.e. pixel-to-pixel) structure in the reset - 600 counts from one pixel to the next. This makes it very hard to see low-signal detail in the raw ramps even if you look at the last one which has the most signal. To alleviate this, we can subtract the first read from every group:
newcube = cube - cube[0,0]
Notice how we didn't need to specify any : or * to tell numpy to use all the groups, rows and columns. Also notice that numpy handles the fact that 'cube' and 'cube[0,0]' have different shapes. This is a feature of numpy known as 'broadcasting'. Numpy looks at the two operands of the '-' operator. cube has a shape of [1,88,2048,2048], while cube[0,0] has a shape of [2048,2048]. To make the operation possible, numpy 'broadcasts' the second operand to be [1, 88, 2048, 2048] by replicating it as needed to make the shapes the same.
newcube.shape
(1, 88, 2048, 2048)
plot(newcube[0,:,99,199])
[<matplotlib.lines.Line2D at 0x4365750>]
That's no good, pylab overplots by default. We can fix this by using the 'hold' keyword argument:
plot(newcube[0,:,99,199],hold=False)
[<matplotlib.lines.Line2D at 0x4639050>]
plot(newcube[0,:,99,199])
plot(newcube[0,:,99,200])
[<matplotlib.lines.Line2D at 0x465bc50>]
Now you can see how well the ramps sit over the top of each other. The bumps and wiggles can actually be taken out for the most part by using the reference pixels.
Most of what you'll do in Matplotlib will be using the pylab interface, which is where the 'plot' command lives. There is a point of confusion in the documentation, however, in that the official name of this interface is pyplot. Most people seem to still use pylab and even the documentation to pyplot refers to pylab all over the place. The point is, to find out more about the plot command, the best thing to google is the pyplot tutorial.
http://matplotlib.org/users/pyplot_tutorial.html
There you will find out that you can add a title and axis labels by doing something like this:
plot(newcube[0,:,99,199])
plot(newcube[0,:,99,200])
xlabel('Read number (1 read every $10\mu{}s$)')
ylabel('Counts')
title('NIRSpec dark ramp')
legend(('Pixel [100,200]','Pixel [200,201]'), 'lower right')
<matplotlib.legend.Legend at 0x3bf9190>
You can change the color and change to plotting symbols using something like this:
plot(newcube[0,:,99,199],'r+',hold=False)
[<matplotlib.lines.Line2D at 0x4c10510>]
This plots each point as a red + symbol. Change the letter to one of the following:
look at the documentation for the pyplot plot command to find out what options exist for choosing the plot style
http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot
If you want to see examples of how to create any plot imaginable, go to the matplot gallery at http://matplotlib.org/gallery.html
numpy allows very flexible specification of indices. The most general specification looks like
array[start:end:stride]
Start is the index of the first item wanted
end is ONE MORE THAN the index of the last item wanted
stride is the increment that is traversed along the array
We call the indices [start:end:stride] a SLICE.
As an example: data[1:100:2] means "take every other element, starting with the second and ending with the 100th (i.e. the one with index 99)"
You can make all of these negative too. In the case of strides, it just means go in the reverse direction. So specifying data[::-1] gives the array reversed.
If the start or end index is negative, it means "starting from the end". So data[-1] means the last element.
More can be found on slicing at http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html
We can do a very simple reference pixel correction using some very basic python and numpy. The reference pixel correction is done by calculating the average signal in the reference pixels for a given amplifier and subtracting this signal from each pixel that is read through this amplifier.
First we need to identify where the reference pixels are. On our NRS1 dark exposure, the reference pixels are on the left and right hand edges, in columns 1-4 and columns 2045-2048. To make this simple, we'll only use the reference pixels on the left hand edge, in columns 1-4. Also, we'll only consider amplifier A, which encompasses rows 1-512.
We specify this region using the following numpy index notation:
ref = newcube[:,:,0:512,0:4]
ref.shape
(1, 88, 512, 4)
Now we need to take the mean of the reference pixels in each plane. We can do this by looping over all the groups and putting the result of each group into a 1-d array of mean reference pixel values. First off, we have to create an empty array of the right size:
meanref = np.zeros(88, dtype=np.float32)
Now we need to loop over the groups. We use the 'range' builtin function in Python to give us a list of all the integers we need, rather like IDL's 'indgen':
for i in range(88):
meanref[i] = newcube[0,i,0:512,0:4].mean()
print meanref
[ 0. 24.65722656 -24.51464844 -2.08642578 -20.03125 -21.64257812 -74.33496094 -94.28515625 -35.43457031 -63.31298828 -77.56152344 -8.43164062 -43.94628906 -37.08154297 -46.32080078 8.484375 -42.23730469 -17.22265625 -69.10791016 -19.24853516 -59.98193359 -40.62988281 -28.23632812 -109.15966797 -15.37841797 -2.18310547 -42.87646484 -40.96875 -77.06445312 -38.67578125 -31.37744141 -35.89941406 48.20117188 -48.14794922 -28.07128906 -66.86328125 -22.19824219 -3.52294922 36.96337891 -41.18994141 -29.60253906 -44.91992188 -57.73974609 -43.97216797 -19.13671875 -22.32470703 -72.12548828 -41.38818359 -0.48681641 -19.09277344 -7.5234375 -9.78564453 -1.14013672 10.33300781 49.76855469 31.59765625 -1.86621094 -9.17626953 52.17382812 21.48583984 43.95703125 38.25048828 66.29443359 83.828125 1.34179688 -11.31884766 36.12792969 -2.66845703 34.83886719 -21.453125 -7.23632812 43.56640625 48.33105469 29.92919922 14.27783203 13.31494141 0.90380859 63.53857422 -5.58007812 77.22021484 15.73339844 21.8125 17.17529297 28.68212891 75.78710938 7.51757812 16.97949219 27.18505859]
Now we can overplot the signal from pixel (200, 100) and the reference pixel signal:
clf()
plot(newcube[0,:,99,199])
plot(meanref)
[<matplotlib.lines.Line2D at 0x7708d50>]
Pretty close, aren't they? Now lets plot the difference:
plot(newcube[0,:,99,199])
plot(meanref)
plot(newcube[0,:,99,199]-meanref)
[<matplotlib.lines.Line2D at 0x78da090>]
PyFITS allows you to do a lot more than just reading the pixel the data. You can examine the headers:
h[0].header
SIMPLE = T / Data conform to FITS standard BITPIX = 8 / Bits per data value NAXIS = 0 / Number of data array dimensions EXTEND = T / File may contain standard extensions NEXTEND = 3 / Number of standard extensions TELESCOP= 'JWST ' / Telescope used to acquire data INSTRUME= 'NIRSPEC ' / Identifier for instrument used to acquire data RADESYS = 'ICRS ' / Reference frame of equatorial or ecliptic coord DATE = '2012-11-07T15:10:05' / [yyyy-mm-ddThh:mm:ss.ssssss] UTC date of file ORIGIN = 'STScI ' / Institution responsible for creating FITS file FILENAME= 'jw00020001001_01101_00001_NRS1_uncal.fits' / Name of file FILETYPE= 'UNCALIBRATED' / Type of data found in file DPSW_VER= '0.1 ' / Data processing software version number Observation identifiers DATE-OBS= '2011-02-28T21:30:03.583' / [yyyy-mm-ddThh:mm:ss.sss] UTC date at star OBS_ID = 'FILL-DET-DARK-SHORT-01' / full programmatic observation identifier VISIT_ID= '' / visit identifier PROGRAM = '00020 ' / program number OBSERVTN= '001 ' / observation number VISIT = '001 ' / visit number VISITGRP= '01 ' / visit group identifier SEQ_ID = '1 ' / parallel sequence identifier ACT_ID = '01 ' / activity identifier EXPOSURE= '00001 ' / exposure request number within a single activit Exposure parameters DETECTOR= 'NRS1 ' / name of detector used to acquire data NINTS = 1 / number of integrations within exposure NGROUPS = 88 / number of groups within integration ZEROFRAM= F / T if a zero frame was read separately READPATT= 'NRSRAPID' / readout pattern DATAPROB= 0 / T if science telemetry indicated any problems w Program information TITLE = 'N/A ' / proposal title PI_NAME = 'UNKNOWN ' / name of principal investigator (last, first, mi CATEGORY= 'N/A ' / program category SUBCAT = 'N/A ' / program sub-category SCICAT = 'N/A ' / science category assigned during the TAC proces CONT_ID = 'N/A ' / continuation of the specified previous program Observation information TEMPLATE= 'N/A ' / proposal instruction template used to specify o OBSLABEL= '#TODO ' / proposer label for observation Visit information VISITYPE= '#TODO ' / type of visit (prime or parallel) VSTSTART= '#TODO ' / UTC visit start time WFSVISIT= '#TODO ' / wavefront sensing and control visit indicator NEXPOSUR= '#TODO ' / total number of exposures in visit INTARGET= '#TODO ' / T if at least one exposure in visit is internal EXTARGET= '#TODO ' / T if at least one exposure in visit is external TARGOOPP= '#TODO ' / visit scheduled as target of opportunity (T/F) Exposure information PNTG_SEQ= '#TODO ' / pointing sequence number EXPCOUNT= '#TODO ' / count of the exposures in visit EXP_TYPE= '#TODO ' / type of data in exposure Target information TARGPROP= 'UNKNOWN ' / proposer's name for the target TARGNAME= 'N/A ' / standard astronomical catalog name for the targ TARGTYPE= 'N/A ' / fixed target, moving target, or generic target TARG_RA = 0.0 / target RA computed at time of exposure TARGURA = 'N/A ' / target RA uncertainty computed at time of expos TARG_DEC= 0.0 / target DEC computed at time of exposure TARRUDEC= 'N/A ' / target Dec uncertainty computed at time of expo PROP_RA = 'N/A ' / proposer specified RA for the target PROP_DEC= 'N/A ' / proposer specified Dec for the target PROPEPOC= 'N/A ' / proposer specified epoch for RA and Dec Exposure times EXPSTART= '2011-02-28T21:30:03.583' / [Modified Julian Date] UTC exposure start EXPMID = '#TODO ' / [Modified Julian Date] UTC exposure mid time EXPEND = '2011-02-28T21:45:48.417' / [Modified Julian Date] UTC exposure end ti Exposure time parameters NSAMPLES= 1 / number of A/D samples per pixel NFRAMES = 1 / number of frames in group GROUPGAP= 0 / number of frames dropped between groups TSAMPLE = 10 / delta time between samples in microseconds NRESET = '#TODO ' / number of resets between integrations in an exp NRSTSTRT= '#TODO ' / number of extra resets at start of exposure TFRAME = 10.73676 / [seconds] time between frames TGROUP = 10.73676 / [seconds] time between groups EFFINTTM= 934.09812 / [seconds] effective integration time EFFEXPTM= '' / [seconds] effective exposure time CHRGTIME= '#TODO ' / [seconds] charge accumulation time DURATION= '#TODO ' / [seconds] total duration of exposure subarray parameters SUBARRAY= F / name of subarray used SUBXSTRT= 1 / starting pixel number in the SIAS x direction SUBXSIZE= 2048 / number of pixels in the SIAS x direction SUBYSTRT= 1 / starting pixel number in the SIAS y direction SUBYSIZE= 2048 / number of pixels in the SIAS y direction NXLIGHT = '#TODO ' / number of light sensitive x values (columns) NIRSpec configuration (NIRSpec only) FILTER = 'OPAQUE ' / name of the filter element used GRATING = 'G140H ' / name of the grating element used FXD_SLIT= '#TODO ' / name of fixed slit aperture used MSASTATE= '#TODO ' / state of MSA: all open, all closed, configured FOCUSPOS= '#TODO ' / [mm] focus position for NIRSpec NIRSpec MSA supporting files (NIRSpec MSA only) MSACONFG= 'N/A ' / MSA configuration file name lamp configuration LAMP = 'CLOSE ' / internal lamp state Guide star information GS_ORDER= 'N/A ' / index of guide star within list of selected gui GSSTRTTM= 'N/A ' / UTC start time of guide star acquisition activi GSENDTIM= 'N/A ' / UTC end time of guide star acquisition activity GDSTARID= 'N/A ' / guide star identifier GS_RA = 'N/A ' / guide star right ascension GS_DEC = 'N/A ' / guide star declination GSURA = 'N/A ' / guide star right ascension uncertainty GSUDEC = 'N/A ' / guide star declination uncertainty GS_MAG = 'N/A ' / guide star magnitude in FGS detector GSUMAG = 'N/A ' / guide star magnitude uncertainty PCS_MODE= 'N/A ' / Pointing Control System mode GSCENTX = 'N/A ' / guide star centroid x position in ideal coordin GSCENTY = 'N/A ' / guide star centroid y position in ideal coordin JITTERMS= 'N/A ' / [arcsec] RMS jitter over the exposure JWST ephemeris information COORDSYS= 'N/A ' / ephemeris coordinate system EPH_TIME= 'N/A ' / [sec] UTC time from ephemeris start time JWST_X = 'N/A ' / [km] X spatial coordinate of JWST JWST_Y = 'N/A ' / [km] Y spatial coordinate of JWST JWST_Z = 'N/A ' / [km] Z spatial coordinate of JWST JWST_DX = 'N/A ' / [km/sec] X component of JWST velocity vector JWST_DY = 'N/A ' / [km/sec] Y component of JWST velocity vector JWST_DZ = 'N/A ' / [km/sec] Z component of JWST velocity vector Spacecraft pointing information PA_V3 = 'N/A ' / [deg] position angle of V3-axis of JWST RA_V1 = 'N/A ' / [deg] RA of telescope V1 axis DEC_V1 = 'N/A ' / [deg] Dec of telescope V1 axis Aperture pointing information APERNAME= '#TODO ' / mnemonic for PDB science aperture used PA_APER = '#TODO ' / [deg] position angle of aperture used with targ WCS parameters WCSAXES = 2 / number of World Coordinate System axes CRPIX1 = 'N/A ' / x-coordinate of the reference pixel CRPIX2 = 'N/A ' / y-coordinate of the reference pixel CRVAL1 = 'N/A ' / [degrees] right ascension value at the referenc CRVAL2 = 'N/A ' / [degrees] declination value at the reference pi CTYPE1 = 'RA---TAN' / first axis coordinate type CTYPE2 = 'DEC--TAN' / second axis coordinate type CDELT1 = 'N/A ' / [degrees] increment per pixel, increasing eastw CDELT2 = 'N/A ' / [degrees] increment per pixel, increasing north PC1_1 = 'N/A ' / linear transformation matrix element cos(theta) PC1_2 = 'N/A ' / linear transformation matrix element -sin(theta PC2_1 = 'N/A ' / linear transformation matrix element sin(theta) PC2_2 = 'N/A ' / linear transformation matrix element cos(theta) S_REGION= 'N/A ' / spatial extent of the observation Velocity aberration correction DVA_RA = 'N/A ' / [radians] velocity aberration correction RA off DVA_DEC = 'N/A ' / [radians] velocity aberration correction Dec of VA_SCALE= 'N/A ' / velocity aberration scale factor Time related keywords BARTDELT= 'N/A ' / [sec] calculated Barycentric time correction fr BSTRTIME= 'N/A ' / [days] Solar System Barycentric exposure start BENDTIME= 'N/A ' / [days] Solar System Barycentric exposure end ti BMIDTIME= 'N/A ' / [days] Solar System Barycentric exposure mid ti HELIDELT= 'N/A ' / [days] calculated Heliocentric time correction HSTRTIME= 'N/A ' / [days] Heliocentric exposure start time in MJD HENDTIME= 'N/A ' / [days] Heliocentric exposure end time in MJD HMIDTIME= 'N/A ' / [days] Heliocentric exposure mid time in MJD
You can change individual header entries:
h[0].header['HMIDTIME'] = 0.0
h[0].header['HMIDTIME']
0.0
But nothing changes until you write out the data
h.writeto('newversion.fits')
Now we have two FITS files, and the only difference between them is the value of the HMIDTIME keyword.
You can also change the data:
h[1].data = newcube
h.writeto('newversion.fits',clobber=True)
Overwriting existing file 'newversion.fits'.
ls
DQInit/ NRS1_dark88.fits dqinit.tar pyfits/ Exposure_Log_FM_Phase5_Latest_29Jul11.csv* Session2_RIJ.ipynb fits_generator/ test/ install_ureka_v1.0beta3* Ureka/ newversion.fits Ureka_linux-rhe5_64_dev.tar DQInit/ NRS1_dark88.fits dqinit.tar pyfits/ Exposure_Log_FM_Phase5_Latest_29Jul11.csv* Session2_RIJ.ipynb fits_generator/ test/ install_ureka_v1.0beta3* Ureka/ newversion.fits Ureka_linux-rhe5_64_dev.tar
If you display this version of the data using ds9 or similar, you'll be able to see much more than is visible in the version without the first group subtracted
You can display images into a ds9 window from within Python using the numdisplay package.
from stsci import numdisplay
numdisplay.display(newcube[0,50],zscale=True)
Image displayed with Z1: -110.0 Z2: 159.234070434
It's a bit like using the IRAF 'display' command, in that your image is scaled into the range 0-199 before being sent to ds9. Personally, I find it much more useful to write the files out and read them directly into ds9, that way you see all the planes of a cube and you can easily scroll through them.
You can work with text tables using the ascii module in the astropy package. This is the same as the asciitables package of Tom Aldcroft that some of you may have used.
Start by importing the package:
import sys
sys.path.insert(1,'/home/robert/py-lib')
from astropy.io import ascii
The ascii table reader can read a large number of different table formats:
First we need some suitable table data. I picked a large .csv file with information about MIRI FM test data.
ls
astropy-2012-12-05/ newversion.fits astropy-2012-12-05.tar NRS1_dark88.fits DQInit/ pyfits/ dqinit.tar Session2_RIJ.ipynb Exposure_Log_FM_Phase5_Latest_29Jul11.csv* test/ fits_generator/ Ureka/ install_ureka_v1.0beta3* Ureka_linux-rhe5_64_dev.tar
This file has a couple of lines of junk, followed by a line with the column names, followed by the data. If we tell the ascii.read function this, it will read the table nicely...
data = ascii.read('Exposure_Log_FM_Phase5_Latest_29Jul11.csv', header_start=2, data_start=3)
That seemed to work. What did we get?
data
OBSID | EXPID | Test Label | ... | PSS LED Current | Pupil X | Pupil Y |
---|---|---|---|---|---|---|
FM1T00010019 | 1 | MIRI_FT_FPE_IM_FULL_TP | ... | 0.0 | 0.0 | 0.0 |
FM1T00010020 | 1 | MIRI_FT_FPE_MRS_FULL_TP | ... | 0.0 | 0.0 | 0.0 |
FM1T00010021 | 1 | MIRI_FT_FPE_COMMS_CHECK | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 1 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 2 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 3 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 4 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 5 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 6 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010023 | 7 | FPS Grounding Test | ... | 0.0 | 0.0 | 0.0 |
FM1T00010024 | 1 | MIRI_FT_FPE_IM_AMB | ... | 0.0 | 0.0 | 0.0 |
... | ... | ... | ... | ... | ... | ... |
FM1T00014002 | 2 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014002 | 3 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014002 | 4 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014002 | 5 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014002 | 6 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014003 | 1 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014003 | 2 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014003 | 3 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014003 | 4 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014003 | 5 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
FM1T00014003 | 6 | IMG_OPT_05_POL | ... | 0.0 | 0.0 | 0.0 |
We clearly got some complex object, but what is it?
type(data)
astropy.table.table.Table
The Table class is part of astropy, and its documentation can be found in the data tables section of astropy's documentation:
http://astropy.readthedocs.org/en/v0.1/table/index.html
We can find out what we can do with this object by tabbing
data.
File "<ipython-input-6-1adda7b3306c>", line 1 data. ^ SyntaxError: invalid syntax
data.colnames
['OBSID', 'EXPID', 'Test Label', 'Exp Start Time', 'FPE Side', 'Detectors', 'Nframes', 'Nints', 'SubArray', 'NColumns', 'NRows', 'Read Mode', 'Det Mode', 'Start Column', 'Start Row', 'Detector Temp', 'Frame_Resets', 'Row_Resets', 'ADC_Delay', 'RPC_Delay', 'Row_Underlap', 'Col_Underlap', 'VDETCOM', 'VDDUC', 'VP Bias', 'VSSOUT', 'FPE Temp', 'ICE Side', 'CCC Pos', 'FWA Pos', 'DGAA Pos', 'DGAB Pos', 'IMG Cal I', 'IMG Cal V', 'IMG Cal Mode', 'MRS Cal I', 'MRS Cal V', 'MRS Cal Mode', 'MTS FW', 'VAS %', 'BB Temp', 'FM4 Temp', 'Source', 'Pinhole', 'FOV X', 'FOV Y', 'Imager X', 'Imager Y', 'SS Focus', 'Focus Offset', 'MRS Alpha', 'MRS Slice', 'MRS Subchan', 'PSS LED Current', 'Pupil X', 'Pupil Y']
data.pprint()
OBSID EXPID Test Label ... PSS LED Current Pupil X Pupil Y ------------ ----- ----------------------- ... --------------- ------- ------- FM1T00010019 1 MIRI_FT_FPE_IM_FULL_TP ... 0.0 0.0 0.0 FM1T00010020 1 MIRI_FT_FPE_MRS_FULL_TP ... 0.0 0.0 0.0 FM1T00010021 1 MIRI_FT_FPE_COMMS_CHECK ... 0.0 0.0 0.0 FM1T00010023 1 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010023 2 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010023 3 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010023 4 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010023 5 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010023 6 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010023 7 FPS Grounding Test ... 0.0 0.0 0.0 FM1T00010024 1 MIRI_FT_FPE_IM_AMB ... 0.0 0.0 0.0 ... ... ... ... ... ... ... FM1T00014002 2 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014002 3 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014002 4 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014002 5 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014002 6 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014003 1 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014003 2 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014003 3 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014003 4 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014003 5 IMG_OPT_05_POL ... 0.0 0.0 0.0 FM1T00014003 6 IMG_OPT_05_POL ... 0.0 0.0 0.0
One thing that isn't obvious from the list of methods is that we can use the column name to refer to the column, rather like a FITS header keyword
data['OBSID']
<Column name='OBSID' units=None format=None description=None> array(['FM1T00010019', 'FM1T00010020', 'FM1T00010021', ..., 'FM1T00014003', 'FM1T00014003', 'FM1T00014003'], dtype='|S12')
data['OBSID"] is another object, looks like a Column object. For some reason the IPython notebook doesn't let you tab the data['OBSID'] object to see what methods it has...
type(data['OBSID'])
astropy.table.table.Column
But if we assign this to a new variable, we can do it...
obsidcol = data['OBSID']
obsidcol.
It looks a lot like a numpy ndarray, but the actual array part is in the 'data' member.
We can use the numpy 'where' method to find all the exposures taken with the imaging array and with at least 5 integrations:
imaging = np.where(data['Detectors'].data == 'IMAGER')
imaging
(array([ 0, 10, 11, ..., 8323, 8324, 8325]),)
The 'where' function of numpy returns a tuple with an associative array as the first element. This array gives the indices (or row numbers) where the condition was found to be True, and can be used as an index for a numpy array.
len(imaging[0])
4454
desired = np.where(data['Nints'].data[imaging] >= 5)
desired
(array([ 128, 129, 130, ..., 4451, 4452, 4453]),)
len(desired[0])
2392
data['OBSID'].data[desired]
array(['FM1T00010235', 'FM1T00010235', 'FM1T00010235', ..., 'FM1T00012149', 'FM1T00012149', 'FM1T00012149'], dtype='|S12')
data['Detector Temp'].data
array([ 293.2877724 , 293.3529064 , 293.4533905 , ..., 6.69992188, 6.69992188, 6.70012893])
data['BB Temp'].data
array([ 294.7900085, 294.7839966, 294.776001 , ..., 800.0009766, 799.9990234, 800.0020142])
plot(data['Detector Temp'].data,data['FM4 Temp'].data,'+')
[<matplotlib.lines.Line2D at 0x8c44dd0>]
data[0].data
('FM1T00010019', 1, 'MIRI_FT_FPE_IM_FULL_TP', '2011-04-26:16:24:53', 'NOMINAL', 'IMAGER', 2, 3, 'FULL', 1032, 1024, 'FAST', 'TEST_PATTERN', 1, 1, 293.2877724, 0, 3, 260, 24, 12, 24, -0.000383, -0.000383, -0.000217, -0.002797, 295.6150632, 'NOMINAL', 'LOCKED', 'F770W', 'LONG', 'LONG', 0.0, 0.0, 'NA', 0.0, 0.0, 'NA', 'UNKNOWN', 0.0, 294.7900085, 292.2000122, 'UNKNOWN', 'NONE', 0.0, 0.0, -999.0, -999.0, 0.0, 0, -999.0, -999.0, 'NONE', 0.0, 0.0, 0.0)