Repository: geopandas/geopandas

Issue:
API (BUG?): GeoDataFrame(gdf) constructor: preserve geometry column name?
When passing a GeoDataFrame object to the GeoDataFrame object, it currently does not preserve the geometry column name if not "geometry":

```
In [1]: df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))

In [2]: df2 = df.rename_geometry("polygons")

In [3]: res = geopandas.GeoDataFrame(df2)

In [4]: res.geometry
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-4-ec92b07bcf43> in <module>
----> 1 res.geometry

~/scipy/pandas/pandas/core/generic.py in __getattr__(self, name)
   5464         ):
   5465             return self[name]
-> 5466         return object.__getattribute__(self, name)
   5467 
   5468     def __setattr__(self, name: str, value) -> None:

~/scipy/geopandas/geopandas/geodataframe.py in _get_geometry(self)
    196     def _get_geometry(self):
    197         if self._geometry_column_name not in self:
--> 198             raise AttributeError(
    199                 "No geometry data set yet (expected in"
    200                 " column '%s'.)" % self._geometry_column_name

AttributeError: No geometry data set yet (expected in column 'geometry'.)

In [5]: res._geometry_column_name
Out[5]: 'geometry'
```

Of course, this somewhat follows from "geometry" being the default column name for the geometry column if not specified by the user. But, it seems we should preserve the name if passing a GeoDataFrame and `geometry=None` in the constructor?

cc @m-richards 

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/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 916bee75..2eb65713 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -117,6 +117,7 @@ Serialization / IO / conversion
    GeoSeries.from_file
    GeoSeries.from_wkb
    GeoSeries.from_wkt
+   GeoSeries.from_xy
    GeoSeries.to_file
    GeoSeries.to_json
    GeoSeries.to_wkb
diff --git a/doc/source/docs/reference/sindex.rst b/doc/source/docs/reference/sindex.rst
index 3201106d..5f4f0f79 100644
--- a/doc/source/docs/reference/sindex.rst
+++ b/doc/source/docs/reference/sindex.rst
@@ -30,6 +30,8 @@ methods:
 
     intersection
     is_empty
+    nearest
+    nearest_all
     query
     query_bulk
     size
