Module wcs¶
Introduction¶
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¶
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 dictionarylike object containing FITSstyle 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 CRPIXrelative 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)¶ 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”).

topixel
(world)¶ Worldtopixel transformation. Similar to
toworld()
, this method can also be called without arguments.

toworld1d
(pixel)¶ 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.

topixel1d
(world)¶ 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.

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()
ortoworld()
is performed. For noncelestial 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 nonoptimal 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 ‘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.
 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 attributealtspecarg
of the new Projection object.  axindex – Index of the spectral axis (0relative). If not specified, the first spectral axis identified by the CTYPE values of the object is assumed.
Returns: a new Projection object.  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

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

grid2pixel
(grid)¶ 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.

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

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

lonpole
¶ The native longitude of the celestial pole.

latpole
¶ The native latitude of the celestial pole.

euler
¶ 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.

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

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

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

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

dateobs
¶ The date of observation (string) as specified by the ‘DATEOBS’ key in the source object or None if not present.

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

epobs
¶ 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.

gridmode
¶ 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.

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 (1relative). None if not defined.

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

specaxnum
¶ Spectral axis number (1relative). 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'), ('VOPTF2W', 'm/s'), ..., ('BETAF2V', '')]
. 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 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
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 thescipy.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'] = '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.
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.
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]  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. 