Iterative implementation of wcs.all_world2pix

This overlaps somewhat with #815.

At the moment there is no full-transformation from world to pixel coordinates. From what I understand, this is because some of the WCS math is not analytically invertible, so only all_pix2world can be computed reliably. However, I think the translation between pixel and world coordinates is still (usually) 1-to-1, and an iterative root-finding all_world2pix would be very useful for these cases.

This is of course subject to the vagaries of what gets done when, but this might be more easily implemented in the future "generalized" WCS, which will probably be based on astropy.modeling. It should be quite a bit more straightforward to use modeling to implement this sort of transform, then.

@mdboom

Maybe @mdboom could comment on how feasible this would be? Also, what milestone do we want for the new generalized WCS?

There are plans for such functionality in modeling. For now it's a method of the Model class, Model.invert, not implemented yet. Note that there's also inverse method which is meant to provide an analytical inverse of a model.

As for the generalized WCS, at the last coordination meeting we came up with a basic class. SHould we put this in 0.3?

The only thing that's non-invertible is the lookup table stuff (SIP and the wcslib stuff are already invertible, assuming the SIP headers have the inverse coefficients). It would just be a matter of writing an iterative reverse transformation on that.

But long term, I'd prefer to see the algorithm for doing an inverse lookup transformation live in astropy.models and just have that used from astropy.wcs.

@nden

@nden - yes, 0.3 makes sense - at least we should make sure that the API is ready for the generalized WCS, even if it's not all implemented, so that future releases are backward-compatible.

I don't mind creating a PR with what we talked about but where should it go? astropy.wcs.gwcs.py or a separate subpackage, astropy.gwcs. For now it's just the basic WCS class but later we will need to add modules handling the coordination systems.

@nden

@nden - I think this is a debate that hasn't been resolved, and has been a matter of debate which it is maybe time to re-ignite. For instance, my personal view is that WCS to me does not mean FITS-WCS, but just world to/from pixel conversions in general, so I would say it belongs in astropy.wcs - but I know that's not the view of others :)

Maybe it's time to raise the issue of generalized WCS on the mailing list? (with some background explaining the motivation for this generalized WCS)

@astrofrog

@astrofrog Definitely not FITS-WCS.

I could write something up this week I think if you want.

On May 29, 2013, at 10:41 AM, Thomas Robitaille wrote:

Maybe it's time to raise the issue of generalized WCS on the mailing list? (with some background explaining the motivation for this generalized WCS)

—
Reply to this email directly or view it on GitHub.

How do you envisage the general WCS working in terms of classes - will there be a single WCS class, or will there be different sub-classes for e.g. FITS-WCS, and other types of WCS?

@perrygreenfield

@perrygreenfield - that would be great! I don't think a full API document is needed, but at least just some background to frame the discussion would be good.

I'll try to make a PR in the API repository this week so we can start a discussion.

@nden

@nden - great, thanks!

