Repository: geopandas/geopandas

Issue:
ENH: Check GeoSeries.intersects(GeoSeries) against all geometries, not only pair-wise
edit (20/03/04) by @martinfleis: _See the discussion below for the context._
<hr>

```
def merge_boundry(image_path, vector_path, clip_size):
    path_img=os.path.dirname(image_path)
    L=os.listdir(image_path)
    for i in L:  
        dname=os.path.join(path_img,i)
        #add
        with rasterio.open(dname) as src:
            transform = src.transform
            crs = src.crs
            img_arr = reshape_as_image(src.read())
            out_transform= _clip_transform(src.meta['transform'], clip_size)
            origin_height, origin_width, _ =img_arr.shape
        
            origin_height=origin_height-2*clip_size
            origin_width=origin_width-2*clip_size
        
            min_x=out_transform[2]
            min_y=out_transform[5]
            max_x=min_x+out_transform[0]*origin_height
            max_y=min_y+out_transform[4]*origin_width
        
           #every location coordinate
            left_top=min_x, min_y
            left_bottom=min_x, max_y 
            right_top=max_x, min_y
            right_bottom=max_x, max_y    
        
            p1, p2, p3, p4, p5= left_top, left_bottom, right_bottom, right_top, left_top
            draw_frame=LineString([p1,p2,p3,p4,p5])
            #invert format
            draw_frame=geopandas.GeoSeries(draw_frame)
            draw_frame.crs=crs
    
    path_vector=os.path.dirname(vector_path)
    vector_name=os.listdir(vector_path)
    
    for vector in vector_name:
        draw_frame_buffer=draw_frame.buffer(1)
        all_path_merge=os.path.join(path_vector, vector)
#         vector_data=geopandas.read_file(all_path_merge)
        vector_data=geopandas.GeoDataFrame.from_file(all_path_merge)
        block_intersect=draw_frame_buffer.intersects(vector_data)
        print(block_intersect)
```
Result:
0      True
1     False
2     False
3     False
4     False
5     False
6     False
7     False
8     False
9     False
10    False
11    False
12    False
13    False
14    False
15    False
16    False
17    False
18    False
19    False
20    False
21    False
22    False
dtype: bool
0     False
1     False
2     False
3     False
4     False
5     False
6     False
7     False
8     False
9     False
10    False
11    False
12    False
13    False
14    False
15    False
16    False
17    False
18    False
19    False
20    False
21    False
22    False
23    False
dtype: bool
0    False
dtype: bool
0     False
1     False
2     False
3     False
4     False
5     False
6     False
7     False
8     False
9     False
10    False
11    False
12    False
13    False
14    False
15    False
16    False
17    False
18    False
19    False
20    False
21    False
22    False
23    False
24    False
25    False
26    False
dtype: bool
0      True
1     False
2     False
3     False
4     False
5     False
6     False
7     False
8     False
9     False
10    False
11    False
12    False
13    False
14    False
15    False
16    False
17    False
18    False
19    False
dtype: bool

reason :
 why not true, only return false ,in fact, exits true

Discussion:
@jorisvandenbossche: @jiezhengxu without more context it is very difficult to understand your issue or question. Can you explain what you are trying to do? And can you also try to provide a reproducible example? (some code that we can run)@b4l: @jorisvandenbossche I think i met the same issue when intersecting two GeoDataFrames after reprojecting using Pandas 1.0.1 and GeoPandas 0.7. Here is an example:

