Forward (and possibly inverse) transforms do not work with arbitrarily shaped inputs

[mcara]

I need to convert pixel coordinates of 3D images (cubes) to world coordinates. The forward transformation crashes when input coordinates are 3D arrays (NOTE: cube is a cube model):

>>> ind = np.indices((372,32,34))[::-1]
>>> cube.meta.wcs(*ind) 
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/.../lib/python2.7/site-packages/gwcs/wcs.py", line 228, in __call__
    result = self.forward_transform(*args)
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 386, in __call__
    __call__, args, [('model_set_axis', None)])
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 382, in __call__
    return super(cls, self).__call__(*inputs, **kwargs)
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 687, in __call__
    outputs = self.evaluate(*chain(inputs, parameters))
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 2459, in evaluate
    return self.__class__.evaluate(*args)
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 1929, in evaluate
    result = cls._evaluate(inputs, params)
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 1755, in <lambda>
    return (lambda inputs, params: (f[0](inputs[:f[1]], params) +
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 1744, in <lambda>
    return (lambda inputs, params: g[0](f[0](inputs, params), params),
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 1744, in <lambda>
    return (lambda inputs, params: g[0](f[0](inputs, params), params),
  File "/.../lib/python2.7/site-packages/astropy/modeling/core.py", line 2374, in <lambda>
    evaluate(*chain(inputs, islice(params, n_params)))
  File "/.../lib/python2.7/site-packages/astropy/modeling/rotations.py", line 206, in evaluate
    phi, theta, psi)
  File "/.../lib/python2.7/site-packages/astropy/modeling/rotations.py", line 170, in _evaluate
    lon_pole, self.axes_order)
  File "/.../lib/python2.7/site-packages/astropy/modeling/rotations.py", line 88, in evaluate
    result = np.dot(matrix, inp)
ValueError: shapes (3,3) and (3,372,32,34) not aligned: 3 (dim 1) != 32 (dim 2)

Cube file will be e-mailed to @nden

[jdavies-st]

Mihai, have you tried using the utility gwcs.wcstools.grid_from_bounding_box instead of np.indices? This orders the slices correctly. Just pass it the bounding box from the wcs object in the cube.

bb = cube.meta.wcs.bounding_box
ind = grid_from_bounding_box(bb, step=(1, 1, 1))

Not sure if this is related to your problem.

[mcara]

@jdavies-st Have you actually tried your solution on a cube and it worked? The reason I am asking this is because I tried it on two cubes (one from @jemorrison and one produced by me by running cube builder) and in both cases I get something like this:

>>> m = IFUCubeModel('model_0_test_single_ch1_s3d.fits')
>>> bb = m.meta.wcs.bounding_box
>>> bb is None
True

So, maybe it is a problem with cube builder? CC: @jemorrison @hbushouse

Also, this is for a piece of code that should work independent of WCS so I do not want to have dependencies on gwcs. Most importantly though, does it matter how I construct indices? The result should be the same.

[jdavies-st]

I have not tried the above solution as I don't have any cubes to play with. That said, if the cubes have wcs objects without bounding boxes, then I would think that is a bug? They should have them.

You could also generate a bounding box from the shape of the cube. I use a function jwst.resample.resample_utils.bounding_box_from_shape() to generate bounding boxes for my output frames when I do resampling. Perhaps this is of more general interest than just for resample? I also realize this is only valid for 2D images, and not for 3D cubes.

[mcara]

@jdavies-st So, I tested your idea. Here is what I get:

>>> bb = resample.resample_utils.bounding_box_from_shape((372 - 1, 32 - 1, 34 - 1))
>>> x1, y1, z1 = gwcs.wcstools.grid_from_bounding_box(bb, step=(1, 1, 1))
>>> x2, y2, z2 = np.indices((372, 32, 34))[::-1]
>>> np.allclose(x1, x2)
True

So, since the results are the same, I do not see how using gwcs.wcstools and/or jwst.resample would change the fact that gwcs's forward transform does not work with multi-dimensional data.

[jdavies-st]

Bummer. I would verify that that WCS object in the cube is not malformed, i.e. that it takes the 3 inputs and does something with each of them for each step in the transform. Perhaps @jemorrison and @nden can have a look at the cube's WCS object.

[mcara]

After some discussion with @jdavies-st and digging some more into the offending WCS of the cube obtained from cube_builder, I get:

<WCS(output_frame=world, input_frame=None, forward_transform=Model: CompoundModel584
Inputs: (u'x0', u'x1', u'x')
Outputs: (u'alpha_C', u'delta_C', u'x')
Model set size: 1
Expression: ([0] & [1] | [2] | [3] & [4] | [5] | [6]) & ([7] | [8] | [9])
Components: 
    [0]: <Shift(offset=-17.0, name=u'crpix1')>
    [1]: <Shift(offset=-16.0, name=u'crpix2')>
    [2]: <AffineTransformation2D(matrix=[[ 1., 0.], [ 0., 1.]], translation=[ 0., 0.], name=u'pc_matrix')>
    [3]: <Scale(factor=3.6111109786563447e-05, name=u'cdelt1')>
    [4]: <Scale(factor=3.6111109786563447e-05, name=u'cdelt2')>
    [5]: <Pix2Sky_Gnomonic(name=u'TAN')>
    [6]: <RotateNative2Celestial(lon=44.99995636409714, lat=0.00017278673100025292, lon_pole=180.0, name=u'crval')>
    [7]: <Shift(offset=-0.0)>
    [8]: <Scale(factor=0.0024999999441206455)>
    [9]: <Shift(offset=4.846315415774628)>

This does not look right, in particular Expression: ([0] & [1] |... indicates that only two coordinates are being used.

I have a feeling it is the cube builder that does not build a correct WCS for the output cubes.

CC: @jemorrison @hbushouse

[jemorrison]

I have a test routine that builds cubes with single = True I then print out the wcsinfo from the returned single mapped IFU cube the WCS information is correct. What are you getting that is not correct ? When I print the data to the FITS file the WCS is also correct. So I am a bit confused about the what you are trying to do.

example of test code all printed out values are correct.

single_IFUCube_result = CubeBuildStep.call(container, config_file=config, channel=ch, single='true')

nsize = single_IFUCube_result[0].data.shape

naxis1 = nsize[2] naxis2 = nsize[1] naxis3 = nsize[0] print(naxis1,naxis2,naxis3)

for j in range(len(single_IFUCube_result)) : modelnum = str(j) single = single_IFUCube_result[j] print('test_modelcontainer',single.meta.filename) print('instrument',single.meta.instrument.name) print('wcsinfo crpix ',single.meta.wcsinfo.crpix1,single.meta.wcsinfo.crpix2, single.meta.wcsinfo.crpix3) print('wcsinfo crval ',single.meta.wcsinfo.crval1,single.meta.wcsinfo.crval2, single.meta.wcsinfo.crval3) print('wcsinfo cdelt ',single.meta.wcsinfo.cdelt1,single.meta.wcsinfo.cdelt2, single.meta.wcsinfo.cdelt3)

On 07/06/2017 12:06 PM, Mihai Cara wrote:

After some discussion with @jdavies-st https://github.com/jdavies-st and digging into the offending WCS of the cube obtained from |cube_builder|, I get:

|<WCS(output_frame=world, input_frame=None, forward_transform=Model: CompoundModel584 Inputs: (u'x0', u'x1', u'x') Outputs: (u'alpha_C', u'delta_C', u'x') Model set size: 1 Expression: ([0] & [1] | [2] | [3] & [4] | [5] | [6]) & ([7] | [8] | [9]) Components: [0]: <Shift(offset=-17.0, name=u'crpix1')> [1]: <Shift(offset=-16.0, name=u'crpix2')> [2]: <AffineTransformation2D(matrix=[[ 1., 0.], [ 0., 1.]], translation=[ 0., 0.], name=u'pc_matrix')> [3]: <Scale(factor=3.6111109786563447e-05, name=u'cdelt1')> [4]: <Scale(factor=3.6111109786563447e-05, name=u'cdelt2')> [5]: <Pix2Sky_Gnomonic(name=u'TAN')> [6]: <RotateNative2Celestial(lon=44.99995636409714, lat=0.00017278673100025292, lon_pole=180.0, name=u'crval')> [7]: <Shift(offset=-0.0)> [8]: <Scale(factor=0.0024999999441206455)> [9]: <Shift(offset=4.846315415774628)> |

This does not look right, in particular |Expression: ([0] & [1] |...| indicates that only two coordinates are being used.

I have a feeling it is the cube builder that does not build a correct WCS for the output cubes.

CC: @jemorrison https://github.com/jemorrison @hbushouse https://github.com/hbushouse

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/spacetelescope/gwcs/issues/89#issuecomment-313490231, or mute the thread https://github.com/notifications/unsubscribe-auth/ARonWf7aG2hahsEnnjd-Ta9A9VTOyIKWks5sLTA2gaJpZM4OAKbX.

[mcara]

@jemorrison Could you, please, run the two commands in https://github.com/spacetelescope/gwcs/issues/89#issue-237346629:

>>> ind = np.indices((372,32,34))[::-1]
>>> cube.meta.wcs(*ind) 

where cube is a CubeModel produced by the cube builder, on one of your cubes? This may show the problem (I hope).

[jemorrison]

Ok one quick thing first. The cubes produced by cube_builder are IFU cubes (IFUCubeModel) defined by the datamodel in ifucube.py not cube.py (CubeModel)

cube.py defines the ramp data - a series of 2-D images. Is that the problem ?

Jane

On 07/06/2017 12:51 PM, Mihai Cara wrote:

@jemorrison https://github.com/jemorrison Could you, please, run the two commands in #89 (comment) https://github.com/spacetelescope/gwcs/issues/89#issue-237346629:

ind= np.indices((372,32,34))[::-1] cube.meta.wcs(*ind)

where |cube| is a |CubeModel| produced by the cube builder, on one of your cubes? This may show the problem (I hope).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/spacetelescope/gwcs/issues/89#issuecomment-313500845, or mute the thread https://github.com/notifications/unsubscribe-auth/ARonWceTj4duPXoz6kCDcKszRxZuR5rMks5sLTqjgaJpZM4OAKbX.

[mcara]

@jemorrison No, I already use IFUCubeModel:

>>> from jwst.datamodels import IFUCubeModel
>>> cube = IFUCubeModel('MIRM103-Q0-SHORT_495_rate_bsub_cal_ch1-short_s3d.fits')
>>> cube.meta.wcs
<WCS(output_frame=world, input_frame=None, forward_transform=Model: CompoundModel602
Inputs: (u'x0', u'x1', u'x')
Outputs: (u'alpha_C', u'delta_C', u'x')
Model set size: 1
Expression: ([0] & [1] | [2] | [3] & [4] | [5] | [6]) & ([7] | [8] | [9])

[jemorrison]

Ok maybe I am not understanding something. But the wcs information about the IFUCube is stored in cube.meta.wcsinfo (not just wcs)

I thought that was the way it was defined ? Jane

On 07/06/2017 01:12 PM, Mihai Cara wrote:

@jemorrison https://github.com/jemorrison No, I already use |IFUCubeModel|:

from jwst.datamodelsimport IFUCubeModel cube= IFUCubeModel('MIRM103-Q0-SHORT_495_rate_bsub_cal_ch1-short_s3d.fits') cube.meta.wcs <WCS(output_frame=world,input_frame=None,forward_transform=Model: CompoundModel602 Inputs: (u'x0',u'x1',u'x') Outputs: (u'alpha_C',u'delta_C',u'x') Modelset size:1 Expression: ([0]& [1]| [2]| [3]& [4]| [5]| [6])& ([7]| [8]| [9])

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/spacetelescope/gwcs/issues/89#issuecomment-313506155, or mute the thread https://github.com/notifications/unsubscribe-auth/ARonWXiXdsE_8h21Cem2uPeq6EEANiKyks5sLT-vgaJpZM4OAKbX.

[nden]

It looks like the WCS is trying to apply a 2D rotation matrix to more than two inputs. I think it's something to look at in cube_build. So I'll close this here and open an issue/PR on jwst.

[jemorrison]

Looking at the core schema - it seems that meta.wcs holds the asdf wcs from the input images ? The IFUCube WCS is stored in meta.wcsinfo

Does this solve the problem Mihai ? Jane

On 07/06/2017 01:12 PM, Mihai Cara wrote:

@jemorrison https://github.com/jemorrison No, I already use |IFUCubeModel|:

from jwst.datamodelsimport IFUCubeModel cube= IFUCubeModel('MIRM103-Q0-SHORT_495_rate_bsub_cal_ch1-short_s3d.fits') cube.meta.wcs <WCS(output_frame=world,input_frame=None,forward_transform=Model: CompoundModel602 Inputs: (u'x0',u'x1',u'x') Outputs: (u'alpha_C',u'delta_C',u'x') Modelset size:1 Expression: ([0]& [1]| [2]| [3]& [4]| [5]| [6])& ([7]| [8]| [9])

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/spacetelescope/gwcs/issues/89#issuecomment-313506155, or mute the thread https://github.com/notifications/unsubscribe-auth/ARonWXiXdsE_8h21Cem2uPeq6EEANiKyks5sLT-vgaJpZM4OAKbX.

[mcara]

• Shouldn't meta.wcsinfo be consistent with meta.wcs?
• wcsinfo is not a WCS object (just WCSLIB-stype WCS info - my understanding). I need a complete WCS object represented by meta.wcs.

[jemorrison]

It was my understanding that I filled in meta.wcsinfo for the IFUCube wcs If that is not the case or if I need to also update meta.wcs I can do that. Howard could you advise ? Jane

On 07/06/2017 03:42 PM, Mihai Cara wrote:

• Shouldn't |meta.wcsinfo| be consistent with |meta.wcs|?
• |wcsinfo| is not a WCS object (just WCSLIB-stype WCS info - my understanding). I need a complete WCS /object/ represented by |meta.wcs|.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/spacetelescope/gwcs/issues/89#issuecomment-313539077, or mute the thread https://github.com/notifications/unsubscribe-auth/ARonWZgYM9vFW9jNdO_uzWJ8qYdJhi5Gks5sLWLWgaJpZM4OAKbX.

[hbushouse]

Not really, since I'm not the WCS expert. But I guess I contribute this much: given that the IFU cubes are constructed on a regular, linearized frame, their WCS should be pretty simple to document and fully describable by the traditional FITS standard WCS keywords (e.g. crpix, crval, cdelt, etc.), which as far as I know is what's contained in model.meta.wcsinfo. So documenting the cube's WCS in the wcsinfo object is valid, but I agree with Mihai that it should also always be documented in the "official" (gwcs-based) WCS object contained in model.meta.wcs. In the case of a nice regular frame like the cubes have, the results from both should be identical. But the point of the gwcs-based wcs object is that it also provides methods that allow for easy transformation back and forth between pixel and sky frames, which is what Mihai apparently needs. I don't believe the simple wcsinfo object allows for that.

Of course to get the correct answer to all of this, we need to pull in the big guns: @nden

[nden]

OK, after some investigation I confirmed that the WCS is correct (meaning it has all necessary components, I don't know about results). The problem (which I think was discussed before elsewhere) is caused by the fact that the inputs are 3D cubes. I will keep this issue open until there's resolution but I don't see this as a high priority. For now I think it's best the calling function flattens the input arrays before passing them to the WCS object. I don't want to do this in GWCS yet because I am not sure I have thought about all use cases. And it's really an issue with models.

[nden]

One other thing, @jemorrison , you may want to add a bounding_box to the WCS. Then indices can be created using wcstools.grid_from_bounding_box.

[stscijgbot]

This ticket is now being tracked at AL-61