Skip to content
Back to top
View this page

# Using the WCS object

This section uses the `imaging_wcs_wdist.asdf` created in Adding distortion to the imaging example to read in a WCS object and demo its methods.

    >>> import asdf
    >>> asdf_file = asdf.open("imaging_wcs_wdist.asdf")
    >>> wcsobj = asdf_file.tree["wcs"]
    >>> print(wcsobj)
             From          Transform
    ----------------- ----------------
             detector       distortion
    undistorted_frame linear_transform
                 icrs             None

## Inspecting Available Coordinate Frames

To see what frames are defined:

    >>> print(wcsobj.available_frames)
    ['detector', 'undistorted_frame', 'icrs']
    >>> wcsobj.input_frame
    <Frame2D(name="detector", unit=(Unit("pix"), Unit("pix")), axes_names=('x', 'y'), axes_order=(0, 1))>
    >>> wcsobj.output_frame
    <CelestialFrame(name="icrs", unit=(Unit("deg"), Unit("deg")), axes_names=('lon', 'lat'), axes_order=(0, 1), reference_frame=<ICRS Frame>)>

Because the `output_frame` is a `CoordinateFrame` object we can get the result of the WCS transform as an `SkyCoord` object and transform them to other standard coordinate frames supported by `astropy.coordinates`.

    >>> skycoord = wcsobj.pixel_to_world(1, 2)
    >>> print(skycoord)
    <SkyCoord (ICRS): (ra, dec) in deg
        (5.50090023, -72.04553535)>
    >>> print(skycoord.transform_to("galactic"))
    <SkyCoord (Galactic): (l, b) in deg
        (306.12713109, -44.8996588)>

## Using Bounding Box

The WCS object has an attribute `bounding_box` (default value of `None`) which describes the range of acceptable values for each input axis.

    >>> wcsobj.bounding_box = ((0, 2048), (0, 1000))
    >>> wcsobj((2,3), (1020, 980))
    [array([       nan, 5.54527989]), array([         nan, -72.06454341])]

The WCS object accepts a boolean flag called `with_bounding_box` with default value of `True`. Output values which are outside the `bounding_box` are set to `NaN`. There are cases when this is not desirable and `with_bounding_box=False` should be passes.

Calling the `footprint()` returns the footprint on the sky.

    >>> wcsobj.footprint()
    array([[  5.5264392 , -72.05149213],
    ...    [  5.54580408, -72.06422321],
    ...    [  5.63010439, -72.05426843],
    ...    [  5.610708  , -72.04173847]])

Warning

GWCS and astropy default to different tuple ordering conventions for representing multi-dimensional bounding boxes.

- GWCS uses the `"F"` ordering convention, where the tuples are ordered `((x0min, x0max), (x1min, x1max), ..., (xnmin, xnmax))` (x,y,z ordering).

- While astropy uses the `"C"` ordering convention, where tuples are ordered `((xnmin, xnmax), ..., (x1min, x1max), (x0min, x0max))` (z, y, x ordering).

This means that given the same tuple of tuples, say `((a, b), (c, d))`, setting the bounding box on the transform prior to creating the GWCS will result in a different bounding box than if one sets the same tuple of tuples on the GWCS object itself. Indeed, in this case the former will assume `(c, d)` is the bounding box for `x` while the latter will assume `(a, b)` is the bounding box for `x`.

It is recommended that when working on GWCS objects that one sets the bounding box on the GWCS object itself, rather than on the transform prior to creating the GWCS object.

Note if one wants to set the bounding box on the transform itself rather than the GWCS object then it should be done with `bind_bounding_box` with the `order` argument properly set.

Note

The GWCS will always convert or assume the bounding box to the `"F"` ordering convention when setting the bounding box on the GWCS object itself and will perform this conversion on the first access to the bounding box through the GWCS object. If conversion occurs on first access, GWCS will issue a warning to alert the user that the bounding box has been converted.

## Manipulating Transforms

Some methods allow managing the transforms in a more detailed manner.