```python
import pandas as pd
import geopandas as gpd
import requests
import io
import zipfile
import tempfile

# some pois
pois = pd.DataFrame(
    {
        "latitude": [47.58, 46.23, 47.46, 46.91, 47.44],
        "longitude": [7.52, 6.10, 8.54, 7.49, 8.23],
    }
)
pois = gpd.GeoDataFrame(
    pois, geometry=gpd.points_from_xy(pois.longitude, pois.latitude)
)
pois.crs = "EPSG:4326"

# some polygons
url = "https://data.geo.admin.ch/ch.swisstopo.swissboundaries3d-land-flaeche.fill/shp/2056/ch.swisstopo.swissboundaries3d-land-flaeche.fill.zip"
r = requests.get(url)
z = zipfile.ZipFile(io.BytesIO(r.content))
tmp = tempfile.mkdtemp()
z.extractall(tmp)
ch = gpd.read_file(
    tmp + "/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp"
)
ch.geometry = ch.buffer(10000).to_crs("EPSG:4326")

# plot
base = ch.plot(color='white', edgecolor='black')
pois.plot(ax=base, marker='o', color='red', markersize=5)

# intersection
pois.intersects(ch)

# what i expect
pois.apply(lambda x: ch.intersects(x.geometry).any(), axis="columns")
```@martinfleis: @b4l `gdf1.intersects(gdf2)` does pair-wise intersection. Rows are matched and only geometries within the same rows are checked against each other. But I agree that it is not what one might expect. The solution might be to add a new keyword to have `gdf.intersects(gdf2, any=True)` which would do what you expected (and could even use spatial index under the hood to do so).@b4l: Thank you for the insights @martinfleis. I think the [spatial predicates](https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates) should be handled by default as set operations, whereas the elementwise variant could be accessed through an optional keyword. What do you think?@martinfleis: I feel the same about it, but to changing this behaviour to default on `.any()` would be backward incompatible and might break a lot of code if implemented. Moreover, this element wise logic is something we inherit from pandas and having some functions working like that and other not is something I'd like to avoid.@b4l: As far as I know, there are many functions in pandas that do not operate elementwise, like all the merge, join, append etc. The 1.0 Milestone would probably be a good opportunity to introduce these changes. It might be painfull, though even more justifying a counterintuitive api to the user.@ljwolf: For a bit of context, the align-then-relate pattern used in geopandas is a common rough edge I think in part because it's not how [PostGIS works](https://postgis.net/docs/ST_Touches.html) for spatial predicates (in addition to being different from the DE9IM semantics). This is why it's all-pairs by default in R's `sf`, too.

I understand how we got here, but I can see how it's a confusing behavior for users, so an all-pairs relates between geoseries by default would definitely make sense for a 1.0 release, but it's a big ask... like @martinfleis, am concerned that this is a huge breakage.
 
If so, a keyword argument is definitely preferred (though I wouldn't want `any()`, I'd want `pairwise=True`. Users should explicitly use `df.intersects(df2, pairwise=True).any(axis=1)` for
> Does `df` intersect __any element__ of `df2`?

Alternatively, a more explicit separation could put fast index-enabled pairwise predicates under a specific namespace in the dataframe, like `df.pairwise.intersects(df2)`. This is used elsewhere [say, sklearn](https://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics.pairwise) to separate the functions that might apply to a single coordinate pair versus applying the function to all pairs between two matrices of coordinates.@b4l: I would prefer `df.intersects(df2, pairwise=True)` with the default value set to `False` for the set operation. Alternatively the default could be set to `True` and changed later to avoid breakage and time the transition.@martinfleis: I think we could do `df.intersects(df2, pairwise=False)` with a deprecation warning that default will change. I would keep it there for a few releases and then change it. It may be 1.0.0 or 0.11.0 or whatever will fit to give enough time for a transition.@ljwolf: If we use `pairwise=True` to mean row-by-row, that's different from the semantics of `sklearn.metrics.pairwise` as well as the implied meaning of `scipy.spatial.distance.pdist`. I don't think that's a good idea... am I understanding your comment about the "default value set to `False` for the set operation." correctly @b4l? 

If so & folks find `pairwise=True` confusing to mean the relation for all pairs, What about `all_pairs=True`?@b4l: What about `paired=True` relatet to `paired_distance` in sklearn? scipy has cdist which seems odd.@martinfleis: Good points. I like `paired` as it is a bit more explanatory than `all_pairs`.@b4l: Perhaps it might be useful to take other comparison methods into account as well in designing the interface, like for example bbox intersection.  The latter could be achieved by another parameter or combined in a string typed mode parameter for example.@martinfleis: I have recently played with this and came up with a possible implementation. See https://gist.github.com/martinfleis/abc7cdbf9f9266bf9ed369080eec7cea 

I am not sure about the API yet, so there are two main questions to resolve (from an implementation point of view):

1. Which options should we offer? Existing + indices or existing + indices + "any"?
2. How should keywords look like?

Note that there might be some temporary workarounds as pygeos does not support all predicates yet (https://github.com/pygeos/pygeos/pull/157).

And then there's the question of change of default...@brendan-ward: Follow up note: the additional predicates that make sense in the context of a spatial index are now available as of [pygeos version 0.8](https://github.com/pygeos/pygeos/blob/master/CHANGELOG.rst#version-08-2020-09-06).@felix-ht: any updates on this?@martinfleis: @felix-ht There is a draft PR and a related discussion in https://github.com/geopandas/geopandas/pull/1674 I need to update whenever I have more time for that. If you need to use it now, you can get the code from that PR.

Solution patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b9d6dacd..4224727e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,10 @@
 
 New features and improvements:
 
+- Added options to return the result of `SpatialIndex.query` in a form of a dense or a
+  sparse boolean array. This adds optional dependency on `scipy` for the sparse output.
+  Note that this also changes the previously undocumented behaviour of the `output_format`
+  keyword (#1674).
 - Add ``grid_size`` parameter to ``union_all`` and ``dissolve`` (#3445).
 - `GeoDataFrame.plot` now supports `pd.Index` as an input for the `column` keyword (#3463).
 - Avoid change of the plot aspect when plotting missing values (#3438).
diff --git a/doc/source/getting_started/install.rst b/doc/source/getting_started/install.rst
index 6246f998..30377507 100644
--- a/doc/source/getting_started/install.rst
+++ b/doc/source/getting_started/install.rst
@@ -158,6 +158,7 @@ Further, optional dependencies are:
 - `GeoAlchemy2`_ (optional; for writing to PostGIS)
 - `geopy`_ (optional; for geocoding)
 - `pointpats`_ (optional; for advanced point sampling)
+- `scipy`_ (optional; for sparse output of spatial indexing)
 
 
 For plotting, these additional packages may be used:
@@ -214,3 +215,5 @@ For plotting, these additional packages may be used:
 .. _pointpats: https://pysal.org/pointpats/
 
 .. _folium: https://python-visualization.github.io/folium/stable/
+
+.. _scipy: https://docs.scipy.org/doc/scipy/
\ No newline at end of file
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 10c11a07..159539c4 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -54,20 +54,46 @@ class SpatialIndex:
         return PREDICATES
 
     def query(
-        self, geometry, predicate=None, sort=False, distance=None, output_format="tuple"
+        self,
+        geometry,
+        predicate=None,
+        sort=False,
+        distance=None,
+        output_format="indices",
     ):
         """
-        Return the integer indices of all combinations of each input geometry
+        Return all combinations of each input geometry
         and tree geometries where the bounding box of each input geometry
         intersects the bounding box of a tree geometry.
 
-        If the input geometry is a scalar, this returns an array of shape (n, ) with
-        the indices of the matching tree geometries.  If the input geometry is an
-        array_like, this returns an array with shape (2,n) where the subarrays
-        correspond to the indices of the input geometries and indices of the
-        tree geometries associated with each.  To generate an array of pairs of
-        input geometry index and tree geometry index, simply transpose the
-        result.
+        The result can be returned as an array of 'indices' or a boolean 'sparse' or
+        'dense' array. This can be controlled using the ``output_format`` keyword.
+        Options are as follows.
+
+        ``'indices'``
+            If the input geometry is a scalar, this returns an array of shape (n, ) with
+            the indices of the matching tree geometries.  If the input geometry is an
+            array_like, this returns an array with shape (2,n) where the subarrays
+            correspond to the indices of the input geometries and indices of the
+            tree geometries associated with each.  To generate an array of pairs of
+            input geometry index and tree geometry index, simply transpose the
+            result.
+        ``'sparse'``
+            If the input geometry is a scalar, this returns a boolean scipy.sparse COO
+            array of shape (len(tree), ) with boolean values marking whether the
+            bounding box of a geometry in the tree intersects a bounding box of a given
+            scalar. If the input geometry is an array_like, this returns a boolean
+            scipy.sparse COO array with shape (len(tree), n) with boolean values marking
+            whether the bounding box of a geometry in the tree intersects a bounding box
+            of a given scalar.
+        ``'dense'``
+            If the input geometry is a scalar, this returns a boolean numpy
+            array of shape (len(tree), ) with boolean values marking whether the
+            bounding box of a geometry in the tree intersects a bounding box of a given
+            scalar. If the input geometry is an array_like, this returns a boolean
+            numpy array with shape (len(tree), n) with boolean values marking
+            whether the bounding box of a geometry in the tree intersects a bounding box
+            of a given scalar.
 
         If a predicate is provided, the tree geometries are first queried based
         on the bounding box of the input geometry and then are further filtered
@@ -106,22 +132,56 @@ class SpatialIndex:
             as the secondary key.
             If False, no additional sorting is applied (results are often
             sorted but there is no guarantee).
+            Applicable only if output_format="indices".
         distance : number or array_like, optional
             Distances around each input geometry within which to query the tree for
             the 'dwithin' predicate. If array_like, shape must be broadcastable to shape
             of geometry. Required if ``predicate='dwithin'``.
+        output_format : {"indices", "sparse", "dense"}, default "indices"
+            Type of the output format representing the result of the query.
 
         Returns
         -------
-        ndarray with shape (n,) if geometry is a scalar
+        `If geometry is a scalar:`
+
+        ndarray with shape (n,)
             Integer indices for matching geometries from the spatial index
-            tree geometries.
+            tree geometries.  If ``output_format="indices"``.
 
         OR
 
-        ndarray with shape (2, n) if geometry is an array_like
+        scipy.sparse COO array with shape (len(tree), )
+            Boolean array aligned with array of geometries in the tree.
+            If ``output_format="sparse"``.
+
+        OR
+
+        ndarray with shape (len(tree), )
+            Boolean array aligned with array of geometries in the tree.
+            If ``output_format="dense"``.
+
+
+        `If geometry is an array_like:`
+
+        ndarray with shape (2, n)
             The first subarray contains input geometry integer indices.
             The second subarray contains tree geometry integer indices.
+            If ``output_format="indices"``.
+
+        OR
+
+        scipy.sparse COO array with shape (len(tree), n)
+            Boolean array aligned with array of geometries in the tree along axis 0 and
+            with ``geometry`` along axis 1.
+            If ``output_format="sparse"``.
+
+        OR
+
+        ndarray with shape (len(tree), n)
+            Boolean array aligned with array of geometries in the tree along axis 0 and
+            with ``geometry`` along axis 1.
+            If ``output_format="dense"``.
+
 
         Examples
         --------
@@ -170,6 +230,32 @@ class SpatialIndex:
         >>> s.sindex.query(box(1, 1, 3, 3), predicate="dwithin", distance=2)
         array([0, 1, 2, 3, 4])
 
+        Returning boolean arrays:
+
+        >>> s.sindex.query(box(1, 1, 3, 3), output_format="sparse")
+        <COOrdinate sparse array of dtype 'bool'
+            with 3 stored elements and shape (10,)>
+
+        >>> s.sindex.query(box(1, 1, 3, 3), output_format="dense")
+        array([False,  True,  True,  True, False, False, False, False, False,
+               False])
+
+        >>> s.sindex.query(s2, output_format="sparse")
+        <COOrdinate sparse array of dtype 'bool'
+            with 5 stored elements and shape (10, 2)>
+
+        >>> s.sindex.query(s2, output_format="dense")
+        array([[False, False],
+               [False, False],
+               [ True, False],
+               [ True, False],
+               [ True, False],
+               [False,  True],
+               [False,  True],
+               [False, False],
+               [False, False],
+               [False, False]])
+
         Notes
         -----
         In the context of a spatial join, input geometries are the "left"
@@ -211,10 +297,7 @@ class SpatialIndex:
 
         indices = self._tree.query(geometry, predicate=predicate, **kwargs)
 
-        if output_format != "tuple":
-            sort = True
-
-        if sort:
+        if output_format == "indices" and sort:
             if indices.ndim == 1:
                 indices = np.sort(indices)
             else:
@@ -224,23 +307,37 @@ class SpatialIndex:
                 indices = np.vstack((geo_idx[sort_indexer], tree_idx[sort_indexer]))
 
         if output_format == "sparse":
-            from scipy.sparse import coo_array
+            scipy = compat.import_optional_dependency("scipy")
 
-            return coo_array(
-                (np.ones(len(indices[0]), dtype=np.bool_), indices),
+            if indices.ndim == 1:
+                return scipy.sparse.coo_array(
+                    (np.ones(len(indices), dtype=np.bool_), indices.reshape(1, -1)),
+                    shape=(len(self.geometries),),
+                    dtype=np.bool_,
+                )
+            return scipy.sparse.coo_array(
+                (np.ones(len(indices[0]), dtype=np.bool_), indices[::-1]),
                 shape=(len(self.geometries), len(geometry)),
                 dtype=np.bool_,
             )
 
         if output_format == "dense":
-            dense = np.zeros((len(self.geometries), len(geometry)), dtype=bool)
-            dense[indices] = True
+            if indices.ndim == 1:
+                dense = np.zeros(len(self.geometries), dtype=bool)
+                dense[indices] = True
+            else:
+                dense = np.zeros((len(self.geometries), len(geometry)), dtype=bool)
+                tree, other = indices[::-1]
+                dense[tree, other] = True
             return dense
 
-        if output_format == "tuple":
+        if output_format == "indices":
             return indices
 
-        raise ValueError(f"Invalid output_format: {output_format}")
+        raise ValueError(
+            f"Invalid output_format: '{output_format}'. "
+            "Use one of 'indices', 'sparse', 'dense'."
+        )
 
     @staticmethod
     def _as_geometry_array(geometry):
diff --git a/pyproject.toml b/pyproject.toml
index 9717a007..d3d964be 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,6 +40,7 @@ all = [
     "folium", # >=0.12 (implied minimum version, but older versions may work)
     "GeoAlchemy2",
     "pyarrow>=10.0.0",
+    "scipy",
 ]
 
 # Minimum supported additional deps, not installed as part of `all` extra


Test patch:
diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index efaf3801..ab7ea414 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -19,6 +19,14 @@ from geopandas import _compat as compat
 import pytest
 from numpy.testing import assert_array_equal
 
+try:
+    from scipy.sparse import coo_array  # noqa: F401
+
+    HAS_SCIPY = True
+except ImportError:
+    HAS_SCIPY = False
+SCIPY_MARK = pytest.mark.skipif(not HAS_SCIPY, reason="scipy not installed")
+
 
 class TestSeriesSindex:
     def test_has_sindex(self):
@@ -261,6 +269,9 @@ class TestShapelyInterface:
             self.df.sindex.intersection(test_geom)
 
     # ------------------------------ `query` tests ------------------------------ #
+    @pytest.mark.parametrize(
+        "output_format", ("indices", pytest.param("sparse", marks=SCIPY_MARK), "dense")
+    )
     @pytest.mark.parametrize(
         "predicate, test_geom, expected",
         (
@@ -359,11 +370,22 @@ class TestShapelyInterface:
             ),  # contains but does not contains_properly
         ),
     )
-    def test_query(self, predicate, test_geom, expected):
+    def test_query(self, predicate, test_geom, expected, output_format):
         """Tests the `query` method with valid inputs and valid predicates."""
         res = self.df.sindex.query(test_geom, predicate=predicate)
         assert_array_equal(res, expected)
 
+        if output_format != "indices":
+            dense = np.zeros(len(self.df), dtype=bool)
+            dense[expected] = True
+
+            res = self.df.sindex.query(
+                test_geom, predicate=predicate, output_format=output_format
+            )
+            if output_format == "sparse":
+                res = res.todense()
+            assert_array_equal(res, dense)
+
     def test_query_invalid_geometry(self):
         """Tests the `query` method with invalid geometry."""
         with pytest.raises(TypeError):
@@ -535,7 +557,15 @@ class TestShapelyInterface:
                 )
             raise e
 
+    def test_unsupported_output(self):
+        with pytest.raises(ValueError, match="Invalid output_format: 'dataarray'."):
+            test_geom = box(-1, -1, -0.5, -0.5)
+            self.df.sindex.query(test_geom, output_format="dataarray")
+
     # ------------------------- `query_bulk` tests -------------------------- #
+    @pytest.mark.parametrize(
+        "output_format", ("indices", pytest.param("sparse", marks=SCIPY_MARK), "dense")
+    )
     @pytest.mark.parametrize(
         "predicate, test_geom, expected",
         (
@@ -635,15 +665,26 @@ class TestShapelyInterface:
             ),  # contains but does not contains_properly
         ),
     )
-    def test_query_bulk(self, predicate, test_geom, expected):
+    def test_query_bulk(self, predicate, test_geom, expected, output_format):
         """Tests the `query` method with valid
         inputs and valid predicates.
         """
-        res = self.df.sindex.query(
-            [box(*geom) for geom in test_geom], predicate=predicate
-        )
+        test_geoms = [box(*geom) for geom in test_geom]
+        res = self.df.sindex.query(test_geoms, predicate=predicate)
         assert_array_equal(res, expected)
 
+        if output_format != "indices":
+            dense = np.zeros((len(self.df), len(test_geoms)), dtype=bool)
+            tree, other = expected[::-1]
+            dense[tree, other] = True
+
+            res = self.df.sindex.query(
+                test_geoms, predicate=predicate, output_format=output_format
+            )
+            if output_format == "sparse":
+                res = res.todense()
+            assert_array_equal(res, dense)
+
     @pytest.mark.parametrize(
         "test_geoms, expected_value",
         [