Windows installer¶

An experimental installer for Microsoft Windows (together with other packages that the Kapteyn Package depends on) is also available. Currently only for Python 2.6 on 32-bit systems. http://www.astro.rug.nl/software/kapteyn_windows/

Scisoft problem¶

If you have Scisoft installed on your computer, it may interfere with the installation of the Kapteyn Package. To install it properly, disable the setup of Scisoft in your startup file (e.g. ~/.cshrc, .profile) by commenting it out.

Mac OS X Compiler problem¶

There is a known problem with Apple’s llvm-gcc-4.2 compiler. This compiler is known to crash with an internal compiler error (Segmentation fault: 11) when WCSLIB routine wcserr.c is compiled. For this reason, setup.py tries to detect this compiler and use the clang compiler instead. If compilation still fails, one could try to prefix a shell variable definition to the install command like this:

export CC=CLANG; python setup.py install ...

Wcs¶

• Given a FITS header or a Python dictionary with header information about a World Coordinate System (WCS), transform between pixel- and world coordinates.
• Different coordinate representations are possible (tuple of scalars, NumPy array etc.)
• Transformations between sky and reference systems.
• Epoch transformations
• Support for ‘alternate’ headers (a header can have more than one description of a WCS)
• Support for mixed coordinate transformations (i.e. pixel- and world coordinates at input are mixed).
• Spectral coordinate translations, e.g. convert a frequency axis to an optical velocity axis.

Celestial¶

• Coordinate transformations between sky and reference systems. Also available via module wcs
• Epoch transformations. Also available via module wcs
• Many utility functions e.g. to convert epochs, to parse strings that define sky- and reference systems, calculate Julian dates, precession angles etc.

Wcsgrat¶

• Most of the functionality in this module is provided via user friendly methods in module maputils.
• Calculate grid lines showing constant latitude as function of varying longitude or vice versa.
• Methods to set the properties of various plot elements like tick marks, tick labels and axis labels.
• Methods to calculate positions of labels inside a plot (e.g. for all sky plots).

Maputils¶

• Easy to combine with Matplotlib
• Convenience methods for methods of modules wcs, celestial, wcsgrat
• Overlays of different graticules (each representing a different sky system),
• Plots of data slices from a data set with more than two axes (e.g. a FITS file with channel maps from a radio interferometer observation)
• Plots with a spectral axis with a ‘spectral translation’ (e.g. Frequency to Radio velocity)
• Rulers with distances in world coordinates, corrected for projections.
• Plots for data that cover the entire sky (allsky plot)
• Mosaics of multiple images (e.g. HI channel maps)
• A simple movie loop program to view ‘channel’ maps.
• Interactive colormap selection and modification.

Positions¶

• Convert strings to positions in pixel- and world coordinates

Rulers¶

• Plot a straight line with markers at constant distance in world coordinates. Its functionality is available in module maputils

Shapes¶

• Advanced plotting with user interaction. A user defines a shape (polygon, ellipse, circle, rectangle, spline) in an image and the shape propagates (in world coordinates) to other images. A shape object keeps track of its area (in pixels) and the sum of the pixels within the shape. From these a flux can be calculated.

Tabarray¶

• Fast I/O for data in ASCII files on disk.

Mplutil¶

• Various advanced utilities for event handling in Matplotlib. Most of its functionality is used in module maputils.

Please cite the Kapteyn Package¶

If you have used the Kapteyn Package in the preparation of a publication, please cite. We need these citations to justify time and resources spent on the software.

You may cite it as follows (BibTeX format):

One coordinate axis¶

For experiments and debug sessions, module wcs allows for very simple and flexible input and output of coordinates. This module interfaces with Mark Calabretta’s WCSLIB and is, because of the flexible I/O, a valuable tool to test this well known library.

Main goal of module wcs is to enable transformations between pixel coordinates and world coordinates The pixel coordinates are defined by the FITS standard. The transformation is defined by meta data which are usually found in FITS headers. So it may be obvious that FITS files play an important role in the use of Module wcs.

However, FITS data processed by wcs can also be FITS keywords that are stored in a Python dictionary. This invites to experiment with WCSLIB even more because one can create a (minimal) FITS header from scratch. In an attempt to create the most simple use of wcs we started to write a minimal FITS header. It defines only one axis. The minimal requirement for FITS keywords are CTYPE, CRVAL, CRPIX and CDELT. A description of these keywords can be found in The FITS standard.

We entered an axis type in CTYPE1 that WCSLIB does not recognize as a known type. With this trick we force the system to do a linear transformation. It shows that you have to be careful with values for CTYPE because you will not be warned when a CTYPE is not recognized.

For the conversions between pixel coordinates and world coordinates we defined methods in a class which we called the wcs.Projection class. An object of this class is created using the header of the FITS file for which we want WCS transformations. It accepts also a user defined Python dictionary with FITS keywords and values. We use this last option in this tutorial to be more flexible when we want to apply changes to the header.

The methods for single axes are called wcs.Projection.toworld1d() and wcs.Projection.topixel1d(). FITS defines CRVAL as the world coordinate that corresponds to the pixel value in CRPIX. Let’s check this with the most basic example we could think of:

#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS'  : 1,
           'CTYPE1' : 'PARAM',
           'CRVAL1' : 5,
           'CRPIX1' : 10,
           'CDELT1' : 1
         }
proj = wcs.Projection(header)
print proj.toworld1d(10)

# Output:
# 5.0

Indeed, at pixel coordinate 10 (=CRPIX), the world coordinate is 5 (=CRVAL). If we want to know which pixel coordinate corresponds to world coordinate 5, then we use proj.topixel1d(5) to get the answer (which is the value of CRPIX: 10). Note that we forced the system to apply linear transformations only.

In many of the examples that we present in this tutorial we included a so called closure test. This is a test which uses the result of a transformation to test the inverse transformation which should result into the original value. Sometimes the result is not exactly what you expect because we work with a limited number precision. A simple closure test is:

proj = wcs.Projection(header)
w = proj.toworld1d(10)
p = proj.topixel1d(w)
print "CRPIX: ", p

# Output:
# CRPIX:  10.0

Coordinate transformations are often done in bulk, so of course the transformation methods accept more than one coordinate to convert. They can be represented as a Python list, a Python tuple or a NumPy array. The representation of the output is the same as that of the input coordinates. The output of the next statements therefore is not a surprise:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

#!/usr/bin/env python
from kapteyn import wcs
import numpy

header = { 'NAXIS'  : 1,
           'CTYPE1' : 'PARAM',
           'CRVAL1' : 5,
           'CRPIX1' : 10,
           'CDELT1' : 1
         }

proj = wcs.Projection(header)

w1 = proj.toworld1d( range(9,12) )
w2 = proj.toworld1d( [9,10,11] )
w3 = proj.toworld1d( (9,10,11) )
w4 = proj.toworld1d( numpy.array([9,10,11]) )
print w1, type(w1)
print w2, type(w2)
print w3, type(w3)
print w4, type(w4)
closure = proj.topixel1d(w4)    # Closure test
print closure, type(closure)

# Output:
# [4.0, 5.0, 6.0] <type 'list'>
# [4.0, 5.0, 6.0] <type 'list'>
# (4.0, 5.0, 6.0) <type 'tuple'>
# [ 4.  5.  6.] <type 'numpy.ndarray'>
# [ 9.  10. 11.] <type 'numpy.ndarray'>

The first two sequences are lists. The third is a tuple and the last is a NumPy array. The pixel coordinates 9, 10 and 11 should give values in the neighbourhood of CRVAL1 and the step size is 1 (CDELT1=1), in arbitrary units.

For representation purposes we often want to print a pixel coordinate and the corresponding world coordinate on one line. Then we often use Pythons built-in function zip to combine two sequences to avoid a call to transformation methods in the print loop:

p = range(5,15)
w = proj.toworld1d(p)
for pix,wor in zip(p,w):
   print "%d: %f" % (pix,wor)

# Output:
# 9: 4.000000
# 10: 5.000000
# 11: 6.000000

wcs.debug = True
proj = wcs.Projection(header)

Next we apply the procedures described above to a real example where we created an artificial header with FITS data. The header describes a single axis of spectral type. Units are standard FITS units and are given in keyword CUNIT1. The example shows that we can access the keywords from the artificial header (or a real FITS header) directly and use their values for example to find the length of the axis in pixels, or to find the units of the world coordinates of that axis:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

#!/usr/bin/env python
from kapteyn import wcs
header  = { 'NAXIS'  : 1,
            'NAXIS1' : 64,
            'CTYPE1' : 'FREQ',
            'CRVAL1' : 1.37835117405e9,
            'CRPIX1' : 32,
            'CUNIT1' : 'Hz',
            'CDELT1' : 9.765625e4
         }
proj = wcs.Projection(header)
n = header['NAXIS1']               # Get the length of the spectral axis
p = range(1, n+1)                  # Set pixel range accordingly
w = proj.toworld1d(p)              # Do the transformation
print "Pixel  %s (%s)" % (header['CTYPE1'],header['CUNIT1'])
for pix,frq in zip(p,w):
   print "%5d: %f" % (pix,frq)

# Output:
# Pixel  FREQ (Hz)
#   1: 1375323830.300000
#   2: 1375421486.550000
#   3: 1375519142.800000
#   4: 1375616799.050000
#   5: 1375714455.300000

In the example we wanted to make a table with pixel coordinates and the corresponding world coordinates. According to the header there are 64 pixels (NAXIS1) along the axis so the first pixel coordinate is 1 and the last is 64. The axis represents frequencies. A start frequency is given by CRVAL1 and a step size is given by CDELT1. Note that the coordinate transformation is linear.

Generic methods toworld() and topixel()¶

The methods wcs.Projection.toworld1d() and wcs.Projection.topixel1d() are special versions of the more general methods wcs.Projection.toworld() and wcs.Projection.topixel(). These methods can be used to convert pixel data for more than one axis at the same time which is necessary for coupled axes, for example in spatial maps where longitude and latitude are not independent axes.

These general methods wcs.Projection.toworld() and wcs.Projection.topixel() accept the same sequences as the ‘1d’ versions. The reason that we introduced the ‘1d’ versions is that for non-experienced Python programmers it usually is confusing that in the one dimensional case the general methods only accept tuples and not scalars and that a tuple with one element (for example 10) needs to be written as (10,).

If you want to replace method toworld1d() by topixel1d() in the first example, then the relevant lines become:

>>> p = proj.toworld( (10,) )
>>> (5.0,)

for one scalar and for a list of values:

>>> p = proj.toworld( (range(9,12),) )
>>> ([4.0, 5.0, 6.0],)

If you want to extract the scalar or the list from the tuple, use element 0 of the tuple.

>>> p = proj.toworld( (range(9,12),) )
>>> print p[0]
>>> [4.0, 5.0, 6.0]

Two coordinate axes¶

As described in the previous section we use wcs.Projection.toworld() and wcs.Projection.topixel() if the number of axes in our data is more than 1. The input and output tuples for projection objects with two coordinate axes consist of two elements. The first element corresponds to the first axis in the projection object and the second element to the second axis. The following Python code constructs an artificial header which describes the world coordinate system of two spatial axes. Then we want to find the world coordinates of the reference pixels (CRPIX1, CRPIX2) and expect the reference values (CRVAL1, CRVAL2) as output tuple:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

#!/usr/bin/env python
from kapteyn import wcs
header  = { 'NAXIS'  : 2,
            'NAXIS1' : 5,
            'CTYPE1' : 'RA---NCP',
            'CRVAL1' : 45,
            'CRPIX1' : 5,
            'CUNIT1' : 'deg',
            'CDELT1' : -0.01,
            'NAXIS2' : 10,
            'CTYPE2' : 'DEC--NCP',
            'CRVAL2' : 30,
            'CRPIX2' : 5,
            'CUNIT2' : 'deg',
            'CDELT2' : +0.01,
         }
proj = wcs.Projection(header)
pixel = (5,5)
world = proj.toworld(pixel)
print world

# Output:
# (45.0, 30.0)

Comments about the composed header: the header is composed from scratch. but it could very well have been copied from an existing FITS header. In either case you should verify items CUNITn and CTYPEn because they are are important. In section 2.1.1 of [Ref1] we read that in WCSLIB:

