1# Load the WCS information from a fits header, and use it
2# to convert pixel coordinates to world coordinates.
3
4import numpy as np
5from astropy import wcs
6from astropy.io import fits
7import sys
8
9
10def load_wcs_from_file(filename):
11    # Load the FITS hdulist using astropy.io.fits
12    hdulist = fits.open(filename)
13
14    # Parse the WCS keywords in the primary HDU
15    w = wcs.WCS(hdulist[0].header)
16
17    # Print out the "name" of the WCS, as defined in the FITS header
18    print(w.wcs.name)
19
20    # Print out all of the settings that were parsed from the header
21    w.wcs.print_contents()
22
23    # Three pixel coordinates of interest.
24    # Note we've silently assumed an NAXIS=2 image here.
25    # The pixel coordinates are pairs of [X, Y].
26    # The "origin" argument indicates whether the input coordinates
27    # are 0-based (as in Numpy arrays) or
28    # 1-based (as in the FITS convention, for example coordinates
29    # coming from DS9).
30    pixcrd = np.array([[0, 0], [24, 38], [45, 98]], dtype=np.float64)
31
32    # Convert pixel coordinates to world coordinates
33    # The second argument is "origin" -- in this case we're declaring we
34    # have 0-based (Numpy-like) coordinates.
35    world = w.wcs_pix2world(pixcrd, 0)
36    print(world)
37
38    # Convert the same coordinates back to pixel coordinates.
39    pixcrd2 = w.wcs_world2pix(world, 0)
40    print(pixcrd2)
41
42    # These should be the same as the original pixel coordinates, modulo
43    # some floating-point error.
44    assert np.max(np.abs(pixcrd - pixcrd2)) < 1e-6
45
46    # The example below illustrates the use of "origin" to convert between
47    # 0- and 1- based coordinates when executing the forward and backward
48    # WCS transform.
49    x = 0
50    y = 0
51    origin = 0
52    assert (w.wcs_pix2world(x, y, origin) ==
53            w.wcs_pix2world(x + 1, y + 1, origin + 1))
54
55
56if __name__ == '__main__':
57    load_wcs_from_file(sys.argv[-1])
58