DMTN-010: WCS and Distortion Requirements and Existing Options

  • John Parejko, Russell Owen

Latest Revision: 2016-05-25

Introduction

LSST-DM currently uses WCSLIB to persist/un-persist and manipulate World Coordinate System (WCS) transformations. WCSLIB is based on the work by Greisen & Calabretta 2002 [3], which were incorporated into the FITS WCS standard. While the FITS WCS standard does not support non-linear distortion corrections, WCSLIB does support some community extensions that work within the constraints of the standard. These extensions provide single-polynomial distortion models, which severely limits the ability to describe complex CCD, focal plane and on-sky distortions. Most previous projects have employed this or a related FITS WCS-based library–usually with some home-built custom functionality–to manage their astrometric results.

There is no standardized method in the astronomical community to improve upon or extend the FITS WCS standard. The original paper describing the FITS distortion standard, Calabretta et al. 2004 (in prep), was not adopted by the FITS community and the paper remains unfinished and unpublished. The Simple Imaging Polynomial convention [5] allows a distortion model represented by a polynomial of up to 9th order. The DECam community pipeline uses the TPV convention which allows a 7th-order polynomial distortion correction. SDSS produced their own asTran model to map the (row,column) coordinates from each field into (mu,nu) great circle spherical coordinates via a 3rd order polynomial.

Two much more flexible, powerful, and extensible systems are Starlink AST and STScI’s GWCS. The Starlink AST package [1], developed by David Berry (East Asian Observatory) in C, with a python interface (PyAST) written by Tim Jenness, provides models that can be combined in a variety of ways. These more advanced models are not widely used outside the Starlink software suite, although ds9 links with AST and so files with those features are viewable in ds9. Nadia Dencheva and Perry Greenfield (STScI) are developing a python-based Generalized World Coordinate System package (GWCS) for JWST, building on top of the generic mathematical modeling system provided by astropy.modeling.

This document discusses the expected requirements for the LSST distortion model and coordinate transform system, the options we have to select from, and provides a recommendation for how we should achieve our requirements.

Requirements

The initial call for discussion of future WCS/distortion model requirements was on this Community posting. Independently, Jim Bosch’s post about fitted models presented some items which may guide our selection of transformation modeling systems.

LSST’s most critical requirements are:

  • Specifying complex parametric and non-parametric models.
  • Arbitrary combinations/compositions of those models.
  • Ability to provide approximate WCSLIB-style FITS WCS standard output, for legacy software use. The particular choice of non-standard convention to produce (TPV, SIP, etc.) could be a user-supplied parameter.
  • Fast pixel-to-pixel performance for image warping and small (~400px) postage stamps in multifit. These transforms are generally not vectorizable, and thus may be difficult to optimize in e.g. numpy. For small arrays, it may be enough to produce an affine transform over that small region.
  • Shared serialization format with GWCS, to allow LSST files to be used in non-LSST code and vice-versa (at a recent joint meeting, LSST and AstroPy representatives discussed leveraging the work the IVOA has done on STC2 as a possible route toward this). This requires that all transforms used by LSST are available within GWCS: LSST could contribute the required code to GWCS for any of our externally-visible transforms that they do not have.

Coordinate and Transformation Requirements

  • The ability to model all pixel distortion effects that are frozen through one exposure (not Brighter-Fatter). These are likely our “exotic” ones–e.g. tree rings, edge roll-off–that are not yet included in any currently existing distortion modeling system.
  • Mappings between on-camera coordinate systems (i.e. the future versions of afw.cameraGeom and afw.geom.XYTransform) must be entirely interoperable with image->sky (i.e. the future afw.image.Wcs) transformations.
    • A method for easily creating simple WCS from e.g. existing files or a handful of on-sky values.
    • A method to easily produce an initial guess WCS by combining the post-ISR CCD geometry and telescope pointing, to feed into our astrometric solver.
  • Compositions of transforms should be simplifiable, for optimal performance and simpler serialization.
    • Compositions should only simplify when it can be done exactly (to machine precision), or when explicitly requested with a bound on accuracy.
  • Transforms should know their domain to ensure that composed transforms are valid between different domains.
  • Transforms should know whether they involve spherical-spherical, spherical-Cartesian, Cartesian-spherical, or Cartesian-Cartesian coordinates, to allow interoperability with geometry libraries that distinguish these coordinate systems.
  • Transforms should be invertible, or if not invertible should allow pre-defined one-way transforms to be identified as each-other’s inverses.
  • Ability to efficiently obtain the exact transform at one point for warping. This needs to be fast and parallel.
  • Ability to efficiently obtain a local linear approximation around a point (must be parallel).
  • A method for obtaining a local TAN WCS approximation.
  • Transforms must be persistable as components of arbitrary objects (e.g. Exposure, Psf).
    • Groups of composed transforms should be persisted efficiently, e.g. for all CCDs in a visit, to reduce data size. It is unclear if this is actually necessary: the only obvious “heavyweight” transform is pixel distortion (e.g. tree rings), which is per-CCD anyway, while other transforms will likely be a very small part of our data volume.
    • An efficient (compressed?) persistence format for pixel-grid-based transformations will be necessary if we end up using full pixel grids to represent pixel distortions.
  • The ability to compute or provide derivatives with respect to pixel coordinates (to compute local affine transformations).
    • If we use the same system for Fitters as we do for Consumers (see Models for Consumers vs. Fitters), some transforms will need the ability to compute or provide (at least) first derivatives with respect to the parameters.
  • The ability to add more transforms in the future as we discover a need for them: polynomial/Chebyshev transforms are not enough.
  • We likely do not need to include wavelength-dependent effects (e.g. Differential Chromatic Refraction) in the WCS, if we define our PSFs with offset centroids.

Models for Consumers vs. Fitters

We may want to have separate model representations for astrometric fitting (Fitters) and for using the result of the fit (Consumers). In the current LSST stack, we have an afw.geom.XYTransform as the interface for the consumer (fitted) model, while the Gtransfo object introduced in jointcal is an interface for a particular fitter.

  • It is useful to separate the fitter from the consumer, as they may have different “best” internal representations.
  • It can be conceptually and programmatically helpful to have the same underlying system (or API) for both, to allow easy transfer between them.
  • Fitters must be mutable. Consumers need not be and may be better as immutable to allow the complex object to be safely shared across threads.
  • Consumers must be persistable. Fitters may not need to be.

As a related point, it could be useful to have the same model description system available for other purposes (e.g. representing galaxy shapes, photometric calibration).

Options

There are essentially 5 options available to us, with varying trade-offs between work required, flexibility, likely performance, ease of calling from C++, and standardization in the broader community. These options are not necessarily mutually exclusive; in particular we could begin with 3. Adopt Starlink AST while developing a new system per 4. Adopt AstroPy GWCS or 5. Work with David Berry to develop modern C++ version of AST. Note that 3. Adopt Starlink AST represents a continuum in which our interface to AST could change as our needs and API design evolves.

1. Develop our own

Following the grand tradition of past astronomy surveys, we could develop our own WCS/distortion software (likely in C++, with a python interface), independent of any existing implementation. This seems like an obviously bad choice, given the work that has already gone into AST and GWCS.

Advantages

  • We have full control over the implementation of and interface to the models.

Disadvantages

  • Significant time investment.
  • FITS WCS standard and extensions (e.g. SIP, TPV) not immediately available to us.
  • Yet Another WCS “Standard”.
  • WCS and distortion models are complex objects: usable interface is challenging to develop.
  • Lessons learned by previous groups would be hard to capture.

2. Build on top of WCSLIB

