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