WCSXY2SPH
Name
WCSXY2SPH
Purpose
Convert x and y (map) coordinates to spherical coordinates
Explanation
To convert x and y (map) coordinates to spherical (longitude and
latitude or sky) coordinates. This procedure is the inverse of
WCSSPH2XY.
This is a lower level procedure -- given a FITS header, the user will
usually use XYAD which will then call WCSXY2SPH with the appropriate
parameters.
Category
Mapping and Auxilary FITS Routine
Calling Sequence
wcsxy2sph, x, y, longitude, latitude, [map_type], [ CTYPE = ,$
FACE = ,PV2 = ,CRVAL =, CRXY =, LONGPOLE=, LATPOLE=]
Input Parameters
x - x coordinate of data, scalar or vector, in degrees, NOTE: x
increases to the left, not the right
y - y coordinate of data, same number of elements as x, in degrees
map_type - optional positional parameter, scalar 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/ pv2_1 = 0
SIN 3 Orthographic pv2_1, pv2_2 optional
STG 4 Stereographic AZP w/ pv2_1 = 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 inverse converges poorly
Spherical Cube
QCS 24 Quadrilateralized
Spherical Cube
TSC 25 Tangential Spherical Cube
SZP 26 Slant Zenithal perspective PV2_1,PV2_2, PV2_3 optional
HPX 27 HealPix
HCT 28 HealCart (Cartesian approximation of Healpix)
OPTIONAL 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.
FACE - a input variable used for spherical cube projections to
designate the face of the cube on which the x and y
coordinates lie. Must contain the same number of elements
as X and Y.
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 of the x
coordinate is assumed to be 0.
LATPOLE - native latitude of the standard system's North Pole
LONGPOLE - native longitude of standard system's North Pole, default
is 180 degrees, numeric scalar
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.
Output Parameters
longitude - longitude of data, same number of elements as x, in degrees
latitude - latitude of data, same number of elements as x, in degrees
Longitude and latitude will be set to NaN, wherever elements of X,Y
have no corresponding longitude, latitude values.
Notes
The conventions followed here are described in more detail in the paper
"Representations of Celestial Coordinates in FITS" by Calabretta &
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 convert x and y coordinates into a
"native" longitude and latitude and then rotate the system into one of
three generally recognized systems (celestial, galactic or ecliptic).
This procedure necessitates two basic sections. The first converts
x and y coordinates to "native" coordinates while the second converts
"native" to "standard" coordinates. The first section contains the
guts of the code in which all of the map projection is done. The
second step is performed by WCS_ROTATE and only involves rotation of
coordinate systems. WCSXY2SPH can be called in a form similar to
AITOFF, EQPOLE, or QDCB by calling wcsxy2sph with a fifth parameter
specifying the map projection by number and by not using any of the
keywords related to the map projection type (eg ctype1 and ctyp2).
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 a map projection
type). All variables are converted into double precision values.
The second task of the procedure is to take x and y coordinates and
convert them into "native" latitude and longitude coordinates.
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). WCS_ROTATE is then called to
convert the "native" coordinates to "standard" coordinates by rotating
the coordinate system. This rotation is governed by the keywords
CRVAL, and LONGPOLE. The transformation is a straightforward
application of euler angles. Finally, longitude values are converted
into the range from 0 to 360 degrees.
Common Blocks
none
Procedures Called
WCS_ROTATE
Author
Rick Balsano
Modifications/revision Level
1.1 8/31/93
1.2 9/12/93 W. Landsman Vectorized CRXY, CRVAL, CTYPE
1.3 29/12/93 I. Freedman Eliminated LU decomposition
1.4 22/09/94 W. Landsman If scalar input, then scalar output
1.5 02/03/05 W. Landsman Change variable name BETA for V4.0 compatibility
1.6 06/07/05 W. Landsman Change loop index from integer to long
1.7 02/18/99 W. Landsman Fixed implementation of ARC algorithm
1.8 June 2003 W. Landsman Update conic projections, add LATPOLE keyword
1.81 Sep 2003 W. Landsman Avoid divide by zero
1.82 Sep 2003 W. Landsman CTYPE keywords need not be 8 characters
1.83 Sep 2003 W. Landsman Preserve input array sizes
1.9 Jan 2004 W. Landsman don't modify scalars, fix PARabolic code
2.0 Feb 2004 W. Landsman Fix AIR and AZP projections
2.1 Feb 2004 W. Landsman Fix tangent projection for matrix input
3.0 May 2004 W. Landsman Support extended SIN (=NCP), slant zenithal
(SZP), and zenithal polynomial (ZPN) projections, use
PV2 keyword vector instead of PROJP1, PROJP2
3.1 May 2004 W. Landsman/J. Ballet Handle NaN values, flag invalid output
for AITOFF projection
3.1.1 Dec 2005 W. Landsman/W. Thompson Fixed problem with Airy projection
centered on 90 degree latitude
3.1.2 May 2006 W. Landsman/Y.Sato Fix problem selecting the correct root
for the ZPN projection
3.2 Aug 2007 W. Landsman Correct treatment of PVi_j parameters
3.3 Oct 2007 Sergey Koposov Support HEALPIX projection
3.4 May 2012 Benjamin Alan Weaver, Add nonstandard HEALCART
projection, Allow map_index to be > 25