Repository: geopandas/geopandas

Issue:
ENH: preserve relative order in `clip`
- [x] I have checked that this issue has not already been reported.

- [x] I have confirmed this bug exists on the latest version of geopandas.

- [ ] (optional) I have confirmed this bug exists on the main branch of geopandas.

---

I wasn't sure if this behavior is intended so I left the BUG label.

#### Code Sample

```python
>>> import geopandas as gpd
>>> import numpy as np
>>> import shapely
>>>
>>> z = np.array([5, 2, 4, 0, 3, 1])
>>> gdf = gpd.GeoSeries.from_xy(z, z).to_frame("geometry")
>>> bbox = shapely.geometry.box(0, 0, 2, 2)
>>> gdf
                  geometry
0  POINT (5.00000 5.00000)
1  POINT (2.00000 2.00000)
2  POINT (4.00000 4.00000)
3  POINT (0.00000 0.00000)
4  POINT (3.00000 3.00000)
5  POINT (1.00000 1.00000)

>>> gdf.clip(bbox)
                  geometry
3  POINT (0.00000 0.00000)
5  POINT (1.00000 1.00000)
1  POINT (2.00000 2.00000)
```

#### Problem description

Clipping a `GeoDataFrame` produces a result that has been shuffled. It looks like the result is sorted based on a spatial index. This is not noted in the documentation and was not what I expected to happen.

#### Expected Output

I expected clip to maintain the relative order in the result or for this behavior to be documented.

#### Output of ``geopandas.show_versions()``

<details>

SYSTEM INFO
-----------
python     : 3.9.16 | packaged by conda-forge | (main, Feb  1 2023, 21:39:03)  [GCC 11.3.0]
executable : /home/fred/anaconda3/envs/rstools/bin/python3.9
machine    : Linux-5.15.0-73-generic-x86_64-with-glibc2.35

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : 3.11.2
GEOS lib   : None
GDAL       : 3.7.0
GDAL data dir: /home/fred/anaconda3/envs/rstools/share/gdal
PROJ       : 9.2.1
PROJ data dir: /home/fred/anaconda3/envs/rstools/share/proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.13.2
numpy      : 1.24.3
pandas     : 2.0.2
pyproj     : 3.6.0
shapely    : 2.0.1
fiona      : 1.9.4
geoalchemy2: None
geopy      : None
matplotlib : 3.7.1
mapclassify: 2.5.0
pygeos     : 0.14
pyogrio    : 0.6.0
psycopg2   : None
pyarrow    : 12.0.0
rtree      : 1.0.1

</details>


Discussion:
martinfleis: > Doesn't 100% address the request as they requested preservation of the input order, but I imagine sorting would achieve the same outcome 9/10 times.

Can you elaborate on that? In which case this does not work as expected?
nicholas-ys-tan: > > Doesn't 100% address the request as they requested preservation of the input order, but I imagine sorting would achieve the same outcome 9/10 times.
> 
> 
> 
> Can you elaborate on that? In which case this does not work as expected?

The poster requested preservation of the input order. In the case where the input dataframe is not sorted, but clip is run with sort=True, it will return a sorted dataframe. If run with sort=False, it may return with a different unsorted order from the input.

I can't think of any cases where a user would want to maintain the order of an unsorted dataframe though.
martinfleis: > The poster requested preservation of the input order. In the case where the input dataframe is not sorted, but clip is run with sort=True, it will return a sorted dataframe. If run with sort=False, it may return with a different unsorted order from the input.
> 
> I can't think of any cases where a user would want to maintain the order of an unsorted dataframe though.

I don't think that is correct. The indices coming from `sindex.query` are integer positions of geometries in the original GeoSeries, it does not matter how their index looks like. So even if the GeoSeries to be clipped has a custom, unsorted index, with `sort=True`, the original order is fully preserved. Taking an example from your tests:

```py
cities = gpd.read_file(naturalearth_cities)
world = gpd.read_file(naturalearth_lowres)
south_america = world[world["continent"] == "South America"]
cities.reset_index().set_index("name").clip(south_america, sort=True)

                index                     geometry
name                                              
Georgetown         55    POINT (-58.16703 6.80197)
Paramaribo         59    POINT (-55.16703 5.83503)
Asunción           62  POINT (-57.62583 -25.29067)
Quito              88     POINT (-78.502 -0.21304)
Valparaíso        101  POINT (-71.61703 -33.04774)
Sucre             114  POINT (-65.25952 -19.04097)
La Paz            122  POINT (-68.15193 -16.49603)
Brasília          169    POINT (-47.918 -15.78139)
Caracas           181   POINT (-66.91898 10.50294)
Lima              189  POINT (-77.05201 -12.04607)
Buenos Aires      210  POINT (-58.43251 -34.61071)
Bogota            230    POINT (-74.08529 4.59837)
Santiago          236   POINT (-70.6505 -33.44021)
Rio de Janeiro    238  POINT (-43.21212 -22.90731)
São Paulo         239  POINT (-46.62697 -23.55673)
```

