Tutorial wcs module

Introduction

This tutorial aims at starters. Experienced users find relevant but compact documentation in the module documentation. In this tutorial we address different practical situations where we need to convert between pixel- and world coordinates. Many examples are working scripts, others are very useful to try in an interactive Python session.

wcs is the core of the Kapteyn Package. An important feature of that package is that it provides a world coordinate system which is easy to incorporate in your own (Python) environment and wcs provides the basic methods to do this. Together with module celestial it allows a user to transform between pixel coordinates and world coordinates for a set of supported projections and sky systems. Module celestial provides a rotation matrix for sky transformations and is more or less embedded in wcs, so (for standard work) there is no need to import it separately.

Module wcs module has a number of important features:

  • Flexible I/O of coordinates
  • Support for spatial and spectral data
  • Support for ‘mixed’ coordinates
  • Support for conversions between different celestial systems
  • Objects have useful attributes
  • Easy to combine with other software written in Python

Coordinate representations

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:

#!/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( list(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.

Note

An advantage of NumPy arrays is that you can use them in mathematical expressions to process the array content. For example: assume you have a sequence of velocities in a numpy array V but want to express the numbers in km/s, then change the content with expression: V /= 1000

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 = list(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

Note

Class wcs has an attribute called debug. If you set its value to True then you get debug information from WCSLIB showing what has been correctly parsed from the given header data. Use it as follows:

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:

#!/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 = list(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( (list(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( (list(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:

#!/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:

Note

any CTYPEi not covered by convention and agreement shall be taken to be linear.

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 = list(range(3,8))
y = list(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.

Note

Note that longitude and latitude are not independent. You need always two pixel coordinates (x,y) to get a world coordinate pair (RA,DEC).

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.

Note

The coordinate representation in methods wcs.Projection.toworld() and wcs.Projection.topixel() of the output is the same as that of the input.

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():

#!/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:

Note

The transformation routines accept sequences of coordinates. Calculations with sequences are more efficient than repetitive calls in a loop.

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 = list(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:

#!/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 = list(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 arrays and matrices

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:

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.        ]]

Note

The rowvec attribute can also be set in the constructor of the projection object as follows: proj = wcs.Projection(header, rowvec=True)

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:

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:

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).

Note

In this module the values in the NumPy arrays and matrices are of type ‘f8’ (64 bit).

Attributes

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.

Note

If WCSLIB cannot find a value of one of the requested attributes, its value is set to None

#!/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 = list(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.

Invalid coordinates

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:

#!/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 data from a FITS file

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:

#!/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:

#!/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:

#!/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 transformations with wcs

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.

Note

If you need transformations of world coordinates between any of the supported input sky-/reference system then you should use objects and methods from class wcs.Transformation.

If you need to convert pixel coordinates in a system defined by (FITS) header information, then set the skyout attribute of a Projection object and use methods wcs.Projection.toworld() and wcs.Projection.topixel()

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:

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:

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:

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:

#!/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.

Note

We are not restricted to the transformation of one coordinate. The input of positions follow the rules of coordinate representations as described for methods wcs.Projection.toworld() and wcs.Projection.topixel().

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:

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.

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)

Spectral transformations

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.):

from kapteyn import wcs
header = { 'NAXIS'  :  1,
           'CTYPE1' : 'FREQ',
           'CRVAL1' :  1.37835117405e9,
           'CRPIX1' :  32,
           'CUNIT1' : 'Hz',
           'CDELT1' :  9.765625e4
         }
proj = wcs.Projection(header)
pixels = list(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.

Note

The first four characters of a spectral CTYPE specify the new coordinate type, the fifth character is ‘-’ and the next three characters specify a predefined algorithm for computing the world coordinates from intermediate physical coordinates ([Ref3] ).

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

Type Name Symbol Units Associated with
FREQ Frequency ν Hz ν
ENER Energy E J ν
WAVN Wavenumber κ 1/m ν
VRAD Radio velocity V m/s ν
WAVE Vacuum wavelength λ m λ
VOPT Optical velocity Z m/s λ
ZOPT Redshift z - λ
AWAV Air wavelength λa m λa
VELO Apparent radial velocity v m/s v
BETA Beta factor (v/c) β - v

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

Code sampled in Expressed as
F2W Frequency Wavelength
F2V Frequency Apparent radial velocity
F2A Frequency Air wavelength
W2F Wavelength Frequency
W2V Wavelength Apparent radial velocity
W2A Wavelength Air wavelength
V2F Apparent radial velocity Frequency
V2W Apparent radial velocity Wavelength
V2A Apparent radial velocity Air wavelength
A2F Air wavelength Frequency
A2W Air wavelength Wavelength
A2V Air wavelength Apparent radial velocity

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:

#!/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 = list(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:

(1)\[ Z = c ( \frac{\lambda - \lambda_0}{\lambda_0}) = c\ (\frac{\nu_0 - \nu}{\nu})\]

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

(2)\[ \nu = \nu_{ref} + (N - N_{\nu_{ref}}) \delta \nu\]

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:

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 = list(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.

Warning

Reference frequencies given in FITS keyword CRVALn refer to a reference system. This system should be given with FITS keyword SPECSYS (e.g. SPECSYS=’TOPOCENT’). Module wcs converts between frequencies and velocities in the same reference system. You should inspect your FITS header to find what this system is.

Warning

Legacy FITS headers often define frequencies in a Topocentric system. Also a reference velocity is given in another reference system. WCSLIB needs instructions how to convert between these systems. If legacy headers are recognized, module wcs tries to convert the frequency system to the reference system of the reference velocity. See also the next section and the background documentation about spectral coordinates

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:

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.

Note

For a given header the attribute wcs.Projection.altspec stores all possible spectral translations.

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:

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:

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 = list(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:

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:

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.

References

[Ref1]Representations of world coordinates in FITS http://www.atnf.csiro.au/people/mcalabre/WCS/wcs.pdf Greisen E.W. and Calabretta M.R.
[Ref2]Representations of celestial coordinates in FITS http://www.atnf.csiro.au/people/mcalabre/WCS/ccs.pdf Calabretta M.R. and Greisen E.W.
[Ref3](1, 2, 3, 4, 5) Representations of spectral coordinates in FITS http://www.atnf.csiro.au/people/mcalabre/WCS/scs.pdf E. W. Greisen, M. R. Calabretta, F. G. Valdes, and S. L. Allen
[FITS]Definition of the Flexible Image Transport System (FITS), FITS Standard Version 3.0 http://fits.gsfc.nasa.gov/fits_standard.html FITS Working Group , Commission 5: Documentation and Astronomical Data, International Astronomical Union