Handle units on input

[jdavies-st]

The gwcs.WCS.__call__ function has a with_units parameter that can be used to return an astropy.coordinates.SkyCoord object with units.

But if one tries to round trip this by taking this astropy.coordinates.SkyCoord object and feed it back into the backward_transform, then it fails, as it's expecting N inputs, where N is the number of dimensions inside the SkyCoord object.

Following the example from the docs, we construct a simplified gwcs.WCS object:

import numpy as np
from astropy.modeling.models import (Shift, Scale, Rotation2D,
      Pix2Sky_TAN, RotateNative2Celestial)
det2sky = (Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) |
           Scale(.01) & Scale(.04) | Pix2Sky_TAN() |
           RotateNative2Celestial(5.6, -72.05, 180)).rename("det2sky")

from astropy import units as u
from astropy import coordinates as coord
from gwcs import coordinate_frames as cf
detector_frame = cf.Frame2D(name="detector", unit=(u.pix, u.pix))
sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs')

from gwcs import wcs
pipeline = [(detector_frame, det2sky),
            (sky_frame, None)
            ]
wcsobj = wcs.WCS(pipeline)

And then test it out:

In [12]: wcsobj(1, 1)
Out[12]: (5.599990008144789, -71.73800308407476)

In [13]: wcsobj(1, 1, with_units=True)
Out[13]: 
<SkyCoord (ICRS): (ra, dec) in deg
    (5.59999001, -71.73800308)>

In [14]: w = wcsobj(1, 1, with_units=True)

In [15]: wcsobj.backward_transform(w, with_units=True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-23-fb3d7869d5a0> in <module>()
----> 1 wcsobj.backward_transform(w, with_units=True)

TypeError: __call__() got an unexpected keyword argument 'with_units'

In [16]: wcsobj.backward_transform(w)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-24-99c2ede38a9c> in <module>()
----> 1 wcsobj.backward_transform(w)

TypeError: __call__() missing 1 required positional argument: 'delta_C'

Now, one can decompose the SkyCoord object by hand, but even here we have to be careful.

In [34]: w.ra
Out[34]: <Longitude 5.28344155 deg>

In [35]: w.dec
Out[35]: <Latitude -72.53775313 deg>

In [36]: wcsobj.backward_transform(w.ra, w.dec)
---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
<ipython-input-36-3d1c306b417f> in <module>()
----> 1 wcsobj.backward_transform(w.ra, w.dec)

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in __call__(self, alpha_C, delta_C, model_set_axis, with_bounding_box, fill_value, equivalencies)
    380                                      ('with_bounding_box', False),
    381                                      ('fill_value', np.nan),
--> 382                                      ('equivalencies', None)])
    383 
    384             # The following makes it look like __call__ was defined in the class

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in __call__(self, *inputs, **kwargs)
    361             def __call__(self, *inputs, **kwargs):
    362                 """Evaluate this model on the supplied inputs."""
--> 363                 return super(cls, self).__call__(*inputs, **kwargs)
    364 
    365             # When called, models can take two optional keyword arguments:

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in __call__(self, *inputs, **kwargs)
    813                     outputs = [np.asarray(r) for r in result]
    814         else:
--> 815             outputs = self.evaluate(*chain(inputs, parameters))
    816         if self.n_outputs == 1:
    817             outputs = (outputs,)

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in evaluate(self, *args)
   2919     @sharedmethod
   2920     def evaluate(self, *args):
-> 2921         return self.__class__.evaluate(*args)
   2922 
   2923     # TODO: The way this works is highly inefficient--the inverse is created by

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in evaluate(cls, *args)
   2360         inputs = args[:cls.n_inputs]
   2361         params = iter(args[cls.n_inputs:])
-> 2362         result = cls._evaluate(inputs, params)
   2363         if cls.n_outputs == 1:
   2364             return result[0]

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in <lambda>(inputs, params)
   2192     #
   2193     # and similarly for g
-> 2194     return (lambda inputs, params: g[0](f[0](inputs, params), params),
   2195             f[1], g[2])
   2196 

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in <lambda>(inputs, params)
   2192     #
   2193     # and similarly for g
-> 2194     return (lambda inputs, params: g[0](f[0](inputs, params), params),
   2195             f[1], g[2])
   2196 

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in <lambda>(inputs, params)
   2192     #
   2193     # and similarly for g
