Author: Hans Terlouw <gipsy@astro.rug.nl>
This Python module interfaces to Mark Calabretta’s WCSLIB and also provides a selfcontained 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, FK4NOE, FK5, ICRS and dynamic J2000, and galactic and supergalactic coordinates.
Coordinates can be represented in a number of different ways:
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).
Parameters: 


Methods:
Pixeltoworld 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”).
Worldtopixel transformation. Similar to toworld(), this method can also be called without arguments.
Simplified method for onedimensional projection objects. Its argument can be a list, a tuple, an array or a scalar. An object of the same class will be returned.
Simplified method for onedimensional projection objects. Its argument can be a list, a tuple, an array or a scalar. An object of the same class will be returned.
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 noncelestial coordinate elements any NaN value will be replaced by a value derived from the corresponding element in the other coordinate.
Returns:  a tuple (world, pixel) containing the resulting coordinates. 

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

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

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.
Pixeltogrid 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 CRPIXrelative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.
Gridtopixel 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 CRPIXrelative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.
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: 


Returns: 
This method returns a tuple with four elements:
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.
WCSLIBrelated 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.
The projection category: one of the strings undefined, zenithal, cylindrical, pseudocylindrical, conventional, conic, polyconic, quadcube, HEALPix.
A tuple with the axes’ types in the axis order of the object.
A tuple with the axes’ physical units in the axis order of the object.
A tuple with the axes’ reference values in the axis order of the object.
A tuple with the axes’ coordinate increments in the axis order of the object.
A tuple with the axes’ reference points in the axis order of the object.
A tuple with the axes’ coordinate rotations, or None if no rotations have been specified.
A NumPy matrix for the linear transformation between pixel axes and intermediate coordinate axes, or None if not specified.
A NumPy matrix for the linear transformation (with scale) between pixel axes and intermediate coordinate axes, or None if not specified.
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.
A list with charactervalued coordinate parameters. Each list element is a tuple consisting of the world coordinate axis number i, the parameter number m and the parameter value.
The native longitude of the celestial pole.
The native latitude of the celestial pole.
A fiveelement 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 XEuler angle.
The equinox (formerly ‘epoch’) of the projection.
Rest frequency in Hz.
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 readonly.
The projection’s ‘native’ sky system. E.g., (equatorial, fk5, 'J2000.0').
Alternative sky system. Can be specified according to the rules of the module celestial. See: Sky definitions. For pixeltoworld transformations, the result in the projection’s ‘native’ system is transformed to the specified one and for worldtopixel transformations, the given coordinates are first transformed to the native system, then to pixels.
Reference frame of equatorial or ecliptic coordinates: one of the (symbolic) values as defined in module celestial. E.g. icrs, fk5 or fk4.
The projection’s epoch string as derived from the attributes equinox and radesys. E.g., “B1950.0” or “J2000.0”.
The date of observation (string) as specified by the ‘DATEOBS’ key in the source object or None if not present.
The date of observation (floating point number) as specified by the ‘MJDOBS’ key in the source object or None if not present.
The date of observation as specified by either the ‘MJDOBS’ or the ‘DATEOBS’ 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.
True or False. If True, the object will use grid coordinates instead of pixel coordinates. Grid coordinates are CRPIXrelative pixel coordinates, e.g. used in GIPSY. If CRPIX is not integer, the nearest integer is used as reference.
If set to True, no exception will be raised for invalid coordinates. Invalid coordinates will be indicated by numpy.NaN (‘not a number’) values.
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.
If set to True, input and output coordinates, when given as NumPy matrices, will be row vectors instead of the standard column vectors.
Indicates whether the date of observation is to be used for the appropriate celestial transformations. True or False.
A tuple with the axes’ coordinate types (‘longitude’, ‘latitude’, ‘spectral’ or None) in the axis order of the object.
A tuple with the axes’ lengths in the axis order of the object. (Convenience attribute not directly related to WCS.)
Longitude axis number (1relative). None if not defined.
Latitude axis number (1relative). None if not defined.
Spectral axis number (1relative). None if not defined.
Convenience attribute. The object from which the Projection object was created.
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'), ('VOPTF2W', 'm/s'), ..., ('BETAF2V', '')]. If there is no spectral axis, the attribute will have the value None.
If the object was created with a call to spectra(), the argument ctype as specified in that call. Otherwise None.
The object was created with the argument minimal=True, using only the NAXIS and NAXISi items.
Example:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25  #!/bin/env python
from kapteyn import wcs
import pyfits
hdulist = pyfits.open('aurora.fits') # open 3dimensional FITS file
proj3 = wcs.Projection(hdulist[0].header) # create Projection object
pixel = ([51, 32], [17, 60], [11, 12]) # two 3dimensional 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 2dimensional 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

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


Method:
Parameters: 


Instead of calling this method, the object itself can also be called in the same way.
Attribute:
If set to True, input and output coordinates, when given as NumPy matrices, will be row vectors instead of the standard column vectors.
Example:
1 2 3 4 5 6 7 8 9 10 11 12  #!/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)

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:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25  #!/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'] = 'GLONTAN', 'GLATTAN'
# 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('ngc6946gal.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.
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.
Flexible epoch parser.
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.
Convert an angle in degrees to degrees, minutes, seconds.
Convert an angle in degrees to hours, minutes, seconds format.
Sky systems (imported from celestial)
Reference systems (imported from celestial)
Physical
Velocity of light
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]  For convenience, a slightly modified version of this module is also available in the Kapteyn Package as kapteyn.interpolation. The modification replaces NaN values in the array to a finite value in case order>1, preventing the result becoming all NaN. 