Repository: geopandas/geopandas

Issue:
ENH: Add groupby's observed kwarg to dissolve
Using `dissolve` on a `GeoDataFrame` with categorical and non-categorical columns produces unwanted results.

Consider this `DataFrame`, which currently has no categorical columns.
```
>>> df = pd.DataFrame({"cat": ['a', 'a', 'b', 'b'],
                       "noncat": [1, 1, 1, 2],
                       "to_agg": [1, 2, 3, 4],
                       "geometry": ['POINT (0 0)', 'POINT (1 1)', 'POINT (2 2)', 'POINT (4 4)']})
```
If we turn this into a `GeoDataFrame` and dissolve it by the columns `cat` and `noncat`, we get the desired results.
```
>>> gdf = geopandas.GeoDataFrame(df, geometry=geopandas.array.from_wkt(df['geometry'].values))
>>> gdf.dissolve(['cat', 'noncat'])
                                                 geometry  to_agg
cat noncat
a   1       MULTIPOINT (0.00000 0.00000, 1.00000 1.00000)       1
b   1                             POINT (2.00000 2.00000)       3
    2                             POINT (4.00000 4.00000)       4
```

However, if we turn one of the groupby columns into a categorical, we do not get the desired results.
```
>>> df2 = df.astype({'cat': 'category'})
>>> gdf2 = geopandas.GeoDataFrame(df2, geometry=geopandas.array.from_wkt(df2['geometry'].values))
>>> gdf2.dissolve(['cat', 'noncat'])
                                                 geometry  to_agg
cat noncat
a   1       MULTIPOINT (0.00000 0.00000, 1.00000 1.00000)     1.0
    2                                                None     NaN
b   1                             POINT (2.00000 2.00000)     3.0
    2                             POINT (4.00000 4.00000)     4.0
```
We have an additional entry in the index for `cat='a'`, `noncat=2` which doesn't exist in the original data.

In [`pandas.groupby`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.groupby.html), this behavior can be controlled by the `observed` parameter.
> **observed** _bool, default False_
This only applies if any of the groupers are Categoricals. If True: only show observed values for categorical groupers. If False: show all values for categorical groupers.

`observed` is not used in `GeoDataFrame.dissolve`'s calls to `groupby`, nor is it controllable with an arg. If we modify `GeoDataFrame.dissolve` to have an `observed=False` kwarg and pass `observed` into both `groupby` operations that it performs internally, then the desired result is obtained.

```python
def new_dissolve(self, by=None, aggfunc="first", as_index=True, observed=False):
    if by is None:
        by = np.zeros(len(self), dtype="int64")

    # Process non-spatial component
    data = self.drop(labels=self.geometry.name, axis=1)
    aggregated_data = data.groupby(by=by, observed=observed).agg(aggfunc)

    # Process spatial component
    def merge_geometries(block):
        merged_geom = block.unary_union
        return merged_geom

    g = self.groupby(by=by, group_keys=False, observed=observed)[self.geometry.name].agg(
        merge_geometries
    )

    # Aggregate
    aggregated_geometry = GeoDataFrame(g, geometry=self.geometry.name, crs=self.crs)
    # Recombine
    aggregated = aggregated_geometry.join(aggregated_data)

    # Reset if requested
    if not as_index:
        aggregated = aggregated.reset_index()

    return aggregated
```
```
>>> GeoDataFrame.new_dissolve = new_dissolve
>>> gdf2.new_dissolve(['cat', 'noncat'], observed=True)
                                                 geometry  to_agg
cat noncat
a   1       MULTIPOINT (0.00000 0.00000, 1.00000 1.00000)       1
b   1                             POINT (2.00000 2.00000)       3
    2                             POINT (4.00000 4.00000)       4
```

#### API breaking implications

None. We simply expose a kwarg on `dissolve` that already exists on the underlying `groupby` function. Giving it the default value `False` preserves the existing behavior.

#### Describe alternatives you've considered

In my application code, we have worked around this by converting all categorical columns into non-categorical columns (with the values being their category codes) before `groupby` and converting them back after.

```python
def categorical_workaround_pre_groupby(df):
    """Workaround (part a) for dissolve not exposing the 'observed' parameter for categorical"""
    categories_map = dict()
    for col in df:
        if isinstance(df[col].dtype, pd.CategoricalDtype):
            categories_map[col] = df[col].cat.categories
            df[col] = df[col].cat.codes
    return categories_map

def categorical_workaround_post_groupby(df, categories_map):
    """Workaround (part b) for dissolve not exposing the 'observed' parameter for categorical"""
    for col in categories_map:
        df[col] = pd.Categorical.from_codes(df[col], categories_map[col])
```

This is possible, it does what we need, but it is not ideal.


Discussion:
johnflavin: > Thanks! Are there any other keywords in pandas groupby which would make sense to support in dissolve?

That's a good question. Potentially several of them could make sense.

Here are [all the args to `pandas.DataFrame.groupby`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.groupby.html) as of this version:
* `by` - Already supported
* `axis` - Could make sense
* `level` - Could make sense
* `as_index` - Already supported
* `sort` - Could make sense
* `group_keys` - Would be tricky to support since we use it differently in the two calls to `groupby` in the `dissolve` implementation
* `squeeze` - Deprecated, so no
* `observed` - This PR
* `dropna` - Could make sense

The tricky part for me would be not be adding the args but writing tests that exercise all the functions of those args and making sure it makes sense to support them. I knew how to test `observed` because I had an example where I could see the need for it. `sort` and `dropna` seem fairly straightforward. But I would have to do some exploration to figure out how to test `axis` and `level`.
martinfleis: Thanks for the overview!

I can imagine a situation when dissolution along rows as opposed to columns (`axis` kwd) can be useful, it just requires a bit of mental shift from traditional GIS thinking geom+attributes.

`level` should be straightforward and makes sense. `sort` as well and may give us tiny performance gain (although compared to the cost of `unary_union` it probably does not matter much).

`group_keys` are probably not useful since we always use `.agg(aggfunc)`, never `apply`.

`dropna` is certainly useful and we even have an issue on that #1814

