NOAO IRAF ZPX World Coordinate System

The ZPX World Coordinate System is a non-standard system for evaluating celestial coordinates from the image pixel coordinates. It follows the the FITS conventions for a zenithal polynomial projection (ZPN) but adds an additional two-dimensional polynomial distortion term to the evaluation. This discussion concentrates on the non-ZPN distortion extension and assumes the reader understands the FITS WCS conventions including applying a zenithal polynomial projection. The reference for the FITS WCS standard for undistorted celestial coordinates systems is Representations of celestial coordinates in FITS Calabretta, M. R., and Greisen, E. W., Astronomy & Astrophysics, 395, 1077-1122, 2002. (PDF, HTML). (These links to the publisher's web site are currently available by subscription only. Reprints are available from the author's web site in PDF format.)

One thing to note is that generally the radial "pin-cushion" distortion included in the ZPN projection is generally fixed and possibly defined by known optical design terms. For NOAO Mosaic data, each WCS calibration exposure has different distortion terms from fitting the refraction and other small effects but the calibration does not change the radial polynomial coefficients.

The ZPN World Coordinate System projection has a FITS keyword representation as illustrated in figure 1.

Figure 1: Sample header with ZPX WCS projection

WCSAXES =                    2 / Number of WCS axes
CTYPE1  = 'RA---ZPX'           / Coordinate type
CTYPE2  = 'DEC--ZPX'           / Coordinate type
CRVAL1  =     320.687374999995 / Coordinate reference value
CRVAL2  =      36.908555555556 / Coordinate reference value
CRPIX1  =     4167.56175625891 / Coordinate reference pixel
CRPIX2  =     4120.25894749731 / Coordinate reference pixel
CD1_1   =  -5.2588308681025E-8 / Coordinate matrix
CD2_1   =  -7.2772379161132E-5 / Coordinate matrix
CD1_2   =  -7.2753930850119E-5 / Coordinate matrix
CD2_2   =  -1.8637632244742E-8 / Coordinate matrix
WAT0_001= 'system=image'       / Coordinate system
WAT1_001= 'wtype=zpx axtype=ra projp0=0. projp1=1. projp2=0. projp3=337.74 proj'
WAT1_002= 'p4=0. projp5=632052. lngcor = "3. 3. 3. 2. 0.001876397956622823 0.29'
WAT1_003= '99113930557312 0.1542460039112511 0.3032873851581314 1.9247409545894'
WAT1_004= '95E-5 -1.348328290485618E-5 1.414186703253352E-4 -1.792784764381400E'
WAT1_005= '-4 -1.276226238774833E-4 4.339217671825231E-4 "'
WAT2_001= 'wtype=zpx axtype=dec projp0=0. projp1=1. projp2=0. projp3=337.74 pro'
WAT2_002= 'jp4=0. projp5=632052. latcor = "3. 3. 3. 2. 0.001876397956622823 0.2'
WAT2_003= '999113930557312 0.1542460039112511 0.3032873851581314 9.963957331149'
WAT2_004= '402E-5 -1.378185066830135E-4 1.559892401479664E-4 -8.280442729203771'
WAT2_005= 'E-4 3.966701903249366E-4 0.001678960379199465 "'

The WCSAXES keyword (possible seen as WCSDIM) will always be 2. That this is a ZPX projection is indicated by the CTYPE keywords. These keywords also indicate that the first image axis corresponds to RA and the second to DEC.

The ZPX projection is evaluated as follows.

Compute the first order standard coordinates xi and eta from the linear part of the solution stored in CRPIX and the CD matrix.
```
         xi = CD1_1 * (x - CRPIX1) + CD1_2 * (y - CRPIX2)
        eta = CD2_1 * (x - CRPIX1) + CD2_2 * (y - CRPIX2)
```
Add the non-linear part of the projection using the coefficients in the WAT keywords as described below.
```
        xi' =  xi + lngcor (xi, eta)
       eta' = eta + latcor (xi, eta)
```
Apply the zenithal polynomial projection to xi' and eta' as described in Calabretta and Greisen where the Pm coefficients in that paper are are given by the projpm coefficients in the WAT keywords. The reference tangent point given by the CRVAL values lead to the final RA and DEC in degrees. Note that the units of xi, eta, lngcor, and latcor are also degrees.

The coefficients for the zenithal polynomial (the Pm) and the non-linear polynomial distortion functions lngcor(xi,eta) and latcor(xi,eta) are stored as FITS keywords under the indexed WATj_nnn keywords. The j refers to the image axis and the nnn give a sequence number. The cards for a particular image axis are sorted by the sequence number and then concatenated together into one long string. Take care not to add spaces between the concatenated strings since the coefficients may be split across strings.

