Initial implementation of CircleSkyRegion.to_pixel in full mode

[astrofrog]

This will likely require some more discussion, but it is a first implementation of the full mode in CircleSkyRegion.to_pixel based on the code in WCSAxes. Basically we start by setting up a spherical circle with the right radius and then rotating it in 3D to the correct coordinates, then converting the polygon to pixel coordinates.

The main point I'd like to discuss is how to set the precision of this conversion. At the moment, the API specifies tolerance but that is ill defined, and requiring a final tolerance in pixel space would require iterating. Instead, I think we could do what I'm doing here and have a parameter that specifies the number of points to use in the polygon since in full mode the shapes always get converted to a polygon. I think in that case we should rename from tolerance to something else. This is much easier to implement.

Here's an example of this code in action:

from regions import CircleSkyRegion

from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

import matplotlib.pyplot as plt

wcs = WCS('2MASS_k.fits')

c = SkyCoord(1 * u.deg, 2 * u.deg, frame='galactic')

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, aspect='equal')

for radius in [1, 5, 10, 20]:

    reg = CircleSkyRegion(c, radius=radius * u.deg)

    pix1 = reg.to_pixel(wcs, mode='local')
    pix2 = reg.to_pixel(wcs, mode='full', tolerance=200)

    ax.add_patch(pix1.as_patch(edgecolor='blue', facecolor='none'))
    ax.add_patch(pix2.as_patch(edgecolor='green', facecolor='none'))

ax.set_xlim(-20000, 20000)
ax.set_ylim(-20000, 20000)
fig.savefig('circles.png')

This also needs some kind of tests of course - not sure what the best way to test this would be?

@cdeil - obviously this doesn't have to go in 0.2, so no need to wait for it.

[astrofrog]

@adonath - you might want to try this out with your examples in #76?

[cdeil]

@astrofrog - Thanks!

This is what you currently do in SkyPolygonRegion.to_pixel(mode='full'):

        elif mode == 'full':

            # TODO: avoid converting to unit spherical or spherical if already
            #       using a spherical representation

            # Extract longitude/latitude, either from a tuple of two quantities, or
            # a single 2-element Quantity.
            rep = self.center.represent_as('unitspherical')
            longitude, latitude = rep.lon, rep.lat

            # Start off by generating the circle around the North pole
            lon = np.linspace(0., 2 * np.pi, tolerance + 1)[:-1] * u.radian
            lat = np.repeat(0.5 * np.pi - self.radius.to(u.radian).value, tolerance) * u.radian

            # Now rotate it to the correct longitude/latitude
            lon, lat = rotate_polygon(lon, lat, longitude, latitude)

            # Make a new SkyCoord
            vertices_sky = SkyCoord(lon, lat, frame=self.center)

            # Convert to PixCoord
            x, y = skycoord_to_pixel(vertices_sky, wcs)
            vertices_pix = PixCoord(x, y)

            # Make polygon
            return PolygonPixelRegion(vertices_pix)

There's clearly two steps:

1. compute sky polygon vertices
2. transform to pixel coordinates

@adonath and I think that computing the sky polygon approximation should be a separate method SkyPolygonRegion.to_poly.

This also needs some kind of tests of course - not sure what the best way to test this would be?

If it's a separate method, testing would be easy: compute the sky separate of the SkyPolygonRegion to the circle center --- they should all be equal to the radius.

---

I'm not convinced that this:

sky_circle.to_pixel(wcs, mode='full', <other pars for precision>)

is better than

sky_circle.to_poly(<pars for precision>).to_pixel(wcs)

We can discuss whether to_pixel(mode='full') is good or bad API in the telcon next week.
But to have more modular and easier to test code, can you please separate out the to_poly method for now?

[cdeil]

FYI: I've put this on the 0.2 milestone.
If we're getting such great pull requests, I think we should finish and merge them and delay the release to later next week or even later.

[astrofrog]

@cdeil - done (I called it to_polygon since it's not much longer, but we can discuss this on hangout)

[cdeil]

@astrofrog - Thanks!

Do you want to add a test or docs example here?
If no - I'd like to merge when travis-ci passes, just in case I have time to continue with regions on the weekend.

As far as I can see for several sky regions (rectangle, ellipse) the contains and to_polygon methods aren't well-defined. My suggestion would be to not add them to the base class like you do here, and try to spend time coming up with their exact definition.
If you agree, maybe remove for now.

Or we merge and discuss next week in the hangout ...

[astrofrog]

@cdeil - no time to add docs/tests here unfortunately. If you'd like to merge anyway, that works for me.

[cdeil]

I'll merge this now, because it's key functionality for debugging sky region plotting.

Just for the record: currently to_pixel has a tolerance option which really is n_points for to_polygon and is undocumented / misnamed. I'd like to remove it completely, but will try to resist until the telcon next week to remove it since you want to keep it, right?

[astrofrog]

Just for the record: currently to_pixel has a tolerance option which really is n_points for to_polygon and is undocumented / misnamed. I'd like to remove it completely, but will try to resist until the telcon next week to remove it since you want to keep it, right?

Yep, let's discuss this on the telecon.

[cdeil]

Thank you!