Notice that the original `index` is sorted even though the actual index is not, hence the order is correctly preserved. I suppose that we could improve the docstring within `sindex` to be a bit clearer about what we mean by "index" there.
nicholas-ys-tan: > > The poster requested preservation of the input order. In the case where the input dataframe is not sorted, but clip is run with sort=True, it will return a sorted dataframe. If run with sort=False, it may return with a different unsorted order from the input.
> > I can't think of any cases where a user would want to maintain the order of an unsorted dataframe though.
> 
> I don't think that is correct. The indices coming from `sindex.query` are integer positions of geometries in the original GeoSeries, it does not matter how their index looks like. So even if the GeoSeries to be clipped has a custom, unsorted index, with `sort=True`, the original order is fully preserved. Taking an example from your tests:
> 
> ```python
> cities = gpd.read_file(naturalearth_cities)
> world = gpd.read_file(naturalearth_lowres)
> south_america = world[world["continent"] == "South America"]
> cities.reset_index().set_index("name").clip(south_america, sort=True)
> 
>                 index                     geometry
> name                                              
> Georgetown         55    POINT (-58.16703 6.80197)
> Paramaribo         59    POINT (-55.16703 5.83503)
> Asunción           62  POINT (-57.62583 -25.29067)
> Quito              88     POINT (-78.502 -0.21304)
> Valparaíso        101  POINT (-71.61703 -33.04774)
> Sucre             114  POINT (-65.25952 -19.04097)
> La Paz            122  POINT (-68.15193 -16.49603)
> Brasília          169    POINT (-47.918 -15.78139)
> Caracas           181   POINT (-66.91898 10.50294)
> Lima              189  POINT (-77.05201 -12.04607)
> Buenos Aires      210  POINT (-58.43251 -34.61071)
> Bogota            230    POINT (-74.08529 4.59837)
> Santiago          236   POINT (-70.6505 -33.44021)
> Rio de Janeiro    238  POINT (-43.21212 -22.90731)
> São Paulo         239  POINT (-46.62697 -23.55673)
> ```
> 
> Notice that the original `index` is sorted even though the actual index is not, hence the order is correctly preserved. I suppose that we could improve the docstring within `sindex` to be a bit clearer about what we mean by "index" there.

Thank you for taking the time to explain, that makes sense 

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 124ed9a9..79bca301 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,7 @@ API changes:
 
 New methods:
 