Instead of starting entirely from scratch, we could continue to build on top of WCSLIB. This has the advantage of not having to re-implement the FITS WCS standard, but may be limiting in what we would be able to build on top of it, in addition to requiring nearly as much effort as option 1, above.

Advantages

  • We have nearly full control over the implementation of and interface to the models.
  • FITS WCS standard immediately available to us.

Disadvantages

  • Significant time investment.
  • Enhancements on top of FITS are Yet Another WCS “Standard.”
  • FITS WCS has inherent limitations in namespace, extensibility, flexibility.
  • WCS and distortion models are complex objects: usable interface is challenging to develop.
  • Lessons learned by previous groups would be hard to capture.

4. Adopt AstroPy GWCS

GWCS is a Generalized World Coordinate System library currently being developed by STScI for use by JWST. It is written in pure python, and built on top of the astropy.modeling framework. Complex models can be built from more simple models via standard mathematical operations, and can be composed and chained in serial and parallel. It is under active development, so LSST could have a hand in shaping its future path.

Advantages

  • FITS WCS standard immediately available to us (not clear if all portions of Greisen & Calabretta 2002 [3], Calabretta & Greisen 2002 [2], or the SIP/TPV conventions are currently implemented).
  • More complicated distortion models immediately available to us.
  • Pure python, allowing easy extension.
  • Clean API for adding additional models.
  • Significant and understandable documentation already exists.
  • Community adoption likely very high.
  • Would share development effort with STScI.
  • Serialization format would be automatically shared with GWCS.

Disadvantages

  • Significant time investment: current code manipulates WCS in C++.
  • Not directly callable from C++: calls to python from C++ may incur significant overhead.
  • Model description framework is pure python: unclear if performance requirements can be met, particularly for warping.
  • Ongoing development work: not all features we may need are available.
  • No effort yet on performance optimizations.
  • Less likely that serialization format would be available outside the python community.

5. Work with David Berry to develop modern C++ version of AST

Section 6 of the AST paper discusses “lessons learned”, including a statement that they would have developed it in C++, if they were starting development now. David Berry is interested in re-implementing AST in a modern language as a legacy to the community. LSST could contract him out and guide the development of a new implementation of AST that we could use from C++, while solving some of the current limitations in AST (e.g. adding quad-double precision for time, better unit support, clearer API).

As part of this process, the astropy.modeling API should be used as a reference for how to create and combine models. Their method of using mathematical operations to combine transforms makes the creation of complicated models from simpler components highly intuitive, and presents a good design to build a C++ transformation system from.

To be of greatest benefit to the community, the new AST should be independent of the LSST stack. Some necessary features of the stack, e.g. lsst.afw.coords, could be pushed up into AST, to make them more widely available to the community. This could also simplify the “astropy integration question” [4], by pushing much of the low-level astropy linkages into the new AST and out of afw.

Advantages

  • Lessons learned from AST development can be directly applied.
  • AST has significant test suite and would be a reference implementation to guide development.
  • LSST has influence on new API.
  • LSST can take long-term ownership of new system.
  • David Berry willing to be contracted out for development.
  • major portions of AST code likely can be copied to new interface with minimal changes (e.g. FITS WCS support).

Disadvantages

  • Significant time investment (shared with David Berry).
  • Details of contract with East Asian Observatory need to be developed.
  • Requires LSST C++ expertise to design new API, and produce idiomatic C++.
  • Unclear how much LSST guidance would be required to make a long-term supportable, well documented API.

Recommendations

There are three clearly viable choices: some variation on 3. Adopt Starlink AST; 4. Adopt AstroPy GWCS; or 5. Work with David Berry to develop modern C++ version of AST. The choice between these is a balance between having a workable solution in a relatively short time (3.) vs. having a modern API and functionality whose details we have more direct control over (4. and 5.). We estimate 2 developer months would be required to implement a usable abstraction layer between AST and the LSST stack, whereas implementing the LSST requirements (not all current features of AST, e.g. no spectra/time mappings) in a new C++-based AST would likely require at least 6 months of David Berry’s time, with a comparable amount of LSST developer time for design and guidance. Similarly, we expect that adapting our C++ warping code into python (and possibly making our Exposure object pure-python) and implementing our required transforms in GWCS would be at least 2 months of developer time, while probably a year would be required to attempt (see below) to bring GWCS up to our required performance levels and make all LSST WCS and transforms operate in pure python.