The long string for each image axis is composed of a set of keyword/value pairs where the value is quoted if it contains whitespace. Figure 2 shows how the WAT keywords in figure 1 would be decomposed into parameters and coefficients.

Figure 2: Decomposing the WAT keywords from figure 1

    AXIS 1			    AXIS 2
    ----------------------------    --------------------------
    wtype=zpx			    wtype=zpx
    axtype=ra			    axtype=dec
    projp0=0.                       projp0=0.
    projp1=1.                       projp1=1.
    projp2=0.                       projp2=0.
    projp3=337.74                   projp3=337.74
    projp4=0.                       projp4=0.
    projp5=632052.                  projp5=632052.
    lngcor=                         latcor=
            3.                              3.
	    3.                              3.
	    3.                              3.
	    2.                              2.
	    0.001876397956622823            0.001876397956622823
	    0.2999113930557312              0.2999113930557312
            0.1542460039112511              0.1542460039112511
	    0.3032873851581314              0.3032873851581314
	    1.924740954589495E-5            9.963957331149402E-5
	   -1.348328290485618E-5           -1.378185066830135E-4
	    1.414186703253352E-4            1.559892401479664E-4
	   -1.792784764381400E-4           -8.280442729203771E-4
	   -1.276226238774833E-4            3.966701903249366E-4
	    4.339217671825231E-4            0.001678960379199465

The list of coefficients are interpreted as follows.

The first number is the function type encoded as 1=chebyshev, 2=legendre, 3=polynomial. The example has a function of type 3 which is the simple polynomial.
The next two numbers represent the "order" of the function in xi and eta. The order is the one less than the highest polynomial power. The powers are represented below by m and n such at m = 0 to xiorder-1 and n = 0 to etaorder-1. In the example the orders are 4 which means cubic polynomials (m=0 to 3 and n=0 to 3).
The next (fourth) number specifies the type of cross-terms encoded as 0=no cross-terms, 1=full cross-terms, 2=half-cross terms. The cross-terms are terms of xi^m*eta^n where m and n are non-zero. Full cross-terms mean that both m and n will go to the their maximum values independently while half-cross terms mean that m + n will only go to the maximum of xiorder-1 and etaorder-1.
The next 4 numbers describe the region of validity of the fits in xi and eta space, e.g. ximin, ximax, etamin, etamax. They are used to compute normalized values for xi and eta used in the chebyshev and legendre polynomial functions:
```
    xin = (2 * xi - (ximax + ximin)) / (ximax - ximin)
    etan = (2 * eta - (etamax + etamin)) / (etamax - etamin)
```

The remaining terms are the coefficients of the polynomial terms. The functions are evaluated by summing polynomial terms Pmn(xi,eta) multiplied by the coefficients Cmn as

     lngcor(xi,eta) = sum (Cmn * Pmn(xi,eta))
     latcor(xi,eta) = sum (Cmn * Pmn(xi,eta))

Representing the coeffients as Cmn for the polynomials Pmn, where m and n are the powers of xi and eta, they are ordered as

In the example with the half cross-terms and orders of 4 the ten coefficients would be C00, C10, C20, C30, C01, C11, C21, C02, C12, and C03.

The polynomials Pmn are defined below. The chebyshev and legendre polynomials are define iteratively as functions of the normalized coordinates defined earlier.

    Pmn = xi ** m * eta ** n    (polynomial)

    Pmn = Pm(xin) * Pn(etan)     (chebyshev)

	P0(xin) = 1.0
	P1(xin) = xin
	Pm+1(xin) = 2.0 * xin * Pm(xin) - Pm-1(xin) 

	P0(etan) = 1.0
	P1(etan) = etan
	Pn+1(etan) = 2.0 * etan * Pn(etan) - Pn-1(etan) 

    Pmn = Pm(xin) * Pn(etan)     (legendgre)

	P0(xin) = 1.0
	P1(xin) = xin
	Pm+1(xin) = ((2m+1) * xin * Pm(xin) - m * Pm-1(xin))/ (m+1)   

	P0(etan) = 1.0
	P1(etan) = etan
	Pn+1(etan) = ((2n+1) * etan * Pn(etan) - n * Pn-1(etan))/ (n+1)

In the example with with a simple polynomial the functions would be given as follows.

    lngcor/latcor = C00
        + C10 * xi     + C20*xi^2       + C30 * xi^3
	+ C01 * eta    + C02*eta^2      + C03 * eta^3
	+ C11 * xi*eta + C21 * xi^2*eta + C12 * xi*eta^2

(copied from NOAO