Add spherical sky region classes

[sedonaprice]

TL;DR: Proposed spherical regions classes (implemented, with demos); if greenlit, need to implement test suite, docs before this would be ready to mark as a "non-draft" pull request. I would be very grateful for initial feedback by ~mid Sept (09/15).

Hello!

First, a huge thank you to regions' developers, for making this extremely useful package available!

I'd like to propose an implementation of spherical sky regions, complementing the planar pixel and (implicitly planar) sky region classes. (I've seen a number of issues / past discussions on this topic and how it might be implemented, but my understanding is this development has not yet been done.)

I've discussed these proposed classes with @eteq, and I'm keen to get feedback from other key folks (gleaned from past issue discussions --- @pllim, @larrybradley, @keflavich, and of course anyone else!)

Some background: I originally developed these classes as part of my work at STScI. Specifically, I needed "helper class" functionality to help demonstrate how to specify a selection region in one coordinate frame, and then transform to a different frame (and from there, construct a database query in new frame). The infrastructure and general nature of this problem led me to first develop these classes within the regions framework -- and I hope these classes could be released to the broader astronomical community as part of the regions package. If these classes could be included in regions, this would be an ideal way of providing this spherical selection & transformation functionality needed as part of the Roman Research Nexus example user notebooks.

I've put together a demo jupyter notebook walking through the core functionality (including plot demos) here: https://github.com/sedonaprice/regions/blob/add-spherical-regions-with-demo-NB/demo_spherical_regions.ipynb

I have marked this as a draft pull request because I have "tested" this code in a few demos (including the one linked above), but some work is outstanding as follows (pending a green light):