-> 2194     return (lambda inputs, params: g[0](f[0](inputs, params), params),
   2195             f[1], g[2])
   2196 

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in <lambda>(inputs, params)
   2192     #
   2193     # and similarly for g
-> 2194     return (lambda inputs, params: g[0](f[0](inputs, params), params),
   2195             f[1], g[2])
   2196 

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in <lambda>(inputs, params)
   2203     #
   2204     # and similarly for g
-> 2205     return (lambda inputs, params: (f[0](inputs[:f[1]], params) +
   2206                                     g[0](inputs[f[1]:], params)),
   2207             f[1] + g[1], f[2] + g[2])

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in f(inputs, params)
   2800             def f(inputs, params):
   2801                 param_values = tuple(islice(params, n_params))
-> 2802                 return evaluate_wrapper(model, inputs, param_values)
   2803         else:
   2804             # Where previously model was a class, now make an instance

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/core.py in evaluate_wrapper(model, inputs, param_values)
   2792         def evaluate_wrapper(model, inputs, param_values):
   2793             inputs = model._validate_input_units(inputs)
-> 2794             outputs = model.evaluate(*inputs, *param_values)
   2795             if n_outputs == 1:
   2796                 outputs = (outputs,)

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/modeling/functional_models.py in evaluate(x, offset)
    479     def evaluate(x, offset):
    480         """One dimensional Shift model function"""
--> 481         return x + offset
    482 
    483     @staticmethod

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/units/quantity.py in __array_ufunc__(self, function, method, *inputs, **kwargs)
    618         # consistent units between two inputs (e.g., in np.add) --
    619         # and the unit of the result (or tuple of units for nout > 1).
--> 620         converters, unit = converters_and_unit(function, method, *inputs)
    621 
    622         out = kwargs.get('out', None)

~/anaconda3/envs/jwst_dev/lib/python3.6/site-packages/astropy/units/quantity_helper.py in converters_and_unit(function, method, *args)
    565                                      "argument is not a quantity (unless the "
    566                                      "latter is all zero/infinity/nan)"
--> 567                                      .format(function.__name__))
    568             except TypeError:
    569                 # _can_have_arbitrary_unit failed: arg could not be compared

UnitsError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)

One can jettison the units altogether, and it eventually works.

In [37]: wcsobj.backward_transform(w.ra.value, w.dec.value)
Out[37]: (0.999999999999897, 0.9999999999997442)

But the with_units parameter only seems defined on the WCS.__call__ method, not WCS.forward_transform or WCS.backward_transform.

In [38]: wcsobj.backward_transform(w.ra.value, w.dec.value, with_units=True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-38-72f326c8ef7a> in <module>()
----> 1 wcsobj.backward_transform(w.ra.value, w.dec.value, with_units=True)

TypeError: __call__() got an unexpected keyword argument 'with_units'

This is mostly to document this behavior to take into account any changes to the calling API that might be made in the future.

[nden]

@jdavies-st wcsobj.forward_transform and wcsobj.backward_transform use the astropy.modeling APi for models while with_units is part of the GWCS API.
wcsobj.invert(skyCoordObj) works.
So yes, this needs to be documented.

[jdavies-st]

Ah, I see!

[pllim]

I encountered this same error at ejeschke/ginga#719 but it is not obvious how I can use invert for my usage. Here is the method that calls your method:

    def radectopix(self, ra_deg, dec_deg, coords='data', naxispath=None):
        # NOTE: origin is always 0, coords unused.

        args = [ra_deg, dec_deg]
        if naxispath:
            args += [0] * len(naxispath)
        skycrd = np.array([args], np.float_)

        try:
            xy = np.squeeze(self.wcs.world_to_pixel_values(skycrd))[:2]
        except Exception as e:
            self.logger.error(
                "Error calculating radectopix: {}".format(str(e)))
            raise common.WCSError(e)

        return xy

p.s. skycrd here is not SkyCoord instance, but rather RA and DEC already in degrees.

[stscijgbot]

This ticket is now being tracked at AL-34

[Cadair]

I think this can be closed?!

[jdavies-st]

Yes, I think so, as APE14 interface solves this. That said, there's still not very good documentation within the gwcs docs that show how to actually take coordinates and transform them from one frame to another via APE14 interface. So that would be nice to document. Also, the bounding box documentation is very out-of-date.

[nden]

Right, the documentation needs a lot of work. Hopefully we can do some soon.
The APE 14 interface cannot be used to transform between intermediate frames. And this should be documented as well.