SkyCircleRegion matplotlib plot issue

[cdeil]

In #74 @joleroi added docs/plot_compound.py which generates this plot: http://astropy-regions.readthedocs.io/en/latest//plot_compound.png

As discussed in #74, there is an issue, the region is not drawn correctly. To @joleroi and me the script looks OK, so we don't know what's wrong.

@astrofrog or @larrybradley or anyone - if you could have a look, that would be great!

[larrybradley]

The plotting script also looks OK to me, but the plot is wrong. From looking at the code and the plot results, I think the problem lies somewhere in the WCSAxes transform machinery. Since the plotted regions are too small, I tried manually changing the radius, but quickly realized that doesn't work (bad points also appear within the shape). It appears that the distorted/transformed region size and shape are wrong. @astrofrog?

[astrofrog]

I need to look into this - no time today but will look asap

[adonath]

@cdeil @joleroi There is a bug in the plotting script as well. The sky coordinates should be computed using the WCS from the example dataset and not from a new regular grid (see here).

From the plot one can see that the regions are not circular any more, because of the projection defined in the wcs of the dataset. I presume using SkyCoord.from_pixel(...wcs=wcs) and passing the right WCS should fix the issue with the pixels contained or not contained in the region.

[cdeil]

@adonath - I don't think so. To me it still looks like the circle is simply drawn incorrectly.

Yes, one could create SkyCoords from PixCoords, but that won't change the SkyCircle's being drawn incorrectly. Note that the circles are really distorted in a weird way, that's not consistent with it being an AIT projection with reference point at the center of the image.

@adonath - Do you still think there's a bug in the script? If yes, can you please fix it and post the new version and output plot somewhere so that I can understand exactly what was wrong?

[joleroi]

@adonath and I discussed about the script on friday.

Yes, one could create SkyCoords from PixCoords, but that won't change the SkyCircle's being drawn incorrectly.

Right, it does not change anything. The only difference is that the coordinates (points in the plot) are also distored and indeed some of the points are not marked incorrectly anymore.

What does actually fix the problem is changing the transformation to 'galactic' in the as_patch method of SkyCircle

https://github.com/astropy/regions/blob/master/regions/shapes/circle.py#L199

I don't know if there is a general way to handle this in the transformations but maybe it not necessary since @astrofrog fixed the issue in #77

[cdeil]

@joleroi - Thanks for pointing out that SkyCircleRegion.as_patch has transform=ax.get_transform('icrs') hard-coded. That means it's completely broken and could be fixed by changing to ax.get_transform('world')?

@astrofrog - In #77 you still have transform=ax.get_transform('icrs') hard-coded. But I think the discussion there is leaning towards removing SkyCircleRegion.as_patch completely?

[adonath]

@cdeil Yes, the change I proposed doesn't fix the distorted circle, but it makes the plotted coordinates aligned with the transformation the circles are plotted with. This fixes the issue with the incorrectly marked points.

Because of the hard-coded transform=ax.get_transform('icrs'), I presume that the circle is actually drawn correctly and the distortion we see is the 'natural' distortion from the icrs system at the galactic center.

Here is a plot for illustration:

The white line marks the galactic plane.

So, I'm :+1: to remove the as_patch() methods and let the users define the matplotlib patches.

[astrofrog]

The hard-coded ICRS is not a problem btw because I always convert the SkyCoord to ICRS first too.

[cdeil]

The hard-coded ICRS is not a problem btw because I always convert the SkyCoord to ICRS first too.

I think it is!? It leads to incorrectly distorted circles because you're using the wrong WCS, no?

[adonath]

@astrofrog Just to clarify my point: I think if you draw a circle at the galactic center using the icrs coordinate system (I mean not only transforming the center of the circle...), the result is exactly the one we saw in http://astropy-regions.readthedocs.io/en/latest//plot_compound.png or the plot I've made. The result will be distorted, because of the intrinsic distortion of the coordinate system. So it's just the confusion of the two coordinate systems that leads to an unexpected result.

[astrofrog]

@adonath - just to be clear, I think the plot you showed in your previous command is actually correct (unlike the plotting issue we saw before). If you change the code in my pull request to have galactic instead of ICRS, I'm pretty sure the result will be the same (it looks like the projection you are using is not centered on the Galactic center, so intrinsic spherical distortion is expected regardless of coordinate system?)

[adonath]

@astrofrog Here's the code I used to produce my plot:

import numpy as np
from matplotlib.patches import Circle
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
from astropy.wcs import WCS