• Test suite for SphericalSkyRegion classes (complete)
• Create documentation pages + demos for the new functionality (and edit other doc pages as necessary) --- based on the below demo, but chunked out into separate topical entries. (complete
• Decision on handling of serialization for spherical regions (keyword switch to use spherical?), or add note to docs that serialization will default read in planar SkyRegion classes only. No serialization support for now / for consideration in a future PR.
• Decision on whether or not (and how) to implement planar -> spherical conversion if boundary distortions are to be included. For consideration in a future PR.

If folks think these classes would be a good candidate for inclusion, I will implement the outstanding test suite and documentation (and other code additions, pending discussion on the points noted above).

Related issues:
These new classes would address the following issues:

• fix Add spherical circles and polygons? #276
• fix Support for three IVOA region types as astropy.regions objects: what to do about RANGE? #303

and also help to achieve the aims in:

• clarification on RectangleSkyRegion? #564

An implementation (to be done) supporting conversion from planar to spherical accounting for boundary distortions would also address #217 (I have ideas on how to implement this, based on the spherical region class discretization).

• (In the future, a solid angle property could be implemented in the SphericalSkyRegion classes, addressing Introduce SkyRegion.solid_angle #434)

Thanks, and looking forward to any feedback! I'd be very grateful for initial feedback by mid Sept (09/15), so I can work on next steps (either on the outstanding tasks, or figuring out some other way of enabling the functionality for the RRN example notebooks).

TODO

• Remove demo notebook.
• Add proper documentation.
• Add tests.
• Squash out demo notebook in the commits or make sure maintainer use "Squash and Merge" button when merging.

[keflavich]

Awesome to see this finally implemented!

The jupyter notebook is great. I agree, it needs to be chunked out into docs.

For serialization, we are defining a new standard here, afaik. I recommend we make something that is still compatible with CARTA/DS9 serializations but, when read with astropy, can be interpreted as SphericalSkyRegions instead of projected sky regions.

The big concern is that, for very large images loaded into ds9/carta, there will be inconsistency between the astropy & viewer interpretation of the same file, and I don't like introducing that. So... yeah, there's an argument that we serialize in a way that ds9/carta will ignore these regions rather than show potentially incorrect regions.

I'd recommend splitting out the serialization and distortion questions into separate issues to handle after the main functionality is merged.

[sedonaprice]

Definitely, will chunk it out!

That sounds very reasonable as far as serialization and distortion (for planar -> spherical).

I'll note for completeness that I did implement distortion for spherical->planar (which I'd argue is good to keep in this PR and not split out; but perhaps I'm reading too much into your comment).

[keflavich]

No need to split that out, imo. I meant this bullet: "Decision on whether or not (and how) to implement planar -> spherical conversion if boundary distortions are to be included." should be a separate issue

[pllim]

I don't see any tests. Do you plan to add them?

I can't really do any meaningful review of the contents but maybe @mcara (maintainer of spherical-geometry) and @perrygreenfield can, or point to people who can.

Thanks!

[sedonaprice]

@pllim -- Thanks!

Re tests: I was holding off on implementing proper tests until there was consensus that this draft pull request is something that the maintainers would want to include in the package. I definitely plan to add tests as soon as there is go-ahead consensus!

(Also before finalizing, I would remove the demo NB and all changes to .pre-commit-config.yaml from the repo -- those were only temporary measures to get the NB into my fork for discussion purposes.)

[pllim]

Sounds good. Thanks! I added a "todo" section in your OP above, please update as needed. I will do a general technical re-review when this PR is at final stage. Good luck!

[sedonaprice]

I've moved the demo notebook into a separate branch on my fork for reference, and to "clean up" the draft pull-request branch. (https://github.com/sedonaprice/regions/blob/add-spherical-regions-with-demo-NB/demo_spherical_regions.ipynb; also updated to this link in the initial comment)

Splitting the demo notebook out into separate topical docs sections is on the to-do list, along with the test suite (pending discussion of the big picture as above).

[sedonaprice]

(Demo notebook + precommit changes squashed out of this branch)

[astrofrog]

Thank you for working on this! I can help review this in September once I am back at work 😊

[sedonaprice]

Quick note: I've gone ahead and implemented tests for the code as it stands right now. (And found a few fairly minimal logic bugs to fix, while I was at it.)

Looking forward to feedback & further discussion!

[astrofrog]

I've had a chance to take a look at this, and while I am not one of the regions maintainers, I think it would be great to include this so I would encourage you to proceed with writing the docs etc. This is very useful functionality! 💯

I agree with @keflavich that we should worry about serialization after this PR though, because I think the main code/functionality is needed irrespective of whether we find ways to serialize it to some of the formats where there might be ambiguities.

At this point, the only thing I am not sure about API wise is the transformation of sky regions to other systems. In particular, while it's cool that this works:

>>> lonlat_range_transf = lonlat_range.transform_to("icrs")
>>> lonlat_range_transf
<RangeSphericalSkyRegion(
frame=<class 'astropy.coordinates.builtin_frames.icrs.ICRS'>,
longitude_bounds=<LuneSphericalSkyRegion(center_gc1=<SkyCoord (ICRS): (ra, dec) in deg
    (288.42757587, 10.72370325)>, center_gc2=<SkyCoord (ICRS): (ra, dec) in deg
    (217.98010947, -60.49575721)>)>,
latitude_bounds=<CircleAnnulusSphericalSkyRegion(center=<SkyCoord (ICRS): (ra, dec) in deg
    (192.85947789, 27.12825241)>, inner_radius=45.0 deg, outer_radius=90.0 deg)>
)>

It looks pretty complex (in particular the use of nested regions) and I'm not sure I understand the benefit of this. If someone accesses the vertices, they can always convert these between different frames, but I'm not sure that it makes sense to transform the whole spherical sky region. I guess fundamentally, what is the value in allowing sky regions to be transformed between frames since they are still fundamentally the same region? For instance, contains should work regardless of frames, and if someone does to_pixel for plotting they have to pass a WCS anyway. What use cases require transform_to?

[sedonaprice]

Hi @astrofrog --

For the specific case above ("range"), the issue with just using the vertices is that a lon/lat "range" is bounded by great circles along the lines of constant longitude, but small or great circles along the lines of constant latitude. A spherical polygon defined with the same vertices would end up "bulging out" past the lines of constant latitude, so isn't equivalent.

I had added the "transforms_to()", including the nested regions, because circles transform neatly (as long as it's a spherical to spherical tranformation). Thus, since all spherical regions region can be represented as the appropriate circles (for polygons, indeed the transformation boils down to "transform the vertices and then remake the polygon, because great circles' centers transform easily; for small circles, again only the center transforms, the radius stays the same), yes the region is the same in all frames --- as long as the appropriate circle centers + radii are kept track of.

As I've implemented it, "contains()" actually does work regardless of frame (because the SkyCoord to be evaluated will be transformed to the region frame, if necessary)

The "transforms_to()" addresses the issue of "a user can define a region in their preferred non-ICRS frame, but needs to query against a database that only knows about ICRS" -- by transforming the region to ICRS, it's then possible to do bounding lon/lat in ICRS for search constraints, or a bounding circle (though of course the bounding circle transforms "easily".

(Definitely plotting requires "to_pixel" and a WCS, but the use case this was implemented to address is database/table searches with mismatches in coordinate frame.)

[sedonaprice]

Hi folks,

I've gone through and put together (draft) documentation for spherical regions (and edits to existing sections where relevant). The changes are now in the PR branch, but are also visible "live" at https://sedonaprice.github.io/regions-docs-tmp/ (temporary; I wasn't sure about how other github actions from various branches in my fork would execute, so I created a separate repo for the build + temporary publication of these draft docs).

Please let me know if you have and comments or suggestions on the docs!

(Also, now that there is a full set of draft docs and a test suite, please let me know if it would be a good time to convert this from a draft PR to a normal PR.)

[pllim]

I don't maintain this package, so I am not qualified to review. Sorry and good luck!

Stale

[astrofrog]

@sedonaprice - I should have time to review this next week, could you rebase though to resolve the conflict?

[astrofrog]

@larrybradley @keflavich - since I'm not a maintainer of the package, one of you should obviously have a look at this too, so once I've reviewed this and am happy with this I'll approve but won't merge.

[astrofrog]

Overall this looks good, but one remaining big picture comment/question in addition to a few in-line comments below: are we expecting that some of the to_spherical_sky operations that are currently NotImplementedError might be implemented in future? My understanding is that some of them will never be implemented because they don't make sense/can't be done. The Python docs say about NotImplementedError:

Note It should not be used to indicate that an operator or method is not meant to be supported at all – in that case either leave the operator / method undefined or, if a subclass, set it to None.

Basically we should not be using NotImplementedError and instead should just not add the method where it doesn't make sense. In addition, if it can't be implemented for include_boundary_distortions=True, I would suggest either removing that keyword argument altogether, or raising a ValueError if it's set to True.

Aside from this, please add an entry to the changelog.

[astrofrog]

Do you plan to do this in future or is it easy to do now?

[sedonaprice]

Right, this slipped through -- I'll implement this after addressing other points.

[sedonaprice]

Implemented in 14d4d23, for your review @astrofrog.

[sedonaprice]

Hi @astrofrog --

Thank you! I'll work on addressing the review now, and add an entry to the changelog.

That's a great point about raising a ValueError for invalid inputs instead of NotImplementedError for cases when the transformation would never make sense/can't be done.

My thoughts on this are aggregated below:

-----------------------------------------------------------------------------------------------------
| Shape         |  Distortions?  |  planar to spherical    |  spherical to planar                   |
-----------------------------------------------------------------------------------------------------
| Circle        |  True          | Implement (eventually)  |   Implemented                          |
| Circle        |  False         | Implemented             |   Implemented                          |
| CircleAnnulus |  True          | Implement (eventually)  |   Implemented                          |
| CircleAnnulus |  False         | Implemented             |   Implemented                          |
| Polygon       |  True          | Implement (eventually)  |   Implemented                          |
| Polygon       |  False         | Implemented             |   Implemented                          |
| Range         |  True          | (N/A, No equiv shape)   |   Implemented                          |
| Range         |  False         | (N/A, No equiv shape)   |   raise ValueError                     |
| Lune          |  True          | (N/A, No equiv shape)   |   Implement (eventually)               |
| Lune          |  False         | (N/A, No equiv shape)   |   raise ValueError                     |
| (MiscAnnuli)  |  True          | Implement (eventually)  |  (Implement shapes+transf eventually)  |
| (MiscAnnuli)  |  False         | Implement (eventually)  |  (Implement shapes+transf eventually)  |
| (Rectangle)   |  True          | Implement (eventually?) |  (Implement shapes+transf eventually?) |
| (Rectangle)   |  False         | Implement (eventually?) |  (Implement shapes+transf eventually?) |
-----------------------------------------------------------------------------------------------------

Would you agree with this categorization on which cases should raise a ValueError to indicate the transformation can't be done/doesn't make sense?

This would involve the following changes from the original code:

• Raise ValueError if transforming Range with to_pixel() or to_sky() with include_boundary_distortions=False
• Potentially support transforming lune to planar so long as boundary distortions are included. Upon reflection, users might want to plot or do some other manipulation with overlaps of a lune on an image, and so while portions of that shape would be badly distorted, this might be something worth including in the future. Raise ValueError for lune spherical to planar, if include_boundary_distortions=False.

[sedonaprice]

@astrofrog I've implemented:

• validation of lon/lat range inputs (and also the bounds inputs) in commit 14d4d23
• An update to spherical <-> planar error messages, to reflect the logical changes noted in the table above. Namely, lune and range now raise ValueError when include_boundary_distortion=False (as these aren't analogous to planar shapes), and lune raises NotImplementedError for include_boundary_distortion=True (as discretizing the boundary into a polygon does have a planar equivalent, even though it would be rather irregular). (Commit 327df12)

I've also added an entry to the changelog. If too verbose, I can reduce it to point users to docs pages with relevant details.

[sedonaprice]

Note: rebased as this had fallen behind the v0.11 release.

Also fixed a bug in validation error raising (incorrectly had as assert) noted from CI/CD check failure (and then fixing an introduced logic mistake).

[sedonaprice]

Thanks folks!

I have:

• Rebased against astropy:main
• Added documentation that spherical polygon only supports convex polygons (0df756e). Suggestions of edits for additional clarity welcome!
• Updated some failing tests and documentation tests (given changes in Improve pixel <-> sky region conversions #650, Add helper function to format floats in repr #655).

@astrofrog --- the rebase looks good, and the test failures that were corrected are understandable as far as other changes that were made.

@keflavich @larrybradley -- please let me know if I should make additional clarifications to document the limitations for PolygonSphericalSkyRegions for now, as far as taking route (1).

[sedonaprice]

As an aside point:

Given #649, would there be interest in changing the spherical regions' SphericalSkyRegion.discretize_boundary() methods to SphericalSkyRegion.to_polygon(), for parity with the new planar pixel and sky to_polygon() methods?

This would be a bit recursive as far as naming for PolygonSphericalSkyRegion, but discretizing/"extra polygonizing" a polygon's edges is key to capturing boundary distortions in planar<->spherical.... retaining the name discretize_boundary() would avoid this confusion, but would do effectively the same thing as to_polygon() with a different name.

(Alternatively an additional method to_polygon() could be added for the spherical shapes where that is logical, providing an alias of discretize_boundary().

Not a critical item, but I wanted to raise the possibility pre-merge.

[sedonaprice]

@keflavich @larrybradley I wanted to ping again about your high-level review, comments on how I addressed changes to take route 1 (document limitations of convex polygons for now), and if we could move ahead with merging the spherical regions --- thanks!

[keflavich]

I have a minor suggestion on something to add to the docs, but in my view, this is good to go and high-value - I'd like to see it merged soon.

[keflavich]

Could we add an example here showing the difference - i.e., showing how a point could go from inside the spherical circle to outside the planar circle, and vice-versa? The more clear we can make the limitations of each approach, the better.

[keflavich]

also I'd be happy to add the to_polygon method, but I don't feel strongly about it. It might be better understood as to_projected_polygon, but then you lose the symmetry advantage. But maybe the mnemonic value is worth adding it?

[larrybradley]

Thanks again, @sedonaprice. It's getting close. I have a few small changes.

For consistency, I do recommend adding a to_polygon method that delegates to discretize_boundary. discretize_boundary uses n_points, but to_polygon uses npoints. Let's go with n_points everywhere for consistency. The to_polygon methods are recent additions and have not yet been in a release, so we don't need to deprecate npoints.

[larrybradley]

@sedonaprice Not needed for this PR, but I think it would be good for a followup PR to add a check for input non-convex polygons in PolygonSphericalSkyRegion that raises an error. I think this can silently return wrong answers for non-convex inputs (users may not read the docs).

[larrybradley]

I'm also wondering if n_points should be 100 everywhere. In a couple places, you have n_points=10.

[sedonaprice]

Hi folks, thanks for the detailed comments!

At glance I can implement basically everything relatively quickly. I'll also try to add the docs example of "spherical vs planar" circle in/outside a circle as @keflavich suggested, as that's not too time intensive.

I'll go through these in the coming days and let folks know when things are ready for re-review!

@larrybradley re checking for non-convexity of input polygons: I'd definitely like to chat with you and others about how to actually implement that check, as I don't know a convenient algorithm for that off the top of my head. (Maybe checking with folks working on spherical-geometry).

(Thanks also for the suggestion to remove discretize_kwargs --- I'd thought about that at one point after the development matured, but didn't follow through!)

[larrybradley]

@sedonaprice For the concave (non-convex) check, perhaps something like comparing the convex hull area to the polygon area would work. For a concave polygon, the convex hull area will be be larger than the polygon area.

The spherical_geometry SphericalPolygon class (https://spherical-geometry.readthedocs.io/en/latest/api/spherical_geometry.polygon.SphericalPolygon.html#spherical_geometry.polygon.SphericalPolygon) has area and convex_hull methods. It's probably a good idea to input an interior point (inside=).
CC: @mcara for spherical_geometry advice

[sedonaprice]

Sure -- but my understanding was there were past decisions to not have spherical_geometry be a (required) dependency of regions. It might be more expedient to just a pure-pythonic spherical polygon fix that works for convex and non-convex (eg, as I prototyped in a separate branch, starting from the same algorithm as in spherical_geometry, but of course lacks quad precision on cross products).

(I'd also argue that it might be good to permit users to not pass interior points, as eg users would not be able to simply pass vertices from an s_region.)

[mcara]

@larrybradley My understanding is that these "spherical sky regions" (or polygons) are fundamentally planar polygons. To check that polygons are convex:

1. detect self-intersecting polygons. By definition they are not convex. AND
2. perform cross-product sign test: https://math.stackexchange.com/a/1745427 (pretty simple to implement in Python using numpy).

For polygons on the sphere, the analog on the sphere of method from #2 above would be the half‑space sign test.

With regard to the inside point: It is used to indicate the inside face of the spherical polygon. We are investigating the use of this in the spherical_geometry. At this moment, my preliminary "feeling" about it is that it introduces more confusion than it helps. So, at this moment, I cannot recommend introducing this parameter

[sedonaprice]

Hi @mcara -- regarding "My understanding is that these "spherical sky regions" (or polygons) are fundamentally planar polygons.". These spherical sky regions are not planar: the spherical polygons are sky coordinates joined by a great circle arcs, the spherical circles have a circumference defined by constant spherical projected distance from a sky coordinate, etc. Could you clarify what you meant by that statement?

[sedonaprice]

And thank you @mcara! A check for (1) self-intersection and (2) leveraging the half-space sign test would be fairly straightforward to implement.

For completeness as far as implementing a convexity check, @larrybradley, I will note that if this test was attached to the creation of a SphericalSkyPolygon, then in some cases when discretizing to transform from a spherical shape, like a circle, to a planar shape while accounting for the WCS transformation boundary distortions, the transformed, discretized boundary can be non-convex, so the full functionality of those transformations wouldn't be supported until arbitrary polygons are supported.

[mcara]

Could you clarify what you meant by that statement?

I am not familiar with regions package and I guess I did not read the getting started rst completely and I was thinking of "planar Sky Regions". In this case use "half‑space sign test".

With regard to include_boundary_distortions - from a quick look it seems this would throw a NotImplemented error. Is that correct? If true, I would suggest not having this parameter at all. Also, documentation about include_boundary_distortions mentions only projection distortions but I could not find any mention of instrument distortions. Are those ignored?

[sedonaprice]

Thanks for the clarification!

For include_boundary_distortions, no, it doesn't throw a NotImplementedError in general. For LuneSphericalSkyRegion and RangeSphericalSkyRegion (aka an lon/lat bounding box), trying to transform to a planar shape with include_boundary_distortions=False will throw NotImplementedError, as these spherical shapes have no direct analogs --- for a lune, it's full "orange slice" and a canonical planar shape can't approximate it, and for range this ensure large lon/lat ranges aren't approximated as a rectangle, which is pretty far off. Note that for LuneSphericalSkyRegion and RangeSphericalSkyRegion a transformation to a planar shape with include_boundary_distortions=True will work, as the polygonization enables it to be described on a plane.

By "projection distortions", I intended to imply "all the distortions happening to a boundary when mapping from spherical geometry to a planar projection as described by the WCS considered, or vice versa". This would indeed include instrument distortion if included in the WCS. Very welcome to suggestions for a more clear way to phrase this --- would saying "WCS projection distortions" make this more clear?

[mcara]

I use "projection distortion" to describe distortions arising from projecting a "straight" line (an arc) from a sphere to a plane. The "distortion" will depend on the projection used in the WCS. It is separate from instrument distortions usually described using SIP, tabular distortions, or TPV.

[sedonaprice]

My understanding is in general the WCS for reduced images/mosaics as parsed by astropy.wcs would include both projection distortion and also components described by SIP.

I'll look back at the docstrings & docs to see if it makes sense to add "WCS distortion" or "distortion described by the WCS transformation" to clarify things!