Repository: geopandas/geopandas
Pull request #2135 (2021-09-25 15:14:44)

PERF: avoid calculating total_bounds in sjoin
Closes https://github.com/geopandas/geopandas/issues/2116

This removes it altogether. Alternative is put it behind a `if left_df.has_sindex:` check.


Patch:
diff --git a/doc/source/docs/reference/sindex.rst b/doc/source/docs/reference/sindex.rst
index 3201106d..8deab80a 100644
--- a/doc/source/docs/reference/sindex.rst
+++ b/doc/source/docs/reference/sindex.rst
@@ -30,6 +30,7 @@ methods:
 
     intersection
     is_empty
+    nearest_all
     query
     query_bulk
     size
@@ -42,4 +43,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 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..39346df5 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,21 @@ 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` but does not support the ``"inner"`` option.
+
+:func:`sjoin_nearest` 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..93709f4f 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),"
@@ -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/io/arrow.py b/geopandas/io/arrow.py
index 9d2ebdcd..9aab7c3e 100644
--- a/geopandas/io/arrow.py
+++ b/geopandas/io/arrow.py
@@ -8,7 +8,7 @@ from geopandas._compat import import_optional_dependency
 from geopandas.array import from_wkb
 from geopandas import GeoDataFrame
 import geopandas
-
+from .file import _expand_user
 
 METADATA_VERSION = "0.1.0"
 # reference: https://github.com/geopandas/geo-arrow-spec
@@ -243,6 +243,7 @@ def _to_parquet(df, path, index=None, compression="snappy", **kwargs):
         "pyarrow.parquet", extra="pyarrow is required for Parquet support."
     )
 
+    path = _expand_user(path)
     table = _geopandas_to_arrow(df, index=index)
     parquet.write_table(table, path, compression=compression, **kwargs)
 
@@ -290,6 +291,7 @@ def _to_feather(df, path, index=None, compression=None, **kwargs):
     if pyarrow.__version__ < LooseVersion("0.17.0"):
         raise ImportError("pyarrow >= 0.17 required for Feather support")
 
+    path = _expand_user(path)
     table = _geopandas_to_arrow(df, index=index)
     feather.write_feather(table, path, compression=compression, **kwargs)
 
@@ -452,6 +454,7 @@ def _read_parquet(path, columns=None, storage_options=None, **kwargs):
         path, filesystem=filesystem, storage_options=storage_options
     )
 
+    path = _expand_user(path)
     kwargs["use_pandas_metadata"] = True
     table = parquet.read_table(path, columns=columns, filesystem=filesystem, **kwargs)
 
@@ -513,5 +516,6 @@ def _read_feather(path, columns=None, **kwargs):
     if pyarrow.__version__ < LooseVersion("0.17.0"):
         raise ImportError("pyarrow >= 0.17 required for Feather support")
 
+    path = _expand_user(path)
     table = feather.read_table(path, columns=columns, **kwargs)
     return _arrow_to_geopandas(table)
diff --git a/geopandas/io/file.py b/geopandas/io/file.py
index acae8b3c..59467818 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -1,7 +1,8 @@
+import os
 from distutils.version import LooseVersion
 from pathlib import Path
-
 import warnings
+
 import numpy as np
 import pandas as pd
 
@@ -13,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
 
@@ -60,6 +66,15 @@ _EXTENSION_TO_DRIVER = {
 }
 
 
+def _expand_user(path):
+    """Expand paths that use ~."""
+    if isinstance(path, str):
+        path = os.path.expanduser(path)
+    elif isinstance(path, Path):
+        path = path.expanduser()
+    return path
+
+
 def _check_fiona(func):
     if fiona is None:
         raise ImportError(
@@ -149,6 +164,8 @@ def _read_file(filename, bbox=None, mask=None, rows=None, **kwargs):
     by using the encoding keyword parameter, e.g. ``encoding='utf-8'``.
     """
     _check_fiona("'read_file' function")
+    filename = _expand_user(filename)
+
     if _is_url(filename):
         req = _urlopen(filename)
         path_or_bytes = req.read()
@@ -336,6 +353,8 @@ def _to_file(
     by using the encoding keyword parameter, e.g. ``encoding='utf-8'``.
     """
     _check_fiona("'to_file' method")
+    filename = _expand_user(filename)
+
     if index is None:
         # Determine if index attribute(s) should be saved to file
         index = list(df.index.names) != [None] or type(df.index) not in (
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 4a7e6fc5..42141190 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -206,6 +206,49 @@ class BaseSpatialIndex:
         """
         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 : {GeoSeries, GeometryArray, numpy.array of PyGEOS geometries}
+            Accepts GeoPandas geometry iterables (GeoSeries, GeometryArray)
+            or a numpy array of PyGEOS geometries.
+        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
+
     @property
     def size(self):
         """Size of the spatial index
@@ -470,6 +513,8 @@ if compat.HAS_PYGEOS:
     from . import array  # noqa
     import pygeos  # noqa
 
+    _PYGEOS_PREDICATES = {p.name for p in pygeos.strtree.BinaryPredicate} | set([None])
+
     class PyGEOSSTRTreeIndex(pygeos.STRtree):
         """A simple wrapper around pygeos's STRTree.
 
@@ -509,7 +554,7 @@ if compat.HAS_PYGEOS:
             {'contains', 'crosses', 'covered_by', None, 'intersects', 'within', \
 'touches', 'overlaps', 'contains_properly', 'covers'}
             """
-            return {p.name for p in pygeos.strtree.BinaryPredicate} | set([None])
+            return _PYGEOS_PREDICATES
 
         @doc(BaseSpatialIndex.query)
         def query(self, geometry, predicate=None, sort=False):
@@ -531,6 +576,30 @@ 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
+                or a GeoPandas GeoSeries/GeometryArray.
+
+            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
+            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:
@@ -539,12 +608,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)
 
@@ -556,6 +621,13 @@ if compat.HAS_PYGEOS:
 
             return res
 
+        @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 ef47fa1e..ca516f81 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -185,8 +185,8 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
     2       2.0       2.0  POLYGON ((4.00000 4.00000, 4.00000 3.00000, 3....
     3       1.0       NaN  POLYGON ((2.00000 0.00000, 0.00000 0.00000, 0....
-    4       2.0       NaN  MULTIPOLYGON (((4.00000 3.00000, 4.00000 2.000...
-    5       NaN       1.0  MULTIPOLYGON (((3.00000 2.00000, 3.00000 1.000...
+    4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
+    5       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     6       NaN       2.0  POLYGON ((3.00000 5.00000, 5.00000 5.00000, 5....
 
     >>> geopandas.overlay(df1, df2, how='intersection')
@@ -198,8 +198,8 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     >>> geopandas.overlay(df1, df2, how='symmetric_difference')
        df1_data  df2_data                                           geometry
     0       1.0       NaN  POLYGON ((2.00000 0.00000, 0.00000 0.00000, 0....
-    1       2.0       NaN  MULTIPOLYGON (((4.00000 3.00000, 4.00000 2.000...
-    2       NaN       1.0  MULTIPOLYGON (((3.00000 2.00000, 3.00000 1.000...
+    1       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
+    2       NaN       1.0  MULTIPOLYGON (((2.00000 2.00000, 3.00000 2.000...
     3       NaN       2.0  POLYGON ((3.00000 5.00000, 5.00000 5.00000, 5....
 
     >>> geopandas.overlay(df1, df2, how='difference')
@@ -213,7 +213,7 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
     1       2.0       1.0  POLYGON ((2.00000 2.00000, 2.00000 3.00000, 3....
     2       2.0       2.0  POLYGON ((4.00000 4.00000, 4.00000 3.00000, 3....
     3       1.0       NaN  POLYGON ((2.00000 0.00000, 0.00000 0.00000, 0....
-    4       2.0       NaN  MULTIPOLYGON (((4.00000 3.00000, 4.00000 2.000...
+    4       2.0       NaN  MULTIPOLYGON (((3.00000 3.00000, 4.00000 3.000...
 
     See also
     --------
@@ -317,32 +317,64 @@ def overlay(df1, df2, how="intersection", keep_geom_type=None, make_valid=True):
             result = dfunion[dfunion["__idx1"].notnull()].copy()
 
     if keep_geom_type:
-        key_order = result.keys()
-        exploded = result.reset_index(drop=True).explode()
-        exploded = exploded.reset_index(level=0)
+        geom_type = df1.geom_type.iloc[0]
+
+        # First we filter the geometry types inside GeometryCollections objects
+        # (e.g. GeometryCollection([polygon, point]) -> polygon)
+        # we do this separately on only the relevant rows, as this is an expensive
+        # operation (an expensive no-op for geometry types other than collections)
+        is_collection = result.geom_type == "GeometryCollection"
+        if is_collection.any():
+            geom_col = result._geometry_column_name
+            collections = result[[geom_col]][is_collection]
+
+            exploded = collections.reset_index(drop=True).explode()
+            exploded = exploded.reset_index(level=0)
+
+            orig_num_geoms_exploded = exploded.shape[0]
+            if geom_type in polys:
+                exploded = exploded.loc[exploded.geom_type.isin(polys)]
+            elif geom_type in lines:
+                exploded = exploded.loc[exploded.geom_type.isin(lines)]
+            elif geom_type in points:
+                exploded = exploded.loc[exploded.geom_type.isin(points)]
+            else:
+                raise TypeError(
+                    "`keep_geom_type` does not support {}.".format(geom_type)
+                )
+            num_dropped_collection = orig_num_geoms_exploded - exploded.shape[0]
+
+            # level_0 created with above reset_index operation
+            # and represents the original geometry collections
+            # TODO avoiding dissolve to call unary_union in this case could further
+            # improve performance (we only need to collect geometries in their
+            # respective Multi version)
+            dissolved = exploded.dissolve(by="level_0")
+            result.loc[is_collection, geom_col] = dissolved[geom_col].values
+        else:
+            num_dropped_collection = 0
 
+        # Now we filter all geometries (in theory we don't need to do this
+        # again for the rows handled above for GeometryCollections, but filtering
+        # them out is probably more expensive as simply including them when this
+        # is typically about only a few rows)
         orig_num_geoms = result.shape[0]
-        geom_type = df1.geom_type.iloc[0]
         if geom_type in polys:
-            exploded = exploded.loc[exploded.geom_type.isin(polys)]
+            result = result.loc[result.geom_type.isin(polys)]
         elif geom_type in lines:
-            exploded = exploded.loc[exploded.geom_type.isin(lines)]
+            result = result.loc[result.geom_type.isin(lines)]
         elif geom_type in points:
-            exploded = exploded.loc[exploded.geom_type.isin(points)]
+            result = result.loc[result.geom_type.isin(points)]
         else:
             raise TypeError("`keep_geom_type` does not support {}.".format(geom_type))
+        num_dropped = orig_num_geoms - result.shape[0]
 
-        # level_0 created with above reset_index operation
-        # and represents the original geometry collections
-        result = exploded.dissolve(by="level_0")[key_order]
-
-        if (result.shape[0] != orig_num_geoms) and keep_geom_type_warning:
-            num_dropped = orig_num_geoms - result.shape[0]
+        if (num_dropped > 0 or num_dropped_collection > 0) and keep_geom_type_warning:
             warnings.warn(
                 "`keep_geom_type=True` in overlay resulted in {} dropped "
                 "geometries of different geometry types than df1 has. "
                 "Set `keep_geom_type=False` to retain all "
-                "geometries".format(num_dropped),
+                "geometries".format(num_dropped + num_dropped_collection),
                 UserWarning,
                 stacklevel=2,
             )
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index ca6f921e..893a403a 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)
@@ -144,7 +127,9 @@ stria    AUT    416600.0
     return joined
 
 
-def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
+def _basic_checks(
+    left_df, right_df, how, lsuffix, rsuffix, allowed_hows=("left", "right", "inner")
+):
     """Checks the validity of join input parameters.
 
     `how` must be one of the valid options.
@@ -155,7 +140,7 @@ def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
     ------------
     left_df : GeoDataFrame
     right_df : GeoData Frame
-    how : str, one of 'left', 'right', 'inner'
+    how : str, one of allowed_hows
         join type
     lsuffix : str
         left index suffix
@@ -172,7 +157,6 @@ def _basic_checks(left_df, right_df, how, lsuffix, rsuffix):
             "'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)
@@ -249,16 +233,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 +293,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 +311,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 +333,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 +350,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 = "left",
+    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 'left'
+        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
+    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, allowed_hows=("left", "right")
+    )
+
+    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/io/tests/test_arrow.py b/geopandas/io/tests/test_arrow.py
index 81d6132b..4732599a 100644
--- a/geopandas/io/tests/test_arrow.py
+++ b/geopandas/io/tests/test_arrow.py
@@ -7,6 +7,7 @@ import pytest
 from pandas import DataFrame, read_parquet as pd_read_parquet
 from pandas.testing import assert_frame_equal
 import numpy as np
+from shapely.geometry import box
 
 import geopandas
 from geopandas import GeoDataFrame, read_file, read_parquet, read_feather
@@ -541,3 +542,21 @@ def test_non_fsspec_url_with_storage_options_raises():
 def test_prefers_pyarrow_fs():
     filesystem, _ = _get_filesystem_path("file:///data.parquet")
     assert isinstance(filesystem, pyarrow.fs.LocalFileSystem)
+
+
+def test_write_read_parquet_expand_user():
+    gdf = geopandas.GeoDataFrame(geometry=[box(0, 0, 10, 10)], crs="epsg:4326")
+    test_file = "~/test_file.parquet"
+    gdf.to_parquet(test_file)
+    pq_df = geopandas.read_parquet(test_file)
+    assert_geodataframe_equal(gdf, pq_df, check_crs=True)
+    os.remove(os.path.expanduser(test_file))
+
+
+def test_write_read_feather_expand_user():
+    gdf = geopandas.GeoDataFrame(geometry=[box(0, 0, 10, 10)], crs="epsg:4326")
+    test_file = "~/test_file.feather"
+    gdf.to_feather(test_file)
+    f_df = geopandas.read_feather(test_file)
+    assert_geodataframe_equal(gdf, f_df, check_crs=True)
+    os.remove(os.path.expanduser(test_file))
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index ee304fc1..3267b202 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -837,3 +837,14 @@ def test_to_file__undetermined_driver(tmp_path, df_nybb):
     df_nybb.to_file(shpdir)
     assert shpdir.is_dir()
     assert list(shpdir.glob("*.shp"))
+
+
+@pytest.mark.parametrize(
+    "test_file", [(pathlib.Path("~/test_file.geojson")), "~/test_file.geojson"]
+)
+def test_write_read_file(test_file):
+    gdf = geopandas.GeoDataFrame(geometry=[box(0, 0, 10, 10)], crs=_CRS)
+    gdf.to_file(test_file, driver="GeoJSON")
+    df_json = geopandas.read_file(test_file)
+    assert_geodataframe_equal(gdf, df_json, check_crs=True)
+    os.remove(os.path.expanduser(test_file))
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index a07bf62b..ec9b6820 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -575,6 +575,11 @@ def test_keep_geom_type_geometry_collection():
     df1 = read_file(os.path.join(DATA, "geom_type", "df1.geojson"))
     df2 = read_file(os.path.join(DATA, "geom_type", "df2.geojson"))
 
+    with pytest.warns(UserWarning, match="`keep_geom_type=True` in overlay"):
+        intersection = overlay(df1, df2, keep_geom_type=None)
+    assert len(intersection) == 1
+    assert (intersection.geom_type == "Polygon").all()
+
     intersection = overlay(df1, df2, keep_geom_type=True)
     assert len(intersection) == 1
     assert (intersection.geom_type == "Polygon").all()
@@ -584,6 +589,44 @@ def test_keep_geom_type_geometry_collection():
     assert (intersection.geom_type == "GeometryCollection").all()
 
 
+def test_keep_geom_type_geometry_collection2():
+    polys1 = [
+        box(0, 0, 1, 1),
+        box(1, 1, 3, 3).union(box(1, 3, 5, 5)),
+    ]
+
+    polys2 = [
+        box(0, 0, 1, 1),
+        box(3, 1, 4, 2).union(box(4, 1, 5, 4)),
+    ]
+    df1 = GeoDataFrame({"left": [0, 1], "geometry": polys1})
+    df2 = GeoDataFrame({"right": [0, 1], "geometry": polys2})
+
+    result1 = overlay(df1, df2, keep_geom_type=True)
+    expected1 = GeoDataFrame(
+        {
+            "left": [0, 1],
+            "right": [0, 1],
+            "geometry": [box(0, 0, 1, 1), box(4, 3, 5, 4)],
+        }
+    )
+    assert_geodataframe_equal(result1, expected1)
+
+    result1 = overlay(df1, df2, keep_geom_type=False)
+    expected1 = GeoDataFrame(
+        {
+            "left": [0, 1, 1],
+            "right": [0, 0, 1],
+            "geometry": [
+                box(0, 0, 1, 1),
+                Point(1, 1),
+                GeometryCollection([box(4, 3, 5, 4), LineString([(3, 1), (3, 2)])]),
+            ],
+        }
+    )
+    assert_geodataframe_equal(result1, expected1)
+
+
 @pytest.mark.parametrize("make_valid", [True, False])
 def test_overlap_make_valid(make_valid):
     bowtie = Polygon([(1, 1), (9, 9), (9, 1), (1, 9), (1, 1)])
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index d3f72e18..48224886 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,323 @@ 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": "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"))
+    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("distance_col", (None, "distance"))
+    def test_empty_left_df_how_left(self, 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,
+        )
+        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()
+
+    def test_empty_join_due_to_max_distance_how_left(self):
+        # 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="left",
+            max_distance=1,
+            distance_col="distances",
+        )
+        expected = left.copy()
+        expected["index_right"] = [np.nan]
+        expected["distances"] = [np.nan]
+        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)
+
+    def test_max_distance_how_left(self):
+        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",
+            max_distance=1,
+            distance_col="distances",
+        )
+        expected = left.copy()
+        expected["index_right"] = [np.nan, 0]
+        expected["distances"] = [np.nan, 0]
+        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(
+        "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],
+    ):
+        left = geopandas.GeoDataFrame({"geometry": geo_left})
+        right = geopandas.GeoDataFrame({"geometry": geo_right})
+        expected_gdf = left.iloc[expected_left]
+        expected_gdf["index_right"] = expected_right
+        # without distance col
+        joined = sjoin_nearest(left, right, how="left")
+        assert_geodataframe_equal(expected_gdf, joined)
+        # 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)
+
+    @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]
+        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)