Module wcs

Introduction

This Python module interfaces to Mark Calabretta’s WCSLIB and also provides a self-contained suite of celestial transformations. The WCSLIB routines “implement the FITS World Coordinate System (WCS) standard which defines methods to be used for computing world coordinates from image pixel coordinates, and vice versa.” The celestial transformations have been implemented in Python, using NumPy, and support equatorial and ecliptic coordinates of any epoch and reference systems FK4, FK4-NO-E, FK5, ICRS and dynamic J2000, and galactic and supergalactic coordinates.

Coordinates

Coordinates can be represented in a number of different ways:

• as a tuple of scalars, e.g. (ra, dec).

• as a tuple of lists or NumPy arrays, e.g. ([ra_1, ra_2, …], [dec_1, dec_2, …], [vel_1, vel_2, …]).

• as a NumPy matrix. The standard representation is a matrix with column vectors, but row vectors are also supported.

• as a NumPy array. This array can have any shape. The individual coordinate components are stored contiguously along the last axis.

• as a list of tuples. Every tuple represents one position, e.g. [(ra_1, dec_1), (ra_2, dec_2), …].

Results delivered by the transformations done by the classes described below will have the same representation as their inputs. NumPy arrays and matrices will always be returned as type ‘f8’ (64 bit).

Class Projection

class kapteyn.wcs.Projection(source[, rowvec=False, skyout=None, usedate=False, gridmode=False, alter='', minimal=False])

Parameters:

• source – a Python dictionary or dictionary-like object containing FITS-style keys and values, e.g. a header object from PyFITS.

• rowvec – indicates whether input and output coordinates, when given as NumPy matrices, will be row vectors instead of the standard column vectors. True or False.

• skyout – can be used to specify a system different from the sky system specified by the projection. This can be given as a string e.g., "equatorial fk4_no_e B1950.0" or as a tuple: (equatorial fk4_no_e 'B1950.0'). For a complete description see: Sky definitions.

• usedate – indicates whether the date of observation is to be used for the appropriate celestial transformations. True or False.

• gridmode – True or False. If True, the object will use grid coordinates instead of pixel coordinates. Grid coordinates are CRPIX-relative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.

• alter – an optional letter from ‘A’ through ‘Z’, indicating an alternative WCS axis description.

• minimal – True or False. If True, the object will be constructed from only the NAXIS and NAXISi items in the source. All other items are ignored. In this way world- and pixel coordinates will have the same values. This can be useful when it is impossible to build an object from all items in the source, e.g., when there is an error in a FITS header.

Methods:

toworld(pixel)

Pixel-to-world transformation. pixel is an object containing one or more pixel coordinates and a similar object with the corresponding world coordinates will be returned. Note that FITS images are indexed from (1,1), not from (0,0) like Python arrays. Coordinates can be specified in a number of different ways. See section Coordinates. When an exception due to invalid coordinates has occurred, this method can be called again without arguments to retrieve the result in which the invalid positions will have the value numpy.NaN (“not a number”).

topixel(world)

World-to-pixel transformation. Similar to toworld(), this method can also be called without arguments.

toworld1d(pixel)

Simplified method for one-dimensional projection objects. Its argument can be a list, a tuple, an array or a scalar. An object of the same class will be returned.

topixel1d(world)

Simplified method for one-dimensional projection objects. Its argument can be a list, a tuple, an array or a scalar. An object of the same class will be returned.

mixed(world, pixel[, span=None, step=0.0, iter=7])

Hybrid transformation.

When either the celestial longitude or latitude plus an element of the pixel coordinate is given, the remaining elements are solved by iteration on the unknown celestial coordinate. Which elements are to be solved, is indicated by assigning NaN to those elements. In case of multiple coordinates, the same elements must be indicated for every coordinate. This operation is only possible for the projection’s “native” sky system. When a different sky system has been specified, an exception will be raised. When either both celestial coordinates or both pixel coordinates are given, an operation equivalent to topixel() or toworld() is performed. For non-celestial coordinate elements any NaN value will be replaced by a value derived from the corresponding element in the other coordinate.

• span – a sequence containing the solution interval for the celestial coordinate, in degrees. The ordering of the two limits is irrelevant. Longitude ranges may be specified with any convenient normalization, for example [-120,+120] is the same as [240,480], except that the solution will be returned with the same normalization, i.e. lie within the interval specified. The default is the appropriate CRVAL value ±15°.

• step – step size for solution search, in degrees. If zero, a sensible, although perhaps non-optimal default will be used.

• iter – if a solution is not found then the step size will be halved and the search recommenced. iter controls how many times the step size is halved. The allowed range is 5 - 10.

Returns:

a tuple (world, pixel) containing the resulting coordinates.

sub([axes=None, nsub=None])

Extract a new Projection object for a subimage from an existing one.

• axes is a sequence of image axis numbers to extract. Order is significant; axes[0] is the axis number of the input image that corresponds to the first axis in the subimage, etc. If not specified, the first nsub axes are extracted.

• nsub is the number of axes to extract when axes is not specified.

Returns:

a new Projection object.

spectra(ctype[, axindex=None])

Create a new Projection object in which the spectral axis is translated. For example, a ‘FREQ’ axis may be translated into ‘ZOPT-F2W’ and vice versa. For non-standard frequency types, e.g. FREQ-OHEL as used by GIPSY, corrections are applied first to obtain barycentric frequencies. For more information, see chapter Background information spectral translations.

• ctype – Required spectral CTYPEi. Wildcarding may be used, i.e. if the final three characters are specified as ‘???’, or if just the eighth character is specified as ‘?’, the correct algorithm code will be substituted and returned. The attribute altspec provides a list of acceptable spectral types. For later reference, the value of ctype is stored in the attribute altspecarg of the new Projection object.

• axindex – Index of the spectral axis (0-relative). If not specified, the first spectral axis identified by the CTYPE values of the object is assumed.

Returns:

a new Projection object.

inside(coords, mode)

Test whether one or more coordinates are inside the area defined by the attribute naxis. This is a convenience method not directly related to WCS. coords is an object containing one or more coordinates which depending on mode can be either world- or pixel coordinates. The argument mode must be ‘world’ or ‘pixel’. The method returns a value True or False or, in the case of multiple coordinates, a list with these values.

pixel2grid(pixel)

Pixel-to-grid conversion. pixel is an object containing one or more pixel coordinates and a similar object with the corresponding grid coordinates will be returned. Grid coordinates are CRPIX-relative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.

grid2pixel(grid)

Grid-to-pixel conversion. grid is an object containing one or more grid coordinates and a similar object with the corresponding pixel coordinates will be returned. Grid coordinates are CRPIX-relative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.

str2pos(postxt[, mixpix=None])

This method accepts a string that represents one or more positions in the projection object’s coordinate system. If the string contains a valid position, the method returns the arrays with the corresponding world- and pixel coordinates. If a position could not be converted, then an error message is returned.

Parameters:

• postxt (string) – one or more positions to be parsed.

• mixpix (integer or None) – for a world coordinate system with one spatial axis, a pixel coordinate for the missing spatial axis is required to be able to convert between world- and pixel coordinates.

Returns:

This method returns a tuple with four elements:

• a NumPy array with the parsed positions in world coordinates

• a NumPy array with the parsed positions in pixel coordinates

• A tuple with the units that correspond to the axes in your world coordinate system.

• An error message when a position could not be parsed

Each position in the input string is returned in the output as an element of a numpy array with parsed positions. A position has the same number of coordinates as there are axes in the data defined by the projection object.

For its implementation, this method uses the function positions.str2pos() from module positions. Please refer to that module’s documentation for a detailed explanation.

WCSLIB-related attributes:

The following attributes contain values which are parameters to WCSLIB, after interpretation. So they can differ from the values in the source object. These attributes should not be modified.

category

The projection category: one of the strings undefined, zenithal, cylindrical, pseudocylindrical, conventional, conic, polyconic, quadcube, HEALPix.

ctype

A tuple with the axes’ types in the axis order of the object.

cunit

A tuple with the axes’ physical units in the axis order of the object.

crval

A tuple with the axes’ reference values in the axis order of the object.

cdelt

A tuple with the axes’ coordinate increments in the axis order of the object.

crpix

A tuple with the axes’ reference points in the axis order of the object.

crota

A tuple with the axes’ coordinate rotations, or None if no rotations have been specified.

pc

A NumPy matrix for the linear transformation between pixel axes and intermediate coordinate axes, or None if not specified.

cd

A NumPy matrix for the linear transformation (with scale) between pixel axes and intermediate coordinate axes, or None if not specified.

pv

A list with numeric coordinate parameters. Each list element is a tuple consisting of the world coordinate axis number i, the parameter number m and the parameter value.

ps

A list with character-valued coordinate parameters. Each list element is a tuple consisting of the world coordinate axis number i, the parameter number m and the parameter value.

lonpole

The native longitude of the celestial pole.

latpole

The native latitude of the celestial pole.

euler

A five-element list: Euler angles and associated intermediaries derived from the coordinate reference values. The first three values are the Z-, X-, and Z’-Euler angles, and the remaining two are the cosine and sine of the X-Euler angle.

equinox

The equinox (formerly ‘epoch’) of the projection.

restfrq

Rest frequency in Hz.

restwav

Vacuum rest wavelength in m.

Other Attributes:

The attributes skyout, allow_invalid, rowvec, epobs, gridmode and usedate can be modified at any time. The others are read-only.

skysys

The projection’s ‘native’ sky system. E.g., (equatorial, fk5, 'J2000.0').

skyout

Alternative sky system. Can be specified according to the rules of the module celestial. See: Sky definitions. For pixel-to-world transformations, the result in the projection’s ‘native’ system is transformed to the specified one and for world-to-pixel transformations, the given coordinates are first transformed to the native system, then to pixels.

radesys

Reference frame of equatorial or ecliptic coordinates: one of the (symbolic) values as defined in module celestial. E.g. icrs, fk5 or fk4.

epoch

The projection’s epoch string as derived from the attributes equinox and radesys. E.g., “B1950.0” or “J2000.0”.

dateobs

The date of observation (string) as specified by the ‘DATE-OBS’ key in the source object or None if not present.

mjdobs

The date of observation (floating point number) as specified by the ‘MJD-OBS’ key in the source object or None if not present.

epobs

The date of observation as specified by either the ‘MJD-OBS’ or the ‘DATE-OBS’ key in the source object or None if both are absent. This attribute is a string with the prefix ‘MJD’ or ‘F’ which can be parsed by the function epochs() in the module ‘celestial’ and consequently be part of the arguments sky_in and sky_out when creating a Transformation object.

gridmode

True or False. If True, the object will use grid coordinates instead of pixel coordinates. Grid coordinates are CRPIX-relative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.

allow_invalid

If set to True, no exception will be raised for invalid coordinates. Invalid coordinates will be indicated by numpy.NaN (‘not a number’) values.

invalid

True or False, indicating whether invalid coordinates were detected in the last transformation. In the output, invalid coordinates are indicated by numpy.NaN (‘not a number’) values.

rowvec

If set to True, input and output coordinates, when given as NumPy matrices, will be row vectors instead of the standard column vectors.

usedate

Indicates whether the date of observation is to be used for the appropriate celestial transformations. True or False.

types

A tuple with the axes’ coordinate types (‘longitude’, ‘latitude’, ‘spectral’ or None) in the axis order of the object.

naxis

A tuple with the axes’ lengths in the axis order of the object. (Convenience attribute not directly related to WCS.)

lonaxnum

Longitude axis number (1-relative). None if not defined.

lataxnum

Latitude axis number (1-relative). None if not defined.

specaxnum

Spectral axis number (1-relative). None if not defined.

source

Convenience attribute. The object from which the Projection object was created.

altspec

A list of tuples with alternative spectral types and units. The first element of such a tuple is a string with an allowed alternative spectral type which can be used as the argument of method spectra() and the second element is a string with the corresponding units. Example: [('FREQ', 'Hz'),  ('ENER', 'J'), ('VOPT-F2W', 'm/s'), ..., ('BETA-F2V', '')]. If there is no spectral axis, the attribute will have the value None.

altspecarg

If the object was created with a call to spectra(), the argument ctype as specified in that call. Otherwise None.

minimal

The object was created with the argument minimal=True, using only the NAXIS and NAXISi items.

Example:

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

hdulist = pyfits.open('aurora.fits')      # open 3-dimensional FITS file

proj3 = wcs.Projection(hdulist[0].header) # create Projection object

pixel = ([51, 32], [17, 60], [11, 12])    # two 3-dimensional pixel coordinates
world = proj3.toworld(pixel)              # transform pixel to world coordinates
print(world)
print(proj3.topixel(world))               # back from world to pixel coordinates

proj2 = proj3.sub([2,1])                  # subimage projection, axes 2 and 1

pixel = ([1, 2, 4, 3], [7, 6, 8, 2])      # four 2-dimensional pixel coordinates
world = proj2.toworld(pixel)              # transform pixel to world coordinates
print(world)

proj2.skyout = (wcs.equatorial, wcs.fk5,
                'J2008')                  # specify alternative sky system

world = proj2.toworld(pixel)              # transform to that sky system
print(world)
print(proj2.topixel(world))               # back to pixel coordinates

Class Transformation

Celestial transformations are handled by objects of the class Transformation. These objects are callable. Currently supported sky systems are equatorial and ecliptic of any epoch and galactic and supergalactic.

class kapteyn.wcs.Transformation(sky_in, sky_out[, rowvec=False])

Parameters:

• sky_out (sky_in,) – the input- and output sky system. Can be specified as e.g., “(equatorial, fk4, ‘B1950.0’)” or “galactic”.

• rowvec – if set to True, input and output coordinates, when given as NumPy matrices, will be row vectors instead of the standard column vectors.

Method:

transform(in[, reverse=False])

Parameters:

• in – an object containing one or more coordinates to be transformed and out will receive a similar object with the transformed coordinates. Coordinates can be specified in a number of different ways. See section Coordinates.

• reverse – if True, the inverse transformation will be performed.

Instead of calling this method, the object itself can also be called in the same way.

Attribute:

rowvec

If set to True, input and output coordinates, when given as NumPy matrices, will be row vectors instead of the standard column vectors.

Example:

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

tran = wcs.Transformation((wcs.equatorial, wcs.fk4, 'B1950.0'), wcs.galactic)

radec = numpy.matrix(([33.3, 177.2, 230.1],
                      [66.2, -11.5,  13.0]))

lbgal = tran(radec)
print(lbgal)
print(tran(lbgal, reverse=True))

Functions

Function coordmap

kapteyn.wcs.coordmap(proj_src, proj_dst[, dst_shape=None, dst_offset=None, src_offset=None])

• proj_src, proj_dst – the source- and destination projection objects.

• dst_shape – the destination image’s shape. Must be compatible with the projections’ dimensionality. The elements are in Python order, i.e., the first element corresponds to the last FITS axis. If dst_shape is None (the default), the shape is derived from the proj_dst.naxis attribute.

• dst_offset – the destination image’s offset. If None, the offset for all axes will be zero. Otherwise it must be compatible with the projections’ dimensionality. The elements are in Python order, i.e., the first element corresponds to the last FITS axis.

• src_offset – the source image’s offset. If None, the offset for all axes will be zero. Otherwise it must be compatible with the projections’ dimensionality. The elements are in Python order, i.e., the first element corresponds to the last FITS axis.

This function returns a coordinate map which can be used as the argument coordinates in calls to the function map_coordinates() from the scipy.ndimage.interpolation module. [1] The resulting coordinate map can be used for reprojecting an image into another image with a different coordinate system.

Example:

#!/bin/env python
from kapteyn import wcs
import numpy, pyfits
from kapteyn.interpolation import map_coordinates

hdulist = pyfits.open('ngc6946.fits')
header = hdulist[0].header

proj1 = wcs.Projection(header)                       # source projection
trans = wcs.Transformation(proj1.skysys, skyout=wcs.galactic)

header['CTYPE1'], header['CTYPE2'] = 'GLON-TAN', 'GLAT-TAN'
                                                     # new axis types
header['CRVAL1'], header['CRVAL2'] = trans((header['CRVAL1'],header['CRVAL2']))
                                                     # new reference point

proj2 = wcs.Projection(header)                       # destination projection

coords = wcs.coordmap(proj1, proj2)

image_in = hdulist[0].data
image_out = map_coordinates(image_in, coords, order=1, cval=numpy.NaN)

hdulist[0].data = image_out
hdulist.writeto('ngc6946-gal.fits')

This example is a complete program and illustrates how a FITS file containing an image with arbitrary coordinates can be reprojected into an image with galactic coordinates. The image can have two or more dimensions.

Utility functions

The following are functions from the module celestial which have been made available within the namespace of this wcs module: For detailed information, refer to celestial’s documentation.

kapteyn.wcs.epochs(spec)

Flexible epoch parser.

kapteyn.wcs.lat2dms(a[, prec=1])

Convert an angle in degrees into the degrees, minutes, seconds format assuming it was a latitude of which the value should be in the range -90 to 90 degrees.

kapteyn.wcs.lon2dms(a[, prec=1])

Convert an angle in degrees to degrees, minutes, seconds.

kapteyn.wcs.lon2hms(a[, prec=1])

Convert an angle in degrees to hours, minutes, seconds format.

Constants

Sky systems (imported from celestial)

kapteyn.wcs.equatorial

kapteyn.wcs.eq

kapteyn.wcs.ecliptic

kapteyn.wcs.ecl

kapteyn.wcs.galactic

kapteyn.wcs.gal

kapteyn.wcs.supergalactic

kapteyn.wcs.sgal

Reference systems (imported from celestial)

kapteyn.wcs.fk4

kapteyn.wcs.fk4_no_e

kapteyn.wcs.fk5

kapteyn.wcs.icrs

kapteyn.wcs.dynj2000

kapteyn.wcs.j2000

Physical

kapteyn.wcs.c

Velocity of light

Error handling

Errors are reported through the exception mechanism. Two exception classes have been defined: WCSerror for unrecoverable errors and WCSinvalid for situations where a partial result may be available.

Footnotes

[1]

In version 3.4 of the Kapteyn Package we only use the plain SciPy routines i.e. without modifications that were applied in 3.0. This implies that the interpolation is more sensitive to NaN’s in the data.