Table 1 Estimated work required
  minimal optimal
effort
(months)
 result
effort
(months)
 result
AST 2
minimal AST wrapper replacing
afw.image.wcs and afw.geom.XYTransform
 12
C++ AST meeting LSST requirement
(i.e. no spectra/time mappings)
GWCS 2
LSST warping code in python using GWCS;
transforms implemented; not performant
 6-12
performant GWCS; LSST transforms in python;
approximation output as e.g. FITS SIP

Although adopting GWCS would be ideal from the perspective of getting involvement from the broader astronomical python community, there are two main reasons we are not recommending that option at this time:

  1. It is unclear whether GWCS would be able to achieve our required performance targets when computing transformations on small pixel regions. Our testing (GWCS/PyAST Performance comparison) found a very significant overhead (10-200 times slower) when using GWCS over small (~100-1000) pixel regions (see appendix pyast/gwcs). Some of this overhead could be removed if LSST put effort into optimizing GWCS, but it is unclear whether optimizations to a python library would be sufficient for our needs. It is even less clear whether we could use a python-based WCS and transform library from within C++ without sustaining a significant performance penalty.
  2. Our current warping code–afw.math.warpExposure–is written purely in C++ and would incur a significant effort to rewrite in python. Warping involves calculations on small patches in a manner that is not easily vectorized. Because of the concerns about performance on small patches described above, it is unclear if the new product would be performant enough to justify the effort.

The LSST warping code is one of our most WCS-related performance intensive calculations. We achieved approximately a 2x performance improvement by computing the WCS on a grid and linearly interpolating across the grid. This suggests that the actual WCS calculation is a more time-intensive part of the warping calculation than the convolution step, implying that any WCS implementation we choose must be equally performant. This comparison becomes worse when corrections for tree rings and other high-order distortions come into play: they vary on the few-pixel level, and thus linear interpolation across dozens of pixels will likely not properly account for them.

Our AST TAN-SIP test found AST’s pixel->sky->pixel performance for a TAN-SIP FITS file is two times faster than the current LSST stack’s implementation, so long as the input points are provided in a single call (e.g. vectorized) rather than as individual calls per point. Thus, using AST for warping would require changing the current warping code to do some vectorization, but the code could remain at the C++ layer.

Given the requirements, options, and caveats listed above, our recommendation is to immediately begin implementing a rewrite of afw.image.Wcs, afw.cameraGeom and afw.geom.XYTransform on top of AST (option 3.), while pursuing a new C++ rewrite of AST (option 5.) that takes into account the lessons learned from the design and API of astropy.modeling and GWCS. We will decide how much to abstract AST as we design the new afw API, and that API can help guide the new C++-based AST rewrite. In addition, so long as we insist on sharing a serialization format with GWCS and work together to ensure we can round-trip data between the projects, we would retain the option of using GWCS in the future.

References

[1]D. S. Berry, R. F. Warren-Smith, and T. Jenness. AST: A library for modelling and manipulating coordinate systems. Astronomy and Computing, 15:33–49, April 2016. arXiv:1602.06681, doi:10.1016/j.ascom.2016.02.003.
[2]M. R. Calabretta and E. W. Greisen. Representations of celestial coordinates in FITS. \aap , 395:1077–1122, December 2002. arXiv:astro-ph/0207413, doi:10.1051/0004-6361:20021327.
[3](1, 2) E. W. Greisen and M. R. Calabretta. Representations of world coordinates in FITS. A&A, 395:1061–1075, December 2002. arXiv:astro-ph/0207407, doi:10.1051/0004-6361:20021326.
[4]T. et al. Jenness. Investigating interoperability of the LSST Data Management software stack with Astropy. SPIE proceedings, in-prep.
[5]D. L. Shupe, M. Moshir, J. Li, D. Makovoz, R. Narron, and R. N. Hook. The SIP Convention for Representing Distortion in FITS Image Headers. In P. Shopbell, M. Britton, and R. Ebert, editors, Astronomical Data Analysis Software and Systems XIV, volume 347 of Astronomical Society of the Pacific Conference Series, 491. December 2005.