The CTYPE consists of a coordinate type (max 4 characters) followed by ‘-‘ followed by a three character code that represents the algorithm to calculate the world coordinates (‘ABCD-XYZ’). Shorter coordinate types are padded with the ‘-‘ character, shorter algorithm codes are padded on the right with blanks (‘RA—NCP’, ‘RA—UV_ ‘). So if we were sloppy and wrote RA–NCP and DEC-NCP then WCSLIB assigns a linear conversion algorithm. It does not complain, but you get unexpected results. If your CTYPEs are correct but the units are not standard and are not recognized by WCSLIB, then you get an Python exception after you try to create the Projection object. For example, if you specified CUNIT1=’Degree’ then the error message displayed by the exception is: “Invalid coordinate transformation parameters”.

If you want to be sure that WCSLIB recognizes your coordinate type and unit, you can print the Projection attributes wcs.Projection.types and wcs.Projection.units as in the example below. Unrecognized types are returned as None.

>>> proj = wcs.Projection(header)
>>> print "WCS units: ",proj.units
    WCS units:  ('deg', 'deg')
>>> print "WCS type: ",proj.types
    WCS type:  ('longitude', 'latitude')

With the same variable header as in the previous script we demonstrate that each element in the coordinate tuple can be a list of scalars. Let’s convert pixel positions (3,3), (4,4), ..., (7,7) etc. to their corresponding world coordinates:

proj = wcs.Projection(header)
x = range(3,8)
y = range(3,8)
pixel = (x,y)
world = proj.toworld(pixel)
print world

# Output:
# ([45.023089356221305, 45.011545841750113, 45.0, 44.988451831142257, 44.97690133535837],
#  [29.979985885372404, 29.989996472289789, 30.0, 30.009996474046854, 30.019985899953429])

The output is a tuple with two elements. Each element is a list. The first list contains the longitude coordinates for input pixel coordinates (3,3), (4,4) etc. The second list contains the latitude coordinates for the input pixel coordinates (3,3), (4,4) etc.

Here input and output coordinates for the methods wcs.Projection.toworld() and wcs.Projection.topixel() are tuples. The dimension of the tuple corresponds to the number of axes in the Projection object, and each element in the tuple can be a list of scalars. In some situations it is more intuitive to start with a list of 2 dimensional positions. The wcs module allows for this type of input. You can get the same coordinate output as the previous script if you replace the body by:

proj = wcs.Projection(header)
pixels = [(3,3), (4,4), (5,5), (6,6), (7,7)]
world = proj.toworld(pixels)
print world

# Output:
# [(45.023089356221305, 29.979985885372404), (45.011545841750113, 29.989996472289789), (45.0, 30.0),
# (44.988451831142257, 30.009996474046854), (44.97690133535837, 30.019985899953429)]

Note that the representation of the output differs from the previous script because the representation of the input differs, i.e.: a list with tuples. The dimension of the tuples being the number of axes in your projection object.

Mixed transformations (pixel- and world coordinates) using method wcs.Projection.mixed()¶

We describe the mixed() method in some detail in the section about data sets with three or more axes. Here we show how to use the method in a simple case. Suppose you want to mark data in a plot at constant declination in pixels (i.e. parallel to the x-axis of the plot) but with equal steps in Right Ascension, then you need method wcs.Projection.mixed():

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#!/usr/bin/env python
from kapteyn import wcs
import numpy
header  = { 'NAXIS'  : 2,
            'NAXIS1' : 5,
            'CTYPE1' : 'RA---TAN',
            'CRVAL1' : 45,
            'CRPIX1' : 5,
            'CUNIT1' : 'deg',
            'CDELT1' : -0.01,
            'NAXIS2' : 10,
            'CTYPE2' : 'DEC--TAN',
            'CRVAL2' : 30,
            'CRPIX2' : 10,
            'CUNIT2' : 'deg',
            'CDELT2' : +0.01,
         }
proj = wcs.Projection(header)
# 1 pixel and 1 world coordinate pair
pixel_in = (numpy.nan, 10)
world_in = (45.0, numpy.nan)
world_out, pixel_out = proj.mixed(world_in, pixel_in)
print world_out
print pixel_out

# Output:
# (45.0, 30.0)
# (5.0, 10.0)

# A loop over a number of Right Ascensions at constant Declination
for ra in range(44, 47):
   world_in = (ra,numpy.nan)
   world_out, pixel_out = proj.mixed(world_in, pixel_in)
   print "World: ", world_out, "Pixel: ", pixel_out

# Output:
# World:  (44.0, 29.99622120337045) Pixel:  (91.61133499750801, 10.000000000096229)
# World:  (45.0, 30.0) Pixel:  (5.0, 10.0)
# World:  (46.0, 29.99622120337045) Pixel:  (-81.61133499750801, 10.000000000096248)

First we have a pixel position of which the x coordinate is set to unknown. We use a special value for this: numpy.nan which is the representation of NumPy’s Not A Number. The y coordinate is set to 10. For the wcs.Projection.mixed(), we need to specify the unknown values in the pixel position with a world coordinate. In the example we entered 45.0 (deg). The mixed() method returns two tuples. One for the pixel position and one for the position in world coordinates. The unknown values are calculated in an iterative process. The second part of the example is a loop over a number of world coordinates in Right Ascension, and a constant pixel coordinate in the y-direction (i.e. 10). The output (as listed as comment in the code) shows two things that need to be addressed. First we notice that the output pixel is not exactly 10. This is related to finite precision of numbers when a solution is calculated in an iterative way. The second observation is more important: the Declination varies while the y coordinate in pixels is constant. But this is exactly what we expect for spatial data when a projection is involved.

A note about efficiency:

So in our example it is more efficient to avoid the loop over the right ascensions. This can be done by creating an input tuple with two lists. The output is the same as in the example above, but the representation is different. As we stated earlier, the representation of the output is the same as the representation of the input (a tuple with two lists):

# As example above but without a loop
ra = range(44, 47)
dec = [numpy.nan]*len(ra)  # NumPy trick to repeat elements in a list.
world_in = (ra, dec)
x = [numpy.nan]*len(ra)
y = [10]*len(ra)
pixel_in = (x, y)
world_out, pixel_out = proj.mixed(world_in, pixel_in)
print world_out
print pixel_out

# Output:
# ([44.0, 45.0, 46.0], [29.99622120337045, 30.0, 29.99622120337045])
# ([91.61133499750801, 5.0, -81.61133499750801], [10.000000000096229, 10.0, 10.000000000096248])

Three or more coordinate axes¶

In this section we discuss method wcs.Projection.sub() which allows us to define coordinate transformations for positions with less dimensions than the dimension of the data structure. In practice we encounter many astronomical measurements based on three or more independent axes. Well known examples are of course the data sets from radio interferometers. Usually these are spatial maps observed at different frequencies and sometimes as function of Stokes parameters (polarization). If we are only interested in spatial maps and don’t bother about the other axes, we can create a Projection object with only the relevant axes. This is done with the wcs.Projection.sub() method from the Projection class.

map = proj.sub(axes=None, nsub=None)

The method has two parameters. You can specify parameter nsub which sets the first nsub axes from the original Projection object to the actual axes. Or you can use the other parameter axes which is a tuple or a list with axis numbers. Axis numbers in WCSLIB follow the FITS standard so they start with 1. The order in the sequence is important. The axis description sequence in a FITS file is not bound to rules and luckily WCSLIB accepts permuted axis number sequences. This can be illustrated with the next example. First we show the code and then explain the output:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

#!/usr/bin/env python
from kapteyn import wcs
import numpy
header  = { 'NAXIS'  : 3,
            # First spatial axis
            'NAXIS1' : 5,
            'CTYPE1' : 'RA---TAN',
            'CRVAL1' : 45,
            'CRPIX1' : 5,
            'CUNIT1' : 'deg',
            'CDELT1' : -0.01,
            # A dummy axis
            'NAXIS2' : 5,
            'CTYPE2' : 'PARAM',
            'CRVAL2' : 444,
            'CRPIX2' : 99,
            'CDELT2' : 1.0,
            'CUNIT2' : 'wprf',
            # Second spatial axis
            'NAXIS3' : 0,
            'CTYPE3' : 'DEC--TAN',
            'CRVAL3' : 30,
            'CRPIX3' : 10,
            'CUNIT3' : 'deg',
            'CDELT3' : +0.01
         }
proj = wcs.Projection(header)
map = proj.sub( [1,3] )
pixel = (header['CRPIX1'], header['CRPIX3'])
world = map.toworld(pixel)
print world

# Output:
# (45.0, 30.0)

map = proj.sub( [3,1] )
pixel = (header['CRPIX3'], header['CRPIX1'])
world = map.toworld(pixel)
print world

# Output:
# (30.0, 45.0)

line  = proj.sub( 2 )
crpix = header['CRPIX2']
pixels = range(crpix-5,crpix+6)
world = line.toworld1d(pixels)
print world

# Output:
# [439.0, 440.0, 441.0, 442.0, 443.0, 444.0, 445.0, 446.0, 447.0, 448.0, 449.0]

We created a header representing a spatial map as function of some parameter along the CTYPE2=’PARAM’ axis. This axis is not recognized by WCSLIB and a linear transformation is applied. Also special is that the spatial axes do not have conventional numbers. First we want to set up a transformation of pixel (x,y) to (R.A., Dec) for the pixel values in (CRPIX1, CRPIX3) -which should transform to (CRVAL1, CRVAL3)-. Then we reverse the spatial axis sequence to set up a transformation from (y,x) to (Dec, R.A.). Finally we want a transformation only for the PARAM axis. Its axis number is 2. With the output we show that for this axis indeed the transformation between pixels and world coordinates is a linear. transformation.

The axis sequence in the wcs.Projection.sub() method sets the axis order with parameter axes. It sets in fact the order of the coordinates in the transformation methods wcs.Projection.toworld(), wcs.Projection.topixel() and wcs.Projection.mixed(). Parameter axes is either a single integer or a list/tuple of integers e.g. sub(2) vs. sub([3,1]).

NumPy matrices¶

In many Python applications programmers use NumPy arrays and matrices because it is easy to manipulate them. First let’s explore what can be done with a NumPy matrix as coordinate representation. A NumPy matrix is a rank 2 array with special properties. The first list in the numpy.matrix() constructor in the next example is the first row in the matrix and the second list is the second row. The first row contains the x coordinate of the pixels and the second row contains the y coordinates. In the next script we want to convert pixel positions (4,5), (5,5) and (6,5) to world coordinates. So the first list in the matrix constructor are the x coordinates [4,5,6] and the second are the y coordinates [5,5,5]. We convert these with:

proj = wcs.Projection(header)
pixel = numpy.matrix( [[4,5,6],[5,5,5]] )
world = proj.toworld(pixel)
print world
# Output:
# [[ 45.01154701  45.          44.98845299]
# [ 29.99999798  30.          29.99999798]]

pixel = proj.topixel(world)
print pixel

# Output:
# [[ 4.00000001  5.          5.99999999]
# [ 5.          5.          5.        ]]

The output is what we expected. It is a NumPy matrix with two rows. The first row contains the longitudes and the second the latitudes. The numbers seem ok (three RA’s at almost constant declination). We added a closure test by using the output world coordinates as input for the wcs.Projection.topixel() method. As you can see, the closure test returns the original input.

There is also a matrix representation that is equivalent to the list of coordinate tuples in the previous section. We want an input matrix to contain the coordinates: [[4,5],[5,5],[6,5]]. For this representation you have to set an attribute of the projection object. The name of the attribute is wcs.Projection.rowvec. Its default value is False. When you set it to True then each row in the matrix represents a position in x and y. Here is an example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

proj = wcs.Projection(header)
proj.rowvec = True
pixel = numpy.matrix( [[4,5],[5,5],[6,5]] )
world = proj.toworld(pixel)
print world

# Output:
# [[ 45.01154701  29.99999798]
# [ 45.          30.        ]
# [ 44.98845299  29.99999798]]

pixel = proj.topixel(world)
print pixel

# Output:
# [[ 4.00000001  5.        ]
# [ 5.          5.        ]
# [ 5.99999999  5.        ]]

NumPy arrays¶

It is possible to build a NumPy array with x coordinates and another for the y coordinates. You can use these arrays in a tuple. Then the elements in the tuple are not lists, as in the previous section, but NumPy arrays. With the same example in mind as the one with the NumPy matrix we demonstrate this option in the following script:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

proj = wcs.Projection(header)
x = numpy.array( [4,5,6] )
y = numpy.array( [5,5,5] )
pixel = (x, y)
world = proj.toworld(pixel)
print world

# Output:
# (array([ 45.01154701,  45. ,  44.98845299]), array([ 29.99999798,  30. ,  29.99999798]))

pixel = proj.topixel(world)
print pixel

# Output:
# (array([ 4.00000001,  5.        ,  5.99999999]), array([ 5.,  5.,  5.]))

