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 that in the processes of exploring additional functionality (discretization for planar shapes, and completing the planar -> spherical with distortions set of transformations), I identified a bug where the internal calculations for discretization of boundaries (to polygons) could end up with some vertices incorrectly flipped and thus the boundary covering the other part of the great circle.

This has been fixed in d95c1a9.

[sedonaprice]

Also @astrofrog --- I have also realized that the implementation of PolygonSphericalSkyRegion.contains() in this PR is only valid for convex polygons (computation simply uses the intersection of the great circles along which the edges lie on).

I have tested a generalized implementation (using astropy and numpy functions -- see branch here). In investigating solutions, I examined the code in spherical_geometry (particularly here), but when coding it up I worked directly from the math in the source references (here and here; though I have not checked to see the behavior when a point lies on the boundary --- which would make it inconsistent with circle contains excluding the boundary. Choices regarding edge contains conventions is probably best deferred to a separate discussion).

Question: Should I merge in this fix to enable general polygon support to this PR, or should that be the subject of a different PR after this one is approved?

Happy to discuss --- including thoughts from other maintainers or other relevant folks!

This of course relates to previous issue discussions regarding "should spherical_geometry be a dependency", but I figured it was worth exploring (a) a fix for the incorrect implementation in this PR, and (b) a no-dependency solution for at least the interim, or for use in case spherical_geometry was considered for addition as an optional dependency.

[astrofrog]

This is looking good - I found a few issues below but once these are addressed, I'm happy with this, and I don't anticipate having time to review in full again.

One of the official maintainers should make a decision on merging this though (and whether or not they want to review themselves), this is beyond my authority.

[sedonaprice]

@astrofrog Thank you for your thorough review, and your suggestions! I've implemented all changes (including other places in the code where equivalent changes could be made, where relevant).

@keflavich and @larrybradley --- as maintainers of regions, I wanted to check in about next steps for considering this for merging.

I wanted to flag again the following (which I noted on 9 Jan):

I have also realized that the implementation of PolygonSphericalSkyRegion.contains() as currently implemented in the branch for this PR is only valid for convex polygons (computation simply uses the intersection of the great circles along which the edges lie on).

Possible solutions include:

1. Consider this PR code as-is, modifying the doc strings (and/or SphericalSkyPolygon class name) to note it's currently valid for complex polygons. Make modifications to ensure generalized spherical polygons are supported in a future PR.
2. Fix this issue, implementing generalized spherical polygons before completing this PR. (I worked out such a fix in a separate branch.) Longer term, it might be worth considering adding an optional dependency to handle generalized polygon contains logic (e.g., `spherical_geometry).

Happy to discuss how best to tackle this issue, and thank you again everyone for your time and extremely useful feedback!