Appendix: PyAst/GWCS Performance Comparison

Our basic performance comparison results between PyAST and GWCS, and the code to reproduce them, are given below.

The code takes a file with a basic FITS WCS, and adds a 2nd order 2D polynomial to convert actual pixels to mean pixels, with a pure TAN WCS for mean pixels to sky. We consider this a minimal complexity test of the performance of the two systems, while also demonstrating their interfaces. The results shown below were run on a mid-2012 Macbook Pro (2.6GHz i7).

Table 2 GWCS/PyAST Performance comparison
# points GWCS PyAST GWCS PyAST ratio
  time (µs) time (µs) time/pt time/pt  
102 15200 62.3 152 0.623 244
1002 18900 2610 1.89 0.261 7.3
10002 346000 269000 0.346 0.269 1.3

Python comparison code is shown below. This requires having recent versions of both PyAST and GWCS installed. Both are available to install via pip.

#!/usr/bin/env ipython
"""Comparison between GWCS and AST.

You must run this with ipython and you must create simple.fits first. Thus:
$ python makeExposure.py
$ ipython compare_gwcs_ast.py
"""
from __future__ import absolute_import, division, print_function
import os.path

import numpy as np
import astropy.io.fits as pyfits

# PyAST imports
from starlink import Ast, Atl

# astropy imports
import gwcs
from gwcs import coordinate_frames as cf
from astropy.modeling import models
from astropy import units as u
from astropy import coordinates as coord


def get_pix_to_meanpix_coeffs():
    """Return Polynomial2D coefficients for the pix to meanpix map

    Return a tuple of three items:
    - polynomial order
    - x coeff
    - y coeff
    the coefficients are for powers: x0y0, x0y1, x1y0, x0y2, x1y1, x2y2, ...
    """
    return [
        3,
        [
            0,
            1.0001, 1e-4,
            1e-6, 1e-6, 1e-6,
        ],
        [
            0,
            1e-4, 1.0002,
            1e-6, 1e-6, 1e-6,
        ],
    ]


def ast_poly_coeff_iter(ind, coeffs):
    """An iterator for ast poly coeffs for the given index: x=1, y=2

    @param[in] ind: output index: 0 for x, 1 for y
    @param[in] coeffs: list of coefficients for that index, as returned by get_pix_to_meanpix_coeffs

    @return a tuple containing the following (all cast to float):
    - coeff
    - ind (the input argument)
    - x power
    - y power
    """
    xpower = 0
    ypower = 0
    for coeff in coeffs:
        if coeff != 0:
            yield (coeff, float(ind+1), float(xpower), float(ypower))
        if ypower == 0:
            ypower = xpower + 1
            xpower = 0
        else:
            ypower -= 1
            xpower += 1