As you can see, the representation of the output is the same as that of the input. The result is a tuple and the elements of the tuple are 1 dimensional (rank 1, shape N) NumPy arrays. The first array contains the RA’s and the second the Dec’s. The closure test also gives the expected result.

Using NumPy arrays to convert an entire map¶

For applications that transform all the positions in a data set (or in a subset of the data) in one run (e.g. for re-projections of images), it is possible to store all the positions in a NumPy array with shape (NAXIS2, NAXIS1, 2) (note the order). The array can be handled by the wcs.Projection.toworld() and wcs.Projection.topixel() in one step. You could say that we have a two-dimensional array of which the elements are coordinate pairs. The example code below could be part of the body of a real application that re-projects an image:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

from kapteyn import wcs
import numpy

header = {  'NAXIS'  : 2,
            'NAXIS1' : 5,
            'CTYPE1' : 'RA---TAN',
            'CRVAL1' : 45,
            'CRPIX1' : 5,
            'CUNIT1' : 'deg',
            'CDELT1' : -0.01,
            'NAXIS2' : 10,
            'CTYPE2' : 'DEC--TAN',
            'CRVAL2' : 30,
            'CRPIX2' : 10,
            'CUNIT2' : 'deg',
            'CDELT2' : +0.01,
         }

proj = wcs.Projection(header)
n1 = 10
n2 = 8
pixel = numpy.zeros(shape=(n2,n1,2))
for y in xrange(n2):
   for x in xrange(n1):
      pixel[y, x] = (x+1, y+1)

world = proj.toworld(pixel)
print world

# Output:
# [[[ 45.04614616  29.90999204]
#   [ 45.03460962  29.90999556]
#   [ 45.02307308  29.90999807]
#   [ 45.01153654  29.90999957]
# etc.

pixel = proj.topixel(world)
print pixel

# Output:
# [[[  1.   1.]
#   [  2.   1.]
#   [  3.   1.]
#   [  4.   1.]
# etc.

In this example we have NAXIS2=10 y values and NAXIS1=5 x values. The indices start at 0, but the FITS pixel indices start at 1. That’s why the coordinate tuple reads as (x+1, y+1).

Attributes lonaxnum, lataxnum and specaxnum¶

In the previous examples we had foreknowledge of the axis numbers that represented a spatial axis or a spectral axis. If you read a header from a FITS file then it is not always obvious what the axes represent and in which order they are stored in the FITS header. In those circumstances the projection attributes wcs.Projection.lonaxnum, wcs.Projection.lataxnum and wcs.Projection.specaxnum are very useful. These attributes are axis numbers, i.e. they start with 1 and the highest number is equal to header item ‘NAXIS’. In the source below we provide a header which shows an unexpected axis order representing a number of spatial maps as function of frequency. For demonstration purposes we create two separate Projection objects. The first, called line, represents the spectral axis. This is a sub projection of the parent projection object and the axis number is that of the spectral axis. We add a spectral translation to get velocities in the output.

The second, called map, is the spatial map with axis longitude first and latitude second. We try to create these objects in a try/except clause. For any header, this results in the requested sub projections for a spatial map and spectral axis or an error message and an exception. The construction with the attributes and the try/except clause saves us tedious work because without, we need to find and inspect the axis numbers ourselves.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS'  : 3,
           'NAXIS3' : 5,
           'CTYPE3' : 'RA---NCP',
           'CRVAL3' : 45,
           'CRPIX3' : 5,
           'CUNIT3' : 'deg',
           'CDELT3' : -0.01,
           'CTYPE2' : 'FREQ',
           'CRVAL2' : 1378471216.4292786,
           'CRPIX2' : 32,
           'CUNIT2' : 'Hz',
           'CDELT2' : 97647.745732,
           'RESTFRQ': 1.420405752e+9,
           'NAXIS1' : 10,
           'CTYPE1' : 'DEC--NCP',
           'CRVAL1' : 30,
           'CRPIX1' : 15,
           'CUNIT1' : 'deg',
           'CDELT1' : +0.01
         }
try:
   proj = wcs.Projection(header)
   line = proj.sub(proj.specaxnum).spectra('VRAD')
   map  = proj.sub( (proj.lonaxnum, proj.lataxnum) )
except:
   print "Could not find a spatial map AND a spectral line!"
   raise

print proj.lonaxnum, proj.lataxnum, proj.specaxnum

# Output:
# 3 1 2

