Repository: geopandas/geopandas

Issue:
ENH: support how="inner" in sjoin_nearest
Followup to #1865 .

We should also discuss if `"inner"` should be the default.

Discussion:

Patch:
diff --git a/benchmarks/plotting.py b/benchmarks/plotting.py
index 3b4132eb..5731bdbd 100644
--- a/benchmarks/plotting.py
+++ b/benchmarks/plotting.py
@@ -7,44 +7,63 @@ import numpy as np
 
 class Bench:
 
-    param_names = ['geom_type']
-    params = [('Point', 'LineString', 'Polygon', 'MultiPolygon', 'mixed')]
+    param_names = ["geom_type"]
+    params = [("Point", "LineString", "Polygon", "MultiPolygon", "mixed")]
 
     def setup(self, geom_type):
 
-        if geom_type == 'Point':
+        if geom_type == "Point":
             geoms = GeoSeries([Point(i, i) for i in range(1000)])
-        elif geom_type == 'LineString':
-            geoms = GeoSeries([LineString([(random.random(), random.random())
-                                           for _ in range(5)])
-                               for _ in range(100)])
-        elif geom_type == 'Polygon':
-            geoms = GeoSeries([Polygon([(random.random(), random.random())
-                                        for _ in range(3)])
-                               for _ in range(100)])
-        elif geom_type == 'MultiPolygon':
+        elif geom_type == "LineString":
             geoms = GeoSeries(
-                [MultiPolygon([Polygon([(random.random(), random.random())
-                                        for _ in range(3)])
-                               for _ in range(3)])
-                 for _ in range(20)])
-        elif geom_type == 'mixed':
+                [
+                    LineString([(random.random(), random.random()) for _ in range(5)])
+                    for _ in range(100)
+                ]
+            )
+        elif geom_type == "Polygon":
+            geoms = GeoSeries(
+                [
+                    Polygon([(random.random(), random.random()) for _ in range(3)])
+                    for _ in range(100)
+                ]
+            )
+        elif geom_type == "MultiPolygon":
+            geoms = GeoSeries(
+                [
+                    MultiPolygon(
+                        [
+                            Polygon(
+                                [(random.random(), random.random()) for _ in range(3)]
+                            )
+                            for _ in range(3)
+                        ]
+                    )
+                    for _ in range(20)
+                ]
+            )
+        elif geom_type == "mixed":
             g1 = GeoSeries([Point(i, i) for i in range(100)])
-            g2 = GeoSeries([LineString([(random.random(), random.random())
-                                        for _ in range(5)])
-                            for _ in range(100)])
-            g3 = GeoSeries([Polygon([(random.random(), random.random())
-                                     for _ in range(3)])
-                            for _ in range(100)])
+            g2 = GeoSeries(
+                [
+                    LineString([(random.random(), random.random()) for _ in range(5)])
+                    for _ in range(100)
+                ]
+            )
+            g3 = GeoSeries(
+                [
+                    Polygon([(random.random(), random.random()) for _ in range(3)])
+                    for _ in range(100)
+                ]
+            )
 
             geoms = g1
-            geoms.iloc[np.random.randint(0, 100, 50)] = g2
-            geoms.iloc[np.random.randint(0, 100, 33)] = g3
+            geoms.iloc[np.random.randint(0, 100, 50)] = g2.iloc[:50]
+            geoms.iloc[np.random.randint(0, 100, 33)] = g3.iloc[:33]
 
             print(geoms.geom_type.value_counts())
 
-        df = GeoDataFrame({'geometry': geoms,
-                           'values': np.random.randn(len(geoms))})
+        df = GeoDataFrame({"geometry": geoms, "values": np.random.randn(len(geoms))})
 
         self.geoms = geoms
         self.df = df
@@ -53,5 +72,4 @@ class Bench:
         self.geoms.plot()
 
     def time_plot_values(self, *args):
-        self.df.plot(column='values')
-
+        self.df.plot(column="values")
diff --git a/doc/source/docs/user_guide/mergingdata.rst b/doc/source/docs/user_guide/mergingdata.rst
index 39346df5..72e2ae12 100644
--- a/doc/source/docs/user_guide/mergingdata.rst
+++ b/doc/source/docs/user_guide/mergingdata.rst
@@ -145,9 +145,8 @@ Nearest Joins
 
 Proximity-based joins can be done via :func:`sjoin_nearest`.
 
