Merge pull request #495 from cds-astro/feat-alternative-contains-on-s…

…patial-constraints

feat: add intersect modes for registry.Spatial constraint

CHANGES.rst

@@ -1,6 +1,8 @@
 1.5 (unreleased)
 ================
 
+- Add intersect modes for the spatial constraint in the registry module ``pyvo.registry.Spatial`` [#495]
+
 - Added ``alt_identifier``, ``created``, ``updated`` and ``rights`` to the
   attributes of ``pyvo.registry.regtap.RegistryResource`` [#492]
 

docs/registry/index.rst

@@ -65,8 +65,8 @@ keyword arguments.  The following constraints are available:
 * :py:class:`pyvo.registry.Ivoid` (``ivoid``): exactly match a single
   IVOA identifier (that is, in effect, the primary key in the VO).
 * :py:class:`pyvo.registry.Spatial` (``spatial``): match resources
-  covering a certain geometry (point, circle, polygon, or MOC).
-  *RegTAP 1.2 Extension*
+  covering, enclosed or overlapping a certain geometry
+  (point, circle, polygon, or MOC). *RegTAP 1.2 Extension*
 * :py:class:`pyvo.registry.Spectral` (``spectral``): match resources
   covering a certain part of the spectrum (usually, but not limited to,
   the electromagnetic spectrum).  *RegTAP 1.2 Extension*
@@ -146,6 +146,25 @@ interactive data discovery, however, it is usually preferable to use the
        Sloan Digital Sky Survey-II Supernova Survey (Sako+, 2018) ... conesearch, tap#aux, web
   ...
 
+And to look for tap resources *in* a specific cone, you would do
+
+.. doctest-remote-data::
+
+  >>> from astropy.coordinates import SkyCoord
+  >>> registry.search(registry.Servicetype("tap"),
+  ...                 registry.Spatial((SkyCoord("23d +3d"), 3), intersect="enclosed"),
+  ...                 includeaux=True) # doctest: +IGNORE_OUTPUT
+  <DALResultsTable length=1>
+              ivoid                   res_type       short_name                   res_title                  ...  intf_types  intf_roles          alt_identifier         
+                                                                                                             ...                                                         
+              object                   object          object                       object                   ...    object      object                object             
+  ------------------------------ ----------------- ------------- ------------------------------------------- ... ------------ ---------- --------------------------------
+  ivo://cds.vizier/j/apj/835/123 vs:catalogservice J/ApJ/835/123 Globular clusters in NGC 474 from CFHT obs. ... vs:paramhttp        std doi:10.26093/cds/vizier.18350123
+
+Where ``intersect`` can take the following values:
+  * 'covers' is the default and returns resources that cover the geometry provided,
+  * 'enclosed' is for services in the given region,
+  * 'overlaps' returns services intersecting with the region.
 
 The idea is that in notebook-like interfaces you can pick resources by
 title, description, and perhaps the access mode (“interface”) offered.

pyvo/registry/rtcons.py

@@ -528,6 +528,14 @@ class Spatial(Constraint):
 
     .. _MOC: https://www.ivoa.net/documents/MOC/
 
+    To find resources which coverage is enclosed in a region,
+
+        >>> enclosed = registry.Spatial("0/0-11", intersect="enclosed")
+
+    To find resources which coverage intersects a region,
+
+        >>> overlaps = registry.Spatial("0/0-11", intersect="overlaps")
+
     When you already have an astropy SkyCoord::
 
         >>> from astropy.coordinates import SkyCoord
@@ -539,12 +547,11 @@ class Spatial(Constraint):
         >>> resources = registry.Spatial((SkyCoord("23d +3d"), 3))
     """
     _keyword = "spatial"
-    _condition = "1 = CONTAINS({geom}, coverage)"
     _extra_tables = ["rr.stc_spatial"]
 
     takes_sequence = True
 
-    def __init__(self, geom_spec, order=6):
+    def __init__(self, geom_spec, intersect="covers", order=6):
         """
 
         Parameters
@@ -555,11 +562,17 @@ def __init__(self, geom_spec, order=6):
             as a DALI polygon.  Additionally, strings are interpreted
             as ASCII MOCs, SkyCoords as points, and a pair of a
             SkyCoord and a float as a circle.  Other types (proper
-            geometries or pymoc objects) might be supported in the
+            geometries or MOCPy objects) might be supported in the
             future.
+        intersect : str, optional
+            Allows to specify the connection between the resource coverage
+            and the *geom_spec*. The possible values are 'covers' for services
+            that completely cover the *geom_spec* region, 'enclosed' for services
+            completely enclosed in the region and 'overlaps' for services which
+            coverage intersect the region.
         order : int, optional
             Non-MOC geometries are converted to MOCs before comparing
-            them to the resource coverage.  By default, this contrains
+            them to the resource coverage.  By default, this constraint
             uses order 6, which corresponds to about a degree of resolution
             and is what RegTAP recommends as a sane default for the
             order actually used for the coverages in the database.
@@ -592,7 +605,15 @@ def tomoc(s):
         else:
             raise ValueError("This constraint needs DALI-style geometries.")
 
-        self._fillers = {"geom": geom}
+        if intersect == "covers":
+            self._condition = f"1 = CONTAINS({geom}, coverage)"
+        elif intersect == "enclosed":
+            self._condition = f"1 = CONTAINS(coverage, {geom})"
+        elif intersect == "overlaps":
+            self._condition = f"1 = INTERSECTS(coverage, {geom})"
+        else:
+            raise ValueError("'intersect' should be one of 'covers', 'enclosed', or 'overlaps' "
+                             f"but its current value is '{intersect}'.")
 
 
 class Spectral(Constraint):

pyvo/registry/tests/test_rtcons.py

@@ -226,6 +226,19 @@ def test_SkyCoord_Circle(self):
         assert cons.get_search_condition() == "1 = CONTAINS(MOC(6, CIRCLE(3.0, -30.0, 3)), coverage)"
         assert cons._extra_tables == ["rr.stc_spatial"]
 
+    def test_enclosed(self):
+        cons = registry.Spatial("0/1-3", intersect="enclosed")
+        assert cons.get_search_condition() == "1 = CONTAINS(coverage, MOC('0/1-3'))"
+
+    def test_overlaps(self):
+        cons = registry.Spatial("0/1-3", intersect="overlaps")
+        assert cons.get_search_condition() == "1 = INTERSECTS(coverage, MOC('0/1-3'))"
+
+    def test_not_an_intersect_mode(self):
+        with pytest.raises(ValueError, match="'intersect' should be one of 'covers', 'enclosed',"
+                           " or 'overlaps' but its current value is 'wrong'."):
+            registry.Spatial("0/1-3", intersect="wrong")
+
 
 class TestSpectralConstraint:
     # These tests might need some float literal fuzziness.  I'm just