wcs_spec =  {'CDELT1': -1.0,
             'CDELT2': 1.0,
             'CRPIX1': 8.5,
             'CRPIX2': 8.5,
             'CRVAL1': 266.416833,
             'CRVAL2': -29.007806,
             'CTYPE1': 'RA---AIT',
             'CTYPE2': 'DEC--AIT',
             'CUNIT1': 'deg',
             'CUNIT2': 'deg'}

wcs = WCS(wcs_spec)

fig = plt.figure(figsize=(6, 6))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs)

c_1 = Circle((265.07333228, -27.0343367), 5, fc='none', lw=2, transform=ax.get_transform('icrs'), ec='r')
c_2 = Circle((261.003259, -30.6748939), 3, fc='none', lw=2, transform=ax.get_transform('icrs'), ec='y')

ax.imshow(np.zeros((16, 16)), cmap='gray')
ax.add_patch(c_1)
ax.add_patch(c_2)
ax.grid('on', color='w')

l = np.linspace(-10, 10, 100)
b = np.zeros_like(l)
gp = SkyCoord(l, b, frame='galactic', unit='deg').transform_to('icrs')

# add galactic plane
ax.plot(gp.ra.deg, gp.dec.deg, c='w', lw=1, transform=ax.get_transform('icrs'))
ax.set_xlim(0, 15)
ax.set_ylim(0, 15)

plt.savefig('circle_gc_icrs.png', dpi=180)

So the projection is centered on the galactic center (only given in icrs coordinates). And it seems like this is the kind of projection that ax.get_transform('icrs') returns?

[astrofrog]

@adonath - ah yes, that code is wrong, sorry for the confusion. You should use SphericalCircle instead of Circle as in https://github.com/astropy/regions/pull/77

[adonath]

@astrofrog Hm, but the code I've posted in my previous comment, gives exactly the result I expect. The circle has a radius of 5 deg and is distorted according to (at that position) the coordinate system I use for plotting. And this is independent whether I plot a Circle or Rectangle or Polygon, so I don't see how SpericalCircle would make a difference in that case.

The same works the other way around: if I add a circle to my plot, drawn in the galactic coordinate system, the circle is not distorted, as I would expect from the coordinate system at that position:

# draw circle in galactic coordinate system
c_1_g = Circle((1, 2), 5, fc='none', lw=2, transform=ax.get_transform('galactic'), ec='r', ls='--')
c_2_g = Circle((-4, 3), 3, fc='none', lw=2, transform=ax.get_transform('galactic'), ec='y', ls='--')

ax.add_patch(c_1_g)
ax.add_patch(c_2_g)

The circles drawn with the galactic coordinate transformation are shown dashed.

[cdeil]

I've started a new PR in #85 and will try to finish it quickly tomorrow (not much time this week). I'll try to review all the discussion and make docs examples out of the nice examples posted here (like the last one by @adonath).

[cdeil]

I just now merged #85, which fixes this issue.

There's a lot of good discussion and examples here, that should result in further code or docs improvements. So for that reason, I'm keeping the issue open and moving it to the 0.3 milestone.

[astrofrog]

@adonath:

The circle has a radius of 5 deg and is distorted according to (at that position) the coordinate system I use for plotting. And this is independent whether I plot a Circle or Rectangle or Polygon, so I don't see how SpericalCircle would make a difference in that case.

The same works the other way around: if I add a circle to my plot, drawn in the galactic coordinate system, the circle is not distorted, as I would expect from the coordinate system at that position:

So the main issue when using Circle/Rectangle/Polygon is that Matplotlib treats the RA/Dec system as a euclidean system and won't take into account (to first order approximation) the cos(lat) factor. Thus, shapes close to the poles will be very squashed along longitude. The reason you don't see this in the Galactic frame is that you are close to the Galactic plane.

SphericalCircle gets around this by drawing a true spherical circle defined as a polygon so that there is no opportunity for Matplotlib to mess it up. Does that make more sense?

[cdeil]

SphericalCircle gets around this by drawing a true spherical circle defined as a polygon so that there is no opportunity for Matplotlib to mess it up.

So the proper long-term solution is to say people should use regions, and to add as_poly methods to SkyRegion or all region objects? This could even start now, with the code from wcsaxes.SphericalCircle moving to regions.SkyCircleRegion.as_poly?

[astrofrog]

So the proper long-term solution is to say people should use regions, and to add as_poly methods to SkyRegion or all region objects?

No, to_pixel(mode='full') does that conversion to polygons, no need for a new method: https://github.com/astropy/regions/pull/91 (and that implements the code from SphericalCircle)

[cdeil]