def make_ast_actpix_to_meanpix(acc=1e-3, maxacc=1e-4, maxorder=3, minxy=(0, 0), maxxy=(512, 512)):
    """Make an Ast.PolyMap model for actual pixels to pixels

    The forward transform is specified by coefficients from get_pix_to_meanpix_coeffs, whereas
    the reverse transform is fit. The arguments control the reverse transform fit.

    @param[in] acc  desired accuracy of reverse transform fit
    @param[in] maxacc  maximum accuracy of reverse transform fit
    @param[in] maxorder  maximum order of reverse transform fit
    @param[in] minxy  minimum x,y over which reverse transform will be fit
    @param[in] maxxy  maximum x,y over which reverse transform will be fit
    """
    # create forward transformation
    order, xcoeffs, ycoeffs = get_pix_to_meanpix_coeffs()
    xastcoeffs = tuple(ast_poly_coeff_iter(0, xcoeffs))
    yastcoeffs = tuple(ast_poly_coeff_iter(1, ycoeffs))
    polymap = Ast.PolyMap(xastcoeffs + yastcoeffs)
    #
    # fit reverse transformation; 2nd argument is 0 to fit reverse
    polymap = polymap.polytran(0, acc, maxacc, maxorder, minxy, maxxy)
    if polymap is None:
        raise RuntimeError("Could not compute inverse Ast.PolyMap")
    return polymap


def make_astropy_poly2d(order, coeffs):
    """Make an astropy models.Polynomial2D

    @param[in] order  order of pollynomial
    @param[in] coeffs  list of coefficients for one output axis;
        one of the two arrays returned by get_pix_to_meanpix_coeffs
    """
    coeffdict = dict()
    xpower = 0
    ypower = 0
    for coeff in coeffs:
        if coeff != 0:
            coeffname = "c%d_%d" % (xpower, ypower)
            coeffdict[coeffname] = coeff
        if ypower == 0:
            ypower = xpower + 1
            xpower = 0
        else:
            ypower -= 1
            xpower += 1
    return models.Polynomial2D(order, **coeffdict)


def make_astropy_actpix_to_meanpix():
    """Make an astropy.model model for actual pixels to mean pixels

    Uses the coefficients returned by get_pix_to_meanpix_coeffs
    """
    order, xcoeffs, ycoeffs = get_pix_to_meanpix_coeffs()
    xmodel, ymodel = [make_astropy_poly2d(order, coeffs) for coeffs in (xcoeffs, ycoeffs)]
    return models.Mapping((0, 1, 0, 1)) | xmodel & ymodel


def build_gwcs(hdu):
    """Make a gWCS objects from a pyfits HDU of an image, which should contain a simple TAN WCS
    """
    # to get all the parts of the FITS WCS transform, so we can insert
    # our distortion model in between.
    hdu.header["WCSAXES"] = 2
    wcs_info = gwcs.utils.read_wcs_from_header(hdu.header)
    fits_linear = gwcs.utils.fitswcs_linear(wcs_info)
    projcode = gwcs.utils.get_projcode(wcs_info['CTYPE'])
    fits_tan = gwcs.utils.create_projection_transform(projcode).rename(projcode)
    # now get the rotation
    phip, lonp = [wcs_info['CRVAL'][i] for i in gwcs.utils.get_axes(wcs_info)[0]]
    thetap = 180
    n2c = models.RotateNative2Celestial(phip, lonp, thetap, name="crval")
    fits_tan_rot = fits_tan | n2c

    actpix_to_pix = make_astropy_actpix_to_meanpix()

    # Create some coordinate frames to map between
    actpix = cf.Frame2D(name='actual_pixels', axes_order=(0, 1), unit=(u.pix, u.pix))
    pix = cf.Frame2D(name='pixels', axes_order=(0, 1), unit=(u.pix, u.pix))
    focal = cf.Frame2D(name='focal', axes_order=(0, 1), unit=(u.pix, u.pix))
    sky = cf.CelestialFrame(name='icrs', reference_frame=coord.ICRS())

    # tuples of frame:mapping
    # The last frame has None for the transform.
    pipeline = [(actpix, actpix_to_pix),
                (pix, fits_linear),
                (focal, fits_tan_rot),
                (sky, None)
                ]
    return gwcs.wcs.WCS(pipeline)