If you'd be up for this, it would be great to implement all these which make sense, instead of the piecemeal addition of one keyword at a time. We can help to design tests if you'll struggle yourself.
johnflavin: > I can imagine a situation when dissolution along rows as opposed to columns (`axis` kwd) can be useful, it just requires a bit of mental shift from traditional GIS thinking geom+attributes.

Maybe though that is a good reason not to include `axis`. The geometry in a `GeoDataFrame` is always stored as a column, so if we are dissolving it must be the case that we want to group by other columns. If a user wants to store their geometries on a row across columns, 1. I don't know how or why they would do that, since everything in `GeoDataFrame` is built around the geometry column, and 2. they are probably comfortable making their own dissolve mechanism and don't need official support.

I can implement `level`, `sort`, and `dropna` without any trouble. I think I'll be able to come up with a test for each and we can discuss if the way I did them makes sense.
martinfleis: Okay, let's leave `axis` out for now.
jorisvandenbossche: @johnflavin can you merge the latest master? (there is a merge conflict)
johnflavin: > @johnflavin can you merge the latest master? (there is a merge conflict)

👍 Rebased on master

Patch:
diff --git a/doc/source/docs/reference/sindex.rst b/doc/source/docs/reference/sindex.rst
index ffed1c98..3201106d 100644
--- a/doc/source/docs/reference/sindex.rst
+++ b/doc/source/docs/reference/sindex.rst
@@ -18,22 +18,13 @@ Constructor
 
     GeoSeries.sindex
 
-PyGEOS STRtree
---------------
-.. currentmodule:: geopandas.sindex.PyGEOSSTRTreeIndex
-.. autosummary::
-   :toctree: api/
+Spatial Index object
+--------------------
 
-    intersection
-    is_empty
-    query
-    query_bulk
-    size
-    valid_query_predicates
+The spatial index object returned from :attr:`GeoSeries.sindex` has the following
+methods:
 
-rtree Rtree
------------
-.. currentmodule:: geopandas.sindex.RTreeIndex
+.. currentmodule:: geopandas.sindex.SpatialIndex
 .. autosummary::
    :toctree: api/
 
@@ -44,7 +35,11 @@ rtree Rtree
     size
     valid_query_predicates
 
-Furthermore, the ``rtree``-based spatial index offers full capability of
+The concrete implementations currently available are
+``geopandas.sindex.PyGEOSSTRTreeIndex`` and ``geopandas.sindex.RTreeIndex``.
+
+In addition to the methods listed above, the ``rtree``-based spatial index
+(``geopandas.sindex.RTreeIndex``) offers the full capability of
 ``rtree.index.Index`` - see the full API in the `rtree documentation`_.
 
 .. _rtree documentation: https://rtree.readthedocs.io/en/stable/class.html
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index d3d18d77..7c66425a 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -13,7 +13,7 @@ from pyproj import CRS
 
 from geopandas.array import GeometryArray, from_shapely, GeometryDtype
 from geopandas.base import GeoPandasBase, is_geometry_type
-from geopandas.geoseries import GeoSeries
+from geopandas.geoseries import GeoSeries, inherit_doc
 import geopandas.io
 from geopandas.plotting import plot_dataframe
 from . import _compat as compat
@@ -1293,6 +1293,16 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
             result.__class__ = DataFrame
         return result
 
+    @inherit_doc(pd.DataFrame)
+    def apply(self, func, axis=0, raw=False, result_type=None, args=(), **kwargs):
+        result = super().apply(
+            func, axis=axis, raw=raw, result_type=result_type, args=args, **kwargs
+        )
+        if isinstance(result, GeoDataFrame):
+            if self.crs is not None and result.crs is None:
+                result.set_crs(self.crs, inplace=True)
+        return result
+
     @property
     def _constructor(self):
         return GeoDataFrame
@@ -1312,7 +1322,16 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
 
         return self
 
-    def dissolve(self, by=None, aggfunc="first", as_index=True):
+    def dissolve(
+        self,
+        by=None,
+        aggfunc="first",
+        as_index=True,
+        level=None,
+        sort=True,
+        observed=False,
+        dropna=True,
+    ):
         """
         Dissolve geometries within `groupby` into single observation.
         This is accomplished by applying the `unary_union` method
@@ -1331,6 +1350,33 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
             with each group. Passed to pandas `groupby.agg` method.
         as_index : boolean, default True
             If true, groupby columns become index of result.
+        level : int or str or sequence of int or sequence of str, default None
+            If the axis is a MultiIndex (hierarchical), group by a
+            particular level or levels.
+
+            .. versionadded:: 0.9.0
+        sort : bool, default True
+            Sort group keys. Get better performance by turning this off.
+            Note this does not influence the order of observations within
+            each group. Groupby preserves the order of rows within each group.
+
+            .. versionadded:: 0.9.0
+        observed : bool, default False
+            This only applies if any of the groupers are Categoricals.
+            If True: only show observed values for categorical groupers.
+            If False: show all values for categorical groupers.
+
+            .. versionadded:: 0.9.0
+        dropna : bool, default True
+            If True, and if group keys contain NA values, NA values
+            together with row/column will be dropped. If False, NA
+            values will also be treated as the key in groups.
+
+            This parameter is not supported for pandas < 1.1.0.
+            A warning will be emitted for earlier pandas versions
+            if a non-default value is given for this parameter.
+
+            .. versionadded:: 0.9.0
 
         Returns
         -------
@@ -1363,19 +1409,28 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
 
         """
 
-        if by is None:
+        if by is None and level is None:
             by = np.zeros(len(self), dtype="int64")
 
+        groupby_kwargs = dict(
+            by=by, level=level, sort=sort, observed=observed, dropna=dropna
+        )
+        if not compat.PANDAS_GE_11:
+            groupby_kwargs.pop("dropna")
+
+            if not dropna:  # If they passed a non-default dropna value
+                warnings.warn("dropna kwarg is not supported for pandas < 1.1.0")
+
         # Process non-spatial component
         data = self.drop(labels=self.geometry.name, axis=1)
