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:
martinfleis: I am going to merge this now. I think we have discussed this over and over enough. If needed, there can be a follow-up PR prior 1.1.

Patch:
diff --git a/.gitignore b/.gitignore
index fa41cb46..435a84f8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -73,6 +73,8 @@ doc/source/getting_started/my_file.geojson
 
 .ruff_cache
 .env
+doc/test_fiona.gpkg
+doc/test_pyogrio.gpkg
 # pixi environments
 .pixi
 *.egg-info
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",
         [