-:func:`sjoin_nearest` shares the ``how`` argument with :func:`sjoin` but does not support the ``"inner"`` option.
-
-:func:`sjoin_nearest` includes two additional arguments: ``max_distance`` and ``distance_col``.
+:func:`sjoin_nearest` shares the ``how`` argument with :func:`sjoin`, and
+includes two additional arguments: ``max_distance`` and ``distance_col``.
 
 **max_distance**
 
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 2a9b9f4a..f4c51790 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1546,7 +1546,7 @@ individually so that features may have different properties
         return aggregated
 
     # overrides the pandas native explode method to break up features geometrically
-    def explode(self, column=None, **kwargs):
+    def explode(self, column=None, ignore_index=False, index_parts=None, **kwargs):
         """
         Explode muti-part geometries into multiple single geometries.
 
@@ -1554,10 +1554,22 @@ individually so that features may have different properties
         multiple rows with single geometries, thereby increasing the vertical
         size of the GeoDataFrame.
 
-        The index of the input geodataframe is no longer unique and is
-        replaced with a multi-index (original index with additional level
-        indicating the multiple geometries: a new zero-based index for each
-        single part geometry per multi-part geometry).
+        .. note:: ignore_index requires pandas 1.1.0 or newer.
+
+        Parameters
+        ----------
+        column : string, default None
+            Column to explode. In the case of a geometry column, multi-part
+            geometries are converted to single-part.
+            If None, the active geometry column is used.
+        ignore_index : bool, default False
+            If True, the resulting index will be labelled 0, 1, …, n - 1,
+            ignoring `index_parts`.
+        index_parts : boolean, default True
+            If True, the resulting index will be a multi-index (original
+            index with an additional level indicating the multiple
+            geometries: a new zero-based index for each single part geometry
+            per multi-part geometry).
 
         Returns
         -------
@@ -1582,7 +1594,7 @@ individually so that features may have different properties
         0  name1  MULTIPOINT (1.00000 2.00000, 3.00000 4.00000)
         1  name2  MULTIPOINT (2.00000 1.00000, 0.00000 0.00000)
 
-        >>> exploded = gdf.explode()
+        >>> exploded = gdf.explode(index_parts=True)
         >>> exploded
               col1                 geometry
         0 0  name1  POINT (1.00000 2.00000)
@@ -1590,6 +1602,22 @@ individually so that features may have different properties
         1 0  name2  POINT (2.00000 1.00000)
           1  name2  POINT (0.00000 0.00000)
 
+        >>> exploded = gdf.explode(index_parts=False)
+        >>> exploded
+            col1                 geometry
+        0  name1  POINT (1.00000 2.00000)
+        0  name1  POINT (3.00000 4.00000)
+        1  name2  POINT (2.00000 1.00000)
+        1  name2  POINT (0.00000 0.00000)
+
+        >>> exploded = gdf.explode(ignore_index=True)
+        >>> exploded
+            col1                 geometry
+        0  name1  POINT (1.00000 2.00000)
+        1  name1  POINT (3.00000 4.00000)
+        2  name2  POINT (2.00000 1.00000)
+        3  name2  POINT (0.00000 0.00000)
+
         See also
         --------
         GeoDataFrame.dissolve : dissolve geometries into a single observation.
@@ -1601,16 +1629,39 @@ individually so that features may have different properties
             column = self.geometry.name
         # If the specified column is not a geometry dtype use pandas explode
         if not isinstance(self[column].dtype, GeometryDtype):
-            return super().explode(column, **kwargs)
-            # TODO: make sure index behaviour is consistent
+            if compat.PANDAS_GE_11:
+                return super().explode(column, ignore_index=ignore_index, **kwargs)
+            else:
+                return super().explode(column, **kwargs)
+
+        if index_parts is None:
+            if not ignore_index:
+                warnings.warn(
+                    "Currently, index_parts defaults to True, but in the future, "
+                    "it will default to False to be consistent with Pandas. "
+                    "Use `index_parts=True` to keep the current behavior and "
+                    "True/False to silence the warning.",
+                    FutureWarning,
+                    stacklevel=2,
+                )
+            index_parts = True
 
         df_copy = self.copy()
 
-        if "level_1" in df_copy.columns:  # GH1393
-            df_copy = df_copy.rename(columns={"level_1": "__level_1"})
+        level_str = f"level_{df_copy.index.nlevels}"
 
-        exploded_geom = df_copy.geometry.explode().reset_index(level=-1)
-        exploded_index = exploded_geom.columns[0]
+        if level_str in df_copy.columns:  # GH1393
+            df_copy = df_copy.rename(columns={level_str: f"__{level_str}"})
+
+        if index_parts:
+            exploded_geom = df_copy.geometry.explode(index_parts=True)
+            exploded_index = exploded_geom.index
+            exploded_geom = exploded_geom.reset_index(level=-1, drop=True)
+        else:
+            exploded_geom = df_copy.geometry.explode(index_parts=True).reset_index(
+                level=-1, drop=True
+            )
+            exploded_index = exploded_geom.index
 
         df = (
             df_copy.drop(df_copy._geometry_column_name, axis=1)
@@ -1618,13 +1669,19 @@ individually so that features may have different properties
             .__finalize__(self)
         )
 
-        # reset to MultiIndex, otherwise df index is only first level of
-        # exploded GeoSeries index.
-        df.set_index(exploded_index, append=True, inplace=True)
-        df.index.names = list(self.index.names) + [None]
+        if ignore_index:
+            df.reset_index(inplace=True, drop=True)
+        elif index_parts:
+            # reset to MultiIndex, otherwise df index is only first level of
+            # exploded GeoSeries index.
+            df.set_index(exploded_index, inplace=True)
+            df.index.names = list(self.index.names) + [None]
+        else:
+            df.set_index(exploded_index, inplace=True)
+            df.index.names = self.index.names
 
-        if "__level_1" in df.columns:
-            df = df.rename(columns={"__level_1": "level_1"})
+        if f"__{level_str}" in df.columns:
+            df = df.rename(columns={f"__{level_str}": level_str})
 
         geo_df = df.set_geometry(self._geometry_column_name)
         return geo_df
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index ef929791..1cb14e1f 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -757,7 +757,7 @@ class GeoSeries(GeoPandasBase, Series):
         """Interactive map based on folium/leaflet.js"""
         return _explore_geoseries(self, *args, **kwargs)
 
-    def explode(self):
+    def explode(self, ignore_index=False, index_parts=None):
         """
         Explode multi-part geometries into multiple single geometries.
 