+- Added `intersection_all` method from shapely to GeoSeries/GeoDataframe (#3228).
 - Added `line_merge` method from shapely to GeoSeries/GeoDataframe (#3214).
 - Added `set_precision` and `get_precision` methods from shapely to GeoSeries/GeoDataframe (#3175).
 - Added `count_coordinates` method from shapely to GeoSeries/GeoDataframe (#3026).
@@ -48,6 +49,8 @@ New features and improvements:
   the specified columns (#3101).
 - Added support to ``read_file`` for the ``mask`` keyword for the pyogrio engine (#3062).
 - Added support to ``read_file`` for the ``columns`` keyword for the fiona engine (#3133).
+- Add `sort` keyword to `clip` method for GeoSeries and GeoDataFrame to allow optional
+  preservation of the original order of observations. (#3233)
 - Added `show_bbox`, `drop_id` and `to_wgs84` arguments to allow further customization of
   `GeoSeries.to_json` (#3226)
 
@@ -80,7 +83,7 @@ Deprecations and compatibility notes:
 - The deprecation of `geopandas.datasets` has been enforced and the module has been
   removed. New sample datasets are now available in the
   [geodatasets](https://geodatasets.readthedocs.io/en/latest/) package (#3084).
-- Many longstanding deprecated functions, methods and properties have been removed (#3174), (#3189)
+- Many longstanding deprecated functions, methods and properties have been removed (#3174), (#3190)
   - Removed deprecated functions
     `geopandas.io.read_file`, `geopandas.io.to_file` and `geopandas.io.sql.read_postgis`.
     `geopandas.read_file`, `geopandas.read_postgis` and the GeoDataFrame/GeoSeries `to_file(..)`
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index ebeec66a..74d61f0a 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -147,6 +147,7 @@ Aggregating and exploding
    :toctree: api/
 
    GeoSeries.union_all
+   GeoSeries.intersection_all
    GeoSeries.explode
 
 Serialization / IO / conversion
diff --git a/geopandas/array.py b/geopandas/array.py
index 7d18a9dd..d1bcdb26 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -845,6 +845,9 @@ class GeometryArray(ExtensionArray):
     def union_all(self):
         return shapely.union_all(self._data)
 
+    def intersection_all(self):
+        return shapely.intersection_all(self._data)
+
     #
     # Affinity operations
     #
diff --git a/geopandas/base.py b/geopandas/base.py
index c19b4d25..02eff5fd 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -1732,18 +1732,43 @@ GeometryCollection
         --------
 
         >>> from shapely.geometry import box
-        >>> s = geopandas.GeoSeries([box(0,0,1,1), box(0,0,2,2)])
+        >>> s = geopandas.GeoSeries([box(0, 0, 1, 1), box(0, 0, 2, 2)])
         >>> s
         0    POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))
         1    POLYGON ((2 0, 2 2, 0 2, 0 0, 2 0))
         dtype: geometry
 
-        >>> union = s.union_all()
-        >>> print(union)
-        POLYGON ((0 1, 0 2, 2 2, 2 0, 1 0, 0 0, 0 1))
+        >>> s.union_all()
+        <POLYGON ((0 1, 0 2, 2 2, 2 0, 1 0, 0 0, 0 1))>
         """
         return self.geometry.values.union_all()
 
+    def intersection_all(self):
+        """Returns a geometry containing the intersection of all geometries in
+        the ``GeoSeries``.
+
+        This method ignores None values when other geometries are present.
+        If all elements of the GeoSeries are None, an empty GeometryCollection is
+        returned.
+
+        Examples
+        --------
+
+        >>> from shapely.geometry import box
+        >>> s = geopandas.GeoSeries(
+        ...     [box(0, 0, 2, 2), box(1, 1, 3, 3), box(0, 0, 1.5, 1.5)]
+        ... )
+        >>> s
+        0              POLYGON ((2 0, 2 2, 0 2, 0 0, 2 0))
+        1              POLYGON ((3 1, 3 3, 1 3, 1 1, 3 1))
+        2    POLYGON ((1.5 0, 1.5 1.5, 0 1.5, 0 0, 1.5 0))
+        dtype: geometry
+
+        >>> s.intersection_all()
+        <POLYGON ((1 1, 1 1.5, 1.5 1.5, 1.5 1, 1 1))>
+        """
+        return self.geometry.values.intersection_all()
+
     #
     # Binary operations that return a pandas Series
     #
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index d7885fae..f25187bf 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -2311,7 +2311,7 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
             exclusive=exclusive,
         )
 
-    def clip(self, mask, keep_geom_type=False):
+    def clip(self, mask, keep_geom_type=False, sort=False):
         """Clip points, lines, or polygon geometries to the mask extent.
 
         Both layers must be in the same Coordinate Reference System (CRS).
@@ -2334,6 +2334,10 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
             If True, return only geometries of original type in case of intersection
             resulting in multiple geometry types or GeometryCollections.
             If False, return all resulting geometries (potentially mixed types).
+        sort : boolean, default False
+            If True, the order of rows in the clipped GeoDataFrame will be preserved at
+            small performance cost. If False the order of rows in the clipped
+            GeoDataFrame will be random.
 
         Returns
         -------
@@ -2364,7 +2368,7 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
         >>> nws_groceries.shape
         (7, 8)
         """
-        return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type)
+        return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type, sort=sort)
 
     def overlay(self, right, how="intersection", keep_geom_type=None, make_valid=True):
         """Perform spatial overlay between GeoDataFrames.
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 51a655a7..60de6990 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -1306,7 +1306,7 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
         """
         return Series(to_wkt(self.array, **kwargs), index=self.index)
 
-    def clip(self, mask, keep_geom_type: bool = False) -> GeoSeries:
+    def clip(self, mask, keep_geom_type: bool = False, sort=False) -> GeoSeries:
         """Clip points, lines, or polygon geometries to the mask extent.
 
         Both layers must be in the same Coordinate Reference System (CRS).
@@ -1329,6 +1329,10 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
             If True, return only geometries of original type in case of intersection
             resulting in multiple geometry types or GeometryCollections.
             If False, return all resulting geometries (potentially mixed-types).
+        sort : boolean, default False
+            If True, the order of rows in the clipped GeoSeries will be preserved
+            at small performance cost.
+            If False the order of rows in the clipped GeoSeries will be random.
 
         Returns
         -------
@@ -1359,4 +1363,4 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
         >>> nws_groceries.shape
         (7,)
         """
-        return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type)
+        return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type, sort=sort)
diff --git a/geopandas/tools/clip.py b/geopandas/tools/clip.py
index 728ff0d7..267f5d51 100644
--- a/geopandas/tools/clip.py
+++ b/geopandas/tools/clip.py
@@ -22,7 +22,7 @@ def _mask_is_list_like_rectangle(mask):
     )
 
 
-def _clip_gdf_with_mask(gdf, mask):
+def _clip_gdf_with_mask(gdf, mask, sort=False):
     """Clip geometry to the polygon/rectangle extent.
 
     Clip an input GeoDataFrame to the polygon extent of the polygon
@@ -36,6 +36,10 @@ def _clip_gdf_with_mask(gdf, mask):
     mask : (Multi)Polygon, list-like
         Reference polygon/rectangle for clipping.
 
+    sort : boolean, default False
+        If True, the results will be sorted in ascending order using the
+        geometries' indexes as the primary key.
+
     Returns
     -------
     GeoDataFrame
@@ -48,7 +52,9 @@ def _clip_gdf_with_mask(gdf, mask):
     else:
         intersection_polygon = mask
 
-    gdf_sub = gdf.iloc[gdf.sindex.query(intersection_polygon, predicate="intersects")]
+    gdf_sub = gdf.iloc[
+        gdf.sindex.query(intersection_polygon, predicate="intersects", sort=sort)
+    ]
 
     # For performance reasons points don't need to be intersected with poly
     non_point_mask = gdf_sub.geom_type != "Point"
@@ -82,7 +88,7 @@ def _clip_gdf_with_mask(gdf, mask):
     return clipped
 
 
-def clip(gdf, mask, keep_geom_type=False):
+def clip(gdf, mask, keep_geom_type=False, sort=False):
     """Clip points, lines, or polygon geometries to the mask extent.
 
     Both layers must be in the same Coordinate Reference System (CRS).
@@ -113,6 +119,9 @@ def clip(gdf, mask, keep_geom_type=False):
         If True, return only geometries of original type in case of intersection
         resulting in multiple geometry types or GeometryCollections.
         If False, return all resulting geometries (potentially mixed-types).
+    sort : boolean, default False
+        If True, the results will be sorted in ascending order using the
+        geometries' indexes as the primary key.
 
     Returns
     -------
@@ -190,7 +199,7 @@ def clip(gdf, mask, keep_geom_type=False):
     else:
         combined_mask = mask
 
-    clipped = _clip_gdf_with_mask(gdf, combined_mask)
+    clipped = _clip_gdf_with_mask(gdf, combined_mask, sort=sort)
 
     if keep_geom_type:
         geomcoll_concat = (clipped.geom_type == "GeometryCollection").any()


Test patch:
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 67754248..4db468bd 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -1066,6 +1066,29 @@ class TestDataFrame:
         result = left.clip(south_america)
         assert_geodataframe_equal(result, expected)
 
+    def test_clip_sorting(self, naturalearth_cities, naturalearth_lowres):
+        """
+        Test sorting of geodataframe when clipping.
+        """
+        cities = read_file(naturalearth_cities)
+        world = read_file(naturalearth_lowres)
+        south_america = world[world["continent"] == "South America"]
+
+        unsorted_clipped_cities = geopandas.clip(cities, south_america, sort=False)
+        sorted_clipped_cities = geopandas.clip(cities, south_america, sort=True)
+
+        expected_sorted_index = pd.Index(
+            [55, 59, 62, 88, 101, 114, 122, 169, 181, 189, 210, 230, 236, 238, 239]
+        )
+
+        assert not (
+            sorted(unsorted_clipped_cities.index) == unsorted_clipped_cities.index
+        ).all()
+        assert (
+            sorted(sorted_clipped_cities.index) == sorted_clipped_cities.index
+        ).all()
+        assert_index_equal(expected_sorted_index, sorted_clipped_cities.index)
+
     def test_overlay(self, dfs, how):
         """
         Basic test for availability of the GeoDataFrame method. Other
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 058594ef..f0742954 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -16,6 +16,7 @@ from shapely.geometry import (
     MultiPolygon,
     Point,
     Polygon,
+    box,
 )
 from shapely.geometry.collection import GeometryCollection
 from shapely.ops import unary_union
@@ -494,6 +495,18 @@ class TestGeomMethods:
             result = g.unary_union
         assert result == g.union_all()
 
+    def test_intersection_all(self):
+        expected = Polygon([(1, 1), (1, 1.5), (1.5, 1.5), (1.5, 1), (1, 1)])
+        g = GeoSeries([box(0, 0, 2, 2), box(1, 1, 3, 3), box(0, 0, 1.5, 1.5)])
+
+        assert g.intersection_all().equals(expected)
+
+        g2 = GeoSeries([box(0, 0, 2, 2), None])
+        assert g2.intersection_all().equals(g2[0])
+
+        g3 = GeoSeries([None, None])
+        assert g3.intersection_all().equals(shapely.GeometryCollection())
+
     def test_contains(self):
         expected = [True, False, True, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.contains(self.t1))
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index 9ee56150..df9d203a 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -436,6 +436,29 @@ class TestSeries:
         result = left.geometry.clip(south_america)
         assert_geoseries_equal(result, expected)
 
+    def test_clip_sorting(self, naturalearth_cities, naturalearth_lowres):
+        """
+        Test sorting of geodseries when clipping.
+        """
+        cities = read_file(naturalearth_cities)
+        world = read_file(naturalearth_lowres)
+        south_america = world[world["continent"] == "South America"]
+
+        unsorted_clipped_cities = clip(cities, south_america, sort=False)
+        sorted_clipped_cities = clip(cities, south_america, sort=True)
+
+        expected_sorted_index = pd.Index(
+            [55, 59, 62, 88, 101, 114, 122, 169, 181, 189, 210, 230, 236, 238, 239]
+        )
+
+        assert not (
+            sorted(unsorted_clipped_cities.index) == unsorted_clipped_cities.index
+        ).all()
+        assert (
+            sorted(sorted_clipped_cities.index) == sorted_clipped_cities.index
+        ).all()
+        assert_index_equal(expected_sorted_index, sorted_clipped_cities.index)
+
     def test_from_xy_points(self):
         x = self.landmarks.x.values
         y = self.landmarks.y.values
diff --git a/geopandas/tools/tests/test_clip.py b/geopandas/tools/tests/test_clip.py
index 4092d67a..36d63f11 100644
--- a/geopandas/tools/tests/test_clip.py
+++ b/geopandas/tools/tests/test_clip.py
@@ -1,6 +1,7 @@
 """Tests for the clip module."""
 
 import numpy as np
+import pandas as pd
 
 import shapely
 from shapely.geometry import (
@@ -18,6 +19,7 @@ from geopandas import GeoDataFrame, GeoSeries, clip
 from geopandas._compat import HAS_PYPROJ
 
 from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
+from pandas.testing import assert_index_equal
 import pytest
 
 from geopandas.tools.clip import _mask_is_list_like_rectangle
@@ -42,6 +44,14 @@ def point_gdf():
     return gdf
 
 
+@pytest.fixture
+def point_gdf2():
+    """Create a point GeoDataFrame."""
+    pts = np.array([[5, 5], [2, 2], [4, 4], [0, 0], [3, 3], [1, 1]])
+    gdf = GeoDataFrame([Point(xy) for xy in pts], columns=["geometry"], crs="EPSG:3857")
+    return gdf
+
+
 @pytest.fixture
 def pointsoutside_nooverlap_gdf():
     """Create a point GeoDataFrame. Its points are all outside the single
@@ -460,3 +470,16 @@ def test_clip_empty_mask(buffered_locations, mask):
     )
     clipped = clip(buffered_locations.geometry, mask)
     assert_geoseries_equal(clipped, GeoSeries([], crs="EPSG:3857"))
+
+
+def test_clip_sorting(point_gdf2):
+    """Test the sorting kwarg in clip"""
+    bbox = shapely.geometry.box(0, 0, 2, 2)
+    unsorted_clipped_gdf = point_gdf2.clip(bbox)
+    sorted_clipped_gdf = point_gdf2.clip(bbox, sort=True)
+
+    expected_sorted_index = pd.Index([1, 3, 5])
+
+    assert not (sorted(unsorted_clipped_gdf.index) == unsorted_clipped_gdf.index).all()
+    assert (sorted(sorted_clipped_gdf.index) == sorted_clipped_gdf.index).all()
+    assert_index_equal(expected_sorted_index, sorted_clipped_gdf.index)