-        aggregated_data = data.groupby(by=by).agg(aggfunc)
+        aggregated_data = data.groupby(**groupby_kwargs).agg(aggfunc)
 
         # Process spatial component
         def merge_geometries(block):
             merged_geom = block.unary_union
             return merged_geom
 
-        g = self.groupby(by=by, group_keys=False)[self.geometry.name].agg(
+        g = self.groupby(group_keys=False, **groupby_kwargs)[self.geometry.name].agg(
             merge_geometries
         )
 
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 03eab7c9..80f00596 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -1,3 +1,6 @@
+from textwrap import dedent
+import warnings
+
 from shapely.geometry.base import BaseGeometry
 import pandas as pd
 import numpy as np
@@ -21,22 +24,288 @@ def _get_sindex_class():
     )
 
 
+class BaseSpatialIndex:
+    @property
+    def valid_query_predicates(self):
+        """Returns valid predicates for this spatial index.
+
+        Returns
+        -------
+        set
+            Set of valid predicates for this spatial index.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> s = geopandas.GeoSeries([Point(0, 0), Point(1, 1)])
+        >>> s.sindex.valid_query_predicates  # doctest: +SKIP
+        {'contains', 'crosses', 'intersects', 'within', 'touches', \
+'overlaps', None, 'covers', 'contains_properly'}
+        """
+        raise NotImplementedError
+
+    def query(self, geometry, predicate=None, sort=False):
+        """Return the index of all geometries in the tree with extents that
+        intersect the envelope of the input geometry.
+
+        When using the ``rtree`` package, this is not a vectorized function.
+        If speed is important, please use PyGEOS.
+
+        Parameters
+        ----------
+        geometry : shapely geometry
+            A single shapely geometry to query against the spatial index.
+        predicate : {None, 'intersects', 'within', 'contains', \
+'overlaps', 'crosses', 'touches'}, optional
+            If predicate is provided, the input geometry is
+            tested using the predicate function against each item
+            in the tree whose extent intersects the envelope of the
+            input geometry: predicate(input_geometry, tree_geometry).
+            If possible, prepared geometries are used to help
+            speed up the predicate operation.
+        sort : bool, default False
+            If True, the results will be sorted in ascending order.
+            If False, results are often sorted but there is no guarantee.
+
+        Returns
+        -------
+        matches : ndarray of shape (n_results, )
+            Integer indices for matching geometries from the spatial index.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point, box
+        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+        >>> s
+        0    POINT (0.00000 0.00000)
+        1    POINT (1.00000 1.00000)
+        2    POINT (2.00000 2.00000)
+        3    POINT (3.00000 3.00000)
+        4    POINT (4.00000 4.00000)
+        5    POINT (5.00000 5.00000)
+        6    POINT (6.00000 6.00000)
+        7    POINT (7.00000 7.00000)
+        8    POINT (8.00000 8.00000)
+        9    POINT (9.00000 9.00000)
+        dtype: geometry
+
+        >>> s.sindex.query(box(1, 1, 3, 3))
+        array([1, 2, 3])
+
+        >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
+        array([2])
+        """
+        raise NotImplementedError
+
+    def query_bulk(self, geometry, predicate=None, sort=False):
+        """
+        Returns all combinations of each input geometry and geometries in
+        the tree where the envelope of each input geometry intersects with
+        the envelope of a tree geometry.
+
+        In the context of a spatial join, input geometries are the “left”
+        geometries that determine the order of the results, and tree geometries
+        are “right” geometries that are joined against the left geometries.
+        This effectively performs an inner join, where only those combinations
+        of geometries that can be joined based on envelope overlap or optional
+        predicate are returned.
+
+        When using the ``rtree`` package, this is not a vectorized function
+        and may be slow. If speed is important, please use PyGEOS.
+
+        Parameters
+        ----------
+        geometry : {GeoSeries, GeometryArray, numpy.array of PyGEOS geometries}
+            Accepts GeoPandas geometry iterables (GeoSeries, GeometryArray)
+            or a numpy array of PyGEOS geometries.
+        predicate : {None, 'intersects', 'within', 'contains', 'overlaps', \
+'crosses', 'touches'}, optional
+            If predicate is provided, the input geometries are tested using
+            the predicate function against each item in the tree whose extent
+            intersects the envelope of the each input geometry:
+            predicate(input_geometry, tree_geometry).  If possible, prepared
+            geometries are used to help speed up the predicate operation.
+        sort : bool, default False
+            If True, results sorted lexicographically using
+            geometry's indexes as the primary key and the sindex's indexes as the
+            secondary key. If False, no additional sorting is applied.
+
+        Returns
+        -------
+        ndarray with shape (2, n)
+            The first subarray contains input geometry integer indexes.
+            The second subarray contains tree geometry integer indexes.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point, box
+        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+        >>> s
+        0    POINT (0.00000 0.00000)
+        1    POINT (1.00000 1.00000)
+        2    POINT (2.00000 2.00000)
+        3    POINT (3.00000 3.00000)
+        4    POINT (4.00000 4.00000)
+        5    POINT (5.00000 5.00000)
+        6    POINT (6.00000 6.00000)
+        7    POINT (7.00000 7.00000)
+        8    POINT (8.00000 8.00000)
+        9    POINT (9.00000 9.00000)
+        dtype: geometry
+        >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
+        >>> s2
+        0    POLYGON ((4.00000 2.00000, 4.00000 4.00000, 2....
+        1    POLYGON ((6.00000 5.00000, 6.00000 6.00000, 5....
+        dtype: geometry
+
+        >>> s.sindex.query_bulk(s2)
+        array([[0, 0, 0, 1, 1],
+                [2, 3, 4, 5, 6]])
+
+        >>> s.sindex.query_bulk(s2, predicate="contains")
+        array([[0],
+                [3]])
+        """
+        raise NotImplementedError
+
+    def intersection(self, coordinates):
+        """Compatibility wrapper for rtree.index.Index.intersection,
+        use ``query`` intead.
+
+        Parameters
+        ----------
+        coordinates : sequence or array
+            Sequence of the form (min_x, min_y, max_x, max_y)
+            to query a rectangle or (x, y) to query a point.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point, box
+        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+        >>> s
+        0    POINT (0.00000 0.00000)
+        1    POINT (1.00000 1.00000)
+        2    POINT (2.00000 2.00000)
+        3    POINT (3.00000 3.00000)
+        4    POINT (4.00000 4.00000)
+        5    POINT (5.00000 5.00000)
+        6    POINT (6.00000 6.00000)
+        7    POINT (7.00000 7.00000)
+        8    POINT (8.00000 8.00000)
+        9    POINT (9.00000 9.00000)
+        dtype: geometry
+
+        >>> s.sindex.intersection(box(1, 1, 3, 3).bounds)
+        array([1, 2, 3])
+
+        Alternatively, you can use ``query``:
+
+        >>> s.sindex.query(box(1, 1, 3, 3))
+        array([1, 2, 3])
+
+        """
+        raise NotImplementedError
+
+    @property
+    def size(self):
+        """Size of the spatial index
+
+        Number of leaves (input geometries) in the index.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+        >>> s
+        0    POINT (0.00000 0.00000)
+        1    POINT (1.00000 1.00000)
+        2    POINT (2.00000 2.00000)
+        3    POINT (3.00000 3.00000)
+        4    POINT (4.00000 4.00000)
+        5    POINT (5.00000 5.00000)
+        6    POINT (6.00000 6.00000)
+        7    POINT (7.00000 7.00000)
+        8    POINT (8.00000 8.00000)
+        9    POINT (9.00000 9.00000)
+        dtype: geometry
+
+        >>> s.sindex.size
+        10
+        """
+        raise NotImplementedError
+
+    @property
+    def is_empty(self):
+        """Check if the spatial index is empty
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
+        >>> s
+        0    POINT (0.00000 0.00000)
+        1    POINT (1.00000 1.00000)
+        2    POINT (2.00000 2.00000)
+        3    POINT (3.00000 3.00000)
+        4    POINT (4.00000 4.00000)
+        5    POINT (5.00000 5.00000)
+        6    POINT (6.00000 6.00000)
+        7    POINT (7.00000 7.00000)
+        8    POINT (8.00000 8.00000)
+        9    POINT (9.00000 9.00000)
+        dtype: geometry
+
+        >>> s.sindex.is_empty
+        False
+
+        >>> s2 = geopandas.GeoSeries()
+        >>> s2.sindex.is_empty
+        True
+        """
+        raise NotImplementedError
+
+
+def doc(docstring):
+    """
+    A decorator take docstring from passed object and it to decorated one.
+    """
+
+    def decorator(decorated):
+        decorated.__doc__ = dedent(docstring.__doc__ or "")
+        return decorated
+
+    return decorator
+
+
 if compat.HAS_RTREE:
 
     import rtree.index  # noqa
     from rtree.core import RTreeError  # noqa
     from shapely.prepared import prep  # noqa
 