def build_ast(hdu):
    """Make an AST-based WCS from a pyfits HDU of an image, which should contain a simple TAN WCS
    """
    fitschan = Ast.FitsChan(Atl.PyFITSAdapter(hdu))
    # Set Iwc to True to leave a place to insert optical distortion;
    # it seems harmless and doesn't affect timing, but does make the frameset more complicated
    fitschan.Iwc = True
    frameset = fitschan.read()
    pix_frame_ind = frameset.Base
    sky_frame_ind = frameset.Current
    # print "TAN WCS frameset:\n", frameset

    # map of actual pixel to nominal pixel represented by a 2-d polynomial
    actpix_to_pix_map = make_ast_actpix_to_meanpix(maxxy=hdu.data.shape)
    actpix_frame = Ast.Frame(2, "Domain=GRID")

    # add this to the frameset
    # temporarily set inverse=true because it is inserted as a mapping from mean pixels to pixels
    # inserting sets current to the new frameset, so record that index
    actpix_to_pix_map.Invert = 1
    try:
        frameset.addframe(pix_frame_ind, actpix_to_pix_map, actpix_frame)
    finally:
        actpix_to_pix_map.Invert = 0
    actpix_frame_ind = frameset.Current

    # set meanpix as the base and sky as the current
    frameset.Base = actpix_frame_ind
    frameset.Current = sky_frame_ind

    return frameset


if __name__ == "__main__":
    currdir = os.path.dirname(__file__)
    infile = os.path.join(currdir, "simple.fits.gz")
    hdu = pyfits.open(infile)[1] # reading in a MaskedImage; image is in HDU 1
    imarr = hdu.data

    gwcs_wcs = build_gwcs(hdu)

    ast_wcs = build_ast(hdu)

    print("*** Verifying that GWCS and AST give the same results")

    nx = imarr.shape[0]
    ny = imarr.shape[1]
    xx = np.arange(0, nx)
    yy = np.arange(0, ny)
    xv, yv = np.meshgrid(xx, yy)
    result_gwcs = gwcs_wcs(xv, yv)
    result_ast = (180/np.pi)*ast_wcs.tran((xv.flatten(), yv.flatten()))
    result_ast = (result_ast[0].reshape(ny, nx), result_ast[1].reshape(ny, nx))
    result_trangrid = (180/np.pi) * ast_wcs.trangrid([0, 0], [511, 511])
    result_trangrid = (result_trangrid[0].reshape(ny, nx), result_trangrid[1].reshape(ny, nx))

    assert np.allclose(result_ast[0], result_gwcs[0])
    assert np.allclose(result_ast[1], result_gwcs[1])
    assert np.allclose(result_ast[1], result_trangrid[1])

    # What does the tolerence parameter actually mean? Once we know this may be a useful test:
    # print("How divergent is trangrid if we specify a tolerance?")
    # result_trangrid = (180/np.pi) * ast_wcs.trangrid([0, 0], [511, 511], 1000)
    # result_trangrid = (result_trangrid[0].reshape(ny, nx), result_trangrid[1].reshape(ny, nx))
    # print(np.allclose(result_ast[1], result_trangrid[1]))

    print("*** Passed!")


    # Timing comparisons

    magic = get_ipython().magic  # flake8: noqa: ipython-specific code

    nx = imarr.shape[0]
    ny = imarr.shape[1]
    for numxpts in (10, 100, 1000):
        numypts = numxpts
        xx = np.linspace(0, nx, numxpts)
        yy = np.linspace(0, ny, numypts)
        xv, yv = np.meshgrid(xx, yy)
        numpts = numxpts*numypts

        print("\n*** Timing pixel->sky for %d points" % (numpts,))

        print("gwcs:")
        magic('%timeit gwcs_wcs(xv, yv)')
        print("ast trans:")
        magic('%timeit ast_wcs.tran((xv.flatten(), yv.flatten()))')
        print("ast trangrid with tol=1e-3:")
        magic('%timeit ast_wcs.trangrid([0, 0], [511, 511], 1e-3)')