Tutorial wcs module
===================
.. highlight:: python
:linenothreshold: 1000
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.
:mod:`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 :mod:`wcs` provides the basic methods to do this.
Together with module :mod:`celestial` it allows a user to transform between pixel coordinates
and world coordinates for a set of supported projections and sky systems.
Module :mod:`celestial` provides a rotation matrix for sky transformations and
is more or less embedded in :mod:`wcs`, so (for standard work) there is no
need to import it separately.
.. index:: Features of the wcs module
Module :mod:`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
.. index:: Coordinate representations
.. index:: I/O structure
Coordinate representations
--------------------------
One coordinate axis
...................
For experiments and debug sessions, module :mod:`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 :mod:`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 :mod:`wcs`.
However, FITS data processed by :mod:`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 :mod:`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 :class:`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 :meth:`wcs.Projection.toworld1d` and
:meth:`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]
# [4.0, 5.0, 6.0]
# (4.0, 5.0, 6.0)
# [ 4. 5. 6.]
# [ 9. 10. 11.]
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 :meth:`wcs.Projection.toworld1d` and :meth:`wcs.Projection.topixel1d`
are special versions of the more general methods :meth:`wcs.Projection.toworld` and :meth:`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 :meth:`wcs.Projection.toworld` and :meth:`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
:meth:`wcs.Projection.toworld` and :meth:`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 :attr:`wcs.Projection.types` and :attr:`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 :meth:`wcs.Projection.toworld`
and :meth:`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 :mod:`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 :meth:`wcs.Projection.toworld` and :meth:`wcs.Projection.topixel`
of the output is the same as that of the input.
.. index:: Mixing pixel- and world coordinates
Mixed transformations (pixel- and world coordinates) using method :meth:`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 :meth:`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 :meth:`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])
.. index:: Projection objects representing data slices
.. index:: Sub-Projections
Three or more coordinate axes
.............................
In this section we discuss method :meth:`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 :meth:`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 :meth:`wcs.Projection.sub` method sets the axis order with parameter *axes*.
It sets in fact the order of the coordinates in the transformation methods :meth:`wcs.Projection.toworld`,
:meth:`wcs.Projection.topixel` and :meth:`wcs.Projection.mixed`.
Parameter *axes* is either a single integer or a list/tuple of integers e.g. sub(2) vs. sub([3,1]).
.. index:: Data in Numpy arrays
.. index:: Data in a Numpy matrix
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 :meth:`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 :attr:`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 :meth:`wcs.Projection.toworld` and :meth:`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).
.. index:: Attributes of a Projection object
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 :attr:`wcs.Projection.lonaxnum`,
:attr:`wcs.Projection.lataxnum` and :attr:`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.
.. index:: Position-Velocity plots
.. index:: XV maps
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 :meth:`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 :mod:`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.
.. index:: Suppressing exceptions in coordinate transformations
.. index:: Exception suppression
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
:attr:`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])
.. index:: Reading headers from FITS files
.. index:: Header data from a FITS file
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 :mod:`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.
.. index:: Reading headers from FITS files (example)
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 :attr:`wcs.Projection.types` and :attr:`wcs.Projection.units`.
The script is a useful tool to inspect the FITS file and to check its parsing by module
:mod:`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 :mod:`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 :attr:`wcs.Projection.lonaxnum`
and :attr:`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.
.. index:: FITS: Creating a Projection object for a spatial map in a FITS file (example)
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 :meth:`wcs.Projection.toworld` and :meth:`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 :mod:`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 :mod:`wcs` so one can use it
in the context of :mod:`wcs`, with the class :class:`wcs.Transformation`.
for conversions of world coordinates between sky-/reference systems
and also, if pixel coordinates are involved, methods
:meth:`wcs.Projection.toworld` and :meth:`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
:attr:`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 :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
:meth:`wcs.Projection.toworld` and :meth:`wcs.Projection.topixel`
The celestial definitions are described in detail in the background information of
module :mod:`celestial`. We list the most important features of a
celestial definition:
Supported Sky systems (detailed information in :ref:`celestial-skysystems`):
1. Equatorial: Equatorial coordinates (:math:`\alpha`, :math:`\delta`), see next list with reference systems
2. Ecliptic: Ecliptic coordinates (:math:`\lambda`, :math:`\beta`) 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 :ref:`celestial-refsystems`):
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 :ref:`celestial-epochs`):
The equinox and epoch of observations are instants of time and are of type string.
These strings are parsed by a function of module :mod:`celestial` called :func:`celestial.epochs`.
The parser rules are described in the documentation for that function.
Each string starts with a prefix. Supported prefixes are:
#. B: Besselian epoch
#. J: Julian epoch
#. JD: Julian date
#. MJD: Modified Julian Day
#. RJD: Reduced Julian Day
#. 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 :func:`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 :meth:`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 :mod:`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 :mod:`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 :meth:`wcs.Projection.toworld` and :meth:`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 :attr:`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.
Attributes of a Projection object related to celestial systems
..............................................................
There are a number of attributes of an object of class :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 :attr:`wcs.Projection.skysys`::
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 :mod:`celestial`.
.. tabularcolumns:: |p{15mm}|p{135mm}|
========== ===============================================================
Attribute Explanation
========== ===============================================================
skysys A single value or tuple which defines the native system.
Tuples can contain the sky system, the reference system,
the equinox and the date of observation.
equinox equinox is a floating point number. It is read from the
FITS header (keyword EQUINOX).
The equinox is a moment in time
used for the definition of an equatorial system.
epoch This attribute is the epoch of the equinox. That is the
value of the equinox with prefix 'J' or 'B'. The context
(a.o. keyword RADESYS) sets the prefix.
dateobs Date of observation. Floating point number given by FITS
keyword DATE-OBS
mjdobs Date of observation. Floating point number given by FITS
keyword MJD-OBS
epobs Date of observation as an epoch, i.e. copied from
mjdobs or dateobs and prefixed by 'F' or 'MJD'
========== ===============================================================
Available functions from :mod:`celestial`
.........................................
Some of the functions defined in the module :mod:`celestial` are also available in the
namespace of :mod:`wcs`. One of these is :func:`celestial.epochs` for which we wrote an example in
the previous section. Others are :func:`celestial.lon2hms`, :func:`celestial.lon2dms`
and :func:`celestial.lat2hms` to format
degrees into hours, minutes, seconds or degrees, minutes and seconds.
Finally the function :func:`celestial.skymatrix` is also available to :mod:`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 :mod:`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 :mod:`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 :class:`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 :class:`wcs.Projection` needs to know
which spectral type we want to convert to. The translation is set then with :meth:`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 | :math:`\nu` | Hz | :math:`\nu` |
+----------+----------------------------+-------------------+-----------+-------------------+
| ENER | Energy | E | J | :math:`\nu` |
+----------+----------------------------+-------------------+-----------+-------------------+
| WAVN | Wavenumber | :math:`\kappa` | 1/m | :math:`\nu` |
+----------+----------------------------+-------------------+-----------+-------------------+
| VRAD | Radio velocity | V | m/s | :math:`\nu` |
+----------+----------------------------+-------------------+-----------+-------------------+
| WAVE | Vacuum wavelength | :math:`\lambda` | m | :math:`\lambda` |
+----------+----------------------------+-------------------+-----------+-------------------+
| VOPT | Optical velocity | Z | m/s | :math:`\lambda` |
+----------+----------------------------+-------------------+-----------+-------------------+
| ZOPT | Redshift | z | \- | :math:`\lambda` |
+----------+----------------------------+-------------------+-----------+-------------------+
| AWAV | Air wavelength | :math:`\lambda_a` | m | :math:`\lambda_a` |
+----------+----------------------------+-------------------+-----------+-------------------+
| VELO | Apparent radial velocity | v | m/s | v |
+----------+----------------------------+-------------------+-----------+-------------------+
| BETA | Beta factor (v/c) | :math:`\beta` | \- | 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 :mod:`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
:meth:`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:
.. math::
:label: eq1
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:
.. math::
:label: WTeq2
\nu = \nu_{ref} + (N - N_{\nu_{ref}}) \delta \nu
where:
* :math:`\nu_{ref}` is the *reference frequency* (CRVAL1)
* :math:`N` is the pixel coordinate (FITS definition) we are interested in,
* :math:`N_{\nu_{ref}}` is the frequency reference pixel (CRPIX1)
* :math:`\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 :mod:`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 :mod:`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
:mod:`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 :mod:`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 :meth:`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 :attr:`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 :mod:`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 :attr:`wcs.Projection.altspec`.
.. note:: For a given header the attribute :attr:`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 :attr:`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`
``_ Greisen E.W. and Calabretta M.R.
.. [Ref2] `Representations of celestial coordinates in FITS`
``_ Calabretta M.R. and Greisen E.W.
.. [Ref3] `Representations of spectral coordinates in FITS`
``_ 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`
``_ FITS Working Group , Commission 5: Documentation and Astronomical Data, International Astronomical Union