Skip to content
Back to top
View this page

# FITS Equivalent WCS Example

The following example shows how to construct a GWCS object equivalent to a FITS imaging WCS without distortion, defined in this FITS imaging header:

    WCSAXES =                    2 / Number of coordinate axes
    WCSNAME = '47 Tuc     '        / Coordinate system title
    CRPIX1  =               2048.0 / Pixel coordinate of reference point
    CRPIX2  =               1024.0 / Pixel coordinate of reference point
    PC1_1   =   1.290551569736E-05 / Coordinate transformation matrix element
    PC1_2   =  5.9525007864732E-06 / Coordinate transformation matrix element
    PC2_1   =  5.0226382102765E-06 / Coordinate transformation matrix element
    PC2_2   = -1.2644844123757E-05 / Coordinate transformation matrix element
    CDELT1  =                  1.0 / [deg] Coordinate increment at reference point
    CDELT2  =                  1.0 / [deg] Coordinate increment at reference point
    CUNIT1  = 'deg'                / Units of coordinate increment and value
    CUNIT2  = 'deg'                / Units of coordinate increment and value
    CTYPE1  = 'RA---TAN'           / TAN (gnomonic) projection + SIP distortions
    CTYPE2  = 'DEC--TAN'           / TAN (gnomonic) projection + SIP distortions
    CRVAL1  =        5.63056810618 / [deg] Coordinate value at reference point
    CRVAL2  =      -72.05457184279 / [deg] Coordinate value at reference point
    LONPOLE =                180.0 / [deg] Native longitude of celestial pole
    LATPOLE =      -72.05457184279 / [deg] Native latitude of celestial pole
    RADESYS = 'ICRS'                / Equatorial coordinate system

For this example the following imports are needed:

    >>> import numpy as np
    >>> from astropy.modeling import models
    >>> from astropy import coordinates as coord
    >>> from astropy import units as u
    >>> from gwcs import wcs
    >>> from gwcs import coordinate_frames as cf

The `forward_transform` is constructed as a combined model using `astropy.modeling`. The `frames` are subclasses of `CoordinateFrame`. Although strings are acceptable as `coordinate_frames` it is recommended this is used only in testing/debugging.

Using the `modeling` package create a combined model to transform detector coordinates to ICRS following the FITS WCS standard convention.

First, create a transform which shifts the input `x` and `y` coordinates by `CRPIX`. We subtract 1 from the CRPIX values because the first pixel is considered pixel `1` in FITS WCS:

    >>> shift_by_crpix = models.Shift(-(2048 - 1)*u.pix) & models.Shift(-(1024 - 1)*u.pix)

Create a transform which rotates the inputs using the `PC matrix`.

    >>> matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
    ...                    [5.0226382102765E-06 , -1.2644844123757E-05]])
    >>> rotation = models.AffineTransformation2D(matrix * u.deg,
    ...                                          translation=[0, 0] * u.deg)
    >>> rotation.input_units_equivalencies = {"x": u.pixel_scale(1*u.deg/u.pix),
    ...                                       "y": u.pixel_scale(1*u.deg/u.pix)}
    >>> rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix,
    ...                                                  translation=[0, 0] * u.pix)
    >>> rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1*u.pix/u.deg),
    ...                                               "y": u.pixel_scale(1*u.pix/u.deg)}

Create a tangent projection and a rotation on the sky using `CRVAL`.

    >>> tan = models.Pix2Sky_TAN()
    >>> celestial_rotation =  models.RotateNative2Celestial(5.63056810618*u.deg, -72.05457184279*u.deg, 180*u.deg)

    >>> det2sky = shift_by_crpix | rotation | tan | celestial_rotation
    >>> det2sky.name = "linear_transform"

Create a `detector` coordinate frame and a `celestial` ICRS frame.

    >>> detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"),
    ...                             unit=(u.pix, u.pix))
    >>> sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs',
    ...                               unit=(u.deg, u.deg))

This WCS pipeline has only one step - from `detector` to `sky`:

    >>> pipeline = [(detector_frame, det2sky),
    ...             (sky_frame, None)
    ...            ]
    >>> wcsobj = wcs.WCS(pipeline)
    >>> print(wcsobj)
      From      Transform
    -------- ----------------
    detector linear_transform
        icrs             None

Now we have a complete WCS object. The next example will use it to convert pixel coordinates(1, 2) to sky coordinates:

    >>> sky = wcsobj(1*u.pix, 2*u.pix)
    >>> print(sky)
    (<Quantity 5.52515954 deg>, <Quantity -72.05190935 deg>)

The `invert()` method evaluates the `backward_transform()` to provide a mapping from sky coordinates to pixel coordinates if available, otherwise it applies an iterative method to calculate the pixel coordinates.

    >>> wcsobj.invert(*sky)
    (<Quantity 1. pix>, <Quantity 2. pix>)