WCSSPH2XY
Name
WCSSPH2XY
Purpose
Convert spherical coordinates to x and y (map) angular coordinates
Explanation
Convert spherical (longitude and latitude -- sky) coordinates to x
and y (map) angular coordinates. This procedure is the inverse of
WCSXY2SPH. See WCS_DEMO for example of use.
This is a lower level procedure -- given a FITS header, the user will
usually use ADXY which will then call WCSSPH2XY with the appropriate
parameters.
Category
Mapping and Auxiliary FITS Routine
Calling Sequence
wcssph2xy, longitude, latitude, x, y, [ map_type , CTYPE = ,
FACE =,PV2= , CRVAL = , CRXY = , LONGPOLE = ,
LATPOLE = , NORTH_OFFSET =, SOUTH_OFFSET =, BADINDEX =]
Input Parameters
longitude - longitude of data, scalar or vector, in degrees
latitude - latitude of data, same number of elements as longitude,
in degrees
map_type - optional positional parameter, numeric scalar (0-26)
corresponding to a particular map projection. This is not a
FITS standard, it is simply put in to allow function similar
to that of less general map projection procedures (eg AITOFF).
The following list gives the map projection types and their
respective numbers.
FITS Number Name Comments
code code
---- ------ ----------------------- -----------------------------------
DEF 0 Default = Cartesian
AZP 1 Zenithal perspective PV2_1 required
TAN 2 Gnomic AZP w/ mu = 0
SIN 3 Orthographic PV2_1,PV2_2 optional
STG 4 Stereographic AZP w/ mu = 1
ARC 5 Zenithal Equidistant
ZPN 6 Zenithal polynomial PV2_0, PV2_1....PV2_20 possible
ZEA 7 Zenithal equal area
AIR 8 Airy PV2_1 required
CYP 9 Cylindrical perspective PV2_1 and PV2_2 required
CAR 10 Cartesian
MER 11 Mercator
CEA 12 Cylindrical equal area PV2_1 required
COP 13 Conical perspective PV2_1 and PV2_2 required
COD 14 Conical equidistant PV2_1 and PV2_2 required
COE 15 Conical equal area PV2_1 and PV2_2 required
COO 16 Conical orthomorphic PV2_1 and PV2_2 required
BON 17 Bonne's equal area PV2_1 required
PCO 18 Polyconic
SFL 19 Sanson-Flamsteed
PAR 20 Parabolic
AIT 21 Hammer-Aitoff
MOL 22 Mollweide
CSC 23 Cobe Quadrilateralized convergence of inverse is poor
Spherical Cube
QSC 24 Quadrilateralized
Spherical Cube
TSC 25 Tangential Spherical Cube
SZP 26 Slant Zenithal Projection PV2_1,PV2_2, PV2_3 optional
HPX 27 HealPix
HCT 28 HealCart (Cartesian approximation of Healpix)
Optional Input Keyword Parameters
CTYPE - One, two, or three element vector containing 8 character
strings corresponding to the CTYPE1, CTYPE2, and CTYPE3
FITS keywords:
CTYPE[0] - first four characters specify standard system
('RA--','GLON' or 'ELON' for right ascension, Galactic
longitude or ecliptic longitude respectively), second four
letters specify the type of map projection (eg '-AIT' for
Aitoff projection)
CTYPE[1] - first four characters specify standard system
('DEC-','GLAT' or 'ELAT' for declination, galactic latitude
or ecliptic latitude respectively; these must match
the appropriate system of ctype1), second four letters of
ctype2 must match second four letters of ctype1.
CTYPE[2] - if present must be the 8 character string,'CUBEFACE',
only used for spherical cube projections to identify an axis
as containing the face on which each x and y pair of
coordinates lie.
PV2 - Vector of projection parameter associated with latitude axis
PV2 will have up to 21 elements for the ZPN projection, up to 3
for the SIN projection and no more than 2 for any other
projection. The first element corresponds to PV2_1, the
second to PV2_2, etc.
CRVAL - 2 element vector containing standard system coordinates (the
longitude and latitude) of the reference point
CRXY - 2 element vector giving the x and y coordinates of the
reference point, if this is not set the offset is [0,0]
This is not a FITS standard -- it is similar to CRPIX but in
angular X,Y coordinates (degrees) rather than pixel coordinates
LATPOLE - native latitude of the standard system's North Pole
LONGPOLE - native longitude of standard system's North Pole, default
is 180 degrees for Zenithal systems
NORTH_OFFSET - offset (radians) added to input points near north pole.
SOUTH_OFFSET - offset (radians) added to input points near south pole.
BADINDEX - vector, list of transformed points too close to poles.
Output Parameters
x - x coordinate of data, same number of elements as longitude, in
degrees; if CRXY is set, then x will be returned offset by
crxy(0). NOTE: x in all map projections increases to the
left, not the right.
y - y coordinate of data, same number of elements as longitude, in
degrees; if CRXY is set, y will be returned offset by crxy[1]
bad - vector returning index to transformed points close to pole.
OPTIONAL OUTPUT KEYWORD PARAMETERS:
FACE - a output variable used for spherical cube projections to
designate the face of the cube on which the x and y
coordinates lie. Will contain the same number of elements as
X and Y. Must contain at least 1 arbitrary element on input
If FACE is NOT defined on input, it is assumed that the
spherical cube projection is laid out over the whole sky
in the "sideways T" configuration.
Notes
The conventions followed here are described in more detail in
"Representations of Celestial Coordinates in FITS" by Calabretta
and Greisen (2002, A&A, 395, 1077; also see
http://fits.gsfc.nasa.gov/fits_wcs.html). The general
scheme outlined in that article is to first use WCS_ROTATE to convert
coordinates in one of three standard systems (celestial, galactic,
or ecliptic) into a "native system" of latitude and longitude. The
latitude and longitude are then converted into x and y coordinates
which depend on the map projection which is performed. The rotation
from standard to native coordinates can be skipped if one so desires.
This procedure necessitates two basic sections. The first converts
"standard" coordinates to "native" coordinates while the second converts
"native" coordinates to x and y coordinates. The first section is
simply a call to WCS_ROTATE, while the second contains the guts of
the code in which all of the map projection is done. This procedure
can be called in a form similar to AITOFF, EQPOLE, or QDCB by calling
wcssph2xy with a fifth parameter specifying the map projection by
number and by not using any of the keywords related to the map
projection type (e.g. CTYPE).
Procedure
The first task of the procedure is to do general error-checking to
make sure the procedure was called correctly and none of the
parameters or keywords conflict. This is particularly important
because the procedure can be called in two ways (either using
FITS-type keywords or using a number corresponding to a map projection
type). All variables are converted into double precision values and
angular measurements are converted from degrees into radians.
If necessary, longitude values are converted into the range -pi to pi.
Any latitude points close to the of the poles are mapped to a specific
latitude of from the pole so that the map transformations become
completely invertible. The magnitude of this correction is given by
the keywords NORTH_OFFSET and SOUTH_OFFSET and a list of affected
points is optionally returned in the "badindex" output parameter.
The next task of the procedure is to convert the "standard"
coordinates to "native" coordinates by rotating the coordinate system.
This rotation is performed by the procedure WCS_ROTATE and is governed
by the keywords CRVAL and LONGPOLE. The final task of the WCSSPH2XY
is to take "native" latitude and longitude coordinates and convert
them into x and y coordinates. Any map specific error-checking is
done at this time. All of the equations were obtained from
"Representations of Celestial Coordinates in FITS" and cases needing
special attention are handled appropriately (see the comments with
individual map projections for more information on special cases).
Note that a further transformation (using the CD matrix) is required
to convert the (x,y) coordinates to pixel coordinates.
Common Blocks
none
Procedures Called
WCS_ROTATE
Author
Rick Balsano
Modifications/revision Level
1.1 8/31/93
2.3 9/15/93 W. Landsman (HSTX) Update quad cube coords, vectorize
Keywords
2.4 12/29/93 I. Freedman (HSTX) Eliminated LU decomposition
2.5 1/5/93 I. Freedman (HSTX) Offset keywords / bad point index
2.6 Dec 94 Compute pole for transformations where the reference
pixel is at the native origin W. Landsman (HSTX)
2.7 May 95 Change internal variable BETA for V4.0 compatibility
2.8 June 95 Change loop indices from integer to long
2.9 3/18/96 Change FACE usage for cube projections to match WCSLIB
C/FORTRAN software library.
2.10 02/18/99 Fixed implementation of ARC algorithm
2.11 June 2003 Update conic projections, add LATPOLE keyword
2.12 Aug 2003, N.Rich - Fix pre-V5.5 bug from previous update
2.13 Sep 2003, W. Landsman CTYPE keywords need not be 8 characters
2.14 Jan 2004, W. Landsman don't modify scalars, fix PARabolic code
2.15 Feb 2004, W. Landsman Fix AZP and AIR algorithms
3.0 May 2004 W. Landsman Support extended SIN (=NCP), slant zenithal
(SZP), and zenithal polynomail (ZPN) projections, use
PV2 keyword vector instead of PROJP1, PROJP2
3.1 Jul 2005 W.Landsman/C. Markwardt Set unprojectable points in
tangent projection to NaN
3.1.1 Jul 2005 Fixed 3.1 mod to work for scalars
3.2 Dec 2005 Fixed Airy projection for latitude centered at 90 deg
3.3 Aug 2007 R. Munoz, W.Landsman Correct treatment of PV1_2 and
PV2_2 parameters
3.4 Oct 2007 Sergey Koposov Support HEALPIX projection
3.4.1 June 2009 Check for range of validity of ZPN polynomial W.L.
3.5 May 2012 Benjamin Alan Weaver, Add nonstandard HEALCART
projection, Allow map_index to be > 25