@@ -765,6 +765,17 @@ class GeoSeries(GeoPandasBase, Series):
         This is analogous to PostGIS's ST_Dump(). The 'path' index is the
         second level of the returned MultiIndex
 
+        Parameters
+        ----------
+        ignore_index : bool, default False
+            If True, the resulting index will be labelled 0, 1, …, n - 1,
+            ignoring `index_parts`.
+        index_parts : boolean, default True
+            If True, the resulting index will be a multi-index (original
+            index with an additional level indicating the multiple
+            geometries: a new zero-based index for each single part geometry
+            per multi-part geometry).
+
         Returns
         -------
         A GeoSeries with a MultiIndex. The levels of the MultiIndex are the
@@ -782,7 +793,7 @@ class GeoSeries(GeoPandasBase, Series):
         1    MULTIPOINT (2.00000 2.00000, 3.00000 3.00000, ...
         dtype: geometry
 
-        >>> s.explode()
+        >>> s.explode(index_parts=True)
         0  0    POINT (0.00000 0.00000)
            1    POINT (1.00000 1.00000)
         1  0    POINT (2.00000 2.00000)
@@ -795,6 +806,16 @@ class GeoSeries(GeoPandasBase, Series):
         GeoDataFrame.explode
 
         """
+        if index_parts is None and not ignore_index:
+            warnings.warn(
+                "Currently, index_parts defaults to True, but in the future, "
+                "it will default to False to be consistent with Pandas. "
+                "Use `index_parts=True` to keep the current behavior and True/False "
+                "to silence the warning.",
+                FutureWarning,
+                stacklevel=2,
+            )
+            index_parts = True
 
         if compat.USE_PYGEOS and compat.PYGEOS_GE_09:
             import pygeos  # noqa
@@ -819,14 +840,23 @@ class GeoSeries(GeoPandasBase, Series):
 
             # extract original index values based on integer index
             outer_index = self.index.take(outer_idx)
+            if ignore_index:
+                index = range(len(geometries))
 
-            nlevels = outer_index.nlevels
-            index_arrays = [outer_index.get_level_values(lvl) for lvl in range(nlevels)]
-            index_arrays.append(inner_index)
+            elif index_parts:
+                nlevels = outer_index.nlevels
+                index_arrays = [
+                    outer_index.get_level_values(lvl) for lvl in range(nlevels)
+                ]
+                index_arrays.append(inner_index)
+
+                index = MultiIndex.from_arrays(
+                    index_arrays, names=self.index.names + [None]
+                )
+
+            else:
+                index = outer_index
 
-            index = MultiIndex.from_arrays(
-                index_arrays, names=self.index.names + [None]
-            )
             return GeoSeries(geometries, index=index, crs=self.crs).__finalize__(self)
 
         # else PyGEOS is not available or version <= 0.8
@@ -843,11 +873,18 @@ class GeoSeries(GeoPandasBase, Series):
             index.extend(idxs)
             geometries.extend(geoms)
 
-        # if self.index is a MultiIndex then index is a list of nested tuples
-        if isinstance(self.index, MultiIndex):
-            index = [tuple(outer) + (inner,) for outer, inner in index]
+        if ignore_index:
+            index = range(len(geometries))
+
+        elif index_parts:
+            # if self.index is a MultiIndex then index is a list of nested tuples
+            if isinstance(self.index, MultiIndex):
+                index = [tuple(outer) + (inner,) for outer, inner in index]
+            index = MultiIndex.from_tuples(index, names=self.index.names + [None])
+
+        else:
+            index = [idx for idx, _ in index]
 
-        index = MultiIndex.from_tuples(index, names=self.index.names + [None])
         return GeoSeries(geometries, index=index, crs=self.crs).__finalize__(self)
 
     #
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index ca516f81..9a57d0ee 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -328,7 +328,7 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
             geom_col = result._geometry_column_name
             collections = result[[geom_col]][is_collection]
 
-            exploded = collections.reset_index(drop=True).explode()
+            exploded = collections.reset_index(drop=True).explode(index_parts=True)
             exploded = exploded.reset_index(level=0)
 
             orig_num_geoms_exploded = exploded.shape[0]
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index 893a403a..d72cd63d 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -127,9 +127,7 @@ stria    AUT    416600.0
     return joined
 
 
-def _basic_checks(
-    left_df, right_df, how, lsuffix, rsuffix, allowed_hows=("left", "right", "inner")
-):
+def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
     """Checks the validity of join input parameters.
 
     `how` must be one of the valid options.
@@ -140,7 +138,7 @@ def _basic_checks(
     ------------
     left_df : GeoDataFrame
     right_df : GeoData Frame
-    how : str, one of allowed_hows
+    how : str, one of 'left', 'right', 'inner'
         join type
     lsuffix : str
         left index suffix
@@ -157,6 +155,7 @@ def _basic_checks(
             "'right_df' should be GeoDataFrame, got {}".format(type(right_df))
         )
 
+    allowed_hows = ["left", "right", "inner"]
     if how not in allowed_hows:
         raise ValueError(
             '`how` was "{}" but is expected to be in {}'.format(how, allowed_hows)
@@ -403,7 +402,7 @@ def _nearest_query(
 def sjoin_nearest(
     left_df: GeoDataFrame,
     right_df: GeoDataFrame,
-    how: str = "left",
+    how: str = "inner",
     max_distance: Optional[float] = None,
     lsuffix: str = "left",
     rsuffix: str = "right",
@@ -422,11 +421,13 @@ def sjoin_nearest(
     Parameters
     ----------
     left_df, right_df : GeoDataFrames
-    how : string, default 'left'
+    how : string, default 'inner'
         The type of join:
 
         * 'left': use keys from left_df; retain only left_df geometry column
         * 'right': use keys from right_df; retain only right_df geometry column
+        * 'inner': use intersection of keys from both dfs; retain only
+          left_df geometry column
     max_distance : float, default None
         Maximum distance within which to query for nearest geometry.
         Must be greater than 0.
@@ -503,9 +504,7 @@ countries_w_city_data[countries_w_city_data["name_left"] == "Italy"]
     Every operation in GeoPandas is planar, i.e. the potential third
     dimension is not taken into account.
     """
-    _basic_checks(
-        left_df, right_df, how, lsuffix, rsuffix, allowed_hows=("left", "right")
-    )
+    _basic_checks(left_df, right_df, how, lsuffix, rsuffix)
 
     left_df.geometry.values.check_geographic_crs(stacklevel=1)
     right_df.geometry.values.check_geographic_crs(stacklevel=1)


Test patch:
diff --git a/geopandas/tests/test_explore.py b/geopandas/tests/test_explore.py
index 5bb9c319..b72e6188 100644
--- a/geopandas/tests/test_explore.py
+++ b/geopandas/tests/test_explore.py
@@ -613,7 +613,7 @@ class TestExplore:
         assert '"maxNativeZoom":19,"maxZoom":19,"minZoom":0' in out_str
 
     def test_linearrings(self):
-        rings = self.nybb.explode().exterior
+        rings = self.nybb.explode(index_parts=True).exterior
         m = rings.explore()
         out_str = self._fetch_map_string(m)
 
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 895693b3..1c70e855 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -2,7 +2,7 @@ import string
 
 import numpy as np
 from numpy.testing import assert_array_equal
-from pandas import DataFrame, MultiIndex, Series
+from pandas import DataFrame, Index, MultiIndex, Series
 
 from shapely.geometry import LinearRing, LineString, MultiPoint, Point, Polygon
 from shapely.geometry.collection import GeometryCollection
@@ -104,7 +104,7 @@ class TestGeomMethods:
         )
 
     def _test_unary_real(self, op, expected, a):
-        """ Tests for 'area', 'length', 'is_valid', etc. """
+        """Tests for 'area', 'length', 'is_valid', etc."""
         fcmp = assert_series_equal
         self._test_unary(op, expected, a, fcmp)
 
@@ -119,7 +119,7 @@ class TestGeomMethods:
         self._test_unary(op, expected, a, fcmp)
 
     def _test_binary_topological(self, op, expected, a, b, *args, **kwargs):
-        """ Tests for 'intersection', 'union', 'symmetric_difference', etc. """
+        """Tests for 'intersection', 'union', 'symmetric_difference', etc."""
         if isinstance(expected, GeoPandasBase):
             fcmp = assert_geoseries_equal
         else:
@@ -861,7 +861,8 @@ class TestGeomMethods:
             index=MultiIndex.from_tuples(index, names=expected_index_name),
             crs=4326,
         )
-        assert_geoseries_equal(expected, s.explode())
+        with pytest.warns(FutureWarning, match="Currently, index_parts defaults"):
+            assert_geoseries_equal(expected, s.explode())
 
     @pytest.mark.parametrize("index_name", [None, "test"])
     def test_explode_geodataframe(self, index_name):
@@ -869,7 +870,8 @@ class TestGeomMethods:
         df = GeoDataFrame({"col": [1, 2], "geometry": s})
         df.index.name = index_name
 
-        test_df = df.explode()
+        with pytest.warns(FutureWarning, match="Currently, index_parts defaults"):
+            test_df = df.explode()
 
         expected_s = GeoSeries([Point(1, 2), Point(2, 3), Point(5, 5)])
         expected_df = GeoDataFrame({"col": [1, 1, 2], "geometry": expected_s})
@@ -888,7 +890,7 @@ class TestGeomMethods:
         df = GeoDataFrame({"level_1": [1, 2], "geometry": s})
         df.index.name = index_name
 
-        test_df = df.explode()
+        test_df = df.explode(index_parts=True)
 
         expected_s = GeoSeries([Point(1, 2), Point(2, 3), Point(5, 5)])
         expected_df = GeoDataFrame({"level_1": [1, 1, 2], "geometry": expected_s})
@@ -900,13 +902,26 @@ class TestGeomMethods:
         expected_df = expected_df.set_index(expected_index)
         assert_frame_equal(test_df, expected_df)
 
+    @pytest.mark.parametrize("index_name", [None, "test"])
+    def test_explode_geodataframe_no_multiindex(self, index_name):
+        # GH1393
+        s = GeoSeries([MultiPoint([Point(1, 2), Point(2, 3)]), Point(5, 5)])
+        df = GeoDataFrame({"level_1": [1, 2], "geometry": s})
+        df.index.name = index_name
+
+        test_df = df.explode(index_parts=False)
+
+        expected_s = GeoSeries([Point(1, 2), Point(2, 3), Point(5, 5)])
+        expected_df = GeoDataFrame({"level_1": [1, 1, 2], "geometry": expected_s})
+
+        expected_index = Index([0, 0, 1], name=index_name)
+        expected_df = expected_df.set_index(expected_index)
+        assert_frame_equal(test_df, expected_df)
+
     def test_explode_pandas_fallback(self):
         d = {
             "col1": [["name1", "name2"], ["name3", "name4"]],
-            "geometry": [
-                MultiPoint([(1, 2), (3, 4)]),
-                MultiPoint([(2, 1), (0, 0)]),
-            ],
+            "geometry": [MultiPoint([(1, 2), (3, 4)]), MultiPoint([(2, 1), (0, 0)])],
         }
         gdf = GeoDataFrame(d, crs=4326)
         expected_df = GeoDataFrame(
@@ -938,10 +953,7 @@ class TestGeomMethods:
     def test_explode_pandas_fallback_ignore_index(self):
         d = {
             "col1": [["name1", "name2"], ["name3", "name4"]],
-            "geometry": [
-                MultiPoint([(1, 2), (3, 4)]),
-                MultiPoint([(2, 1), (0, 0)]),
-            ],
+            "geometry": [MultiPoint([(1, 2), (3, 4)]), MultiPoint([(2, 1), (0, 0)])],
         }
         gdf = GeoDataFrame(d, crs=4326)
         expected_df = GeoDataFrame(
@@ -977,7 +989,7 @@ class TestGeomMethods:
             index=index,
         )
 
-        test_df = df.explode()
+        test_df = df.explode(index_parts=True)
 
         expected_s = GeoSeries(
             [
@@ -1000,6 +1012,78 @@ class TestGeomMethods:
         expected_df = expected_df.set_index(expected_index)
         assert_frame_equal(test_df, expected_df)
 
+    @pytest.mark.parametrize("outer_index", [1, (1, 2), "1"])
+    def test_explode_pandas_multi_index_false(self, outer_index):
+        index = MultiIndex.from_arrays(
+            [[outer_index, outer_index, outer_index], [1, 2, 3]],
+            names=("first", "second"),
+        )
+        df = GeoDataFrame(
+            {"vals": [1, 2, 3]},
+            geometry=[MultiPoint([(x, x), (x, 0)]) for x in range(3)],
+            index=index,
+        )
+
+        test_df = df.explode(index_parts=False)
+
+        expected_s = GeoSeries(
+            [
+                Point(0, 0),
+                Point(0, 0),
+                Point(1, 1),
+                Point(1, 0),
+                Point(2, 2),
+                Point(2, 0),
+            ]
+        )
+        expected_df = GeoDataFrame({"vals": [1, 1, 2, 2, 3, 3], "geometry": expected_s})
+        expected_index = MultiIndex.from_tuples(
+            [
+                (outer_index, 1),
+                (outer_index, 1),
+                (outer_index, 2),
+                (outer_index, 2),
+                (outer_index, 3),
+                (outer_index, 3),
+            ],
+            names=["first", "second"],
+        )
+        expected_df = expected_df.set_index(expected_index)
+        assert_frame_equal(test_df, expected_df)
+
+    @pytest.mark.parametrize("outer_index", [1, (1, 2), "1"])
+    def test_explode_pandas_multi_index_ignore_index(self, outer_index):
+        index = MultiIndex.from_arrays(
+            [[outer_index, outer_index, outer_index], [1, 2, 3]],
+            names=("first", "second"),
+        )
+        df = GeoDataFrame(
+            {"vals": [1, 2, 3]},
+            geometry=[MultiPoint([(x, x), (x, 0)]) for x in range(3)],
+            index=index,
+        )
+
+        test_df = df.explode(ignore_index=True)
+
+        expected_s = GeoSeries(
+            [
+                Point(0, 0),
+                Point(0, 0),
+                Point(1, 1),
+                Point(1, 0),
+                Point(2, 2),
+                Point(2, 0),
+            ]
+        )
+        expected_df = GeoDataFrame({"vals": [1, 1, 2, 2, 3, 3], "geometry": expected_s})
+        expected_index = Index(range(len(expected_df)))
+        expected_df = expected_df.set_index(expected_index)
+        assert_frame_equal(test_df, expected_df)
+
+        # index_parts is ignored if ignore_index=True
+        test_df = df.explode(ignore_index=True, index_parts=True)
+        assert_frame_equal(test_df, expected_df)
+
     #
     # Test '&', '|', '^', and '-'
     #
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index 773c718a..2444bf25 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -7,6 +7,7 @@ import tempfile
 import numpy as np
 from numpy.testing import assert_array_equal
 import pandas as pd
+from pandas.util.testing import assert_index_equal
 
 from pyproj import CRS
 from shapely.geometry import (
@@ -172,7 +173,7 @@ class TestSeries:
         assert_series_equal(res, exp)
 
     def test_to_file(self):
-        """ Test to_file and from_file """
+        """Test to_file and from_file"""
         tempfilename = os.path.join(self.tempdir, "test.shp")
         self.g3.to_file(tempfilename)
         # Read layer back in?
@@ -239,7 +240,7 @@ class TestSeries:
         # self.na_none.fillna(method='backfill')
 
     def test_coord_slice(self):
-        """ Test CoordinateSlicer """
+        """Test CoordinateSlicer"""
         # need some better test cases
         assert geom_equals(self.g3, self.g3.cx[:, :])
         assert geom_equals(self.g3[[True, False]], self.g3.cx[0.9:, :0.1])
@@ -492,6 +493,27 @@ class TestConstructor:
         s = GeoSeries(
             [MultiPoint([(0, 0), (1, 1)]), MultiPoint([(2, 2), (3, 3), (4, 4)])]
         )
-        s = s.explode()
+        s = s.explode(index_parts=True)
         df = s.reset_index()
         assert type(df) == GeoDataFrame
+
+    def test_explode_without_multiindex(self):
+        s = GeoSeries(
+            [MultiPoint([(0, 0), (1, 1)]), MultiPoint([(2, 2), (3, 3), (4, 4)])]
+        )
+        s = s.explode(index_parts=False)
+        expected_index = pd.Index([0, 0, 1, 1, 1])
+        assert_index_equal(s.index, expected_index)
+
+    def test_explode_ignore_index(self):
+        s = GeoSeries(
+            [MultiPoint([(0, 0), (1, 1)]), MultiPoint([(2, 2), (3, 3), (4, 4)])]
+        )
+        s = s.explode(ignore_index=True)
+        expected_index = pd.Index(range(len(s)))
+        print(expected_index)
+        assert_index_equal(s.index, expected_index)
+
+        # index_parts is ignored if ignore_index=True
+        s = s.explode(index_parts=True, ignore_index=True)
+        assert_index_equal(s.index, expected_index)
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 8bff9e2e..4d18bad2 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -619,7 +619,7 @@ def test_preserve_attrs(df):
     assert df2.attrs == attrs
 
     # https://github.com/geopandas/geopandas/issues/1875
-    df3 = df2.explode()
+    df3 = df2.explode(index_parts=True)
     assert df3.attrs == attrs
 
 
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index 48224886..7dde9083 100644
--- a/geopandas/tools/tests/test_sjoin.py
+++ b/geopandas/tools/tests/test_sjoin.py
@@ -585,13 +585,15 @@ def test_no_nearest_all():
     ),
 )
 class TestNearest:
-    @pytest.mark.parametrize("how_kwargs", ({}, {"how": "left"}, {"how": "right"}))
+    @pytest.mark.parametrize(
+        "how_kwargs", ({}, {"how": "inner"}, {"how": "left"}, {"how": "right"})
+    )
     def test_allowed_hows(self, how_kwargs):
         left = geopandas.GeoDataFrame({"geometry": []})
         right = geopandas.GeoDataFrame({"geometry": []})
         sjoin_nearest(left, right, **how_kwargs)  # no error
 
-    @pytest.mark.parametrize("how", ("inner", "abcde"))
+    @pytest.mark.parametrize("how", ("outer", "abcde"))
     def test_invalid_hows(self, how: str):
         left = geopandas.GeoDataFrame({"geometry": []})
         right = geopandas.GeoDataFrame({"geometry": []})
@@ -629,16 +631,12 @@ class TestNearest:
         if distance_col is not None:
             assert distance_col in joined
 
+    @pytest.mark.parametrize("how", ["inner", "left"])
     @pytest.mark.parametrize("distance_col", (None, "distance"))
-    def test_empty_left_df_how_left(self, distance_col: str):
+    def test_empty_left_df_how_left(self, how, distance_col: str):
         right = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
         left = geopandas.GeoDataFrame({"geometry": []})
-        joined = sjoin_nearest(
-            left,
-            right,
-            how="left",
-            distance_col=distance_col,
-        )
+        joined = sjoin_nearest(left, right, how=how, distance_col=distance_col)
         assert joined.empty
         if distance_col is not None:
             assert distance_col in joined
@@ -658,7 +656,8 @@ class TestNearest:
         if distance_col is not None:
             assert joined[distance_col].isna().all()
 
-    def test_empty_join_due_to_max_distance_how_left(self):
+    @pytest.mark.parametrize("how", ["inner", "left"])
+    def test_empty_join_due_to_max_distance_how_left(self, how):
         # after applying max_distance the join comes back empty
         # (as in NaN in the joined columns)
         left = geopandas.GeoDataFrame({"geometry": [Point(0, 0)]})
@@ -666,13 +665,16 @@ class TestNearest:
         joined = sjoin_nearest(
             left,
             right,
-            how="left",
+            how=how,
             max_distance=1,
             distance_col="distances",
         )
         expected = left.copy()
         expected["index_right"] = [np.nan]
         expected["distances"] = [np.nan]
+        if how == "inner":
+            expected = expected.dropna()
+            expected["index_right"] = expected["index_right"].astype("int64")
         assert_geodataframe_equal(joined, expected)
 
     def test_empty_join_due_to_max_distance_how_right(self):
@@ -693,19 +695,23 @@ class TestNearest:
         expected = expected[["index_left", "geometry", "distances"]]
         assert_geodataframe_equal(joined, expected)
 
-    def test_max_distance_how_left(self):
+    @pytest.mark.parametrize("how", ["inner", "left"])
+    def test_max_distance_how_left(self, how):
         left = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
         right = geopandas.GeoDataFrame({"geometry": [Point(1, 1), Point(2, 2)]})
         joined = sjoin_nearest(
             left,
             right,
-            how="left",
+            how=how,
             max_distance=1,
             distance_col="distances",
         )
         expected = left.copy()
         expected["index_right"] = [np.nan, 0]
         expected["distances"] = [np.nan, 0]
+        if how == "inner":
+            expected = expected.dropna()
+            expected["index_right"] = expected["index_right"].astype("int64")
         assert_geodataframe_equal(joined, expected)
 
     def test_max_distance_how_right(self):
@@ -724,6 +730,7 @@ class TestNearest:
         expected = expected[["index_left", "geometry", "distances"]]
         assert_geodataframe_equal(joined, expected)
 
+    @pytest.mark.parametrize("how", ["inner", "left"])
     @pytest.mark.parametrize(
         "geo_left, geo_right, expected_left, expected_right, distances",
         [
@@ -792,18 +799,21 @@ class TestNearest:
         expected_left: Sequence[int],
         expected_right: Sequence[int],
         distances: Sequence[float],
+        how,
     ):
         left = geopandas.GeoDataFrame({"geometry": geo_left})
         right = geopandas.GeoDataFrame({"geometry": geo_right})
-        expected_gdf = left.iloc[expected_left]
+        expected_gdf = left.iloc[expected_left].copy()
         expected_gdf["index_right"] = expected_right
         # without distance col
-        joined = sjoin_nearest(left, right, how="left")
-        assert_geodataframe_equal(expected_gdf, joined)
+        joined = sjoin_nearest(left, right, how=how)
+        # inner / left join give a different row order
+        check_like = how == "inner"
+        assert_geodataframe_equal(expected_gdf, joined, check_like=check_like)
         # with distance col
         expected_gdf["distance_col"] = np.array(distances, dtype=float)
-        joined = sjoin_nearest(left, right, how="left", distance_col="distance_col")
-        assert_geodataframe_equal(expected_gdf, joined)
+        joined = sjoin_nearest(left, right, how=how, distance_col="distance_col")
+        assert_geodataframe_equal(expected_gdf, joined, check_like=check_like)
 
     @pytest.mark.parametrize(
         "geo_left, geo_right, expected_left, expected_right, distances",
@@ -870,7 +880,7 @@ class TestNearest:
     ):
         left = geopandas.GeoDataFrame({"geometry": geo_left})
         right = geopandas.GeoDataFrame({"geometry": geo_right})
-        expected_gdf = right.iloc[expected_right]
+        expected_gdf = right.iloc[expected_right].copy()
         expected_gdf["index_left"] = expected_left
         expected_gdf = expected_gdf[["index_left", "geometry"]]
         # without distance col
@@ -880,3 +890,31 @@ class TestNearest:
         expected_gdf["distance_col"] = np.array(distances, dtype=float)
         joined = sjoin_nearest(left, right, how="right", distance_col="distance_col")
         assert_geodataframe_equal(expected_gdf, joined)
+
+    @pytest.mark.filterwarnings("ignore:Geometry is in a geographic CRS")
+    def test_sjoin_nearest_inner(self):
+        # check equivalency of left and inner join
+        countries = read_file(geopandas.datasets.get_path("naturalearth_lowres"))
+        cities = read_file(geopandas.datasets.get_path("naturalearth_cities"))
+        countries = countries[["geometry", "name"]].rename(columns={"name": "country"})
+
+        # default: inner and left give the same result
+        result1 = sjoin_nearest(cities, countries, distance_col="dist")
+        assert result1.shape[0] == cities.shape[0]
+        result2 = sjoin_nearest(cities, countries, distance_col="dist", how="inner")
+        assert_geodataframe_equal(result2, result1)
+        result3 = sjoin_nearest(cities, countries, distance_col="dist", how="left")
+        assert_geodataframe_equal(result3, result1, check_like=True)
+
+        # with max_distance: rows that go above are dropped in case of inner
+        result4 = sjoin_nearest(cities, countries, distance_col="dist", max_distance=1)
+        assert_geodataframe_equal(
+            result4, result1[result1["dist"] < 1], check_like=True
+        )
+        result5 = sjoin_nearest(
+            cities, countries, distance_col="dist", max_distance=1, how="left"
+        )
+        assert result5.shape[0] == cities.shape[0]
+        result5 = result5.dropna()
+        result5["index_right"] = result5["index_right"].astype("int64")
+        assert_geodataframe_equal(result5, result4, check_like=True)