Transforms between frames can be retrieved and evaluated separately.

    >>> dist = wcsobj.get_transform('detector', 'undistorted_frame')
    >>> dist(1, 2)
    (41.62325692108675, -12.68101006210054)

Transforms in the pipeline can be replaced by new transforms.

    >>> new_transform = models.Shift(1) & models.Shift(1.5) | distortion
    >>> wcsobj.set_transform('detector', 'undistorted_frame', new_transform)
    >>> wcsobj(1, 2)
    (5.501064280097802, -72.04557376712566)

A transform can be inserted before or after a frame in the pipeline.

    >>> scale = models.Scale(2) & models.Scale(1)
    >>> wcsobj.insert_transform('icrs', scale, after=False)
    >>> wcsobj(1, 2)
    (11.002128560195604, -72.04557376712566)

## Inverse Transformations

Often, it is useful to be able to compute inverse transformation that converts coordinates from the output frame back to the coordinates in the input frame.

Note. the `backward_transform` attribute is equivalent to `forward_transform.inverse`.

In this section, for illustration purpose, we will be using the same 2D imaging WCS from `imaging_wcs_wdist.asdf` created in Adding distortion to the imaging example whose forward transformation converts image coordinates to world coordinates and inverse transformation converts world coordinates back to image coordinates.

    >>> import asdf
    >>> from astropy.utils.data import get_pkg_data_filename
    >>> wcsobj = asdf.open(get_pkg_data_filename('imaging_wcs_wdist.asdf')).tree['wcs']

The most general method available for computing inverse coordinate transformation is the `WCS.invert()` method. This method uses automatic or user-supplied analytical inverses whenever available to convert coordinates from the output frame to the input frame. When analytical inverse is not available as is the case for the `wcsobj` above, a numerical solution will be attempted using `WCS.numerical_inverse()`.

Default parameters used by `WCS.numerical_inverse()` or `WCS.invert()` methods should be acceptable in most situations:

    >>> world = wcsobj(350, 200)
    >>> print(wcsobj.invert(*world))  # convert a single point
    (349.9999994163172, 200.00000017679295)
    >>> world = wcsobj([2, 350, -5000], [2, 200, 6000])
    >>> print(wcsobj.invert(*world))  # convert multiple points at once
    (array([ 1.99999752,  3.49999999e+02, -5.00000000e+03]), array([1.99999972e+00, 2.00000002e+02, 6.00000000e+03])

By default, parameter `quiet` is set to `True` in `WCS.numerical_inverse()` and so it will return results “as is” without warning us about possible loss of accuracy or about divergence of the iterative process.

In order to catch these kind of errors that can occur during numerical inversion, we need to turn off `quiet` mode and be prepared to catch `gwcs.wcs.NoConvergence` exceptions. In the next example, let’s also add a point far away from the image for which numerical inverse fails.

    >>> from gwcs import NoConvergence
    >>> world = wcsobj([-85000, 2, 350, 3333, -5000], [-55000, 2, 200, 1111, 6000],
    ...                with_bounding_box=False)
    >>> try:
    ...     x, y = wcsobj.invert(*world, quiet=False, maxiter=40,
    ...                          detect_divergence=True, with_bounding_box=False)
    ... except NoConvergence as e:
    ...     print(f"Indices of diverging points: {e.divergent}")
    ...     print(f"Indices of poorly converging points: {e.slow_conv}")
    ...     print(f"Best solution:\n{e.best_solution}")
    ...     print(f"Achieved accuracy:\n{e.accuracy}")
    Indices of diverging points: [0]
    Indices of poorly converging points: [4]
    Best solution:
    [[ 1.38600585e+11  6.77595594e+11]
     [ 2.00000000e+00  1.99999972e+00]
     [ 3.49999999e+02  2.00000002e+02]
     [ 3.33300000e+03  1.11100000e+03]
     [-4.99999985e+03  5.99999985e+03]]
    Achieved accuracy:
    [[8.56497375e+02 5.09216089e+03]
     [6.57962988e-06 3.70445289e-07]
     [5.31656943e-06 2.72052603e-10]
     [6.81557583e-06 1.06560533e-06]
     [3.96365344e-04 6.41822468e-05]]