Add spherical circles and polygons?

[cdeil]

This is a reminder issue that we might want to bring back spherical circles and polygons in the regions package.

We had at least spherical circle functionality at some point, but for consistency with DS9 and to avoid confusion of having one class represent a planar and a spherical sky circle sometimes, it was removed in #132 .

Now in #275 I'm removing the last trace of that functionality, an unused and untested helper function. But really, I think we probably should bring that functionality back (with a minor change, not call numpy.matrix, but instead use arrays and numpy.dot or numpy.matmul). I think it was discussed before, and we thought a separate class like CircleSphericalSkyRegion in addition to CircleSkyRegion would be the way to go?
@astrofrog - Do you remember?

[gpdf]

See #303; it would be helpful both in building client-side interfaces to IVOA services, but especially in implementing IVOA protocols in Python, if these spherical geometry region types existed, as they map directly to concepts in the IVOA interfaces.

[cdeil]

The question how to represent and handle "true sky" regions that work without a WCS has come up many times, see e.g. #13, #295, #303, #76, #221

Ideas how to handle this included a flag (e.g. mode={"planar", "spherical"}) either on the region class, or on the methods that require this choice (namely contains and plotting), or separate classes.

I don't think there's a perfect solution, because the most-used serialisation format (DS9) doesn't support the distinction, and there is no great default behaviour. E.g. I think most users would expect CircleSkyRegion to be spherical, but it's planar. Making it spherical would mean that one has to go to polygon on sky_circle.to_pixel(wcs) with some numerical precision (e.g. 100 points) and that wouldn't round-trip or be consistent with DS9 behaviour. So for those reasons we so far chose to remove the "spherical" variants completely and just do "planar" behaviour.

My suggestion would be to start by adding new classes for the spherical cases:

• SphericalCircleSkyRegion
• SphericalPolygonSkyRegion
• SphericalRangeSkyRegion (see Support for three IVOA region types as astropy.regions objects: what to do about RANGE? #303)

It's contains would look like this, and it could or couldn't sub-class CircleSkyRegion, and e.g. to_pixel or plot would polygonise:

class SphericalCircleSkyRegion(CircleSkyRegion):
    """Spherical circle sky region."""

    def contains(self, skycoord, wcs=None):
        """Defined by spherical distance."""
        separation = self.center.separation(skycoord)
        return separation < self.radius

I'm not sure how DS9 or other serialisation should behave.

I'll send a PR today for the SphericalCircleSkyRegion case to get going, but ideally we'd have discussion and find the best possible solution here very soon. Please comment!

[gpdf]

I think it'll be necessary to have some docstring and/or actual documentation to help users understand the difference between circles and spherical circles. It's not at all intuitive. I think the most useful piece of imagery is to make sure that the reader thinks about a planar circle that's away from the point of tangency of a tangent plane. That helps (I think) break the "but how can they be different - they're both rotationally symmetric, so they must be the same" conceptual block.

[cdeil]

I started to implement separate classes in #306 . But now I've left astronomy and am starting a new job, so I'll close that PR now. @gpdf @larrybradley @keflavich @astrofrog or anyone - if you're interested to add spherical regions, and if you want to pursue the separate classes approach, you could start a new PR from that branch. If you prefer another way, e.g. adding methods and possibly a flag to existing classes, it's probably better if you start from master.

[caseyjlaw]

I came across this package recently with an interest in finding sources within irregular regions. It seems like this issue defines the solution most clearly. Since there is not a lot of activity on this issue, I'm wondering if people here can recommend a workaround.
My use case is:

• create catalog of points from list of (RA, Dec)
• define polygon region from (RA, Dec) pairs drawn from a second catalog
• use contains method to find all points in the polygon.

I would think that I could do this well enough with PolygonPixelRegion if my vertices are finely spaced and avoid the poles. Then my coordinates would all be defined as PixelCoords. Does that sound right?