# A transformation along the spectral axis:
pixels = range(30, 35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
   print p, v/1000

# Output:
# 30 8891.97019336
# 31 8871.36054878
# 32 8850.75090419
# 33 8830.14125961
# 34 8809.53161503

# A transformation of a coordinate in a spatial map:
ra  = header['CRVAL'+str(proj.lonaxnum)]
dec = header['CRVAL'+str(proj.lataxnum)]
print map.topixel( (ra,dec) )

# Output:
# (5.0, 15.0)

# Are these indeed the CRPIXn?
ax1 = "CRPIX"+str(proj.lonaxnum)
ax2 = "CRPIX"+str(proj.lataxnum)
print map.topixel( (ra,dec) ) == (header[ax1], header[ax2])

# Output:
# True

Note the check at the end of the code. It should return True (i.e. within some limited precision). We started with world coordinates equal to the values of CRVALn from the header and we assert that these correspond to pixel values equal to the corresponding CRPIXn.

Two dimensional data slices with only one spatial axis¶

Suppose we have a 3D data set with CTYPE’s: (RA—NCP, DEC–NCP, VOPT-F2W) and we want to write coordinate labels in a plot that represents the data as function of one spatial axis and the spectral axis (usually called a position-velocity plot or XV map)? It is obvious that we need extra information about the spatial axis that is left out. Usually this is a pixel position that corresponds to the position on the missing axis along which a data slice is taken. These data slices are fixed on pixel coordinates and not on world coordinates.

Assume the XV data we want to plot has axis types DEC–NCP and VOPT-F2W, then we need to specify at which pixel coordinate in Right Ascension the data is extracted.

What we need is a sub-projection (i.e. a Projection object which is modified by method sub()) which represents the WCSLIB types: (‘latitude’, ‘spectral’, ‘longitude’). Given the CTYPE’s from the header, the axis permutation sequence that is needed for the sub projection is (2,3,1). Now we require a method that for instance calculates for a given world coordinate in Declination (e.g. 60.1538880206 deg) and a velocity (e.g. -243000.0 m/s) and a fixed pixel for R.A. (e.g. 51) the corresponding pixel coordinates.

The required method is called wcs.Projection.mixed(). In a previous section we discussed its use. Method mixed() has for a Projection object p the following syntax and parameters.

world, pixel = p.mixed(world, pixel, span=None, step=0.0, iter=7)

It is a hybrid transformation suited for celestial coordinates. It uses an iterative method to find an unknown pixel- or world coordinate. The iteration is controlled by parameters span, step and iter. They have reasonable defaults which usually give good results. The method needs knowledge about elements that need to be solved. Unknown values that need to be solved are initially set to NaN (i.e. numpy.nan).

With the numbers we listed, the input world coordinate tuple will be world_in = (60.1538880206, -243000.0, numpy.nan). The input pixel tuple will be: pixel_in = (numpy.nan, numpy.nan, 51) then we find the missing coordinates after applying the lines:

subproj = proj.sub([2,3,1])
world_in = (60.1538880206, -243000.0, numpy.nan)
pixel_in = (numpy.nan, numpy.nan, 51)
world_out, pixel_out = subproj.mixed(world_in, pixel_in)
print "world_out = ", world_out
# world_out = (60.1538880206, -243000.0, -51.282084795900005)
print "pixel_out = ", pixel_out
# pixel_out = (51.0, -20.0, 51.0)

The mixed() method in wcs is more powerful than its equivalent in the C-version of WCSLIB. It accepts the same coordinate representations as for topixel() and toworld() whereas the library version accepts only one coordinate pair per call.

Suppress exceptions for invalid coordinates¶

We introduced matrices and arrays as coordinate representations to facilitate the input and output of many coordinates in one call. This is in many practical situations the most efficient way to process those coordinates. However if there is a pixel coordinate in a sequence that could not be converted to a world coordinate then an exception will be raised and your script will stop. One can suppress the exception and flag the unknown coordinate. You need to set the wcs.Projection.allow_invalid attribute of the projection object. Invalid coordinates then are flagged in the output with a NaN (i.e. numpy.nan). On the other hand, if the input contains a NaN, the corresponding converted coordinate will also be a NaN. You can test whether a value is a NaN with function numpy.isnan(). NaN’s cannot be compared so a simple test as in:

>>> x = numpy.nan
>>> if x == numpy.nan:         # ... fails

will fail because the result is always False The test x != x will give True if x is NaN.

In practice it will be difficult to get into problems if you convert from world coordinates to pixel coordinates, but when you start with pixel coordinates then it is possible that a corresponding world coordinate is not available. For a projection like Aitoff’s projection it is obvious that the rectangle in which an all sky map in this the projection is enclosed, contains such pixels.

Here is an example how one can deal with invalid transformations:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

#!/usr/bin/env python
from kapteyn import wcs
import numpy
header = { 'NAXIS'  : 2,
           'NAXIS1' : 5,
           'CTYPE1' : 'RA---AIT',
           'CRVAL1' : 45,
           'CRPIX1' : 5,
           'CUNIT1' : 'deg',
           'CDELT1' : -0.01,
           'NAXIS2' : 10,
           'CTYPE2' : 'DEC--AIT',
           'CRVAL2' : 30,
           'CRPIX2' : 5,
           'CUNIT2' : 'deg',
           'CDELT2' : +0.01,
         }
proj = wcs.Projection(header)
proj.allow_invalid = True
pixel_in = numpy.matrix( [[4000,5000,6000],[5000,5000,7580]] )
world = proj.toworld(pixel_in)
print "World coordinates:\n",world
pixel_out = proj.topixel(world)
print "Back to pixels:\n", pixel_out

if numpy.isnan(pixel_out).any():
   print "Some pixels could not be converted"

indices = numpy.where(numpy.isnan(pixel_out))
print "Index of NaNs: ", indices
print pixel_in[indices]

Reading a FITS header¶

Until now, we created our own header as a Python dictionary. But usually our starting point is a FITS file. A FITS file can contain more than one header. Header data is read from a FITS file with methods from module pyfits. Select the unit you want and store it in a variable (like header) so that it can be parsed by wcs. Below we demonstrate how to read the first header from a FITS file.

A flag is set to enter WCSLIB’s debug mode:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

#!/usr/bin/env python
from kapteyn import wcs
import pyfits

wcs.debug = True
f = raw_input('Enter name of FITS file: ')
hdulist = pyfits.open(f)
header = hdulist[0].header
proj = wcs.Projection(header)

# Part of the output of arbitrary FITS file:
# Output:
#      flag: 137
#      naxis: 3
#      crpix: 0x99b53d8
#               51           51          -20
#         pc: 0x99adf10
#   pc[0][]:   1            0            0
#   pc[1][]:   0            1            0
#   pc[2][]:   0            0            1
#      cdelt: 0x99b71c8
#            -0.007166     0.007166     4200
#      crval: 0x992bd30
#            -51.282       60.154      -2.43e+05
#      cunit: 0x99ad768
#            "deg"
#            "deg"
#            "m/s"
#      ctype: 0x999a7f8
#            "RA---SIN"
#            "DEC--SIN"
#            "VELO"

For testing and debugging one often wants to inspect the items in a FITS header. PyFITS has a nice method to make a list with all the FITS cards. In the next example we added a little filter, using list comprehension, to filter all items that start with ‘HISTORY’. Also we added output for the two projection attributes wcs.Projection.types and wcs.Projection.units. The script is a useful tool to inspect the FITS file and to check its parsing by module wcs:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

#!/usr/bin/env python
from kapteyn import wcs
import pyfits

f = raw_input('Enter name of FITS file: ')
hdulist = pyfits.open(f)
header = hdulist[0].header
clist = header.ascardlist()
c2 = [str(k) for k in header.ascardlist() if not str(k).startswith('HISTORY')]
for i in c2:
   print i

proj = wcs.Projection(header)
print "WCS found types: ", proj.types
print "WCS found units: ", proj.units

Reading WCS data for a spatial map¶

For some world coordinate related applications we want to force the input to represent a spatial map. A spatial map has axes of type longitude and latitude. For example if you need to re-project a map from one projection system to another, then you need a matching axis pair, representing a spatial system. If you don’t know beforehand what the numbers are of the axes in your FITS file that represent these types, you need a way of checking this. There are some rules. First, we must be able to create a Projection object according to the WCSLIB rules (i.e. the axes must have a valid name and extension). For spatial axes, WCSLIB also requires a matching axis pair. So if you have a FITS file with a R.A. axis and not a Dec axis then module wcs will generate an exception with the message Inconsistent or unrecognized coordinate axis types.

Finally, if you have a valid header and made a Projection object, then you still have to find the axis numbers that represent a ‘longitude’ axis and a ‘latitude’ axis (remember: the number of axes in your data could be more than 2) and the latitude axis could be defined earlier than the longitude axis so the order is also important.

In a previous section we discussed the attributes wcs.Projection.lonaxnum and wcs.Projection.lataxnum. They can be used to find the requested spatial axis numbers (remember their value is None if the requested axis is not available). In the following script we try to create the Projection and sub Projection objects with Python’s try/except mechanism.

If we have a valid projection and the right axes, then we check the axes types (and order) with attribute wcs.Projection.types:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

#!/usr/bin/env python
from kapteyn import wcs
import pyfits

f = raw_input('Enter name of FITS file: ')
hdulist = pyfits.open(f)
header = hdulist[0].header
try:
   proj = wcs.Projection(header)
   map = proj.sub((proj.lonaxnum, proj.lataxnum))
except:
   print "Aborting program. Could not find (valid) spatial map."
   raise

# Just a check:
print map.types

Celestial systems¶

Methods wcs.Projection.toworld() and wcs.Projection.topixel() convert between pixel coordinates and world coordinates. If these world coordinates are spatial, they are calculated for the sky- and reference system as defined in the header (FITS header, GIPSY header, header dictionary). To compare positions one must therefore ensure that these positions are all defined in the same sky- and reference system. If such a position is given in another system (e.g. galactic instead of equatorial), then you have to transform the position to the other sky- and/or reference system. Sometimes you might find a so called alternate header in the header information of a FITS file. In an alternate header the WCS related keywords end on a letter A-Z (e.g. CRVAL1A).

Usually these alternate headers describe a world coordinate system for another sky system. But because there could also be different epochs involved, it is worthwhile to have a system that can transform world coordinates between sky- and reference systems and that can do epoch transformations as well.

For the Kapteyn Package we wrote module celestial. This module can be used as stand alone module if one is interested in celestial transformations of world coordinates only. But the module is well integrated in module wcs so one can use it in the context of wcs, with the class wcs.Transformation. for conversions of world coordinates between sky-/reference systems and also, if pixel coordinates are involved, methods wcs.Projection.toworld() and wcs.Projection.topixel() can interpret an alternative sky-/reference system as the system for which a coordinate has to be calculated. The alternative sky-/reference system is stored in attribute wcs.projection.skyout.

The celestial definitions are described in detail in the background information of module celestial. We list the most important features of a celestial definition:

Supported Sky systems (detailed information in Sky systems):

1. Equatorial: Equatorial coordinates (α, δ), see next list with reference systems
2. Ecliptic: Ecliptic coordinates (λ, β) referred to the ecliptic and mean equinox
3. Galactic: Galactic coordinates (lII, bII)
4. Supergalactic: De Vaucouleurs Supergalactic coordinates (sgl, sgb)

Supported Reference systems (detailed information in Reference systems):

1. FK4: Mean place pre-IAU 1976 system.
2. FK4_NO_E: The old FK4 (barycentric) equatorial system but without the “E-terms of aberration”
3. FK5: Mean place post IAU 1976 system.
4. ICRS: The International Celestial Reference System.
5. J2000: This is an equatorial coordinate system based on the mean dynamical equator and equinox at epoch J2000.

Epochs (detailed information in Epochs for the equinox and epoch of observation):

The equinox and epoch of observations are instants of time and are of type string. These strings are parsed by a function of module celestial called celestial.epochs(). The parser rules are described in the documentation for that function. Each string starts with a prefix. Supported prefixes are:

1. B: Besselian epoch
2. J: Julian epoch
3. JD: Julian date
4. MJD: Modified Julian Day
5. RJD: Reduced Julian Day
6. F: Old and new FITS format (old: DD/MM/YY new: YYYY-MM-DD or YYYY-MM-DDTHH:MM:SS)

Example: Next example is a simple test program for epoch specifications. The function celestial.epochs() returns a tuple with three elements:

• the Besselian epoch
• the Julian epoch
• the Julian date.

#!/usr/bin/env python
from kapteyn import wcs

ep = ['J2000', 'j2000', 'j 2000.5', 'B 2000', 'JD2450123.7',
      'mJD 24034', 'MJD50123.2', 'rJD50123.2', 'Rjd 23433',
      'F29/11/57', 'F2000-01-01', 'F2002-04-04T09:42:42.1']

for epoch in ep:
   B, J, JD = wcs.epochs(epoch)
   print "%24s = B%f, J%f, JD %f" % (epoch, B, J, JD)

The output is:

#                  J2000 = B2000.001278, J2000.000000, JD 2451545.000000
#                  j2000 = B2000.001278, J2000.000000, JD 2451545.000000
#               j 2000.5 = B2000.501288, J2000.500000, JD 2451727.625000
#                 B 2000 = B2000.000000, J1999.998723, JD 2451544.533398
#            JD2450123.7 = B1996.109887, J1996.108693, JD 2450123.700000
#              mJD 24034 = B1924.680025, J1924.680356, JD 2424034.500000
#             MJD50123.2 = B1996.109887, J1996.108693, JD 2450123.700000
#             rJD50123.2 = B1996.108518, J1996.107324, JD 2450123.200000
#              Rjd 23433 = B1923.033172, J1923.033539, JD 2423433.000000
#              F29/11/57 = B1957.910029, J1957.909651, JD 2436171.500000
#            F2000-01-01 = B1999.999909, J1999.998631, JD 2451544.500000
# F2002-04-04T09:42:42.1 = B2002.257054, J2002.255728, JD 2452368.904654

The strings that start with prefix ‘F’ are strings read from FITS keywords that represent the date of observation.

The sky definition¶

Given an arbitrary celestial position and a sky system specification you can transform to any of the other sky system specifications. Module wcs recognizes the following built-in sky specifications:

wcs.equatorial - wcs.ecliptic - wcs.galactic - wcs.supergalactic

Reference systems are:

wcs.fk4 - wcs.fk4_no_e - wcs.fk5 - wcs.icrs - wcs.j2000

The syntax for an equatorial sky specification is either a tuple (order of the elements is arbitrary):

(sky system, equinox, reference system, epoch of observation)
e.g.: obj.skyout = (wcs.equatorial, "J1983.5", wcs.fk4, "B1960_OBS")

or a string with minimal match:

(equatorial, equinox, referencesystem, epoch of observation"
e.g.: obj.skyout = "equa J1983.5 FK4 B1960_OBS"

Celestial transformations¶

In this section we check some basic celestial coordinate transformations. Background information can be found in [Ref2] or in the background information for module celestial.

Two parameters instantiate an object from class Transformation. The first is a definition of the input celestial system and the second is a definition for the celestial output system. Method wcs.Transformation.transform() transforms coordinates associated with the celestial input system to coordinates connected to the celestial output system.

The galactic pole has FK4 coordinates (192.25,27.4) in degrees. If we want to verify this, we need to convert this FK4 coordinate to the corresponding galactic coordinate, which should be (0,90) within the limits of precision of the used numbers. The following script shows that this could be true:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

from kapteyn import wcs

world_eq = (192.25, 27.4)   # FK4 coordinates of galactic pole
tran = wcs.Transformation("EQ,fk4,B1950.0", "GAL")
world_gal = tran.transform(world_eq)
print world_gal

# Output:
# (120.8656324107187, 89.999949251695512)

# Closure test:
world_eq = tran.transform(world_gal, reverse=True)
print world_eq

# Output:
# (192.25, 27.400000000000002)

We added a closure test (parameter reverse=True) to give you some feeling about the accuracy. Closure tests usually show errors < 1e-12. We expected the pole at 90 deg., but the difference is about 5e-05 deg. That is too much so there must be another reason for the difference. The reason is described in the background information of module celestial. The galactic pole is not a star and the so called elliptic terms of aberration (only for FK4) are not apply to its position. So in fact the pole is given in FK4-NO-E coordinates. If we repeat the exercise with the appropriate input celestial definition, we get:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

from kapteyn import wcs

world_eq = (192.25, 27.4)   # FK4 coordinates of galactic pole
tran = wcs.Transformation("EQUATORIAL, fk4_no_e, B1950.0", "galactic")
world_gal = tran.transform(world_eq)
print world_gal

# Output:
# (0.0, 90.0)

world_eq = tran.transform(world_gal, reverse=True)
print world_eq

# Output:
# (192.25, 27.400000000000002)

which gives the result as expected. Note that we used attribute reverse of the Transformation class. The two previous examples show that the transformation class is very useful to check basic celestial transformations.

As another test of a standard celestial transformation, let’s check the transformation between galactic and supergalactic coordinates. The supergalactic pole (0,90) deg. has galactic(II) world coordinates (47.37,6.32) deg. The conversion program becomes then:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

from kapteyn import wcs

world_gal = (47.37, 6.32)   # Galactic l,b (II) of supergalactic pole
tran = wcs.Transformation(wcs.galactic, wcs.supergalactic)
world_sgal = tran.transform(world_gal)
print world_sgal

# Output:
# (0.0, 90.0)

world_eq = trans.transform(world_sgal, reverse=True)
print world_gal

# Output:
# (47.369999999999997, 6.3200000000000003)

which agrees with the theory.

The sky system specifications allow for defaults. So if one wants coordinates in the equatorial system with reference system FK5 and equinox J2000 then the specification wcs.fk5 will suffice. Below we demonstrate how to transform a coordinate from the FK4 system to FK5. In fact we want to demonstrate that FK4 is slowly rotating with respect to the inertial FK5 system. We do that by varying the assumed time of observation and convert the position (R.A.,Dec) = (0,0). This behaviour is explained in the background documentation of module celestial:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

#!/usr/bin/env python
from kapteyn import wcs

world_eq1 = (0,0)
s_out = wcs.fk5
epochs = range(1950,2010,10)
for ep in epochs:
   s_in = "EQUATORIAL B1950 fk4 " + 'B'+str(ep)
   tran = wcs.Transformation(s_in, s_out)
   world_eq2 = tran.transform(world_eq1)
   print 'B'+str(ep), world_eq2

# Output:
# B1950 (0.64069100057541584, 0.27840943507737015)
# B1960 (0.64069761256120361, 0.2783973383470032)
# B1970 (0.64070422454697784, 0.27838524161663253)
# B1980 (0.64071083653273853, 0.27837314488625808)
# B1990 (0.64071744851848544, 0.27836104815588009)
# B2000 (0.64072406050421915, 0.27834895142549831)

Usually FK4 catalog values are in equinox and epoch B1950.0, so this program shows an exceptional case.

Combining projections and celestial transformations¶

In previous sections we showed examples how to use methods of an object of class Projection to convert between pixel coordinates and world coordinates. We added the option to change the celestial definition. If your data is a spatial map and its sky system is FK5, then we can convert pixel positions to world coordinates in for example galactic coordinates by specifying a value for attribute wcs.Projection.skyout. In our case this would be for a projection object called proj:

>>> proj.skyout = wcs.galactic

In the next example we test (like in one of the previous examples) a conversion between an equatorial system and the galactic system. The FK4-NO-E coordinates of the galactic pole are the values (CRVAL1, CRVAL2) from the header. First we calculate a couple of world coordinates in the native celestial definition. Then we verify that that native system is indeed FK4-NO-E and the equinox is B1950. That can be verified with:

>>> proj.skyout = (wcs.equatorial, wcs.fk4_no_e, 'B1950')

Finally we test the conversion to galactic coordinates with:

>>> proj.skyout = wcs.galactic

With the output sky set to galactic, we find the galactic pole in galactic coordinates i.e. (90,0) deg. Finally we want to know what the values of the input pixel coordinates are if the output sky system is supergalactic. The galactic pole is (90, 6.32) deg. in supergalactic coordinates. Within the limits of the precision of the used numbers we find the expected output with this script:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

from kapteyn import wcs
header = { 'NAXIS'  : 2,
           'NAXIS1' : 5,
           'CTYPE1' : 'RA---TAN',
           'CRVAL1' : 192.25,
           'CRPIX1' : 5,
           'CUNIT1' : 'degree',
           'CDELT1' : -0.01,
           'NAXIS2' : 10,
           'CTYPE2' : 'DEC--TAN',
           'CRVAL2' : 27.4,
           'CRPIX2' : 5,
           'CUNIT2' : 'degree',
           'CDELT2' : +0.01,
           'RADESYS': 'FK4-NO-E',
           'EQUINOX': 1950.0
         }

proj = wcs.Projection(header)

pixel = [(4,5),(5,5),(6,5)]   # List with coordinate tuples
world = proj.toworld(pixel)
print world
# [(192.26126360281495, 27.399999547653639), (192.25, 27.399999999999999), ...

proj.skyout = "Equatorial FK4-NO-E B1950"
world = proj.toworld(pixel)
print world
# [(192.26126360281495, 27.399999547653639), (192.24999999999997, 27.400000000000002),...

proj.skyout = "galactic"
world = proj.toworld(pixel)
print world
# [(33.00000000001878,  89.990000000101531), (0.0, 90.0), ...

proj.skyout = wcs.supergalactic
world = proj.toworld(pixel)
print world
# [(90.002497049104363, 6.3296871263660073), (90.000000000000014, 6.319999999999995), ...

Note that the second tuple on each line of the output represents the world coordinates at CRPIX. Also important is the observation that the longitude for galactic coordinates shows erratic behaviour. The reason is that close to a pole, the longitudes are less well defined (and undefined on the pole) and the errors in longitudes become important because we are calculating with numbers with a limited precision.

Attributes of a Projection object related to celestial systems¶

There are a number of attributes of an object of class wcs.Projection, related to celestial systems, that can be used to inspect the parsed FITS header. The native system in the previous example could be derived from attribute wcs.Projection.skysys:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

from kapteyn import wcs
header = { 'NAXIS'  : 2,
           'NAXIS1' : 5,
           'CTYPE1' : 'RA---TAN',
           'CRVAL1' : 192.25,
           'CRPIX1' : 5,
           'CUNIT1' : 'degree',
           'CDELT1' : -0.01,
           'NAXIS2' : 10,
           'CTYPE2' : 'DEC--TAN',
           'CRVAL2' : 27.4,
           'CRPIX2' : 5,
           'CUNIT2' : 'degree',
           'CDELT2' : +0.01,
           'RADESYS': 'FK4-NO-E',
           'EQUINOX': 1950.0,
           'MJD-OBS': 36010.2
         }

proj = wcs.Projection(header)
print "Attributes of 'proj':"
print "skysys:    ", proj.skysys
print "equinox:   ", proj.equinox
print "epoch:     ", proj.epoch
print "dateobs:   ", proj.dateobs
print "mjdobs:    ", proj.mjdobs
print "epobs:     ", proj.epobs

# Attributes of 'proj':
# skysys:     (0, 5, 'B1950.0')
# equinox:    1950.0
# epoch:      B1950.0
# dateobs:    None
# mjdobs:     36010.2
# epobs:      MJD36010.2

Below a table with a short explanation of the attributes. More information about epochs and equinoxes can be found in the documentation of celestial.

Available functions from celestial¶

Some of the functions defined in the module celestial are also available in the namespace of wcs. One of these is celestial.epochs() for which we wrote an example in the previous section. Others are celestial.lon2hms(), celestial.lon2dms() and celestial.lat2hms() to format degrees into hours, minutes, seconds or degrees, minutes and seconds. Finally the function celestial.skymatrix() is also available to wcs; it calculates the rotation matrix to convert a coordinate from one sky system to another and it calculates the E-terms (see background documentation for celestial) if appropriate. Usually you will only use this function to compare rotation matrices with matrices from the literature or to do some debugging. Some examples on the Python command line:

Formatting spatial coordinates:

>>> wcs.lon2hms(45.0)
'03h 00m  0.0s'
>>> wcs.lon2hms(23.453839, 4)
'01h 33m 48.9214s'
>>> wcs.lon2dms(245.0, 4)
Out[10]: ' 245d  0m  0.0000s'
>>> wcs.lat2dms(45.0)
'+45d 00m  0.0s'
>>> help(wcs.lon2hms)

Calculate a rotation matrix:

>>> wcs.skymatrix(wcs.galactic, wcs.supergalactic)
(matrix([[ -7.35742575e-01,   6.77261296e-01,  -6.08581960e-17],
        [ -7.45537784e-02,  -8.09914713e-02,   9.93922590e-01],
        [  6.73145302e-01,   7.31271166e-01,   1.10081262e-01]]), None, None)

Introduction¶

The most important documentation about conversions of spectral coordinates in WCSLIB is found paper “Representations of spectral coordinates in FITS” (paper III, [Ref3] ) In the next sections we show how wcs/WCSLIB can deal with spectral conversions with the focus on conversions between frequencies and velocities. We discuss conversion examples shown in the paper in detail and try to illustrate how wcs deals with FITS data from (legacy) AIPS and GIPSY sources. In many of those files the reference frequencies and reference velocities are not given in the same reference system (e.g. topocentric vs. barycentric). It is estimated that there are many of these FITS files and that their headers generate wrong results when they are used to create an object the constructor of wcs.Projection class unmodified. For FITS files generated with legacy software some extra interpretation of the FITS header is applied. This procedure is described in more detail in the background information related to spectral coordinates.

Transformations between frequencies and velocities¶

We built applications that use WCSLIB to convert grid positions, in an image or a spectrum, to world coordinates. For spectral axes with frequency as the primary type (e.g. in the FITS header we read CTYPE3=’FREQ’), it is possible to convert between pixel coordinates and frequencies, but also, if the header provides the correct information, between pixel coordinates and velocities. WCSLIB expects that in a FITS header the given frequencies are bound to the same standard of rest (i.e. reference system) as the given reference velocity. In practice however there are many FITS files that list the frequencies in the topocentric system and a reference velocity in an inertial system (barycentric, lsrk). In those FITS files the inertial systems are usually abbreviated with ‘HEL’ or ‘LSR’ (Heliocentric, Local Standard of Rest) and the velocities are usually not the true velocities but are either the so called radio or optical velocities (of which we give the definitions in the background information about spectral coordinates).

Basic spectral line header example¶

In “Representations of spectral coordinates in FITS” ([Ref3] ) section 10.1 deals with an example of a VLA spectral line cube which is regularly sampled in frequency (CTYPE3=’FREQ’). The section describes how one can define alternative FITS headers to deal with different velocity definitions. We want to examine this exercise in more detail than provided in the article to illustrate how a FITS header can be modified. In the background information you find a more elaborate discussion. Here we summarize some results.

The topocentric spectral properties in the FITS header from the paper are:

CTYPE3= 'FREQ'
CRVAL3=  1.37835117405e9
CDELT3=  9.765625e4
CRPIX3=  32
CUNIT3= 'Hz'
RESTFRQ= 1.420405752e+9
SPECSYS='TOPOCENT'

Usually such descriptions are part of a header that describes a three dimensional data structure where the first two axes represent a spatial map as function of the third axis which is a spectral axis. This example tells us that the spatial data corresponding with channel 32 was observed at a topocentric frequency (SPECSYS=’TOPOCENT’) of 1.37835117405 GHz. The step size in frequency is 97.65625 kHz. A rest frequency (1.420405752e+9 Hz) is needed to convert frequencies to velocities. The description of standard FITS keywords can be found in [FITS]

The topocentric frequency (for the receiver) was derived from a barycentric optical velocity of 9120 km/s that was requested by an observer.

We prepared a minimal header to simulate this FITS header and calculate world coordinates for the spectral axis. The numbers are frequencies. The units are Hz and the central frequency is CRVAL3. The step in frequency is CDELT3. Our minimal header (here presented as a Python dictionary) shows only one axis so our header items got axis number 1 (e.g. CRVAL1, CDELT1, etc.):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

from kapteyn import wcs
header = { 'NAXIS'  :  1,
           'CTYPE1' : 'FREQ',
           'CRVAL1' :  1.37835117405e9,
           'CRPIX1' :  32,
           'CUNIT1' : 'Hz',
           'CDELT1' :  9.765625e4
         }
proj = wcs.Projection(header)
pixels = range(30,35)
Fwcs = proj.toworld1d(pixels)
for p,f in zip(pixels, Fwcs):
   print p, f

# Output:
30 1378155861.55
31 1378253517.8
32 1378351174.05
33 1378448830.3
34 1378546486.55

The output shows frequency as function of pixel coordinate. Pixel coordinate 32 (=*CRPIX1*) shows the value of CRVAL1. Now we have a method to find at which frequency a spatial map in the data cube was observed.

WCSLIB velocities from frequency data¶

Usually similar FITS headers provide information about a velocity. Velocities is what we need for the analysis of the kinematics and dynamics of the observed objects. But there are several definitions for velocities (radio, optical, apparent radial).

For the radio interferometer, like the WSRT, an observer requesting for an observation, needs to specify:

• A rest frequency
• A velocity or Doppler shift
• A frame definition (bary or lsrk)
• A conversion type (z, radio, optical)
• A time of observation. This time is needed (together with the location of the observatory) to calculate the topocentric frequencies needed for the receivers

The observer requests that an observation must correspond to a velocity or Doppler shift (see list below) and a reference system. Only then topocentric frequencies for the receivers can be calculated.

To convert to another spectral type the constructor from class wcs.Projection needs to know which spectral type we want to convert to. The translation is set then with wcs.Projection.spectra(). which stands for spectral translation.

The parameter that we need to set the translation is ctype. Its syntax follows the FITS convention, see note below.

The following spectral types are supported (from [Ref3]):

The non-linear algorithm codes are (from [Ref3]):

If we want to convert pixel coordinates to optical velocities for our example header, then module wcs needs to create a new projection object with ctype = VOPT-F2W because VOPT represents an optical velocity and F2W sets the non linear algorithm which converts from the domain where the step size is constant (frequency) to a velocity associated with wavelength (see table above). The following script shows how to use the method wcs.Projection.spectra() to create this new object and how to convert the pixel coordinates:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS'  : 1,
           'CTYPE1' : 'FREQ',
           'CRVAL1' : 1.37835117405e9,
           'CRPIX1' : 32,
           'CUNIT1' : 'Hz',
           'CDELT1' : 9.765625e4,
           'RESTFRQ': 1.420405752e+9
         }
proj = wcs.Projection(header)
spec = proj.spectra('VOPT-F2W')
pixels = range(30,35)
Vwcs = spec.toworld1d(pixels)
print "Pixel, velocity (%s)" % spec.units
for p,v in zip(pixels, Vwcs):
   print p, v/1000.0

# Output:
# Pixel, velocity (m/s)
# 30 9190.68652655
# 31 9168.7935041
# 32 9146.90358389
# 33 9125.01676527
# 34 9103.13304757

Some comments about this example:

• It shows how to add the spectral translation to the projection object. For a conversion from frequency to optical velocity one can derive a new object with spec = proj.spectra(‘VOPT-F2W’) or proj = wcs.Projection(header).spectra(‘VOPT-F2W’).
• The output is a list with pixel coordinates and topocentric velocities. This explains why we don’t see the requested velocity (9120 km/s) at CRPIX because that velocity was barycentric.
• When we enter an invalid algorithm code for the velocity, the script will raise an exception.

Why do we need a rest frequency?

To get a velocity, the rest frequency needs to be added (RESTFRQ=) to our minimal header. What you get then is a list of velocities according to:

We adopted variable Z for velocities following the optical definition. The frequency as (linear) function of pixel coordinate is:

where:

• \(\nu_{ref}\) is the reference frequency (CRVAL1)
• \(N\) is the pixel coordinate (FITS definition) we are interested in,
• \(N_{\nu_{ref}}\) is the frequency reference pixel (CRPIX1)
• \(\delta \nu\) is the frequency increment (CDELT1)

Let’s check this with a small script:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

from kapteyn import wcs

header = { 'NAXIS'  : 1,
           'CTYPE1' : 'FREQ',
           'CRVAL1' : 1.37835117405e9,
           'CRPIX1' : 32,
           'CUNIT1' : 'Hz',
           'CDELT1' : 9.765625e4,
           'RESTFRQ': 1.420405752e+9
         }
proj = wcs.Projection(header)
spec = proj.spectra(ctype='VOPT-F2W')
pixels = range(30,35)
Vopt = spec.toworld1d(pixels)

print "Pixel coordinate and velocity (%s) with wcs module:" % spec.units
for p,Z in zip(pixels, Vopt):
   print p, Z/1000.0

print "\nPixel coordinate and velocity (%s) with documented formulas:" % spec.units
for p in pixels:
   nu = header['CRVAL1'] + (p-header['CRPIX1'])*header['CDELT1']
   Z = wcs.c*(header['RESTFRQ']-nu)/nu     # wcs.c is speed of light in m/s
   print p, Z/1000.0

# Pixel coordinate and velocity (m/s) with wcs module:
# 30 9190.68652655
# 31 9168.7935041
# 32 9146.90358389
# 33 9125.01676527
# 34 9103.13304757

# Pixel coordinate and velocity (m/s) with documented formulas:
# 30 9190.68652655
# 31 9168.7935041
# 32 9146.90358389
# 33 9125.01676527
# 34 9103.13304757

More checks are documented in the background information for spectral coordinates. This one should give you some idea how WCSLIB transforms spectral coordinates. But we still didn’t address the question about the reference systems. In our code example, this velocity Z is topocentric (defined in the reference system of the observatory) and is not suitable for comparisons because the Earth is moving around its axis and around the Sun. Other reference systems are the barycenter of the Solar system and the Local Standard of Rest. During observations one knows the location of the source, the time of observation and the location of the observatory on Earth. Software then can calculate the (true) velocity of the Earth with respect to a selected inertial reference system and we can transform from topocentric velocities to velocities in another system. Usually these correction velocities (called topocentric correction) are not recorded in the FITS file of the data set. The keyword to look for is VELOSYS=

In the background information about spectral coordinates we give a recipe how one can change the value of the reference frequency in CRVAL1 to a barycentric value. The result is CRVAL1=1.37847121643e+9 If you substitute this value for CRVAL1 in the previous script, the output is:

Pixel coordinate and velocity (m/s) with wcs module:
30 9163.77531673
31 9141.88610757
32 9119.99999984
33 9098.11699288
34 9076.23708605

At pixel coordinate 32 (CRPIX1) the velocity is 9120 km/s as we required. So wcs always returns velocities in the same system as the system of reference frequency.

Spectral CTYPE’s with special extensions¶

There are many (old) FITS headers which describe a system where the reference frequency is topocentric and the required reference velocity is given for another reference system. These velocities are given with keywords like VELR or DRVALn and the reference system for the velocities is given as an extension in CTYPEn (e.g.: CTYPE3=’FREQ-OHEL’). Image processing systems like AIPS and GIPSY have their own tools to deal with this. If wcs recognizes a legacy header, it tries to convert the reference frequency to the system of the required velocity:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

from kapteyn import wcs

header = { 'NAXIS'  : 1,
           'CTYPE1' : 'FREQ-OHEL',
           'CRVAL1' : 1.415418199417E+09,
           'CRPIX1' : 32,
           'CUNIT1' : 'HZ',
           'CDELT1' : -7.812500000000E+04,
           'VELR'   : 1.050000000000E+06,
           'RESTFRQ': 0.14204057520000E+10
         }

proj = wcs.Projection(header)
ctype = 'FREQ-???'
if ctype != None:
   spec = proj.spectra(ctype)
   print "\nSelected spectral translation with algorithm code:", spec.ctype[0]
else:
   spec = proj

crpix = header['CRPIX1']
print "CRVAL from header=%f, CRVAL modified=%f" % (header['CRVAL1'], spec.crval[0])
print "CDELT from header=%f, CDELT modified=%f" % (header['CDELT1'], spec.cdelt[0])
for i in range(-2,+3):
   px = crpix + i
   world = spec.toworld1d(px)
   print "%d %f" % (px, world)

# Output:
# Selected spectral translation with algorithm code: FREQ
# CRVAL from header=1415418199.417000, CRVAL modified=1415448253.482287
# CDELT from header=-78125.000000, CDELT modified=-78123.341180
# 30 1415604500.164647
# 31 1415526376.823467
# 32 1415448253.482287
# 33 1415370130.141107
# 34 1415292006.799927

As spectral translation we selected ‘FREQ’. If you inspect the output list with frequencies then you will see that the list doesn’t show the topocentric frequencies (with CRVAL1 at CRPIX1) but frequencies in the reference system of the given (helocentric) velocity. The attributes spec.crval[0] and spec.cdelt[0] show new values unequal to the header values.

If you want a list with topocentric frequencies then just omit to apply the wcs.Projection.spectra() method (i.e. use ctype = None in example). The output is what we expect:

# Output:
# CRVAL from header=1415418199.417000, CRVAL modified=1415418199.417000
# CDELT from header=-78125.000000, CDELT modified=-78125.000000
# 30 1415574449.417000
# 31 1415496324.417000
# 32 1415418199.417000
# 33 1415340074.417000
# 34 1415261949.417000

A note about algorithm codes¶

It is not always easy to figure out what the algorithm code should be if you want to convert to another spectral type. Therefore WCSLIB allows wildcard characters for the last or the last three characters in CTYPEn. In our example valid entries are:

• spec = proj.spectra(ctype=’VOPT-F2W’)
• spec = proj.spectra(ctype=’VOPT-F2?’)
• spec = proj.spectra(ctype=’VOPT-???’)

The missing algorithm code is returned in wcs.Projection.ctype as in:

>>> spec = proj.spectra(ctype='VOPT-???')
>>> print "Spectral translation with algorithm code:", spec.ctype[0]
    Spectral translation with algorithm code: VOPT-F2W

Module wcs uses this feature to build a list with all spectral translations that are allowed for a given Projection object. For each type in the table with spectral types, the wildcards are used to find the algorithm code (assuming that for the given Projection objects and the spectral type only one algorithm is possible). A tuple is created with the allowed spectral translation as first element and its associated unit as second element) and the tuple is added to the list wcs.Projection.altspec.