@@ -42,4 +44,10 @@ 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`_.
 
+Similarly, the ``pygeos``-based spatial index
+(``geopandas.sindex.PyGEOSSTRTreeIndex``) offers the full capability of
+``pygeos.STRtree``, including nearest-neighbor queries.
+See the full API in the `PyGEOS STRTree documentation`_.
+
 .. _rtree documentation: https://rtree.readthedocs.io/en/stable/class.html
+.. _PyGEOS STRTree documentation: https://pygeos.readthedocs.io/en/latest/strtree.html
diff --git a/doc/source/docs/reference/tools.rst b/doc/source/docs/reference/tools.rst
index ee31e4f8..618fe99b 100644
--- a/doc/source/docs/reference/tools.rst
+++ b/doc/source/docs/reference/tools.rst
@@ -7,6 +7,7 @@ Tools
    :toctree: api/
 
    sjoin
+   sjoin_nearest
    overlay
    clip
    tools.geocode
diff --git a/doc/source/docs/user_guide/mergingdata.rst b/doc/source/docs/user_guide/mergingdata.rst
index 2343571f..72e2ae12 100644
--- a/doc/source/docs/user_guide/mergingdata.rst
+++ b/doc/source/docs/user_guide/mergingdata.rst
@@ -95,8 +95,15 @@ In a Spatial Join, two geometry objects are merged based on their spatial relati
    cities_with_country.head()
 
 
-Sjoin Arguments
-~~~~~~~~~~~~~~~~
+GeoPandas provides two spatial-join functions:
+
+- :func:`sjoin`: joins based on binary predicates (intersects, contains, etc.)
+- :func:`sjoin_nearest`: joins based on proximity, with the ability to set a maximum search radius.
+
+Binary Predicate Joins
+~~~~~~~~~~~~~~~~~~~~~~
+
+Binary predicate joins are available via :func:`sjoin`.
 
 :func:`sjoin` has two core arguments: ``how`` and ``predicate``.
 
@@ -132,3 +139,20 @@ The `how` argument specifies the type of join that will occur and which geometry
 Note more complicated spatial relationships can be studied by combining geometric operations with spatial join.
 To find all polygons within a given distance of a point, for example, one can first use the :meth:`~geopandas.GeoSeries.buffer` method to expand each
 point into a circle of appropriate radius, then intersect those buffered circles with the polygons in question.
+
+Nearest Joins
+~~~~~~~~~~~~~
+
+Proximity-based joins can be done via :func:`sjoin_nearest`.
+
+:func:`sjoin_nearest` shares the ``how`` argument with :func:`sjoin`, and
+includes two additional arguments: ``max_distance`` and ``distance_col``.
+
+**max_distance**
+
+The ``max_distance`` argument specifies a maximum search radius for matching geometries. This can have a considerable performance impact in some cases.
+If you can, it is highly recommended that you use this parameter.
+
+**distance_col**
+
+If set, the resultant GeoDataFrame will include a column with this name containing the computed distances between an input geometry and the nearest geometry.
diff --git a/geopandas/__init__.py b/geopandas/__init__.py
index 79bffdd7..63345b32 100644
--- a/geopandas/__init__.py
+++ b/geopandas/__init__.py
@@ -8,7 +8,7 @@ from geopandas.io.file import _read_file as read_file  # noqa
 from geopandas.io.arrow import _read_parquet as read_parquet  # noqa
 from geopandas.io.arrow import _read_feather as read_feather  # noqa
 from geopandas.io.sql import _read_postgis as read_postgis  # noqa
-from geopandas.tools import sjoin  # noqa
+from geopandas.tools import sjoin, sjoin_nearest  # noqa
 from geopandas.tools import overlay  # noqa
 from geopandas.tools._show_versions import show_versions  # noqa
 from geopandas.tools import clip  # noqa
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index dbe5311b..e5ac045d 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -38,6 +38,10 @@ USE_PYGEOS = None
 PYGEOS_SHAPELY_COMPAT = None
 
 PYGEOS_GE_09 = None
+PYGEOS_GE_010 = None
+
+INSTALL_PYGEOS_ERROR = "To use PyGEOS within GeoPandas, you need to install PyGEOS: \
+'conda install pygeos' or 'pip install pygeos'"
 
 try:
     import pygeos  # noqa
@@ -46,6 +50,7 @@ try:
     if str(pygeos.__version__) >= LooseVersion("0.8"):
         HAS_PYGEOS = True
         PYGEOS_GE_09 = str(pygeos.__version__) >= LooseVersion("0.9")
+        PYGEOS_GE_010 = str(pygeos.__version__) >= LooseVersion("0.10")
     else:
         warnings.warn(
             "The installed version of PyGEOS is too old ({0} installed, 0.8 required),"
@@ -89,7 +94,7 @@ def set_use_pygeos(val=None):
             # validate the pygeos version
             if not str(pygeos.__version__) >= LooseVersion("0.8"):
                 raise ImportError(
-                    "PyGEOS >= 0.6 is required, version {0} is installed".format(
+                    "PyGEOS >= 0.8 is required, version {0} is installed".format(
                         pygeos.__version__
                     )
                 )
@@ -115,10 +120,7 @@ def set_use_pygeos(val=None):
                 PYGEOS_SHAPELY_COMPAT = True
 
         except ImportError:
-            raise ImportError(
-                "To use the PyGEOS speed-ups within GeoPandas, you need to install "
-                "PyGEOS: 'conda install pygeos' or 'pip install pygeos'"
-            )
+            raise ImportError(INSTALL_PYGEOS_ERROR)
 
 
 set_use_pygeos()
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 2a9b9f4a..6b57e037 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -47,6 +47,20 @@ def _ensure_geometry(data, crs=None):
             return out
 
 
+def _crs_mismatch_warning():
+    # TODO: raise error in 0.9 or 0.10.
+    warnings.warn(
+        "CRS mismatch between CRS of the passed geometries "
+        "and 'crs'. Use 'GeoDataFrame.set_crs(crs, "
+        "allow_override=True)' to overwrite CRS or "
+        "'GeoDataFrame.to_crs(crs)' to reproject geometries. "
+        "CRS mismatch will raise an error in the future versions "
+        "of GeoPandas.",
+        FutureWarning,
+        stacklevel=3,
+    )
+
+
 class GeoDataFrame(GeoPandasBase, DataFrame):
     """
     A GeoDataFrame object is a pandas.DataFrame that has a column
@@ -103,9 +117,9 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
     _geometry_column_name = DEFAULT_GEO_COLUMN_NAME
 
-    def __init__(self, *args, geometry=None, crs=None, **kwargs):
+    def __init__(self, data=None, *args, geometry=None, crs=None, **kwargs):
         with compat.ignore_shapely2_warnings():
-            super().__init__(*args, **kwargs)
+            super().__init__(data, *args, **kwargs)
 
         # need to set this before calling self['geometry'], because
         # getitem accesses crs
@@ -117,6 +131,15 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         # but within a try/except because currently non-geometries are
         # allowed in that case
         # TODO do we want to raise / return normal DataFrame in this case?
+
+        # if gdf passed in and geo_col is set, we use that for geometry
+        if geometry is None and isinstance(data, GeoDataFrame):
+            self._geometry_column_name = data._geometry_column_name
+            if crs is not None and data.crs != crs:
+                _crs_mismatch_warning()
+                # TODO: raise error in 0.9 or 0.10.
+            return
+
         if geometry is None and "geometry" in self.columns:
             # Check for multiple columns with name "geometry". If there are,
             # self["geometry"] is a gdf and constructor gets recursively recalled
@@ -136,16 +159,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
                     and crs
                     and not self["geometry"].values.crs == crs
                 ):
-                    warnings.warn(
-                        "CRS mismatch between CRS of the passed geometries "
-                        "and 'crs'. Use 'GeoDataFrame.set_crs(crs, "
-                        "allow_override=True)' to overwrite CRS or "
-                        "'GeoDataFrame.to_crs(crs)' to reproject geometries. "
-                        "CRS mismatch will raise an error in the future versions "
-                        "of GeoPandas.",
-                        FutureWarning,
-                        stacklevel=2,
-                    )
+                    _crs_mismatch_warning()
                     # TODO: raise error in 0.9 or 0.10.
                 self["geometry"] = _ensure_geometry(self["geometry"].values, crs)
             except TypeError:
@@ -165,16 +179,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
                 and crs
                 and not geometry.crs == crs
             ):
-                warnings.warn(
-                    "CRS mismatch between CRS of the passed geometries "
-                    "and 'crs'. Use 'GeoDataFrame.set_crs(crs, "
-                    "allow_override=True)' to overwrite CRS or "
-                    "'GeoDataFrame.to_crs(crs)' to reproject geometries. "
-                    "CRS mismatch will raise an error in the future versions "
-                    "of GeoPandas.",
-                    FutureWarning,
-                    stacklevel=2,
-                )
+                _crs_mismatch_warning()
                 # TODO: raise error in 0.9 or 0.10.
             self.set_geometry(geometry, inplace=True)
 
@@ -1546,7 +1551,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 +1559,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 +1599,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 +1607,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 +1634,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 +1674,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..e02b6abf 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -20,6 +20,7 @@ from .array import (
     from_shapely,
     from_wkb,
     from_wkt,
+    points_from_xy,
     to_wkb,
     to_wkt,
 )
@@ -436,6 +437,52 @@ class GeoSeries(GeoPandasBase, Series):
         """
         return cls._from_wkb_or_wkb(from_wkt, data, index=index, crs=crs, **kwargs)
 
+    @classmethod
+    def from_xy(cls, x, y, z=None, index=None, crs=None, **kwargs):
+        """
+        Alternate constructor to create a :class:`~geopandas.GeoSeries` of Point
+        geometries from lists or arrays of x, y(, z) coordinates
+
+        In case of geographic coordinates, it is assumed that longitude is captured
+        by ``x`` coordinates and latitude by ``y``.
+
+        Parameters
+        ----------
+        x, y, z : iterable
+        index : array-like or Index, optional
+            The index for the GeoSeries.
+        crs : value, optional
+            Coordinate Reference System of the geometry objects. Can be anything
+            accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:4326") or a WKT string.
+        **kwargs
+            Additional arguments passed to the Series constructor,
+            e.g. ``name``.
+
+        Returns
+        -------
+        GeoSeries
+
+        See Also
+        --------
+        GeoSeries.from_wkt
+        points_from_xy
+
+        Examples
+        --------
+
+        >>> x = [2.5, 5, -3.0]
+        >>> y = [0.5, 1, 1.5]
+        >>> s = geopandas.GeoSeries.from_xy(x, y, crs="EPSG:4326")
+        >>> s
+        0    POINT (2.50000 0.50000)
+        1    POINT (5.00000 1.00000)
+        2    POINT (-3.00000 1.50000)
+        dtype: geometry
+        """
+        return cls(points_from_xy(x, y, z, crs=crs), index=index, crs=crs, **kwargs)
+
     @classmethod
     def _from_wkb_or_wkb(
         cls, from_wkb_or_wkt_function, data, index=None, crs=None, **kwargs
@@ -757,7 +804,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 +812,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 +840,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 +853,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 +887,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 +920,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/io/file.py b/geopandas/io/file.py
index 27095bcd..59467818 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -14,17 +14,22 @@ try:
     import fiona
 
     fiona_import_error = None
+
+    # only try to import fiona.Env if the main fiona import succeeded (otherwise you
+    # can get confusing "AttributeError: module 'fiona' has no attribute '_loading'"
+    # / partially initialized module errors)
+    try:
+        from fiona import Env as fiona_env
+    except ImportError:
+        try:
+            from fiona import drivers as fiona_env
+        except ImportError:
+            fiona_env = None
+
 except ImportError as err:
     fiona = None
     fiona_import_error = str(err)
 
-try:
-    from fiona import Env as fiona_env
-except ImportError:
-    try:
-        from fiona import drivers as fiona_env
-    except ImportError:
-        fiona_env = None
 
 from geopandas import GeoDataFrame, GeoSeries
 
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 9970627c..5526372a 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -168,6 +168,110 @@ class BaseSpatialIndex:
         """
         raise NotImplementedError
 
+    def nearest(self, geometry):
+        """
+        Returns the nearest geometry to the input geometries.
+
+        .. note::
+            ``nearest`` currently only works with PyGEOS >= 0.10.
+
+        Note that if pygeos is not available, geopandas
+        will use rtree for the spatial index, where nearest temporarily has a different
+        function signature to temporarily preserve existing functionality. See the
+        documentation of :meth:`rtree.index.Index.nearest` for the details on the
+        ``rtree``-based implementation.
+
+        Parameters
+        ----------
+        geometry : {shapely.geometry, GeoSeries, GeometryArray, numpy.array of PyGEOS \
+geometries}
+            A single shapely geometry, one of the GeoPandas geometry iterables
+            (GeoSeries, GeometryArray), or a numpy array of PyGEOS geometries to query
+            against the spatial index.
+
+        Returns
+        -------
+        matches : 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.head()
+        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)
+        dtype: geometry
+
+        >>> s.sindex.nearest(Point(1, 1))
+        array([[0],
+               [1]])
+
+        >>> s.sindex.nearest([box(4.9, 4.9, 5.1, 5.1)])
+        array([[0],
+               [5]])
+
+        >>> s2 = geopandas.GeoSeries(geopandas.points_from_xy([7.6, 10], [7.6, 10]))
+        >>> s2
+        0    POINT (7.60000 7.60000)
+        1    POINT (10.00000 10.00000)
+        dtype: geometry
+
+        >>> s.sindex.nearest(s2)
+        array([[0, 1],
+               [8, 9]])
+        """
+        raise NotImplementedError
+
+    def nearest_all(self, geometry, max_distance=None, return_distance=False):
+        """Return the nearest geometry in the tree for each input geometry in
+        ``geometry``.
+
+        If multiple tree geometries have the same distance from an input geometry,
+        multiple results will be returned for that input geometry.
+
+        .. note::
+            ``nearest_all`` currently only works with PyGEOS >= 0.10.
+
+        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.
+        If max_distance is not set, this will effectively be a left join
+        because every geometry in ``geometry`` will have a nearest geometry in
+        the tree. However, if ``max_distance`` is used, this becomes an
+        inner join, since some geometries in ``geometry`` may not have a match
+        in the tree.
+
+        For performance reasons, it is highly recommended that you set
+        the ``max_distance`` parameter.
+
+        Parameters
+        ----------
+        geometry : {shapely.geometry, GeoSeries, GeometryArray, numpy.array of PyGEOS \
+geometries}
+            A single shapely geometry, one of the GeoPandas geometry iterables
+            (GeoSeries, GeometryArray), or a numpy array of PyGEOS geometries to query
+            against the spatial index.
+        max_distance : float, optional
+            Maximum distance within which to query for nearest items in tree.
+            Must be greater than 0. By default None, indicating no distance limit.
+        return_distance : bool, optional
+            If True, will return distances in addition to indexes. By default False
+
+        Returns
+        -------
+        Indices or tuple of (indices, distances)
+            Indices is an ndarray of shape (2,n) and distances (if present) an
+            ndarray of shape (n).
+            The first subarray of indices contains input geometry indices.
+            The second subarray of indices contains tree geometry indices.
+        """
+        raise NotImplementedError
+
     def intersection(self, coordinates):
         """Compatibility wrapper for rtree.index.Index.intersection,
         use ``query`` instead.
@@ -288,6 +392,10 @@ if compat.HAS_RTREE:
         def intersection(self, coordinates, *args, **kwargs):
             return super().intersection(coordinates, *args, **kwargs)
 
+        @doc(BaseSpatialIndex.nearest)
+        def nearest(self, *args, **kwargs):
+            return super().nearest(*args, **kwargs)
+
         @property
         @doc(BaseSpatialIndex.size)
         def size(self):
@@ -438,6 +546,72 @@ if compat.HAS_RTREE:
                 input_geometry_index.extend([i] * len(res))
             return np.vstack([input_geometry_index, tree_index])
 
+        def nearest(self, coordinates, num_results=1, objects=False):
+            """
+            Returns the nearest object or objects to the given coordinates.
+
+            Requires rtree, and passes parameters directly to
+            :meth:`rtree.index.Index.nearest`.
+
+            This behaviour is deprecated and will be updated to be consistent
+            with the pygeos PyGEOSSTRTreeIndex in a future release.
+
+            If longer-term compatibility is required, use
+            :meth:`rtree.index.Index.nearest` directly instead.
+
+            Examples
+            --------
+            >>> s = geopandas.GeoSeries(geopandas.points_from_xy(range(3), range(3)))
+            >>> s
+            0    POINT (0.00000 0.00000)
+            1    POINT (1.00000 1.00000)
+            2    POINT (2.00000 2.00000)
+            dtype: geometry
+
+            >>> list(s.sindex.nearest((0, 0)))  # doctest: +SKIP
+            [0]
+
+            >>> list(s.sindex.nearest((0.5, 0.5)))  # doctest: +SKIP
+            [0, 1]
+
+            >>> list(s.sindex.nearest((3, 3), num_results=2))  # doctest: +SKIP
+            [2, 1]
+
+            >>> list(super(type(s.sindex), s.sindex).nearest((0, 0),
+            ... num_results=2))  # doctest: +SKIP
+            [0, 1]
+
+            Parameters
+            ----------
+            coordinates : sequence or array
+                This may be an object that satisfies the numpy array protocol,
+                providing the index’s dimension * 2 coordinate pairs
+                representing the mink and maxk coordinates in each dimension
+                defining the bounds of the query window.
+            num_results : integer
+                The number of results to return nearest to the given
+                coordinates. If two index entries are equidistant, both are
+                returned. This property means that num_results may return more
+                items than specified
+            objects : True / False / ‘raw’
+                If True, the nearest method will return index objects that were
+                pickled when they were stored with each index entry, as well as
+                the id and bounds of the index entries. If ‘raw’, it will
+                return the object as entered into the database without the
+                rtree.index.Item wrapper.
+            """
+            warnings.warn(
+                "sindex.nearest using the rtree backend was not previously documented "
+                "and this behavior is deprecated in favor of matching the function "
+                "signature provided by the pygeos backend (see "
+                "PyGEOSSTRTreeIndex.nearest for details). This behavior will be "
+                "updated in a future release.",
+                FutureWarning,
+            )
+            return super().nearest(
+                coordinates, num_results=num_results, objects=objects
+            )
+
         @doc(BaseSpatialIndex.intersection)
         def intersection(self, coordinates):
             return super().intersection(coordinates, objects=False)
@@ -533,6 +707,41 @@ if compat.HAS_PYGEOS:
 
             return matches
 
+        @staticmethod
+        def _as_geometry_array(geometry):
+            """Convert geometry into a numpy array of PyGEOS geometries.
+
+            Parameters
+            ----------
+            geometry
+                An array-like of PyGEOS geometries, a GeoPandas GeoSeries/GeometryArray,
+                shapely.geometry or list of shapely geometries.
+
+            Returns
+            -------
+            np.ndarray
+                A numpy array of pygeos geometries.
+            """
+            if isinstance(geometry, np.ndarray):
+                return geometry
+            elif isinstance(geometry, geoseries.GeoSeries):
+                return geometry.values.data
+            elif isinstance(geometry, array.GeometryArray):
+                return geometry.data
+            elif isinstance(geometry, BaseGeometry):
+                return array._shapely_to_geom(geometry)
+            elif isinstance(geometry, list):
+                return np.asarray(
+                    [
+                        array._shapely_to_geom(el)
+                        if isinstance(el, BaseGeometry)
+                        else el
+                        for el in geometry
+                    ]
+                )
+            else:
+                return np.asarray(geometry)
+
         @doc(BaseSpatialIndex.query_bulk)
         def query_bulk(self, geometry, predicate=None, sort=False):
             if predicate not in self.valid_query_predicates:
@@ -541,12 +750,8 @@ if compat.HAS_PYGEOS:
                         predicate, self.valid_query_predicates
                     )
                 )
-            if isinstance(geometry, geoseries.GeoSeries):
-                geometry = geometry.values.data
-            elif isinstance(geometry, array.GeometryArray):
-                geometry = geometry.data
-            elif not isinstance(geometry, np.ndarray):
-                geometry = np.asarray(geometry)
+
+            geometry = self._as_geometry_array(geometry)
 
             res = super().query_bulk(geometry, predicate)
 
@@ -558,6 +763,21 @@ if compat.HAS_PYGEOS:
 
             return res
 
+        @doc(BaseSpatialIndex.nearest)
+        def nearest(self, geometry):
+            if not compat.PYGEOS_GE_010:
+                raise NotImplementedError("sindex.nearest requires pygeos >= 0.10")
+
+            geometry = self._as_geometry_array(geometry)
+            return super().nearest(geometry)
+
+        @doc(BaseSpatialIndex.nearest_all)
+        def nearest_all(self, geometry, max_distance=None, return_distance=False):
+            geometry = self._as_geometry_array(geometry)
+            return super().nearest_all(
+                geometry, max_distance=max_distance, return_distance=return_distance
+            )
+
         @doc(BaseSpatialIndex.intersection)
         def intersection(self, coordinates):
             # convert bounds to geometry
diff --git a/geopandas/tools/__init__.py b/geopandas/tools/__init__.py
index 6c5fa610..82e6e9e5 100644
--- a/geopandas/tools/__init__.py
+++ b/geopandas/tools/__init__.py
@@ -1,7 +1,7 @@
 from .crs import explicit_crs_from_epsg
 from .geocoding import geocode, reverse_geocode
 from .overlay import overlay
-from .sjoin import sjoin
+from .sjoin import sjoin, sjoin_nearest
 from .util import collect
 from .clip import clip
 
@@ -12,5 +12,6 @@ __all__ = [
     "overlay",
     "reverse_geocode",
     "sjoin",
+    "sjoin_nearest",
     "clip",
 ]
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 ca6f921e..d72cd63d 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -1,8 +1,11 @@
+from typing import Optional
 import warnings
 
+import numpy as np
 import pandas as pd
 
 from geopandas import GeoDataFrame
+from geopandas import _compat as compat
 from geopandas.array import _check_crs, _crs_mismatch_warn
 
 
@@ -117,26 +120,6 @@ stria    AUT    416600.0
 
     _basic_checks(left_df, right_df, how, lsuffix, rsuffix)
 
-    box_left_gdf = left_df.total_bounds
-    box_right_gdf = right_df.total_bounds
-
-    if not (
-        (
-            (box_left_gdf[0] <= box_right_gdf[2])
-            and (box_right_gdf[0] <= box_left_gdf[2])
-        )
-        and (
-            (box_left_gdf[1] <= box_right_gdf[3])
-            and (box_right_gdf[1] <= box_left_gdf[3])
-        )
-    ):
-        copy_df = left_df.copy()
-        copy_df["index_left"] = 0
-        copy_df["index_right"] = 0
-        indices = pd.DataFrame(columns=["_key_left", "_key_right"], dtype=float)
-        copy_df = _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix)
-        return copy_df.iloc[:0]
-
     indices = _geom_predicate_query(left_df, right_df, predicate)
 
     joined = _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix)
@@ -249,16 +232,18 @@ def _geom_predicate_query(left_df, right_df, predicate):
     return indices
 
 
-def _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix):
+def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
     """Join the GeoDataFrames at the DataFrame level.
 
     Parameters
     ----------
-    indices : DataFrame
-        Indexes returned by the geometric join.
+    join_df : DataFrame
+        Indices and join data returned by the geometric join.
         Must have columns `_key_left` and `_key_right`
         with integer indices representing the matches
         from `left_df` and `right_df` respectively.
+        Additional columns may be included and will be copied to
+        the resultant GeoDataFrame.
     left_df : GeoDataFrame
     right_df : GeoDataFrame
     lsuffix : string
@@ -307,9 +292,9 @@ def _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix):
 
     # perform join on the dataframes
     if how == "inner":
-        indices = indices.set_index("_key_left")
+        join_df = join_df.set_index("_key_left")
         joined = (
-            left_df.merge(indices, left_index=True, right_index=True)
+            left_df.merge(join_df, left_index=True, right_index=True)
             .merge(
                 right_df.drop(right_df.geometry.name, axis=1),
                 left_on="_key_right",
@@ -325,9 +310,9 @@ def _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix):
             joined.index.name = left_index_name
 
     elif how == "left":
-        indices = indices.set_index("_key_left")
+        join_df = join_df.set_index("_key_left")
         joined = (
-            left_df.merge(indices, left_index=True, right_index=True, how="left")
+            left_df.merge(join_df, left_index=True, right_index=True, how="left")
             .merge(
                 right_df.drop(right_df.geometry.name, axis=1),
                 how="left",
@@ -347,7 +332,7 @@ def _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix):
         joined = (
             left_df.drop(left_df.geometry.name, axis=1)
             .merge(
-                indices.merge(
+                join_df.merge(
                     right_df, left_on="_key_right", right_index=True, how="right"
                 ),
                 left_index=True,
@@ -364,3 +349,179 @@ def _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix):
             joined.index.name = right_index_name
 
     return joined
+
+
+def _nearest_query(
+    left_df: GeoDataFrame,
+    right_df: GeoDataFrame,
+    max_distance: float,
+    how: str,
+    return_distance: bool,
+):
+    if not (compat.PYGEOS_GE_010 and compat.USE_PYGEOS):
+        raise NotImplementedError(
+            "Currently, only PyGEOS >= 0.10.0 supports `nearest_all`. "
+            + compat.INSTALL_PYGEOS_ERROR
+        )
+    # use the opposite of the join direction for the index
+    use_left_as_sindex = how == "right"
+    if use_left_as_sindex:
+        sindex = left_df.sindex
+        query = right_df.geometry
+    else:
+        sindex = right_df.sindex
+        query = left_df.geometry
+    if sindex:
+        res = sindex.nearest_all(
+            query, max_distance=max_distance, return_distance=return_distance
+        )
+        if return_distance:
+            (input_idx, tree_idx), distances = res
+        else:
+            (input_idx, tree_idx) = res
+            distances = None
+        if use_left_as_sindex:
+            l_idx, r_idx = tree_idx, input_idx
+            sort_order = np.argsort(l_idx, kind="stable")
+            l_idx, r_idx = l_idx[sort_order], r_idx[sort_order]
+            if distances is not None:
+                distances = distances[sort_order]
+        else:
+            l_idx, r_idx = input_idx, tree_idx
+        join_df = pd.DataFrame(
+            {"_key_left": l_idx, "_key_right": r_idx, "distances": distances}
+        )
+    else:
+        # when sindex is empty / has no valid geometries
+        join_df = pd.DataFrame(
+            columns=["_key_left", "_key_right", "distances"], dtype=float
+        )
+    return join_df
+
+
+def sjoin_nearest(
+    left_df: GeoDataFrame,
+    right_df: GeoDataFrame,
+    how: str = "inner",
+    max_distance: Optional[float] = None,
+    lsuffix: str = "left",
+    rsuffix: str = "right",
+    distance_col: Optional[str] = None,
+) -> GeoDataFrame:
+    """Spatial join of two GeoDataFrames based on the distance between their geometries.
+
+    Results will include multiple output records for a single input record
+    where there are multiple equidistant nearest or intersected neighbors.
+
+    See the User Guide page
+    https://geopandas.readthedocs.io/en/latest/docs/user_guide/mergingdata.html
+    for more details.
+
+
+    Parameters
+    ----------
+    left_df, right_df : GeoDataFrames
+    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.
+        The max_distance used to search for nearest items in the tree may have a
+        significant impact on performance by reducing the number of input
+        geometries that are evaluated for nearest items in the tree.
+    lsuffix : string, default 'left'
+        Suffix to apply to overlapping column names (left GeoDataFrame).
+    rsuffix : string, default 'right'
+        Suffix to apply to overlapping column names (right GeoDataFrame).
+    distance_col : string, default None
+        If set, save the distances computed between matching geometries under a
+        column of this name in the joined GeoDataFrame.
+
+    Examples
+    --------
+    >>> countries = geopandas.read_file(geopandas.datasets.get_\
+path("naturalearth_lowres"))
+    >>> cities = geopandas.read_file(geopandas.datasets.get_path("naturalearth_cities"))
+    >>> countries.head(2).name  # doctest: +SKIP
+        pop_est      continent                      name \
+iso_a3  gdp_md_est                                           geometry
+    0     920938        Oceania                      Fiji    FJI      8374.0  MULTIPOLY\
+GON (((180.00000 -16.06713, 180.00000...
+    1   53950935         Africa                  Tanzania    TZA    150600.0  POLYGON (\
+(33.90371 -0.95000, 34.07262 -1.05982...
+    >>> cities.head(2).name  # doctest: +SKIP
+            name                   geometry
+    0  Vatican City  POINT (12.45339 41.90328)
+    1    San Marino  POINT (12.44177 43.93610)
+
+    >>> cities_w_country_data = geopandas.sjoin_nearest(cities, countries)
+    >>> cities_w_country_data[['name_left', 'name_right']].head(2)  # doctest: +SKIP
+            name_left                   geometry  index_right   pop_est continent name_\
+right iso_a3  gdp_md_est
+    0    Vatican City  POINT (12.45339 41.90328)          141  62137802    Europe      \
+Italy    ITA   2221000.0
+    1      San Marino  POINT (12.44177 43.93610)          141  62137802    Europe      \
+Italy    ITA   2221000.0
+
+    To include the distances:
+
+    >>> cities_w_country_data = geopandas.sjoin_nearest\
+(cities, countries, distance_col="distances")
+    >>> cities_w_country_data[["name_left", "name_right", \
+"distances"]].head(2)  # doctest: +SKIP
+            name_left name_right distances
+    0    Vatican City      Italy       0.0
+    1      San Marino      Italy       0.0
+
+    In the following example, we get multiple cities for Italy because all results are
+    equidistant (in this case zero because they intersect).
+    In fact, we get 3 results in total:
+
+    >>> countries_w_city_data = geopandas.sjoin_nearest\
+(cities, countries, distance_col="distances", how="right")
+    >>> italy_results = \
+countries_w_city_data[countries_w_city_data["name_left"] == "Italy"]
+    >>> italy_results  # doctest: +SKIP
+         name_x        name_y
+    141  Vatican City  Italy
+    141    San Marino  Italy
+    141          Rome  Italy
+
+    See also
+    --------
+    sjoin : binary predicate joins
+
+    Notes
+    -----
+    Since this join relies on distances, results will be innaccurate
+    if your geometries are in a geographic CRS.
+
+    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)
+
+    left_df.geometry.values.check_geographic_crs(stacklevel=1)
+    right_df.geometry.values.check_geographic_crs(stacklevel=1)
+
+    return_distance = distance_col is not None
+
+    join_df = _nearest_query(left_df, right_df, max_distance, how, return_distance)
+
+    if return_distance:
+        join_df = join_df.rename(columns={"distances": distance_col})
+    else:
+        join_df.pop("distances")
+
+    joined = _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix)
+
+    if return_distance:
+        columns = [c for c in joined.columns if c != distance_col] + [distance_col]
+        joined = joined[columns]
+
+    return joined


Test patch:
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index c8aa321b..ca237b0c 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -220,16 +220,16 @@ class TestGeometryArrayCRS:
         assert df.geometry.values.crs == self.osgb
 
         # different passed CRS than array CRS is ignored
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             df = GeoDataFrame(geometry=s, crs=4326)
         assert df.crs == self.osgb
         assert df.geometry.crs == self.osgb
         assert df.geometry.values.crs == self.osgb
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             GeoDataFrame(geometry=s, crs=4326)
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             GeoDataFrame({"data": [1, 2], "geometry": s}, crs=4326)
-        with pytest.warns(FutureWarning):
+        with pytest.warns(FutureWarning, match="CRS mismatch"):
             GeoDataFrame(df, crs=4326).crs
 
         # manually change CRS
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_geodataframe.py b/geopandas/tests/test_geodataframe.py
index cbfebbd5..ea98f478 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -930,11 +930,9 @@ class TestConstructor:
             res = GeoDataFrame(df, geometry="other_geom")
             check_geodataframe(res, "other_geom")
 
-        # when passing GeoDataFrame with custom geometry name to constructor
-        # an invalid geodataframe is the result TODO is this desired ?
+        # gdf from gdf should preserve active geometry column name
         df = GeoDataFrame(gpdf)
-        with pytest.raises(AttributeError):
-            df.geometry
+        check_geodataframe(df, "other_geom")
 
     def test_only_geometry(self):
         exp = GeoDataFrame(
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..0c5c63d4 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])
@@ -322,6 +323,22 @@ class TestSeries:
     def test_to_wkt(self):
         assert_series_equal(pd.Series([self.t1.wkt, self.sq.wkt]), self.g1.to_wkt())
 
+    def test_from_xy_points(self):
+        x = self.landmarks.x.values
+        y = self.landmarks.y.values
+        index = self.landmarks.index.tolist()
+        crs = self.landmarks.crs
+        assert_geoseries_equal(
+            self.landmarks, GeoSeries.from_xy(x, y, index=index, crs=crs)
+        )
+
+    def test_from_xy_points_indexless(self):
+        x = np.array([0.0, 3.0])
+        y = np.array([2.0, 5.0])
+        z = np.array([-1.0, 4.0])
+        expected = GeoSeries([Point(0, 2, -1), Point(3, 5, 4)])
+        assert_geoseries_equal(expected, GeoSeries.from_xy(x, y, z))
+
 
 def test_missing_values_empty_warning():
     s = GeoSeries([Point(1, 1), None, np.nan, BaseGeometry(), Polygon()])
@@ -492,6 +509,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/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index d96d6ee6..806c4a9f 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -15,6 +15,9 @@ from geopandas import GeoDataFrame, GeoSeries, read_file, datasets
 import pytest
 import numpy as np
 
+if compat.USE_PYGEOS:
+    import pygeos
+
 
 @pytest.mark.skip_no_sindex
 class TestSeriesSindex:
@@ -666,6 +669,100 @@ class TestPygeosInterface:
                 )
             raise e
 
+    # ------------------------- `nearest` tests ------------------------- #
+    @pytest.mark.skipif(
+        compat.USE_PYGEOS,
+        reason=("RTree supports sindex.nearest with different behaviour"),
+    )
+    def test_rtree_nearest_warns(self):
+        df = geopandas.GeoDataFrame({"geometry": []})
+        with pytest.warns(
+            FutureWarning, match="sindex.nearest using the rtree backend"
+        ):
+            df.sindex.nearest((0, 0, 1, 1), num_results=2)
+
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS and not compat.PYGEOS_GE_010),
+        reason=("PyGEOS < 0.10 does not support sindex.nearest"),
+    )
+    def test_pygeos_error(self):
+        df = geopandas.GeoDataFrame({"geometry": []})
+        with pytest.raises(NotImplementedError, match="requires pygeos >= 0.10"):
+            df.sindex.nearest(None)
+
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS and compat.PYGEOS_GE_010),
+        reason=("PyGEOS >= 0.10 is required to test sindex.nearest"),
+    )
+    @pytest.mark.parametrize(
+        "geometry,expected",
+        [
+            ([0.25, 0.25], [[0], [0]]),
+            ([0.75, 0.75], [[0], [1]]),
+        ],
+    )
+    def test_nearest_single(self, geometry, expected):
+        geoms = pygeos.points(np.arange(10), np.arange(10))
+        df = geopandas.GeoDataFrame({"geometry": geoms})
+
+        p = Point(geometry)
+        res = df.sindex.nearest(p)
+        assert_array_equal(res, expected)
+
+        p = pygeos.points(geometry)
+        res = df.sindex.nearest(p)
+        assert_array_equal(res, expected)
+
+    @pytest.mark.skipif(
+        not compat.USE_PYGEOS or not compat.PYGEOS_GE_010,
+        reason=("PyGEOS >= 0.10 is required to test sindex.nearest"),
+    )
+    @pytest.mark.parametrize(
+        "geometry,expected",
+        [
+            ([(1, 1), (0, 0)], [[0, 1], [1, 0]]),
+            ([(1, 1), (0.25, 1)], [[0, 1], [1, 1]]),
+        ],
+    )
+    def test_nearest_multi(self, geometry, expected):
+        geoms = pygeos.points(np.arange(10), np.arange(10))
+        df = geopandas.GeoDataFrame({"geometry": geoms})
+
+        ps = [Point(p) for p in geometry]
+        res = df.sindex.nearest(ps)
+        assert_array_equal(res, expected)
+
+        ps = pygeos.points(geometry)
+        res = df.sindex.nearest(ps)
+        assert_array_equal(res, expected)
+
+        s = geopandas.GeoSeries(ps)
+        res = df.sindex.nearest(s)
+        assert_array_equal(res, expected)
+
+        x, y = zip(*geometry)
+        ga = geopandas.points_from_xy(x, y)
+        res = df.sindex.nearest(ga)
+        assert_array_equal(res, expected)
+
+    @pytest.mark.skipif(
+        not compat.USE_PYGEOS or not compat.PYGEOS_GE_010,
+        reason=("PyGEOS >= 0.10 is required to test sindex.nearest"),
+    )
+    @pytest.mark.parametrize(
+        "geometry,expected",
+        [
+            (None, [[], []]),
+            ([None], [[], []]),
+        ],
+    )
+    def test_nearest_none(self, geometry, expected):
+        geoms = pygeos.points(np.arange(10), np.arange(10))
+        df = geopandas.GeoDataFrame({"geometry": geoms})
+
+        res = df.sindex.nearest(geometry)
+        assert_array_equal(res, expected)
+
     # --------------------------- misc tests ---------------------------- #
 
     def test_empty_tree_geometries(self):
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index d3f72e18..828da8ce 100644
--- a/geopandas/tools/tests/test_sjoin.py
+++ b/geopandas/tools/tests/test_sjoin.py
@@ -1,4 +1,7 @@
 from distutils.version import LooseVersion
+import math
+from typing import Sequence
+from geopandas.testing import assert_geodataframe_equal
 
 import numpy as np
 import pandas as pd
@@ -6,12 +9,17 @@ import pandas as pd
 from shapely.geometry import Point, Polygon, GeometryCollection
 
 import geopandas
-from geopandas import GeoDataFrame, GeoSeries, read_file, sjoin
+import geopandas._compat as compat
+from geopandas import GeoDataFrame, GeoSeries, read_file, sjoin, sjoin_nearest
+from geopandas.testing import assert_geoseries_equal
 
 from pandas.testing import assert_frame_equal
 import pytest
 
 
+TEST_NEAREST = compat.PYGEOS_GE_010 and compat.USE_PYGEOS
+
+
 pytestmark = pytest.mark.skip_no_sindex
 
 
@@ -552,3 +560,361 @@ class TestSpatialJoinNaturalEarth:
             self.cities, countries, how="inner", predicate="intersects"
         )
         assert cities_with_country.shape == (172, 4)
+
+
+@pytest.mark.skipif(
+    TEST_NEAREST,
+    reason=("This test can only be run _without_ PyGEOS >= 0.10 installed"),
+)
+def test_no_nearest_all():
+    df1 = geopandas.GeoDataFrame({"geometry": []})
+    df2 = geopandas.GeoDataFrame({"geometry": []})
+    with pytest.raises(
+        NotImplementedError,
+        match="Currently, only PyGEOS >= 0.10.0 supports `nearest_all`",
+    ):
+        sjoin_nearest(df1, df2)
+
+
+@pytest.mark.skipif(
+    not TEST_NEAREST,
+    reason=(
+        "PyGEOS >= 0.10.0"
+        " must be installed and activated via the geopandas.compat module to"
+        " test sjoin_nearest"
+    ),
+)
+class TestNearest:
+    @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", ("outer", "abcde"))
+    def test_invalid_hows(self, how: str):
+        left = geopandas.GeoDataFrame({"geometry": []})
+        right = geopandas.GeoDataFrame({"geometry": []})
+        with pytest.raises(ValueError, match="`how` was"):
+            sjoin_nearest(left, right, how=how)
+
+    @pytest.mark.parametrize("distance_col", (None, "distance"))
+    def test_empty_right_df_how_left(self, distance_col: str):
+        # all records from left and no results from right
+        left = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
+        right = geopandas.GeoDataFrame({"geometry": []})
+        joined = sjoin_nearest(
+            left,
+            right,
+            how="left",
+            distance_col=distance_col,
+        )
+        assert_geoseries_equal(joined["geometry"], left["geometry"])
+        assert joined["index_right"].isna().all()
+        if distance_col is not None:
+            assert joined[distance_col].isna().all()
+
+    @pytest.mark.parametrize("distance_col", (None, "distance"))
+    def test_empty_right_df_how_right(self, distance_col: str):
+        # no records in joined
+        left = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
+        right = geopandas.GeoDataFrame({"geometry": []})
+        joined = sjoin_nearest(
+            left,
+            right,
+            how="right",
+            distance_col=distance_col,
+        )
+        assert joined.empty
+        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(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=how, distance_col=distance_col)
+        assert joined.empty
+        if distance_col is not None:
+            assert distance_col in joined
+
+    @pytest.mark.parametrize("distance_col", (None, "distance"))
+    def test_empty_left_df_how_right(self, distance_col: str):
+        right = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
+        left = geopandas.GeoDataFrame({"geometry": []})
+        joined = sjoin_nearest(
+            left,
+            right,
+            how="right",
+            distance_col=distance_col,
+        )
+        assert_geoseries_equal(joined["geometry"], right["geometry"])
+        assert joined["index_left"].isna().all()
+        if distance_col is not None:
+            assert joined[distance_col].isna().all()
+
+    @pytest.mark.parametrize("how", ["inner", "left"])
+    def test_empty_join_due_to_max_distance(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)]})
+        right = geopandas.GeoDataFrame({"geometry": [Point(1, 1), Point(2, 2)]})
+        joined = sjoin_nearest(
+            left,
+            right,
+            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):
+        # after applying max_distance the join comes back empty
+        # (as in NaN in the joined columns)
+        left = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
+        right = geopandas.GeoDataFrame({"geometry": [Point(2, 2)]})
+        joined = sjoin_nearest(
+            left,
+            right,
+            how="right",
+            max_distance=1,
+            distance_col="distances",
+        )
+        expected = right.copy()
+        expected["index_left"] = [np.nan]
+        expected["distances"] = [np.nan]
+        expected = expected[["index_left", "geometry", "distances"]]
+        assert_geodataframe_equal(joined, expected)
+
+    @pytest.mark.parametrize("how", ["inner", "left"])
+    def test_max_distance(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=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):
+        left = geopandas.GeoDataFrame({"geometry": [Point(1, 1), Point(2, 2)]})
+        right = geopandas.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1)]})
+        joined = sjoin_nearest(
+            left,
+            right,
+            how="right",
+            max_distance=1,
+            distance_col="distances",
+        )
+        expected = right.copy()
+        expected["index_left"] = [np.nan, 0]
+        expected["distances"] = [np.nan, 0]
+        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",
+        [
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1)],
+                [0, 1],
+                [0, 0],
+                [math.sqrt(2), 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0, 0)],
+                [0, 1],
+                [1, 0],
+                [0, 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0, 0), Point(0, 0)],
+                [0, 0, 1],
+                [1, 2, 0],
+                [0, 0, 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0, 0), Point(2, 2)],
+                [0, 1],
+                [1, 0],
+                [0, 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0.25, 1)],
+                [0, 1],
+                [1, 0],
+                [math.sqrt(0.25 ** 2 + 1), 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(-10, -10), Point(100, 100)],
+                [0, 1],
+                [0, 0],
+                [math.sqrt(10 ** 2 + 10 ** 2), math.sqrt(11 ** 2 + 11 ** 2)],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(x, y) for x, y in zip(np.arange(10), np.arange(10))],
+                [0, 1],
+                [0, 1],
+                [0, 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1), Point(0, 0)],
+                [Point(1.1, 1.1), Point(0, 0)],
+                [0, 1, 2],
+                [1, 0, 1],
+                [0, np.sqrt(0.1 ** 2 + 0.1 ** 2), 0],
+            ),
+        ],
+    )
+    def test_sjoin_nearest_left(
+        self,
+        geo_left,
+        geo_right,
+        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].copy()
+        expected_gdf["index_right"] = expected_right
+        # without distance col
+        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=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",
+        [
+            ([Point(0, 0), Point(1, 1)], [Point(1, 1)], [1], [0], [0]),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0, 0)],
+                [1, 0],
+                [0, 1],
+                [0, 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0, 0), Point(0, 0)],
+                [1, 0, 0],
+                [0, 1, 2],
+                [0, 0, 0],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0, 0), Point(2, 2)],
+                [1, 0, 1],
+                [0, 1, 2],
+                [0, 0, math.sqrt(2)],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(1, 1), Point(0.25, 1)],
+                [1, 1],
+                [0, 1],
+                [0, 0.75],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(-10, -10), Point(100, 100)],
+                [0, 1],
+                [0, 1],
+                [math.sqrt(10 ** 2 + 10 ** 2), math.sqrt(99 ** 2 + 99 ** 2)],
+            ),
+            (
+                [Point(0, 0), Point(1, 1)],
+                [Point(x, y) for x, y in zip(np.arange(10), np.arange(10))],
+                [0, 1] + [1] * 8,
+                list(range(10)),
+                [0, 0] + [np.sqrt(x ** 2 + x ** 2) for x in np.arange(1, 9)],
+            ),
+            (
+                [Point(0, 0), Point(1, 1), Point(0, 0)],
+                [Point(1.1, 1.1), Point(0, 0)],
+                [1, 0, 2],
+                [0, 1, 1],
+                [np.sqrt(0.1 ** 2 + 0.1 ** 2), 0, 0],
+            ),
+        ],
+    )
+    def test_sjoin_nearest_right(
+        self,
+        geo_left,
+        geo_right,
+        expected_left: Sequence[int],
+        expected_right: Sequence[int],
+        distances: Sequence[float],
+    ):
+        left = geopandas.GeoDataFrame({"geometry": geo_left})
+        right = geopandas.GeoDataFrame({"geometry": geo_right})
+        expected_gdf = right.iloc[expected_right].copy()
+        expected_gdf["index_left"] = expected_left
+        expected_gdf = expected_gdf[["index_left", "geometry"]]
+        # without distance col
+        joined = sjoin_nearest(left, right, how="right")
+        assert_geodataframe_equal(expected_gdf, joined)
+        # with distance col
+        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)