-    class SpatialIndex(rtree.index.Index):
+    class SpatialIndex(rtree.index.Index, BaseSpatialIndex):
         """Original rtree wrapper, kept for backwards compatibility."""
 
         def __init__(self, *args):
-            super().__init__(self, *args)
+            warnings.warn(
+                "Directly using SpatialIndex is deprecated, and the class will be "
+                "removed in a future version. Access the spatial index through the "
+                "`GeoSeries.sindex` attribute, or use `rtree.index.Index` directly.",
+                FutureWarning,
+                stacklevel=2,
+            )
+            super().__init__(*args)
 
+        @doc(BaseSpatialIndex.intersection)
+        def intersection(self, coordinates, *args, **kwargs):
+            return super().intersection(coordinates, *args, **kwargs)
+
+        @doc(BaseSpatialIndex.size)
         @property
         def size(self):
             return len(self.leaves()[0][1])
 
+        @doc(BaseSpatialIndex.is_empty)
         @property
         def is_empty(self):
             if len(self.leaves()) > 1:
@@ -74,23 +343,9 @@ if compat.HAS_RTREE:
                 [None] * self.geometries.size, dtype=object
             )
 
+        @doc(BaseSpatialIndex.valid_query_predicates)
         @property
         def valid_query_predicates(self):
-            """Returns valid predicates for this spatial index.
-
-            Returns
-            -------
-            set
-                Set of valid predicates for this spatial index.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point
-            >>> s = geopandas.GeoSeries([Point(0, 0), Point(1, 1)])
-            >>> s.sindex.valid_query_predicates  # doctest: +SKIP
-            {'contains', 'crosses', 'intersects', 'within', 'touches', \
-'overlaps', None, 'covers', 'contains_properly'}
-            """
             return {
                 None,
                 "intersects",
@@ -103,60 +358,8 @@ if compat.HAS_RTREE:
                 "contains_properly",
             }
 
+        @doc(BaseSpatialIndex.query)
         def query(self, geometry, predicate=None, sort=False):