The attribute is useful if you want to write code that prompts a user to enter a spectral translation from a list of allowed translations. It can be used as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

from kapteyn import wcs

header = { 'NAXIS'  : 1,
           'CTYPE1' : 'VOPT',
           'CRVAL1' : 9120,
           'CRPIX1' : 32,
           'CUNIT1' : 'km/s',
           'CDELT1' : -21.882651442,
           'RESTFRQ': 1.420405752e+9
         }

proj = wcs.Projection(header)
print "Allowed spectral translations:"
for as in proj.altspec:
   print as
spec = proj.spectra(ctype='FREQ-???')
print "\nSelected spectral translation with algorithm code:", spec.ctype[0]

# Output:
# Allowed spectral translations:
# ('FREQ-W2F', 'Hz')
# ('ENER-W2F', 'J')
# ('VOPT', 'm/s')
# ('VRAD-W2F', 'm/s')
# ('VELO-W2V', 'm/s')
# ('WAVE', 'm')
# ('ZOPT', '')
# ('BETA-W2V', '')

# Selected spectral translation with algorithm code: FREQ-W2F

From velocities to frequencies¶

In the background information about spectral coordinates we calculated that for a barycentric system the step size in barycentric velocity is -21.882651442 km/s. Then we are able to setup a header with velocities and use a spectral translation that converts to frequencies, as in the next example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

