Bug: Circle appear as an ellipses

[hamogu]

Following the docs, I'm overplotting patches onto a wcs axis image. In this image the pixels are square and a circle looks like a circle when I plot it in ds9. However, a circle overplotted with wcsaxis looks like an ellipse. I suspect that there is a cos(delta) missing somewhere.

# Note:
# This example is taken from a larger script and will not run as is.
# I can supply the data by email, if needed to reproduce,
# but I think that any small scale image will do.
# (Maybe images on the equator would work, but I don't have anything like that at hand).
t = F621M[i]

fig = plt.figure()
ax = WCSAxes(fig, [0.1, 0.1, 0.8, 0.8], wcs=t.wcs)
fig.add_axes(ax) 
im = ax.imshow(F621Marr[..., i], cmap=plt.cm.gist_heat, origin='lower', interpolation='nearest')
plt.colorbar(im, ax=ax)
ra = ax.coords['ra']
dec = ax.coords['dec']
ra.grid(color='b', alpha=0.5, linestyle='solid', linewidth=3)
dec.grid(color='b', alpha=0.5, linestyle='solid', linewidth=3)
ra.set_major_formatter('hh:mm:ss.s')

circ = mpl.patches.Circle(posindeg, radius=1./3600., edgecolor='green', facecolor='none',
              transform=ax.get_transform('icrs'))
ax.add_patch(circ)

Note: The image is slightly confusing. Look carefully at the labels to see which lines in the grid are Ra and Dec. The left axis corresponds to the lines that go up and down and the bottom axis to the lines that go left to right (only one of which intersects the bottom axis, thus only one label is plotted).

[hamogu]

I observe that the radius in delta (almost left to right in the image) is 1 arcsec as requested.

[astrofrog]

Hmm, that's an interesting point - we're actually not explicitly computing the new coordinates of the circle but we're letting matplotlib place the points in world coordinates and then we convert that to pixel coordinates. All get_transform does is provide a simple transformation function from ICRS to pixel coordinates, the rest is all matplotlib (not WCSAxes).

If you imagine using the ax.plot command to plot lines between coordinates, this will work correctly, but it's for the special case where we are talking about shapes with radii/sizes that we need to be careful.

In the mean time, a workaround is to add an ellipse where you apply the cos(dec) factor yourself to the RA axis.

[hamogu]

It turn out that in my application it's accentually more convenient to plot the shape in pixel space, it's just that I realized that only after I tried to plot in WCS. That's the upside: I discovered this problem, but it turn out it's not actually a problem for me ;-)

The example in the docs is a rectangle placed very close to lat, lon = 0, but I suspect that in general the aspect ratio of the plotted rectangle will not be exactly the number it should be either.

[hamogu]

If you decide that there is no easy fix, I suggest to add a section to the docs that explains what you said above (we let matplotlib calculate the coordinates and just transform that to pixels) with a list of work-arounds (although that section might be called "Examples where user input isn eed to get coordianates right" or something like that).

Do you want me to take a stab at that for the docs or do you think there is a code fix that will make it "just work"?

Related: Imagine a WCS for a long-slit spectrum with x = distance in arcsec and y = wavelength in Ang. Matplotlib would happily plot a circle into that, although it does not make sense physically. Should that raise an error? Should it just work, because "if the user asks for it, he/she hopefully knows what they want"?

[astrofrog]

Just for my records, here's an example that reproduces the issue with one of the built-in datasets:

from wcsaxes.datasets import fetch_twoMASS_k_hdu

from matplotlib import pyplot as plt
from matplotlib.patches import Circle

from astropy.wcs import WCS

hdu = fetch_twoMASS_k_hdu()
wcs = WCS(hdu.header)

fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs)
fig.add_axes(ax) 
im = ax.imshow(hdu.data, cmap=plt.cm.viridis, origin='lower', interpolation='nearest')
plt.colorbar(im, ax=ax)
ra = ax.coords['ra']
dec = ax.coords['dec']
ra.set_major_formatter('d.dddd')
ra.set_major_formatter('d.dddd')

circ = Circle((266.5, -29), radius=100./3600., edgecolor='green', facecolor='none',
              transform=ax.get_transform('icrs'))
ax.add_patch(circ)

fig.savefig('test_circle.png')

[astrofrog]

And here's a mention of an issue someone had with plotting a circle on a polar plot:

http://stackoverflow.com/questions/19827792/matplotlib-drawing-a-smooth-circle-in-a-polar-plot

[astrofrog]

Fixed by #199