Plotting onto a WCSAxes using get_transform fails when wcs has >2 axes

[astrofrog]

Hi all,

I am constructing a WCSAxes object like this:

ax_lv = WCSAxes(fig, ax_lv_limits, wcs=self.dendrogram.wcs, slices=('x', 0, 'y')) # `dendrogram` is some object that contains a wcs derived from an (l, b, v) datacube.
fig.add_axes(ax_lv)

and then calling "plot" (hoping to plot in coordinate space rather than pixel space) like this:

ax_lv.plot( np.arange(10), transform=ax_lv.get_transform('world') )
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-10-8a6422a25f67> in <module>()
----> 1 iv.ax_lv.plot(np.arange(10), transform=iv.ax_lv.get_transform('world')
      2 )

/Users/tsrice/Documents/Code/wcsaxes/wcsaxes/core.pyc in get_transform(self, frame, equinox, obstime)
    232                   the specified frame to the pixel/data coordinates.
    233         """
--> 234         return self._get_transform_no_transdata(frame, equinox=equinox, obstime=obstime).inverted() + self.transData
    235
    236     def _get_transform_no_transdata(self, frame, equinox=None, obstime=None):

/Users/tsrice/Documents/Code/wcsaxes/wcsaxes/transforms.pyc in inverted(self)
    170         Return the inverse of the transform
    171         """
--> 172         return WCSWorld2PixelTransform(self.wcs, slice=self.slice)
    173
    174

/Users/tsrice/Documents/Code/wcsaxes/wcsaxes/transforms.pyc in __init__(self, wcs, slice)
     69         if self.wcs.wcs.naxis > 2:
     70             if slice is None:
---> 71                 raise ValueError("WCS has more than 2 dimensions, so ``slice`` should be set")
     72             elif len(slice) != self.wcs.wcs.naxis:
     73                 raise ValueError("slice should have as many elements as WCS "

ValueError: WCS has more than 2 dimensions, so ``slice`` should be set

I then get the above error. I have tried passing slice parameters to get_transform to no avail. Can I ask if I am using this functionality incorrectly, or if this is a bug?

[astrofrog]

@tomr-stargazer - could you check if the attached code fixes it for you? (no need to specify slices, it should used the one defined when initializing WCSAxes)

[tomr-stargazer]

Hi @astrofrog,
I'm still hitting that error. I created a file that contains exactly the code that's erroring (to simplify testing) here: https://github.com/tomr-stargazer/astrodendro_analysis/blob/master/wcs_transform_testing.py

the action code is here:

fig = plt.figure()
ax = wcsaxes.WCSAxes(fig, [0.1 ,0.1 ,0.8, 0.8], wcs=wcs_object, slices=('x', 0, 'y'))
ax.plot(np.arange(10), np.arange(10), transform=ax.get_transform('world'))

and produces this error output:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/Users/tsrice/Documents/Code/astrodendro_analysis/wcs_transform_testing.py in <module>()
     49 ax = wcsaxes.WCSAxes(fig, [0.1 ,0.1 ,0.8, 0.8], wcs=wcs_object, slices=('x', 0, 'y'))
     50
---> 51 ax.plot(np.arange(10), np.arange(10), transform=ax.get_transform('world'))
     52

/Users/tsrice/anaconda/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in plot(self, *args, **kwargs)
   1372
   1373         for line in self._get_lines(*args, **kwargs):
-> 1374             self.add_line(line)
   1375             lines.append(line)
   1376

/Users/tsrice/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in add_line(self, line)
   1484             line.set_clip_path(self.patch)
   1485
-> 1486         self._update_line_limits(line)
   1487         if not line.get_label():
   1488             line.set_label('_line%d' % len(self.lines))

/Users/tsrice/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _update_line_limits(self, line)
   1514             if self.transData.is_affine:
   1515                 line_trans_path = line._get_transformed_path()
-> 1516                 na_path, _ = line_trans_path.get_transformed_path_and_affine()
   1517                 data_path = trans_to_data.transform_path_affine(na_path)
   1518             else:

/Users/tsrice/anaconda/lib/python2.7/site-packages/matplotlib/transforms.pyc in get_transformed_path_and_affine(self)
   2600         the path necessary to complete the transformation.
   2601         """
-> 2602         self._revalidate()
   2603         return self._transformed_path, self.get_affine()
   2604

/Users/tsrice/anaconda/lib/python2.7/site-packages/matplotlib/transforms.pyc in _revalidate(self)
   2577             or self._transformed_path is None):
   2578             self._transformed_path = \
-> 2579                 self._transform.transform_path_non_affine(self._path)
   2580             self._transformed_points = \
   2581                 Path(self._transform.transform_non_affine(self._path.vertices),

/Users/tsrice/anaconda/lib/python2.7/site-packages/matplotlib/transforms.pyc in transform_path_non_affine(self, path)
   2265             return path
   2266         elif not self._a.is_affine and self._b.is_affine:
-> 2267             return self._a.transform_path_non_affine(path)
   2268         else:
   2269             return self._b.transform_path_non_affine(

/Users/tsrice/Documents/Code/astrodendro_analysis/wcsaxes/transforms.pyc in transform_path(self, path)
     46             The resulting path
     47         """
---> 48         return Path(self.transform(path.vertices), path.codes)
     49
     50     transform_path_non_affine = transform_path

/Users/tsrice/Documents/Code/astrodendro_analysis/wcsaxes/transforms.pyc in transform(self, world)
     89
     90         if world.shape[1] != self.wcs.wcs.naxis:
---> 91             raise ValueError("Second dimension of input values should match number of WCS coordinates")
     92
     93         pixel = self.wcs.wcs_world2pix(world, 1) - 1

ValueError: Second dimension of input values should match number of WCS coordinates```

[astrofrog]

Hmm, it turns out that this is a limitation of the matplotlib and WCS model. The issue is that in order to convert world coordinates to pixel coordinates, one needs three world coordinates in the case of a 3-d WCS, otherwise there is a degeneracy along the third dimension. However, matplotlib's plot can only take two values. I will raise an issue with the matplotlib folk to see if they have any ideas, but in the mean time you can do:

tr = ax.get_transform('world')
pixel = tr.transform(np.array([np.arange(10),
                               np.arange(10),
                               np.repeat(-90, 10)]).transpose())
x, y = pixel[:,0], pixel[:,1]
ax.plot(x, y)

[astrofrog]

Matplotlib issue: matplotlib/matplotlib#3983

[tomr-stargazer]

Thanks very much Tom for following up on this! I'll see about using the
temporary fix.

On Thu, Jan 8, 2015 at 3:17 PM, Thomas Robitaille notifications@github.com
wrote:

Matplotlib issue: matplotlib/matplotlib#3983
matplotlib/matplotlib#3983

—
Reply to this email directly or view it on GitHub
#148 (comment).

[tomr-stargazer]

Hi Tom -

it turns out even the temporary code doesn't work:

----> 1 ax.get_transform('world')

/Users/tsrice/anaconda/lib/python2.7/site-packages/wcsaxes-0.2.dev328-py2.7.egg/wcsaxes/core.pyc in get_transform(self, frame, equinox, obstime)
    181                   the specified frame to the pixel/data coordinates.
    182         """
--> 183         return self._get_transform_no_transdata(frame, equinox=equinox, obstime=obstime).inverted() + self.transData
    184
    185     def _get_transform_no_transdata(self, frame, equinox=None, obstime=None):

/Users/tsrice/anaconda/lib/python2.7/site-packages/wcsaxes-0.2.dev328-py2.7.egg/wcsaxes/transforms.pyc in inverted(self)
    167         Return the inverse of the transform
    168         """
--> 169         return WCSWorld2PixelTransform(self.wcs, slice=self.slice)
    170
    171

/Users/tsrice/anaconda/lib/python2.7/site-packages/wcsaxes-0.2.dev328-py2.7.egg/wcsaxes/transforms.pyc in __init__(self, wcs, slice)
     64         if self.wcs.wcs.naxis > 2:
     65             if slice is None:
---> 66                 raise ValueError("WCS has more than 2 dimensions, so ``slice`` should be set")
     67             elif len(slice) != self.wcs.wcs.naxis:
     68                 raise ValueError("slice should have as many elements as WCS "

ValueError: WCS has more than 2 dimensions, so ``slice`` should be set

My current workaround is to just grab the actual WCS object and call its wcs_world2pix function:

    world_coordinates = np.vstack([lcen_column, np.repeat(0, len(lcen_column)), vcen_column]).T

    lv_pixels = iv.dendrogram.wcs.wcs_world2pix(world_coordinates, 0)

    l_lv_pixels = lv_pixels[:,0]
    v_lv_pixels = lv_pixels[:,2]

which finally works.

(this is a note for posterity, I suppose!)

[astrofrog]

@tomr-stargazer - are you using this branch with the workaround? It looks like you are using an old version of wcsaxes (it says 0.2 in the path of the package). If you try with this branch I think the workaround should work.

[tomr-stargazer]

Ack - embarrassing on my part, I'd been working on a machine with an old
wcsaxes. I'll update in the morning and let you know. Sorry for the
(likely) false alarm!
On Jan 30, 2015 4:46 AM, "Thomas Robitaille" notifications@github.com
wrote:

@tomr-stargazer https://github.com/tomr-stargazer - are you using this
branch with the workaround? It looks like you are using an old version of
wcsaxes (it says 0.2 in the path of the package). If you try with this
branch I think the workaround should work.

—
Reply to this email directly or view it on GitHub
#148 (comment).

[astrofrog]

@tomr-stargazer - just to check, did you have a chance to try this?

[tomr-stargazer]

@astrofrog
Hi Tom,
sorry this took me two(!) months to address!
The code no longer crashes, and seems to provide sensible output when I make the following adjustment:

pixel = tr.transform(np.array([np.arange(10),
                               np.repeat(-90, 10),
                               np.arange(10)]).transpose())

However, it seems that the math in the transform is incorrect! (or I am using the output very incorrectly.)

To demonstrate, I have updated my script "wcs_transform_testing.py" from before.

For the above WCS, this method transforms the world coordinates into pixel coordinates sensibly

world_coordinates = np.vstack([np.arange(10), np.repeat(0, 10), np.arange(10)]).T

pixel_from_wcs = wcs_object.wcs_world2pix(world_coordinates, 0)

x_from_wcs = pixel_from_wcs[:,0]
y_from_wcs = pixel_from_wcs[:,2]

ax.plot(x_from_wcs, y_from_wcs, 'r')

# produces:
# [ 300.  296.  292.  288.  284.  280.  276.  272.  268.  264.]
# [ 69.41909288  70.18809886  70.95710485  71.72611083  72.49511681 73.26412279  74.03312878  74.80213476  75.57114074  76.34014673]

whereas the suggested code from above gives wildly inaccurate values for the pixel coordinates:

tr = ax.get_transform('world')
pixel_from_tr = tr.transform(np.array([np.arange(10),
                               np.repeat(-90, 10),
                               np.arange(10)]).transpose())
x_from_tr, y_from_tr = pixel_from_tr[:,0], pixel_from_tr[:,1]

ax.plot(x_from_tr, y_from_tr, 'b')

# produces
# [ 153664.  151616.  149568.  147520.  145472.  143424.  141376.  139328.   137280.  135232.] 
# [ 26704.93166613  27000.22996355  27295.52826097  27590.82655839  27886.12485581  28181.42315323  28476.72145065  28772.01974807  29067.31804549  29362.61634292]

(in the above image, the pixel values are so far away that the WCS fails at converting them into world coordinates).

Can you spot an error in the transform I am doing in the second case?

[astrofrog]

@tomr-stargazer - I'm majorly confused because I can't even run your example, and I think it could be due to a matplotlib update... but in the mean time your two examples are not exactly the same - in the second case the second column of your array has values of -90 instead of 0. Could that be it?

[astrofrog]

Will look into all this further later today...

[astrofrog]

@tomr-stargazer - ah, I'm being silly - get_transform returns a transform that then gives display coordinates, not data coordinates, hence the difference. Will figure out a better solution.

[tomr-stargazer]

@astrofrog Thanks Tom. Might it be something as simple as replacing tr = ax.get_transform('world') -> tr = ax.get_transform('data') or something equivalent? Or is the transform not implemented in wcsaxes?

[astrofrog]

@tomr-stargazer - I think what you want is:

tr = ax._get_transform_no_transdata('world').inverted()

It's a private method, so clearly a hack. The inverted is because that actually returns pixel to world by default. Let me know if that works.

[tomr-stargazer]

@astrofrog
That works as desired, actually! Now the outputs are identical and sensible.

I don't think I'll use this in my production code, though, since everywhere I (personally) use a wcsaxes, I also have access to the underlying wcs object, and so

world_coordinates = np.vstack([np.arange(10), np.repeat(0, 10), np.arange(10)]).T

pixel_from_wcs = wcs_object.wcs_world2pix(world_coordinates, 0)

x_from_wcs, y_from_wcs = pixel_from_wcs[:,0], pixel_from_wcs[:,2]

is (in my eyes) slightly cleaner than resorting to a private function in this case:

tr = ax._get_transform_no_transdata('world').inverted()
pixel_from_tr = tr.transform(np.array([np.arange(10),
                               np.repeat(0, 10),
                               np.arange(10)]).transpose())
x_from_tr, y_from_tr = pixel_from_tr[:,0], pixel_from_tr[:,1]

If the wcsaxes API changes such that the ax._get_transform_no_transdata('world').inverted() function is made publicly accessible, then I might use it.

In any case I suspect this pull request is complete and could be closed.

[astrofrog]

I absolutely agree that the use of the private method should make its way into production code :)

I still want to try and find a way that avoids all this - the issue is that I think it's going to be hard for MPL to change plot to accept n coordinates.

[astrofrog]

I'm going to go ahead and merge this since it does have a needed fix, then will open an issue to keep track of what was discussed here.