from kapteyn import wcs

header = { 'NAXIS'  : 1,
           'CTYPE1' : 'VOPT-F2W',
           'CRVAL1' : 9120,
           'CRPIX1' : 32,
           'CUNIT1' : 'km/s',
           'CDELT1' : -21.882651442,
           'RESTFRQ': 1.420405752e+9
         }

proj = wcs.Projection(header)
spec = proj.spectra(ctype='FREQ-???')
print "Spectral translation with algorithm code:", spec.ctype[0]
pixels = range(30,35)
Freq = spec.toworld1d(pixels)

print "Pixel coordinate and frequency (%s)" % spec.units
for p,f in zip(pixels, Freq):
   print p, f

# Output:
# Pixel coordinate and frequency (Hz):
# 30 1378275920.94
# 31 1378373568.68
# 32 1378471216.43
# 33 1378568864.18
# 34 1378666511.92

The reference frequency is at pixel coordinate 32 and its value (1378471216.43) is exactly the barycentric reference frequency that we used before. What happens if we left out the algorithm code in the header? The output differs (except for the reference frequency at pixel 32). That is because it is assumed that the increments in wavelength are constant and not those in frequency. This is confirmed by the returned algorithm code which is FREQ-W2F if CTYPE1=’VOPT’

Processing real FITS data¶

With the knowledge we have at this moment, it is easy to make a small utility that looks for a spectral axis in a FITS file and if it can find one, it converts 5 pixel coordinates in the neighbourhood of CRPIX to world coordinates for all allowed spectral translations:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

from kapteyn import wcs
import pyfits

f = raw_input("Enter name of FITS file: ")
hdulist = pyfits.open(f)
header = hdulist[0].header
proj = wcs.Projection(header)
ax = proj.specaxnum
if ax == None:
   print "No spectral axis available"
else:
   print "Spectral type from header:", proj.ctype[ax-1]
   crpix = header['CRPIX'+str(ax)]
   for alt in proj.altspec:
      line = proj.sub((ax,)).spectra(alt[0])
      print "Pixel, world for translation %s" % alt[0]
      for i in range(-2,+3):
         px = crpix + i
         world = line.toworld1d(px)   #  to world coordinates
         print "%d %.10g (%s)" % (px, world, alt[1])

The projection object reads its header data from the first hdu of the FITS file (hdulist[0].hdr) and is set to only convert the spectral axis of the data set: proj.sub((ax,)). Remember that the argument is a Python tuple but we have only one axis so the tuple has an extra comma. Header items can be read from the header directly (e.g. header[‘CRPIX3’]). That’s how we find the value of CRPIX for the spectral axis. The allowed spectral translations are read from attribute wcs.Projection.altspec.

We ran the example for a fits file called mclean.fits which is a HI data cube and the third axis is the spectral axis:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

Enter name of FITS file: mclean.fits
Spectral type from header: FREQ
Pixel, world for translation FREQ
28 1415604500 (Hz)
29 1415526377 (Hz)
30 1415448253 (Hz)
31 1415370130 (Hz)
32 1415292007 (Hz)
Pixel, world for translation ENER
28 9.379902296e-25 (J)
29 9.379384645e-25 (J)
30 9.378866994e-25 (J)
31 9.378349343e-25 (J)
32 9.377831692e-25 (J)
Pixel, world for translation VOPT-F2W
28 1016794.655 (m/s)
29 1033396.411 (m/s)
30 1050000 (m/s)
31 1066605.422 (m/s)
32 1083212.677 (m/s)
etc. etc.

A simple utility to analyze a FITS file¶

With module maputils one can extract image data from a FITS file, modify it and write it to another FITS file on disk. The methods we use for these purposes are based on package PyFITS (or the equivalent in Astropy, astropy.io.fits) but are adapted to function in the environment of the Kapteyn Package. Note that PyFITS is not part of the Kapteyn Package.

With maputils one can also extract data slices from data described by more than two axes (e.g. images as function of velocity or polarization). For data with only two axes, it can swap those axes (e.g. to swap R.A. and Declination). Also the limits of the data can be set to extract part of 2-dimensional data. To be able to create plots of unfamiliar data without any user interaction, you need to know some of the characteristics of this data before you can extract the right slice. Module maputils provides routines that can display this relevant information.

First you need to create an object from class maputils.FITSimage. Some information is extracted from the FITS header directly. Other information is extracted from attributes of a wcs.Projection object defined in module wcs. Method maputils.FITSimage.str_axisinfo() gets its information from a header and its associated Projection object. It provides information about things like the sky system and the spectral system, as strings, so the method is suitable to get verbose information for display on terminals and in gui’s.

Next we show a simple script which prints meta information of the FITS file ngc6946.fits:

Example: mu_fitsutils.py - Print meta data from FITS header or Python dictionary

from kapteyn import maputils
from matplotlib import pylab as plt

fitsobject = maputils.FITSimage('ngc6946.fits')

print("HEADER:\n")
print fitsobject.str_header()

print("\nAXES INFO:\n")
print fitsobject.str_axisinfo()

print("\nEXTENDED AXES INFO:\n")
print fitsobject.str_axisinfo(long=True)

print("\nAXES INFO for image axes only:\n")
print fitsobject.str_axisinfo(axnum=fitsobject.axperm)

print("\nAXES INFO for non existing axis:\n")
print fitsobject.str_axisinfo(axnum=4)

print("SPECTRAL INFO:\n")
fitsobject.set_imageaxes(axnr1=1, axnr2=3)
print fitsobject.str_spectrans()

This code generates the following output:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