No, to_pixel(mode='full') does that conversion to polygons, no need for a new method: #91 (and that implements the code from SphericalCircle)

I don't think that's good API. Polygonisation should be a separate step from sky-pixel conversion, and be reflected in a separate method. Of course it might be useful to re-expose it in one step, but the basic two steps should be separately available. I'll try to find time to look at #91 to see what you put.

[astrofrog]

@ceil - the conversion to polygon is a necessary side effect though. From the user's perspective, it's just that they are converting the polygon from sky to pixel coordinates using a better algorithm, so that was the idea with the to_pixel(mode='full') API. But I agree that internally we could have a separate as_poly method to do the polygon conversion. I guess I'm just saying I don't want users to have to do reg.as_poly().to_pixel().

[adonath]

@astrofrog Unfortunately it still doesn't make completely sense to me. As far as I can tell the circle is drawn completely correct into to the coordinate system on the sky (as shown by the grid lines in the figures above) in any case.

As far as I know, internally when drawing matplotlib treats patches as a Path (i.e. rendered as a polygon) and applies any defined transformation to this path. When the patch is declared in the icrs sky coordinate system it means the circle will be squashed along longitude when getting closer to the poles, the same way the coordinate grid lines are squashed. Which seems to be the correct behaviour to me.

What I realize now is that there is the use case to draw an actual circle in the plot e.g. to highlight a source or region, independent of the projection (or transformation) of the image. Which is equivalent to only transforming the circle's centre and converting the circle radius using the local pixel size. And this can be achieved by the SphericalCircle. Is that correct?

The same could be achieved by introducing something like as.get_transform(mode='local'), which would return the local euclidean coordinate system to e.g. the image centre. This would be the same idea as the different modes in .to_pixels(mode=['full', 'local', ...]), but only applied to plotting.

As the different transformation modes are needed anyway for regions (e.g. to compute masks), I think the solution proposed in #91 is totally fine for plotting as well. So apply the transformations in SkyRegion.to_pixel() and then only plot the pixel region. But still I don't see where the matplotlib behaviour would be incorrect.

[cdeil]

I guess I'm just saying I don't want users to have to do reg.as_poly().to_pixel().

OK. But as a dev I want it available as two steps, to have code that's simpler to test and maintain. Of course we can discuss what should be private or public, but I think reg.as_poly().to_pixel() should be the internal implementation.

There could even be different algorithms for sky and pixel regions to polygonise. MPL has some polygonisation code, so does sphere I think.

[cdeil]

Should we schedule a Google hangout tomorrow or next week to discuss this, and generally status / next steps for astropy.regions?

[astrofrog]

+1 to a hangout today/tomorrow/next week - I think these issues will be easier to discuss that way

[adonath]

+1 for the hangout

[cdeil]

@astropy/regions-developers - Let's do a Google hangout early next week to discuss Astropy regions (specifically sky region plotting / connection to wcsaxes, and if we're ready to release 0.2 or not). Anyone interested - please fill this doodle to find a time: http://doodle.com/poll/yi9unxmbqdfhhfib I'll close it tomorrow.

[cdeil]

I've pair-coded and discussed an example with @adonath : http://nbviewer.jupyter.org/gist/cdeil/1df42de70326d577e7964be15b2a7396

This resolves the point from @adonath earlier:

But still I don't see where the matplotlib behaviour would be incorrect.

Matplotlib can't polygonise the sky circle correctly. Having the WCS is not sufficient. One needs to know that it's a "sky" circle.

This is what @astrofrog implemented in #91 : correct polygonisation of sky circles.

So I think the conclusion here is still: no concrete action item, we only keep this (very long) issue open as a reminder that we need to write good docs to explain how our regions are exactly plotted with MPL.

[cdeil]

We will have the Astropy regions hangout next week Tue 12/6/16 3:00 PM CEST.

@astrofrog - Should I "book" a room or should we just connect ad-hoc? I've had technical problems either way, not sure what is best.

[astrofrog]

@cdeil - just create a hangout at that time and invite those attending - less hassle

[cdeil]

Astropy regions hangout now: https://hangouts.google.com/call/fif7cot7srabtbxduatoiedthme

[cdeil]

Some of the discussion here led to changes implemented by @astrofrog in #132 . The latest version of the docs is at http://astropy-regions.readthedocs.io/en/latest/

Should we close this issue? Or alternatively - can someone split out concrete, actionable suggestions for further changes to other issues? (the discussion above is very long and hard to re-read and get back to, so I'd prefer to close the issue)

[astrofrog]

@cdeil - I'll try and review this and open better issues later today or tomorrow