-            """Return the index of all geometries in the tree with extents that
-            intersect the envelope of the input geometry.
-
-            Compatibility layer for pygeos-based ``sindex.query``.
-
-            This is not a vectorized function, if speed is important,
-            please use PyGEOS.
-
-            Parameters
-            ----------
-            geometry : shapely geometry
-                A single shapely geometry to query against the spatial index.
-            predicate : {None, 'intersects', 'within', 'contains', \
-'overlaps', 'crosses', 'touches'}, optional
-                If predicate is provided, the input geometry is
-                tested using the predicate function against each item
-                in the tree whose extent intersects the envelope of the
-                input geometry: predicate(input_geometry, tree_geometry).
-                If possible, prepared geometries are used to help
-                speed up the predicate operation.
-            sort : bool, default False
-                If True, the results will be sorted in ascending order.
-                If False, results are often sorted but there is no guarantee.
-
-            Returns
-            -------
-            matches : ndarray of shape (n_results, )
-                Integer indices for matching geometries from the spatial index.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point, box
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.query(box(1, 1, 3, 3))
-            array([1, 2, 3])
-
-            >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
-            array([2])
-            """
-
             # handle invalid predicates
             if predicate not in self.valid_query_predicates:
                 raise ValueError(
@@ -235,71 +438,8 @@ if compat.HAS_RTREE:
             # unsorted
             return np.array(tree_idx, dtype=np.intp)
 
+        @doc(BaseSpatialIndex.query_bulk)
         def query_bulk(self, geometry, predicate=None, sort=False):
-            """
-            Returns all combinations of each input geometry and geometries in
-            the tree where the envelope of each input geometry intersects with
-            the envelope of a tree geometry.
-
-            Compatibility layer for pygeos-based ``sindex.query_bulk``.
-
-            Iterates over ``geometry`` and queries index.
-            This operation is not vectorized and may be slow.
-            Use PyGEOS with ``query_bulk`` for speed.
-
-            Parameters
-            ----------
-            geometry : {GeoSeries, GeometryArray, numpy.array of PyGEOS geometries}
-                Accepts GeoPandas geometry iterables (GeoSeries, GeometryArray)
-                or a numpy array of PyGEOS geometries.
-            predicate : {None, 'intersects', 'within', 'contains', 'overlaps', \
-'crosses', 'touches'}, optional
-                If predicate is provided, the input geometries are tested using
-                the predicate function against each item in the tree whose extent
-                intersects the envelope of the each input geometry:
-                predicate(input_geometry, tree_geometry).  If possible, prepared
-                geometries are used to help speed up the predicate operation.
-            sort : bool, default False
-                If True, results sorted lexicographically using
-                geometry's indexes as the primary key and the sindex's indexes as the
-                secondary key. If False, no additional sorting is applied.
-
-            Returns
-            -------
-            ndarray with shape (2, n)
-                The first subarray contains input geometry integer indexes.
-                The second subarray contains tree geometry integer indexes.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point, box
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-            >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
-            >>> s2
-            0    POLYGON ((4.00000 2.00000, 4.00000 4.00000, 2....
-            1    POLYGON ((6.00000 5.00000, 6.00000 6.00000, 5....
-            dtype: geometry
-
-            >>> s.sindex.query_bulk(s2)
-            array([[0, 0, 0, 1, 1],
-                   [2, 3, 4, 5, 6]])
-
-            >>> s.sindex.query_bulk(s2, predicate="contains")
-            array([[0],
-                   [3]])
-            """
             # Iterates over geometry, applying func.
             tree_index = []
             input_geometry_index = []
@@ -310,69 +450,13 @@ if compat.HAS_RTREE:
                 input_geometry_index.extend([i] * len(res))
             return np.vstack([input_geometry_index, tree_index])
 
+        @doc(BaseSpatialIndex.intersection)
         def intersection(self, coordinates):
-            """Wrapper for rtree.index.Index.intersection.
-
-            Parameters
-            ----------
-            coordinates : sequence or array
-                Sequence of the form (min_x, min_y, max_x, max_y)
-                to query a rectangle or (x, y) to query a point.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point, box
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.intersection(box(1, 1, 3, 3).bounds)
-            array([1, 2, 3])
-
-            Alternatively, you can use ``query``:
-
-            >>> s.sindex.query(box(1, 1, 3, 3))
-            array([1, 2, 3])
-
-            """
             return super().intersection(coordinates, objects=False)
 
+        @doc(BaseSpatialIndex.size)
         @property
         def size(self):
-            """Size of the spatial index
-
-            Number of leaves (input geometries) in the index.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.size
-            10
-            """
             if hasattr(self, "_size"):
                 size = self._size
             else:
@@ -383,34 +467,9 @@ if compat.HAS_RTREE:
                 self._size = size
             return size
 
+        @doc(BaseSpatialIndex.is_empty)
         @property
         def is_empty(self):
-            """Check if the spatial index is empty
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.is_empty
-            False
-
-            >>> s2 = geopandas.GeoSeries()
-            >>> s2.sindex.is_empty
-            True
-            """
             return self.geometries.size == 0 or self.size == 0
 
         def __len__(self):
@@ -464,61 +523,8 @@ if compat.HAS_PYGEOS:
             """
             return pygeos.strtree.VALID_PREDICATES | set([None])
 
+        @doc(BaseSpatialIndex.query)
         def query(self, geometry, predicate=None, sort=False):
-            """
-            Return the index of all geometries in the tree with extents
-            that intersect the envelope of the input geometry.
-
-            Wrapper for pygeos.STRtree.query.
-
-            This also ensures a deterministic (sorted) order for the results.
-
-            Parameters
-            ----------
-            geometry : single PyGEOS or shapely geometry
-            predicate : {None, 'intersects', 'within', 'contains', \
-'overlaps', 'crosses', 'touches'}, optional
-                If predicate is provided, the input geometry is tested
-                using the predicate function against each item in the
-                tree whose extent intersects the envelope of the input
-                geometry: predicate(input_geometry, tree_geometry).
-            sort : bool, default False
-                If True, the results will be sorted in ascending order.
-                If False, results are often sorted but there is no guarantee.
-
-            Returns
-            -------
-            matches : ndarray of shape (n_results, )
-                Integer indices for matching geometries from the spatial index.
-
-            Notes
-            -----
-            See PyGEOS.strtree documentation for more information.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point, box
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.query(box(1, 1, 3, 3))
-            array([1, 2, 3])
-
-            >>> s.sindex.query(box(1, 1, 3, 3), predicate="contains")
-            array([2])
-            """
-
             if predicate not in self.valid_query_predicates:
                 raise ValueError(
                     "Got `predicate` = `{}`; ".format(predicate)
@@ -537,81 +543,8 @@ if compat.HAS_PYGEOS:
 
             return matches
 
+        @doc(BaseSpatialIndex.query_bulk)
         def query_bulk(self, geometry, predicate=None, sort=False):
-            """
-            Returns all combinations of each input geometry and geometries in
-            the tree where the envelope of each input geometry intersects with
-            the envelope of a tree geometry.
-
-            Wrapper to expose underlaying pygeos objects to pygeos.query_bulk.
-
-            In the context of a spatial join, input geometries are the “left”
-            geometries that determine the order of the results, and tree geometries
-            are “right” geometries that are joined against the left geometries.
-            This effectively performs an inner join, where only those combinations
-            of geometries that can be joined based on envelope overlap or optional
-            predicate are returned.
-
-            This also allows a deterministic (sorted) order for the results.
-
-
-            Parameters
-            ----------
-            geometry : {GeoSeries, GeometryArray, numpy.array of PyGEOS geometries}
-                Accepts GeoPandas geometry iterables (GeoSeries, GeometryArray)
-                or a numpy array of PyGEOS geometries.
-            predicate : {None, 'intersects', 'within', 'contains', \
-'overlaps', 'crosses', 'touches'}, optional
-                If predicate is provided, the input geometry is tested
-                using the predicate function against each item in the
-                index whose extent intersects the envelope of the input geometry:
-                predicate(input_geometry, tree_geometry).
-            sort : bool, default False
-                If True, results sorted lexicographically using
-                geometry's indexes as the primary key and the sindex's indexes as the
-                secondary key. If False, no additional sorting is applied.
-
-            Returns
-            -------
-            ndarray with shape (2, n)
-                The first subarray contains input geometry integer indexes.
-                The second subarray contains tree geometry integer indexes.
-
-            Notes
-            -----
-            See PyGEOS.strtree documentation for more information.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point, box
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-            >>> s2 = geopandas.GeoSeries([box(2, 2, 4, 4), box(5, 5, 6, 6)])
-            >>> s2
-            0    POLYGON ((4.00000 2.00000, 4.00000 4.00000, 2....
-            1    POLYGON ((6.00000 5.00000, 6.00000 6.00000, 5....
-            dtype: geometry
-
-            >>> s.sindex.query_bulk(s2)
-            array([[0, 0, 0, 1, 1],
-                   [2, 3, 4, 5, 6]])
-
-            >>> s.sindex.query_bulk(s2, predicate="contains")
-            array([[0],
-                   [3]])
-            """
-
             if predicate not in self.valid_query_predicates:
                 raise ValueError(
                     "Got `predicate` = `{}`, `predicate` must be one of {}".format(
@@ -635,41 +568,8 @@ if compat.HAS_PYGEOS:
 
             return res
 
+        @doc(BaseSpatialIndex.intersection)
         def intersection(self, coordinates):
-            """Wrapper for pygeos.query that uses the RTree API.
-
-            Compatibility wrapper, use ``query`` instead.
-
-            Parameters
-            ----------
-            coordinates : sequence or array
-                Sequence of the form (min_x, min_y, max_x, max_y)
-                to query a rectangle or (x, y) to query a point.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point, box
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.query(box(1, 1, 3, 3))
-            array([1, 2, 3])
-
-            >>> s.sindex.intersection(box(1, 1, 3, 3).bounds)
-            array([1, 2, 3])
-
-            """
             # convert bounds to geometry
             # the old API uses tuples of bound, but pygeos uses geometries
             try:
@@ -698,60 +598,12 @@ if compat.HAS_PYGEOS:
 
             return indexes
 
+        @doc(BaseSpatialIndex.size)
         @property
         def size(self):
-            """Size of the spatial index
-
-            Number of leaves (input geometries) in the index.
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.size
-            10
-            """
             return len(self)
 
+        @doc(BaseSpatialIndex.is_empty)
         @property
         def is_empty(self):
-            """Check if the spatial index is empty
-
-            Examples
-            --------
-            >>> from shapely.geometry import Point
-            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(10), range(10)))
-            >>> s
-            0    POINT (0.00000 0.00000)
-            1    POINT (1.00000 1.00000)
-            2    POINT (2.00000 2.00000)
-            3    POINT (3.00000 3.00000)
-            4    POINT (4.00000 4.00000)
-            5    POINT (5.00000 5.00000)
-            6    POINT (6.00000 6.00000)
-            7    POINT (7.00000 7.00000)
-            8    POINT (8.00000 8.00000)
-            9    POINT (9.00000 9.00000)
-            dtype: geometry
-
-            >>> s.sindex.is_empty
-            False
-
-            >>> s2 = geopandas.GeoSeries()
-            >>> s2.sindex.is_empty
-            True
-            """
             return len(self) == 0
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index 9aebbb39..5b95c922 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -137,7 +137,7 @@ def _overlay_union(df1, df2):
     return dfunion.reindex(columns=columns)
 
 
-def overlay(df1, df2, how="intersection", keep_geom_type=None):
+def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     """Perform spatial overlay between two GeoDataFrames.
 
     Currently only supports data GeoDataFrames with uniform geometry types,
@@ -159,6 +159,9 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
         if False, return all resulting geometries. Default is None,
         which will set keep_geom_type to True but warn upon dropping
         geometries.
+    make_valid : bool, default True
+        If True, any invalid input geometries are corrected with a call to `buffer(0)`,
+        if False, a `ValueError` is raised if any input geometries are invalid.
 
     Returns
     -------
@@ -179,9 +182,9 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
     >>> geopandas.overlay(df1, df2, how='union')
        df1_data  df2_data                                           geometry
     0       1.0       1.0  POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....
-    1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
+    1       2.0       1.0  POLYGON ((3.00000 2.00000, 2.00000 2.00000, 2....
     2       2.0       2.0  POLYGON ((3.00000 4.00000, 4.00000 4.00000, 4....
-    3       1.0       NaN  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....
+    3       1.0       NaN  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....
     4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
     5       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     6       NaN       2.0  POLYGON ((3.00000 4.00000, 3.00000 5.00000, 5....
@@ -189,27 +192,27 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
     >>> geopandas.overlay(df1, df2, how='intersection')
        df1_data  df2_data                                           geometry
     0         1         1  POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....
-    1         2         1  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
+    1         2         1  POLYGON ((3.00000 2.00000, 2.00000 2.00000, 2....
     2         2         2  POLYGON ((3.00000 4.00000, 4.00000 4.00000, 4....
 
     >>> geopandas.overlay(df1, df2, how='symmetric_difference')
        df1_data  df2_data                                           geometry
-    0       1.0       NaN  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....
+    0       1.0       NaN  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....
     1       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
     2       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     3       NaN       2.0  POLYGON ((3.00000 4.00000, 3.00000 5.00000, 5....
 
     >>> geopandas.overlay(df1, df2, how='difference')
                                                 geometry  df1_data
-    0  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....         1
+    0  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....         1
     1  MULTIPOLYGON (((2.00000 3.00000, 2.00000 4.000...         2
 
     >>> geopandas.overlay(df1, df2, how='identity')
        df1_data  df2_data                                           geometry
     0       1.0       1.0  POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....
-    1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
+    1       2.0       1.0  POLYGON ((3.00000 2.00000, 2.00000 2.00000, 2....
     2       2.0       2.0  POLYGON ((3.00000 4.00000, 4.00000 4.00000, 4....
-    3       1.0       NaN  POLYGON ((0.00000 0.00000, 0.00000 2.00000, 1....
+    3       1.0       NaN  POLYGON ((2.00000 1.00000, 2.00000 0.00000, 0....
     4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
 
     See also
@@ -262,12 +265,24 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None):
             )
 
     # Computations
-    df1 = df1.copy()
-    df2 = df2.copy()
-    if df1.geom_type.isin(polys).all():
-        df1[df1._geometry_column_name] = df1.geometry.buffer(0)
-    if df2.geom_type.isin(polys).all():
-        df2[df2._geometry_column_name] = df2.geometry.buffer(0)
+    def _make_valid(df):
+        df = df.copy()
+        if df.geom_type.isin(polys).all():
+            mask = ~df.geometry.is_valid
+            col = df._geometry_column_name
+            if make_valid:
+                df.loc[mask, col] = df.loc[mask, col].buffer(0)
+            elif mask.any():
+                raise ValueError(
+                    "You have passed make_valid=False along with "
+                    f"{mask.sum()} invalid input geometries. "
+                    "Use make_valid=True or make sure that all geometries "
+                    "are valid before using overlay."
+                )
+        return df
+
+    df1 = _make_valid(df1)
+    df2 = _make_valid(df2)
 
     with warnings.catch_warnings():  # CRS checked above, supress array-level warning
         warnings.filterwarnings("ignore", message="CRS mismatch between the CRS")


Test patch:
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index 14f4f719..ea28f070 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -581,6 +581,16 @@ class TestGeometryArrayCRS:
         result = s.apply(lambda x: x.centroid)
         assert result.crs == 27700
 
+    def test_apply_geodataframe(self):
+        df = GeoDataFrame({"col1": [0, 1]}, geometry=self.geoms, crs=27700)
+        assert df.crs == 27700
+
+        # apply preserves the CRS if the result is a GeoDataFrame
+        result = df.apply(lambda col: col, axis=0)
+        assert result.crs == 27700
+        result = df.apply(lambda row: row, axis=1)
+        assert result.crs == 27700
+
 
 class TestSetCRS:
     @pytest.mark.parametrize(
diff --git a/geopandas/tests/test_dissolve.py b/geopandas/tests/test_dissolve.py
index cc5c6eef..d50d25fe 100644
--- a/geopandas/tests/test_dissolve.py
+++ b/geopandas/tests/test_dissolve.py
@@ -3,6 +3,7 @@ import pandas as pd
 
 import geopandas
 from geopandas import GeoDataFrame, read_file
+from geopandas import _compat as compat
 
 from pandas.testing import assert_frame_equal
 import pytest
@@ -128,3 +129,175 @@ def test_dissolve_none_mean(nybb_polydf):
         crs=nybb_polydf.crs,
     )
     assert_frame_equal(expected, test, check_column_type=False)
+
+
+def test_dissolve_level():
+    gdf = geopandas.GeoDataFrame(
+        {
+            "a": [1, 1, 2, 2],
+            "b": [3, 4, 4, 4],
+            "c": [3, 4, 5, 6],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "POINT (1 1)", "POINT (2 2)", "POINT (3 3)"]
+            ),
+        }
+    ).set_index(["a", "b", "c"])
+
+    expected_a = geopandas.GeoDataFrame(
+        {
+            "a": [1, 2],
+            "geometry": geopandas.array.from_wkt(
+                ["MULTIPOINT (0 0, 1 1)", "MULTIPOINT (2 2, 3 3)"]
+            ),
+        }
+    ).set_index("a")
+    expected_b = geopandas.GeoDataFrame(
+        {
+            "b": [3, 4],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "MULTIPOINT (1 1, 2 2, 3 3)"]
+            ),
+        }
+    ).set_index("b")
+    expected_ab = geopandas.GeoDataFrame(
+        {
+            "a": [1, 1, 2],
+            "b": [3, 4, 4],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "POINT (1 1)", "MULTIPOINT (2 2, 3 3)"]
+            ),
+        }
+    ).set_index(["a", "b"])
+
+    assert_frame_equal(expected_a, gdf.dissolve(level=0))
+    assert_frame_equal(expected_a, gdf.dissolve(level="a"))
+    assert_frame_equal(expected_b, gdf.dissolve(level=1))
+    assert_frame_equal(expected_b, gdf.dissolve(level="b"))
+    assert_frame_equal(expected_ab, gdf.dissolve(level=[0, 1]))
+    assert_frame_equal(expected_ab, gdf.dissolve(level=["a", "b"]))
+
+
+def test_dissolve_sort():
+    gdf = geopandas.GeoDataFrame(
+        {
+            "a": [2, 1, 1],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "POINT (1 1)", "POINT (2 2)"]
+            ),
+        }
+    )
+
+    expected_unsorted = geopandas.GeoDataFrame(
+        {
+            "a": [2, 1],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "MULTIPOINT (1 1, 2 2)"]
+            ),
+        }
+    ).set_index("a")
+    expected_sorted = expected_unsorted.sort_index()
+
+    assert_frame_equal(expected_sorted, gdf.dissolve("a"))
+    assert_frame_equal(expected_unsorted, gdf.dissolve("a", sort=False))
+
+
+@pytest.mark.skipif(
+    not compat.PANDAS_GE_025,
+    reason="'observed' param behavior changed in pandas 0.25.0",
+)
+def test_dissolve_categorical():
+    gdf = geopandas.GeoDataFrame(
+        {
+            "cat": pd.Categorical(["a", "a", "b", "b"]),
+            "noncat": [1, 1, 1, 2],
+            "to_agg": [1, 2, 3, 4],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "POINT (1 1)", "POINT (2 2)", "POINT (3 3)"]
+            ),
+        }
+    )
+
+    # when observed=False we get an additional observation
+    # that wasn't in the original data
+    expected_gdf_observed_false = geopandas.GeoDataFrame(
+        {
+            "cat": pd.Categorical(["a", "a", "b", "b"]),
+            "noncat": [1, 2, 1, 2],
+            "geometry": geopandas.array.from_wkt(
+                [
+                    "MULTIPOINT (0 0, 1 1)",
+                    None,
+                    "POINT (2 2)",
+                    "POINT (3 3)",
+                ]
+            ),
+            "to_agg": [1, None, 3, 4],
+        }
+    ).set_index(["cat", "noncat"])
+
+    # when observed=True we do not get any additional observations
+    expected_gdf_observed_true = geopandas.GeoDataFrame(
+        {
+            "cat": pd.Categorical(["a", "b", "b"]),
+            "noncat": [1, 1, 2],
+            "geometry": geopandas.array.from_wkt(
+                ["MULTIPOINT (0 0, 1 1)", "POINT (2 2)", "POINT (3 3)"]
+            ),
+            "to_agg": [1, 3, 4],
+        }
+    ).set_index(["cat", "noncat"])
+
+    assert_frame_equal(expected_gdf_observed_false, gdf.dissolve(["cat", "noncat"]))
+    assert_frame_equal(
+        expected_gdf_observed_true, gdf.dissolve(["cat", "noncat"], observed=True)
+    )
+
+
+@pytest.mark.skipif(
+    not compat.PANDAS_GE_11, reason="dropna groupby kwarg added in pandas 1.1.0"
+)
+def test_dissolve_dropna():
+    gdf = geopandas.GeoDataFrame(
+        {
+            "a": [1, 1, None],
+            "geometry": geopandas.array.from_wkt(
+                ["POINT (0 0)", "POINT (1 1)", "POINT (2 2)"]
+            ),
+        }
+    )
+
+    expected_with_na = geopandas.GeoDataFrame(
+        {
+            "a": [1.0, np.nan],
+            "geometry": geopandas.array.from_wkt(
+                ["MULTIPOINT (0 0, 1 1)", "POINT (2 2)"]
+            ),
+        }
+    ).set_index("a")
+    expected_no_na = geopandas.GeoDataFrame(
+        {
+            "a": [1.0],
+            "geometry": geopandas.array.from_wkt(["MULTIPOINT (0 0, 1 1)"]),
+        }
+    ).set_index("a")
+
+    assert_frame_equal(expected_with_na, gdf.dissolve("a", dropna=False))
+    assert_frame_equal(expected_no_na, gdf.dissolve("a"))
+
+
+@pytest.mark.skipif(
+    compat.PANDAS_GE_11, reason="dropna warning is only emitted if pandas < 1.1.0"
+)
+def test_dissolve_dropna_warn(nybb_polydf):
+    # No warning with default params
+    with pytest.warns(None) as record:
+        nybb_polydf.dissolve()
+
+    for r in record:
+        assert "dropna kwarg is not supported" not in str(r.message)
+
+    # Warning is emitted with non-default dropna value
+    with pytest.warns(
+        UserWarning, match="dropna kwarg is not supported for pandas < 1.1.0"
+    ):
+        nybb_polydf.dissolve(dropna=False)
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index b8281122..f0c095a6 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -2,7 +2,7 @@ import os
 
 import pandas as pd
 
-from shapely.geometry import Point, Polygon, LineString, GeometryCollection
+from shapely.geometry import Point, Polygon, LineString, GeometryCollection, box
 from fiona.errors import DriverError
 
 import geopandas
@@ -571,3 +571,24 @@ def test_keep_geom_type_geometry_collection():
     intersection = overlay(df1, df2, keep_geom_type=False)
     assert len(intersection) == 1
     assert (intersection.geom_type == "GeometryCollection").all()
+
+
+@pytest.mark.parametrize("make_valid", [True, False])
+def test_overlap_make_valid(make_valid):
+    bowtie = Polygon([(1, 1), (9, 9), (9, 1), (1, 9), (1, 1)])
+    assert not bowtie.is_valid
+    fixed_bowtie = bowtie.buffer(0)
+    assert fixed_bowtie.is_valid
+
+    df1 = GeoDataFrame({"col1": ["region"], "geometry": GeoSeries([box(0, 0, 10, 10)])})
+    df_bowtie = GeoDataFrame(
+        {"col1": ["invalid", "valid"], "geometry": GeoSeries([bowtie, fixed_bowtie])}
+    )
+
+    if make_valid:
+        df_overlay_bowtie = overlay(df1, df_bowtie, make_valid=make_valid)
+        assert df_overlay_bowtie.at[0, "geometry"].equals(fixed_bowtie)
+        assert df_overlay_bowtie.at[1, "geometry"].equals(fixed_bowtie)
+    else:
+        with pytest.raises(ValueError, match="1 invalid input geometries"):
+            overlay(df1, df_bowtie, make_valid=make_valid)
diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index 9ef56961..7d8d35c8 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -718,3 +718,16 @@ class TestPygeosInterface:
 
         res = world.sindex.query_bulk(capitals.geometry, predicate)
         assert res.shape == expected_shape
+
+
+@pytest.mark.skipif(not compat.HAS_RTREE, reason="no rtree installed")
+def test_old_spatial_index_deprecated():
+    t1 = Polygon([(0, 0), (1, 0), (1, 1)])
+    t2 = Polygon([(0, 0), (1, 1), (0, 1)])
+
+    stream = ((i, item.bounds, None) for i, item in enumerate([t1, t2]))
+
+    with pytest.warns(FutureWarning):
+        idx = geopandas.sindex.SpatialIndex(stream)
+
+    assert list(idx.intersection((0, 0, 1, 1))) == [0, 1]