print """
HEADER:

SIMPLE  =                    T / SIMPLE FITS FORMAT
BITPIX  =                  -32 / NUMBER OF BITS PER PIXEL
NAXIS   =                    3 / NUMBER OF AXES
NAXIS1  =                  100 / LENGTH OF AXIS
etc.

AXES INFO:

Axis 1: RA---NCP  from pixel 1 to   100
{crpix=51 crval=-51.2821 cdelt=-0.007166 (DEGREE)}
{wcs type=longitude, wcs unit=deg}
etc.

EXTENDED AXES INFO:

axisnr     - Axis number:  1
axlen      - Length of axis in pixels (NAXIS):  100
ctype      - Type of axis (CTYPE):  RA---NCP
axnamelong - Long axis name:  RA---NCP
axname     - Short axis name:  RA
etc.

WCS INFO:

Current sky system:                 Equatorial
reference system:                   ICRS
Output sky system:                  Equatorial
Output reference system:            ICRS
etc.

SPECTRAL INFO:

0   FREQ-V2F (Hz)
1   ENER-V2F (J)
2   WAVN-V2F (1/m)
3   VOPT-V2W (m/s)
etc.
"""

The example extracts data from a FITS file on disk as given in the example code. To make the script a real utility one should allow the user to enter a file name. This can be done with Python’s raw-input function but to make it robust one should check the existence of a file, and if a FITS file has more than one header, one should prompt the user to specify the header. We also have to deal with alternate headers for world coordinate systems (WCS) etc. etc.

The projection class from wcs interprests and stores the header information. It serves as the interface between maputils and the library WCSLIB.

Example: mu_projection.py - Get data from attributes of the projection system

from kapteyn import maputils

print "Projection object from FITSimage:" 
fitsobj = maputils.FITSimage("mclean.fits")
print "crvals:", fitsobj.convproj.crval
fitsobj.set_imageaxes(1,3)
print "crvals after axes specification:", fitsobj.convproj.crval
fitsobj.set_spectrans("VOPT-???")
print "crvals after setting spectral translation:", fitsobj.convproj.crval

print "Projection object from Annotatedimage:"
annim = fitsobj.Annotatedimage()
print "crvals:", annim.projection.crval

The output:

Projection object from FITSimage:
crvals: (178.7792, 53.655000000000001)
crvals after axes specification: (178.7792, 1415418199.4170001, 53.655000000000001)
crvals after setting spectral translation: (178.7792, 1050000.0000000042, 53.655000000000001)
Projection object from Annotatedimage:
crvals: (178.7792, 1050000.0000000042, 53.655000000000001)

Explanation:

As soon as the FITSimage object is created, it is assumed that the first two axes are the axes of the data you want to display. After setting alternative axes, a slice is taken and the projection system is changed. Now it shows attributes in the order of the slice and we see a third value in the tuple with term:CRVAL‘s. That’s because the last value represents the necessary matching spatial axis which is needed to do conversions between pixels and world coordinates.

As a next step we set a spectral translation for the second axis which is a frequency axis. Again the projection system changes. We did set the translation to an optical velocity and the printed CRVAL is indeed the reference optical velocity from the header of this FITS file (1050 km/s).

When an object from class maputils.Annotatedimage is created, the projection data is copied from the maputils.FITSimage object. It can be accessed through the maputils.Annotatedimage.projection attribute (last line of the above example).

Specification of a map¶

Class maputils.FITSimage extracts data from a FITS file so that a map can be plotted with its associated world coordinate system. So we have to specify a number of parameters to get the required image data. This is done with the following methods:

• Header - The constructor maputils.FITSimage needs name and path of the FITS file. It can be a file on disk or an URL. The file can be zipped. A FITS file can contain more than one header data unit. If this is an issue you need to enter the number of the unit that you want to use. A case insensitive name of the hdu is also allowed. A FITS header can also contain one or more alternate headers. Usually these describe another sky or spectral system. We list three examples. The first is a complete description of the FITS header. The second get its parameters from an external ‘prompt’ function (see next section) and the third uses a prompt function with a pre specfication of parameter alter which sets the alternate header.

  >>> fitsobject = maputils.FITSimage('alt2.fits', hdunr=0, alter='A', memmap=1)
  >>> fitsobject = maputils.FITSimage(promptfie=maputils.prompt_fitsfile)
  >>> fitsobject = maputils.FITSimage(promptfie=maputils.prompt_fitsfile, alter='A')
  >>> fitsobject = maputils.FITSimage('NICMOSn4hk12010_mos.fits', hdunr='err')
  

  Later we will also discuss examples where we use external headers (i.e. not from a FITS file, but as a Python dictionary) and external data (i.e. from another source or from processed data). The user/programmer is responsible for the right shape of the data. Here is a small example of processed data:

  >>> f = maputils.FITSimage("m101.fits", externaldata=log(abs(fftA)+1.0))
  

• Axis numbers - Method maputils.FITSimage.set_imageaxes() sets the axis numbers. These numbers follow the FITS standard, i.e. they start with 1. If the data has only two axes then it is possible to swap the axes. This method can be used in combination with an external prompt function. If the data has more than two axes, then the default value for the axis numbers axnr1=1 and axnr2=2. One can also enter the names of the axes. The input is case insensitive and a minimal match is applied to the axis names found in the CTYPE keywords in the header. Examples are:

  >>> fitsobject.set_imageaxes(promptfie=maputils.prompt_imageaxes)
  >>> fitsobject.set_imageaxes(axnr1=2, axnr2=1)
  >>> fitsobject.set_imageaxes(2,1)         # swap
  >>> fitsobject.set_imageaxes(1,3)         # XV map
  >>> fitsobject.set_imageaxes('R','D')     # RA-DEC map
  >>> fitsobject.set_imageaxes('ra','freq') # RA-FREQ map
  

• Data slices from data sets with more than two axes - For an artificial set called ‘manyaxes.fits’, we want to extract one spatial map. The axes order is [frequency, declination, right ascension, stokes]. We extract a data slice at FREQ=2 and STOKES=1. This spatial map is obtained with the following lines:

  >>> fitsobject = maputils.FITSimage('manyaxes.fits') # FREQ-DEC-RA-STOKES
  >>> fitsobject.set_imageaxes(axnr1=3, axnr2=2, slicepos=(2,1))
  

• Coordinate limits - If you want to extract only a part of the image then you need to set limits for the pixel coordinates. This is set with maputils.FITSimage.set_limits(). The limits can be set manually or with a prompt function. Here are examples of both:

  >>> fitsobject.set_limits(pxlim=(20,40), pylim=(22,38))
  >>> fitsobject.set_limits((20,40), (22,38))
  >>> fitsobject.set_limits(promptfie=maputils.prompt_box)
  

• Output sky definition - For conversions between pixel- and world coordinates one can define to which output sky definition the world coordinates are related. The sky parameters are set with maputils.FITSimage.set_skyout(). The syntax for a sky definition (sky system, reference system, equinox, epoch of observation) is documented in celestial.skymatrix().

  >>> fitsobject = maputils.FITSimage('m101.fits')
  >>> fitsobject.set_skyout("Equatorial, J1952, FK4_no_e, J1980")
      or:
  >>> fitsobject.set_skyout(promptfie=maputils.prompt_skyout)
  

Writing data to a FITS file or to append data to a FITS file is also possible. The method written for these purposes is called writetofits(). It has parameters to scale data and it is possible to skip writing history and comment keywords. Have a look at the examples:

>>> # Write data with scaling
>>> fitsobject.writetofits(history=True, comment=True,
                           bitpix=fitsobject.bitpix,
                           bzero=fitsobject.bzero,
                           bscale=fitsobject.bscale,
                           blank=fitsobject.blank)

>>> # Append data to existing FITS file
>>> fitsobject.writetofits("standard.fits", append=True, history=False, comment=False)

Basic image¶

If one is interested in displaying image data only (i.e. without any wcs information) then we need very few lines of code as we show in the next example.

Example: mu_simple.py - Minimal script to display image data

from kapteyn import maputils
from matplotlib import pyplot as plt

f = maputils.FITSimage("m101.fits")
mplim = f.Annotatedimage()
im = mplim.Image()
mplim.plot()

plt.show()

(Source code)

Explanation:

This is a simple script that displays an image using the defaults for the axes, the limits, the color map and many other properties. From an object from class maputils.FITSimage an object from class maputils.FITSimage.Annotatedimage is derived.

This object has methods to create other objects (image, contours, graticule, ruler etc.) that can be plotted with Matplotlib. To plot these objects we need to call method maputils.Annotatedimage.plot()

If you are only interested in displaying the image and don’t want any white space around the plot then you need to specify a Matplotlib frame as parameter for maputils.Annotatedimage(). This frame is created with Matplotlib’s add_subplot() or add_axes(). The latter has options to specify origin and size of the frame in so called normalized coordinates [0,1]. Note that images are displayed while preserving the pixel aspect ratio. Therefore images are scaled to fit either the given width or the given height.

In the next example we reduce whitespace and use keyword parameters cmap, clipmin and clipmax to set a color map and the clip levels between which the color mapping is applied.

Example: mu_withimage.py - Display image data using keyword parameters

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

from kapteyn import maputils
from matplotlib import pyplot as plt

fitsobj = maputils.FITSimage("m101.fits")
fitsobj.set_limits((180,344), (180,344))

fig = plt.figure(figsize=(4,4))
frame = fig.add_axes([0,0,1,1])

annim = fitsobj.Annotatedimage(frame, cmap="spectral", clipmin=10000, clipmax=15500)
annim.Image(interpolation='spline36')
print "clip min, max:", annim.clipmin, annim.clipmax
annim.plot()

plt.show()

(Source code)

Fig.: mu_withimage.py - Image with non default plot frame and parameters for color map and clip levels.

RGB image¶

It is possible to compose an image from three separate images which represent a red, green and blue component. In this case you need to create an maputils.Annotatedimage object first. The data associated with this image can be used to draw e.g. contours, while the parameters of the method maputils.Annotatedimage.RGBimage() compose the RGB image. The maputils.Annotatedimage.RGBimage() method creates a new data array and inserts your R, G & B images in the right place. Then it scales the composed data to a range between 0 and 1. For an RGB image we don’t apply interactive color bar editing, but allow for the use of a function or lambda expression to rescale the data to enhance features. In the example we used keyword parameter fun to square the data to enhance the image a bit.

Example: mu_rgbdemo.py - display RGB image

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

from kapteyn import maputils
from matplotlib import pyplot as plt

# In the comments we show how to set a smaller box t display
f_red = maputils.FITSimage('m101_red.fits')
#f_red.set_limits((200,300),(200,300))
f_green = maputils.FITSimage('m101_green.fits')
#f_green.set_limits((200,300),(200,300))
f_blue = maputils.FITSimage('m101_blue.fits')
#f_blue.set_limits((200,300),(200,300))

# Show the three components R,G & B separately
# Show Z values when moving the mouse
fig = plt.figure()
frame = fig.add_subplot(2,2,1); frame.set_title("Red with noise")
a = f_red.Annotatedimage(frame); a.Image()
a.interact_toolbarinfo(wcsfmt=None, zfmt="%g")
frame = fig.add_subplot(2,2,2); frame.set_title("Greens are 1")
a = f_green.Annotatedimage(frame); a.Image()
a.interact_toolbarinfo(wcsfmt=None, zfmt="%g")
frame = fig.add_subplot(2,2,3); frame.set_title("Blues are 1")
a = f_blue.Annotatedimage(frame); a.Image()
a.interact_toolbarinfo(wcsfmt=None, zfmt="%g")

# Plot the composed RGB image
frame = fig.add_subplot(2,2,4); frame.set_title("RGB composed of previous")
annim = f_red.Annotatedimage(frame)
annim.RGBimage(f_red, f_green, f_blue, fun=lambda x:x*x, alpha=1)


# Note: color interaction not possible (RGB is fixed)
annim.interact_toolbarinfo(wcsfmt=None, zfmt="%g")
# Write RGB values to terminal after clicking left mouse button
annim.interact_writepos(pixfmt=None, wcsfmt="%.12f", zfmt="%.3e", hmsdms=False)
maputils.showall()

(Source code)

Explanation:

Three FITS files contain data in rectangular shapes in different positions. If you want to use only a part of the images you need to set limits (with method set_limits()) for each image separately. The shapes in these example data files have values 1 (or near 1) and do overlap to illustrate the composition of a new color. The region where three shapes overlap is white in the composed output image. For each RGB component a maputils.FITSimage is created. One of these is used to make a maputils.FITSimage.Annotatedimage object. The three FITSimage objects are used as parameters for method maputils.Annotatedimage.RGBimage() to set the individual components of a RGB image. The red component is a bit special because we added some Gaussian noise to it. A second parameter used in the method is alpha. This is an alpha factor which applies to all the pixels in the composed map.

This script also displays a message in the message tool bar with information about mouse positions and the corresponding image value. For an RGB image, all three image values (z values) are displayed. The format of the message is changed with parameters in maputils.Annotatedimage.interact_toolbarinfo() as in:

>>> annim.interact_toolbarinfo(wcsfmt=None, zfmt="%g")

Figure size¶

In the previous example we specified a size in inches for the figure. We provide a method maputils.FITSimage.get_figsize() to set the figure size in cm. Enter only the direction for which you want to set the size. The other size is calculated using the pixel aspect ratio, but then it is not garanteed that all labels will fit. The method is used as follows:

>>> fig = plt.figure(figsize=f.get_figsize(xsize=15, cm=True))

Introduction¶

Module maputils can create graticule objects with method maputils.Annotatedimage.Graticule(). But in fact the method that does all the work is defined in module wcsgrat So in the next sections we often refer to module wcsgrat.

Module wcsgrat creates a graticule for a given header with WCS information. That implies that it finds positions on a curve in 2 dimensions in image data for which one of the world coordinates is a constant value. These positions are stored in a graticule object. The positions at which these lines cross one of the sides of the rectangle (made up by the limits in pixels in both x- and y-direction), are stored in a list together with a text label showing the world coordinate of the crossing.

Simple example¶

Example: mu_axnumdemosimple.py - Simple plot using defaults

from kapteyn import maputils
from matplotlib import pyplot as plt

fitsobj = maputils.FITSimage("m101.fits")
mplim = fitsobj.Annotatedimage()
graticule = mplim.Graticule()
mplim.plot()

plt.show()

(Source code)

Explanation:

The script opens an existing FITS file. Its header is parsed by methods in module wcs and methods from classes in module wcsgrat calculate the graticule data. A plot is made with Matplotlib. Note the small rotation of the graticules.

The recipe:

• Given a FITS file on disk (m101.fits) we want to plot a graticule for the spatial axes in the FITS file.
• The necessary information is retrieved from the FITS header with PyFITS through class maputils.FITSimage.
• To plot something we need to tell method maputils.FITSimage.Annotatedimage() in which frame it must plot. Therefore we need a Matplotlib figure instance and a Matplotlib Axes instance (which we call a frame in the context of maputils).
• A graticule representation is calculated by maputils.Annotatedimage.Graticule() and stored in object grat. The maximum number of defaults are used.
• Finally we tell the Annotated image object mplim to plot itself and display the result with Matplotlib’s function show(). This last step can be compressed to one statement: maputils.showall() which plots all the annotated images in your script at once and then call Matplotlib’s function show().

The wcsgrat module estimates the ranges in world coordinates in the coordinate system defined in your FITS file. The method is just brute force. This is the only way to get good estimates for large maps, rotated maps maps with weird projections (e.g. Bonne) etc. It is also possible to enter your own world coordinates to set limits. Methods in this module calculate ‘nice’ numbers to annotate the plot axes and to set default plot attributes.

Hint: Matplotlib versions older than 0.98 use module pylab instead of pyplot. You need to change the import statement to: from matplotlib import pylab as plt

Probably you already have many questions about what wcsgrat can do more:

• Is it possible to draw labels only and no graticule lines?
• Can I change the starting point and step size for the coordinate labels?
• Is it possible to change the default tick label format?
• Can I change the default titles along the axes?
• Is it possible to highlight (e.g. by changing color) just one graticule line?
• Can I plot graticules in maps with one spatial- and one spectral coordinate?
• Can I control the aspect ratio of the plot?
• Is it possible to set limits on pixel coordinates?

In the following sections we will give a number of examples to answer most of these questions.

Selecting axes for image or graticule¶

For data sets with more than 2 axes or data sets with swapped axes (e.g. Declination as first axis and Right Ascension as second), we need to make a choice of the axes and axes order. To demonstrate this we created a FITS file with four axes. The order of the axes is uncommon and should only demonstrate the flexibility of the maputils module. We list the data for these axes in this ‘artificial’ FITS file:

Filename: manyaxes.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      44  (10, 50, 50, 4)  int32
Axis  1 is FREQ   runs from pixel 1 to    10  (crpix=5 crval,cdelt=1.37835, 9.76563e-05 GHZ)
Axis  2 is DEC    runs from pixel 1 to    50  (crpix=30 crval,cdelt=45, -0.01 DEGREE)
Axis  3 is RA     runs from pixel 1 to    50  (crpix=25 crval,cdelt=30, -0.01 DEGREE)
Axis  4 is POL    runs from pixel 1 to     4  (crpix=1 crval,cdelt=1000, 10 STOKES)

You can download the file manyaxes.fits for testing. The world coordinate system is arbitrary.

Example: mu_manyaxes.py - Selecting WCS axes from a FITS file with NAXIS > 2

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

from kapteyn import maputils
from matplotlib import pyplot as plt

# 1. Read the header
fitsobj = maputils.FITSimage("manyaxes.fits")

# 2. Create a Matplotlib Figure and Axes instance
figsize=fitsobj.get_figsize(ysize=12, xsize=11, cm=True)
fig = plt.figure(figsize=figsize)
frame = fig.add_subplot(1,1,1)

# 3. Create a graticule
fitsobj.set_imageaxes('freq','pol')
mplim = fitsobj.Annotatedimage(frame)
grat = mplim.Graticule(starty=1000, deltay=10)

# 4. Show the calculated world coordinates along y-axis
print "The world coordinates along the y-axis:", grat.ystarts

# 5. Show header information in attributes of the Projection object
#    The projection object of a graticule is attribute 'gmap'
print "CRVAL, CDELT from header:", grat.gmap.crval, grat.gmap.cdelt

# 6. Set a number of properties of the graticules and plot axes
grat.setp_tick(plotaxis="bottom",
               fun=lambda x: x/1.0e9, fmt="%.4f",
               rotation=-30 )

grat.setp_axislabel("bottom", label="Frequency (GHz)")
grat.setp_gratline(wcsaxis=0, position=grat.gmap.crval[0],
                   tol=0.5*grat.gmap.cdelt[0], color='r')
grat.setp_ticklabel(plotaxis="left", position=1000, color='m', fmt="I")
grat.setp_ticklabel(plotaxis="left", position=1010, color='b', fmt="Q")
grat.setp_ticklabel(plotaxis="left", position=1020, color='r', fmt="U")
grat.setp_ticklabel(plotaxis="left", position=1030, color='g', fmt="V")
grat.setp_axislabel("left", label="Stokes parameters")


# 7. Set a title for this frame
title = r"""Polarization as function of frequency at:
            $(\alpha_0,\delta_0) = (121^o,53^o)$"""
t = frame.set_title(title, color='#006400', y=1.01, linespacing=1.4)

# 8. Add labels inside plot
inlabs = grat.Insidelabels(wcsaxis=0, constval=1015,
                           deltapx=-0.15, rotation=90,
                           fontsize=10, color='r',
                           fun=lambda x: x*1e-9, fmt="%.4f.10^9")

w = grat.gmap.crval[0] + 0.2*grat.gmap.cdelt[0]
cv = grat.gmap.crval[1]
# Print without any formatting
inlab2 = grat.Insidelabels(wcsaxis=0, world=w, constval=cv,
                           deltapy=0.1, rotation=20,
                           fontsize=10, color='c')

pixel = grat.gmap.topixel((w,grat.gmap.crval[1]))
frame.plot( (pixel[0],), (pixel[1],), 'o', color='red' )

# 9. Plot the objects
maputils.showall()

(Source code)

The plot shows a system of grid lines that correspond to non spatial axes and it will be no surprise that the graticule is a rectangular system. The example follows the same recipe as the previous ones and it shows how one selects the required plot axes in a FITS file. The axes are set with maputils.FITSimage.set_imageaxes() with two numbers. The first axis of a set is axis 1, the second 2, etc. (i.e. FITS standard). The default is (1,2) i.e. the first two axes in a FITS header.

For a R.A.-Dec. graticule one should enter for this FITS file:

>>> f.set_imageaxes(3,2)

Let’s study the plot in more detail:

• The header shows a Stokes axes with an uncommon value for CRVAL and CDELT. We want to label four graticule lines with the familiar Stokes parameters. With the knowledge we have about this CRVAL and CDELT we tell the Graticule constructor to create 4 graticule lines (starty=1000, deltay=10).

• The four positions are stored in attribute ystarts as in grat.ystarts. We use these numbers to change the coordinate labels into Stokes parameters with method wcsgrat.Graticule.setp_ticklabel()

  >>> grat.setp_ticklabel(plotaxis="left", position=1000, color='m', fmt="I")
  

• We use wcsgrat.Insidelabels() to add coordinate labels inside the plot. We mark a position near CRVAL and plot a label and with the same method we added a single label at that position.

This example shows an important feature of the underlying module wcsgrat and that is its possibility to change properties of graticules, ticks and labels. We summarize:

• Graticule line properties are set with wcsgrat.Graticule.setp_gratline() or the equivalent wcsgrat.Graticule.setp_lineswcs1() or wcsgrat.Graticule.setp_lineswcs1(). The properties are all Matplotlib properties given as keyword arguments. One can apply these to all graticule lines, to one of the wcs types or to one graticule line (identified by its position in world coordinates).
• Graticule ticks (the intersections with the borders) are modified by method wcsgrat.Graticule.setp_tick(). Ticks are identified by either the wcs axis (e.g. longitude or latitude) or by one of the four rectangular plot axes or by a position in world coordinates. Combinations of these are also possible. Plot properties are given as Matplotlib keyword arguments. The labels can be scaled and formatted with parameters fun and fmt. Usually one uses method wcsgrat.Graticule.setp_ticklabel() to change tick labels and wcsgrat.Graticule.setp_tickmark() to change the tick markers.
• Ticks can be native to a plot axis (e.g. an pixel X axis which corresponds to a R.A. axis in world coordinates. But sometimes you can have ticks from two world coordinate axes along the same pixel axis (e.g. for a rotated plot). Then it is possible to control which ticks are plotted and which not. A tick mode for one or more plot/pixel axes is set with wcsgrat.Graticule.set_tickmode().
• The titles along one of the rectangular plot axes can be modified with wcsgrat.Graticule.setp_axislabel() which is a specialization of method wcsgrat.Graticule.setp_plotaxis(). A label text is set with parameter label and the plot properties are given as Matplotlib keyword arguments.
• Properties of labels inside a plot are set in the constructor wcsgrat.Insidelabels.setp_label().
• Properties of labels along a ruler are set with method rulers.Ruler.setp_label(). Properties of the ruler line can be changed with rulers.Ruler.setp_line()
• For labels along the plotaxes which correspond to pixel positions one can change the properties of the labels with maputils.Pixellabels.setp_label() while the properties of the markers can be changed with: maputils.Pixellabels.setp_marker()

Let’s summarize these methods in a table:

-Table- Objects related to graticules and their methods to set properties.

In the following sections we show some examples for changing the graticule properties. Note that for some methods we can identify objects either with the graticule line type (i.e. 0 or 1), the number or name of the plot axis ([0..4] or one of ‘left’, ‘bottom’, ‘right’, ‘top’ (or a minimal match of these strings). Some objects (e.g. tick labels) can also be identified by a position in world coordinates. Often also a combination of these identifiers can be used.

Graticule axis labels¶

Example: mu_labeldemo.py - Properties of axis labels

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

from kapteyn import maputils
from matplotlib import pylab as plt

header = {
'NAXIS':                     2 ,
'NAXIS1':                  100 ,
'NAXIS2':                  100 ,
'CDELT1':  -7.165998823000E-03 ,
'CRPIX1':   5.100000000000E+01 ,
'CRVAL1':  -5.128208479590E+01 ,
'CTYPE1': 'RA---NCP        ' ,
'CUNIT1': 'DEGREE          ' ,
'CDELT2':   7.165998823000E-03 ,
'CRPIX2':   5.100000000000E+01 ,
'CRVAL2':   6.015388802060E+01 ,
'CTYPE2': 'DEC--NCP        ' ,
'CUNIT2': 'DEGREE          ' ,
}

fig = plt.figure(figsize=(6,5.2))
frame = fig.add_axes([0.15,0.15,0.8,0.8])
f = maputils.FITSimage(externalheader=header)
annim = f.Annotatedimage(frame)
grat = annim.Graticule()
grat.setp_axislabel(fontstyle='italic')      # Apply to all
grat.setp_axislabel("top", visible=True, xpos=0.0, ypos=1.0, rotation=180)
grat.setp_axislabel("left",
                    backgroundcolor='y',
                    color='b',
                    style='oblique',
                    weight='bold',
                    ypos=0.3)
grat.setp_axislabel("bottom",                # Label in LaTeX
                    label=r"$\mathrm{Right\ Ascension\ (2000)}$",
                    fontsize=14)
annim.plot()
plt.show()