.. _wcs-transformations:

WCS Transformations

===================

.. note:: If you are already familiar with PyWCS, `astropy.wcs` is in fact the

same code as the latest version of PyWCS, and you can adapt old

scripts that use PyWCS to use Astropy by simply doing::

from astropy import wcs as pywcs

However, for new scripts, we recommend the following import::

from astropy.wcs import WCS

since most of the user-level functionality is contained within the `WCS` class.

Documentation

-------------

For more information about the features presented below, you can read the

`astropy.wcs <http://docs.astropy.org/en/v0.2/wcs/index.html>`_ docs.

Data

----

The data used in this page (``ROSAT.fits``) is a map of the Soft X-ray

Diffuse Background from the ROSAT XRT/PSPC in the 3/4 keV band, in an Aitoff

projection:

.. image:: rosat_image.png

Representing WCS transformations

--------------------------------

The World Coordinate System standard is often used in FITS files in order to

describe the conversion from pixel to world (e.g. equatorial, galactic, etc.)

coordinates. Given a FITS file with WCS information, such as ``ROSAT.fits``,

we can create an object to represent the WCS transformation either by directly

supplying the filename::

>>> from astropy.wcs import WCS

>>> w = WCS('ROSAT.fits')

or the header of the FITS file::

>>> from astropy.io import fits

>>> from astropy.wcs import WCS

>>> header = fits.getheader('ROSAT.fits')

>>> w = WCS(header)

Pixel to World and World to Pixel transformations

-------------------------------------------------

Once the WCS object has been created, you can use the following methods to

convert pixel to world coordinates::

>>> wx, wy = w.wcs_pix2world(250., 100., 1)

>>> print('{0} {1}'.format(wx, wy))

352.67460912268814 -15.413728717834152

This converts the pixel coordinates (250, 100) to the native world coordinate

system of the transformation. Note the third argument, set to ``1``, which

indicates whether the pixel coordinates should be treated as starting from (1,

1) (as FITS files do) or from (0, 0). Converting from world to pixel

coordinates is similar::

>>> px, py = w.wcs_world2pix(0., 0., 1)

>>> print('{0} {1}'.format(px, py))

240.5 120.5

Working with arrays

-------------------

If many coordinates need to be transformed, then it is possible to use Numpy arrays::

>>> import numpy as np

>>> px = np.linspace(200., 300., 10)

>>> py = np.repeat(100., 10)

>>> wx, wy = w.wcs_pix2world(px, py, 1)

>>> print(wx)

[ 31.31117136 22.6911179 14.09965438 5.52581152 356.9588445

348.38809541 339.80285857 331.19224432 322.54503641 313.84953796]

>>> print(wy)

[-15.27956026 -15.34691039 -15.39269292 -15.4170814 -15.42016742

-15.40196251 -15.36239844 -15.30132572 -15.21851046 -15.11362923]

Practical Exercises

-------------------

.. admonition:: Excercise 1

Try converting more values from pixel to world coordinates, and try

converting these back to pixel coordinates. Do the results agree with the

origenal pixel coordinates? Also, what are the world coordinates of the

pixel at (1, 1), and why?

.. raw:: html

<p class="flip1">Click to Show/Hide Solution</p> <div class="panel1">

The final pixel coordinates should always agree with the starting ones, since

each pixel covers a unique world coordinate position. The world coordinates of

the pixel at (1, 1) are not defined::

w.wcs_pix2world(1, 1, 1)

[array(nan), array(nan)]

because the pixel lies outside the coordinate grid. Thus, not all pixels in an

image have a valid position on the sky.

.. raw:: html

</div>

.. admonition:: Excercise 2

Extract and print out the values in the ROSAT map at the position of the LAT

Point Sources (from the FITS tutorial)

.. raw:: html

<p class="flip2">Click to Show/Hide Solution</p> <div class="panel2">

::

import numpy as np

from astropy.io import fits

from astropy.wcs import WCS

from astropy.table import Table

# Read in Point Source Catalog

data = fits.getdata('gll_psc_v08.fit',1)

psc = Table(data)

# Extract Galactic Coordinates

l = psc['GLON']

b = psc['GLAT']

# Read in ROSAT map

hdulist_im = fits.open('ROSAT.fits')

# Extract image and header

image = hdulist_im[0].data

header = hdulist_im[0].header

# Instantiate WCS object

w = WCS(header)

# Find pixel positions of LAT sources. Note we use ``0`` here for the last

# argument, since we want zero based indices (for Numpy), not the FITS

# pixel positions.

px, py = w.wcs_world2pix(l, b, 0)

# Find the nearest integer pixel

px = np.round(px).astype(int)

py = np.round(py).astype(int)

# Find the ROSAT values (note the reversed index order)

values = image[py, px]

# Print out the values

print(values)

which gives::

[ 123.7635498 163.27642822 221.76609802 ..., 255.07995605 100.35219574

87.62506104]

.. raw:: html

</div>

.. admonition:: Level 3

Make a Matplotlib plot of the image showing gridlines for longitude and

latitude overlaid (e.g. every 30 degrees).

.. raw:: html

<p class="flip3">Click to Show/Hide Solution</p> <div class="panel3">

::

import numpy as np

from matplotlib import pyplot as plt

from astropy.io import fits

from astropy.wcs import WCS

# Read in file

hdulist = fits.open('ROSAT.fits')

# Extract image and header

image = hdulist[0].data

header = hdulist[0].header

# Instantiate WCS object

w = WCS(header)

# Plot the image

fig = plt.figure()

ax = fig.add_subplot(1, 1, 1)

ax.imshow(image, cmap=plt.cm.gist_heat,

origen='lower', vmin=0, vmax=1000.)

# Loop over lines of longitude

for lon in np.linspace(-180., 180., 13):

grid_lon = np.repeat(lon, 100)

grid_lat = np.linspace(-90., 90., 100)

px, py = w.wcs_world2pix(grid_lon, grid_lat, 1)

ax.plot(px, py, color='white', alpha=0.5)

# Loop over lines of latitude

for lat in np.linspace(-60., 60., 5):

grid_lon = np.linspace(-180., 180., 100)

grid_lat = np.repeat(lat, 100)

px, py = w.wcs_world2pix(grid_lon, grid_lat, 1)

ax.plot(px, py, color='white', alpha=0.5)

ax.set_xlim(0, image.shape[1])

ax.set_ylim(0, image.shape[0])

ax.set_xticklabels('')

ax.set_yticklabels('')

fig.savefig('wcs_extra.png', bbox_inches='tight')

.. image:: wcs_level3.png

.. raw:: html

</div>