[It was suggested that we start to raise this issue on the list. I don't know if this is the right level of detail to start with (a teaser?). I thought it outlines things at a fairly high level. We can go into detail in many areas if necessary.]

Perry

It was suggested earlier this week that it was time to raise the issue of wcs models for astropy.

Some may think: "Why? We already have pywcs." Unfortunately the FITS WCS system is inadequate for many problems, and trying to patch it to make it work simply isn't worth the trouble (I can elaborate on that if anyone is interested). At STScI we've been working on a more flexible and generalizable WCS system. This is going to be based on the model/fitting framework that has recently been added to astropy (for 0.3).

This wcs system is centered around coordinate transformations that can be much more flexible than what is provided by FITS. For example, it will be possible to construct more complex transforms combinations of simpler ones, either additively (e.g., a basic distortion transformation with an additive correction model), or in series (e.g., applying a coordinate transformation to the output of the previous transformation). In this way, it will be possible to keep WCS transformations partitioned the way they are determined rather then forcing it into one transformation model. It will have no limits applied on WCS models having to do with 8 character keyword sizes (e.g., like the FITS limitations on polynomial models). It will have many of the advantages properties of AST in this regard.

The new scheme does raise issues of what to do about saving the model in a FITS file. There is no universally good solution to this. There are three basic approaches one can take; not mutually exclusive.

Save the exact model in an opaque way in a FITS file (in a non-standard FITS extension). Our libraries can read and write these, but no standard FITS library (and correspondingly, wcslib) will know about this information in any direct way.
Use an existing FITS WCS model to approximate the astropy model. Information is lost in doing so, and depending on the nature of the transformation, some degree of accuracy. But for many cases (e.g., viewing in DS9) the loss of accuracy is not serious.
Use a pixel map transform (images mapping pixel locations for each of the output WCS coordinates) in FITS to keep full accuracy. However this is cumbersome space wise, and by merging transformations into one, information is still lost, but little accuracy.

A clear byproduct of all this is that some new way of storing this information is needed. We are working on that as well. It will not be python-specific, but a well defined format (i.e., having a specifications document), that we hope will also be suitable for all kinds of data (yes, an alternative to FITS). This is not the place to go into the details. (And no, it won't be HDF5, though it could be made isomorphic with some HDF5 representation, though I doubt I would encourage that). So for option 1) above, we would essentially embed this format within a FITS file for just the WCS component. In most WCS cases, this format would be simple text, with organizational properties visually obvious.

So when we consider the wcs attribute of standard astropy data, this is what we have in mind, not FITS wcs objects (though in general, all FITS WCS models will have a representation within this new scheme).

Just an off-topic side comment on that:

And no, it won't be HDF5, though it could be made isomorphic with some HDF5 representation, though I doubt I would encourage that

I wouldn't necessarily encourage it, but I wouldn't discourage it either. The nice thing about Astropy's I/O framework is we've made it easy (at least for some cases, though I expect that to improve) to open some data in FITS or in New Format and dump it out to HDF5 (if it turns out to be better at handling your 30 GB table, which would be no surprise) and then dump it back to the original format. So in other words, less fussing about what format a file is will be great. And implementing WCS across file formats includes a need for a transferable WCS format (even if it's one that's otherwise opaque to, say, existing FITS readers).

I've been using the following awful little subclass in my own code:

class WCSWithInverse(astropy.wcs.WCS):
    """Patched astropy.wcs.WCS, full inverse with SIP distortion parameters"""
    def _all_world2pix(self, world, origin, **kwargs):
        pixes = np.zeros_like(world)
        for i in range(len(world)):
            soln = scipy.optimize.root(
                (lambda pix: (self.all_pix2world(np.atleast_2d(pix), origin, **kwargs) - world[i]).flatten()),
                self.wcs_world2pix(np.atleast_2d(world[i]), origin, **kwargs).flatten()
            )
            pixes[i] = soln.x.flatten()
        return pixes
    def all_world2pix(self, *args, **kwargs):
        return self._array_converter(self._all_world2pix, 'input', *args, **kwargs)

This is at least enough to (slowly) render pixel-perfect images with APLPy and pywcsgrid2.

I think it would be nice to have a functional (even slow) implementation of this as soon as possible (for 0.3) as the generalized WCS may take longer to sort out.

Is Scipy permitted as a dependency of astropy? If so, you could use scipy.optimize.root like I did above.

@lpsinger

Scipy is valid as an optional dependency -- i.e. functions that aren't critical are allowed to use it as long as importing the package still passes.

This probably fits into that category, and we could use @lpsinger's suggestion for the time being, with the assumption that astropy.models will probably grow some of this stuff in the future, at which time we could switch to using that.

Add iterative wcs.all_world2pix

Navigation Menu

Repository navigation

Iterative implementation of wcs.all_world2pix #1066

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Iterative implementation of wcs.all_world2pix #1066

Description

Activity

eteq commented on May 29, 2013

astrofrog commented on May 29, 2013

nden commented on May 29, 2013

mdboom commented on May 29, 2013

astrofrog commented on May 29, 2013

nden commented on May 29, 2013

astrofrog commented on May 29, 2013

astrofrog commented on May 29, 2013

nden commented on May 29, 2013

perrygreenfield commented on May 29, 2013

astrofrog commented on May 29, 2013

astrofrog commented on May 29, 2013

nden commented on May 29, 2013

astrofrog commented on May 29, 2013

perrygreenfield commented on May 31, 2013

embray commented on May 31, 2013

lpsinger commented on Jul 16, 2013

astrofrog commented on Jul 16, 2013

lpsinger commented on Jul 16, 2013

mdboom commented on Jul 16, 2013

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions