Repository: geopandas/geopandas

Issue:
API: default behaviour of alignment in binary ops
#1668 introduces `align` keyword to control whether two GeoSeries gets aligned in binary, row by row operations or not. However, the default behaviour may require a bit more discussion.

Moving from the original PR here (from https://github.com/geopandas/geopandas/pull/1668#discussion_r578692627).

> @jorisvandenbossche:  Idea: make the default None, which is treated as True but gives a warning when alignment is actually done. This gives a way for the user to indicate they are OK with alignment without getting warnings about it.
(small disadvantage is that the signature in the function is less explicit than with align=True)

> @martinfleis: I do understand your point, but considering how troublesome this behaviour is (I responded to countless questions on this), I would prefer having an explicit align=True in the signature. Having None/True/False is a nice trick for advanced users to silence warnings, but in this specific case, I would prefer understanding for basic users (I really don't want to answer the same question again :D).

> @jorisvandenbossche: Then we should maybe rather make it an error by default? To make it really clear for users they are doing something wrong (with a good explanation then)

> @martinfleis: That'd be a breaking change.



Discussion:
jorisvandenbossche: I am certainly fine with not exhaustively testing this for every single method

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ab8cbb04..69c084d6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,12 +11,15 @@ Notes on dependencies:
 API changes:
 
 - `unary_union` is now deprecated and replaced by the `union_all` method (#3007).
+- `align` keyword in binary methods now defaults to `None`, treated as True. Explicit True
+  will silence the warning about mismachted indices. (#3212)
 - The `sjoin` method will now preserve the name of the index of the right
   GeoDataFrame, if it has one, instead of always using `"index_right"` as the
   name for the resulting column in the return value (#846, #2144).
 
 New methods:
 
+- Added `line_merge` method from shapely to GeoSeries/GeoDataframe (#3214).
 - Added `set_precision` and `get_precision` methods from shapely to GeoSeries/GeoDataframe (#3175).
 - Added `count_coordinates` method from shapely to GeoSeries/GeoDataframe (#3026).
 - Added `minimum_clearance` method from shapely to GeoSeries/GeoDataframe (#2989).
@@ -41,7 +44,7 @@ New features and improvements:
   `buffer(0)` (#3113).
 - Passing `"geometry"` as `dtype` to `pd.read_csv` will now return a GeoSeries for
   the specified columns (#3101).
-- 
+-
 API changes:
 - Added support for ``mask`` keyword for pyogrio engine for pyogrio >= 0.7.0 (#3062).
 
@@ -67,16 +70,16 @@ Deprecations and compatibility notes:
   removed. New sample datasets are now available in the
   [geodatasets](https://geodatasets.readthedocs.io/en/latest/) package (#3084).
 - Many longstanding deprecated functions, methods and properties have been removed (#3174)
-  - Removed deprecated functions  
-    `geopandas.io.read_file`, `geopandas.io.to_file` and `geopandas.io.sql.read_postgis`. 
-    `geopandas.read_file`, `geopandas.read_postgis` and the GeoDataFrame/GeoSeries `to_file(..)` 
+  - Removed deprecated functions
+    `geopandas.io.read_file`, `geopandas.io.to_file` and `geopandas.io.sql.read_postgis`.
+    `geopandas.read_file`, `geopandas.read_postgis` and the GeoDataFrame/GeoSeries `to_file(..)`
     method should be used instead.
   - Removed deprecated `GeometryArray.data` property, `np.asarray(..)` or the `to_numpy()`
     method should be used instead.
   - Removed deprecated `sindex.query_bulk` method, using `sindex.query` instead.
   - Removed deprecated `sjoin` parameter `op`, `predicate` should be supplied instead.
-  - Removed deprecated GeoSeries/ GeoDataFrame methods `__xor__`, `__or__`, `__and__` and 
-    `__sub__`. Instead use methods `symmetric_difference`, `union`, `intersection` and 
+  - Removed deprecated GeoSeries/ GeoDataFrame methods `__xor__`, `__or__`, `__and__` and
+    `__sub__`. Instead use methods `symmetric_difference`, `union`, `intersection` and
     `difference` respectively.
 - Fixes for compatibility with psycopg (#3167).
 
diff --git a/ci/envs/310-pd20-conda-forge.yaml b/ci/envs/310-pd20-conda-forge.yaml
index 4854fdb8..91e55bc9 100644
--- a/ci/envs/310-pd20-conda-forge.yaml
+++ b/ci/envs/310-pd20-conda-forge.yaml
@@ -11,8 +11,7 @@ dependencies:
   - pyproj
   - packaging
   # testing
-  # TEMP pin, 8.1 breaks doctestplus, see https://github.com/scientific-python/pytest-doctestplus/issues/239
-  - pytest==8.0.*
+  - pytest
   - pytest-cov
   - pytest-xdist
   - fsspec
diff --git a/ci/envs/311-dev.yaml b/ci/envs/311-dev.yaml
index 93c8c534..313e8758 100644
--- a/ci/envs/311-dev.yaml
+++ b/ci/envs/311-dev.yaml
@@ -5,7 +5,6 @@ dependencies:
   - python=3.11
   - cython
   # required
-  - numpy
   - pyproj
   - geos
   - packaging
@@ -24,8 +23,8 @@ dependencies:
     - mapclassify>=2.4.0
     # dev versions of packages
     - --pre --index-url https://pypi.anaconda.org/scientific-python-nightly-wheels/simple --extra-index-url https://pypi.fury.io/arrow-nightlies/ --extra-index-url https://pypi.org/simple
-    # TEMP pin, see https://github.com/geopandas/geopandas/issues/3207
-    - fiona==1.9.*
+    - numpy>=2.0.0.dev
+    - fiona
     - pandas
     - matplotlib
     - pyarrow
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 6a8fbbee..ebeec66a 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -137,6 +137,7 @@ Linestring operations
    :toctree: api/
 
    GeoSeries.interpolate
+   GeoSeries.line_merge
    GeoSeries.project
 
 Aggregating and exploding
diff --git a/geopandas/array.py b/geopandas/array.py
index 8a8d4caf..5ce017f6 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -665,6 +665,11 @@ class GeometryArray(ExtensionArray):
             shapely.transform(self._data, transformation, include_z), crs=self.crs
         )
 
+    def line_merge(self, directed=False):
+        return GeometryArray(
+            shapely.line_merge(self._data, directed=directed), crs=self.crs
+        )
+
     def set_precision(self, grid_size, mode="valid_output"):
         return GeometryArray(
             shapely.set_precision(self._data, grid_size=grid_size, mode=mode),
@@ -1181,7 +1186,7 @@ class GeometryArray(ExtensionArray):
 
         result = take(self._data, indices, allow_fill=allow_fill, fill_value=fill_value)
         if allow_fill and fill_value is None:
-            result[pd.isna(result)] = None
+            result[~shapely.is_valid_input(result)] = None
         return GeometryArray(result, crs=self.crs)
 
     def _fill(self, idx, value):
diff --git a/geopandas/base.py b/geopandas/base.py
index dfa59bc4..c19b4d25 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -27,13 +27,24 @@ def is_geometry_type(data):
 
 def _delegate_binary_method(op, this, other, align, *args, **kwargs):
     # type: (str, GeoSeries, GeoSeries) -> GeoSeries/Series
+    if align is None:
+        align = True
+        maybe_warn = True
+    else:
+        maybe_warn = False
     this = this.geometry
     if isinstance(other, GeoPandasBase):
         if align and not this.index.equals(other.index):
-            warn(
-                "The indices of the two GeoSeries are different.",
-                stacklevel=4,
-            )
+            if maybe_warn:
+                warn(
+                    "The indices of the left and right GeoSeries' are not equal, and "
+                    "therefore they will be aligned (reordering and/or introducing "
+                    "missing values) before executing the operation. If this alignment "
+                    "is the desired behaviour, you can silence this warning by passing "
+                    "'align=True'. If you don't want alignment and protect yourself of "
+                    "accidentally aligning, you can pass 'align=False'.",
+                    stacklevel=4,
+                )
             this, other = this.align(other.geometry)
         else:
             other = other.geometry
@@ -1606,6 +1617,73 @@ GeometryCollection
         """
         return _delegate_geo_method("force_3d", self, z=z)
 
+    def line_merge(self, directed=False):
+        """Returns (Multi)LineStrings formed by combining the lines in a
+        MultiLineString.
+
+        Lines are joined together at their endpoints in case two lines are intersecting.
+        Lines are not joined when 3 or more lines are intersecting at the endpoints.
+        Line elements that cannot be joined are kept as is in the resulting
+        MultiLineString.
+
+        The direction of each merged LineString will be that of the majority of the
+        LineStrings from which it was derived. Except if ``directed=True`` is specified,
+        then the operation will not change the order of points within lines and so only
+        lines which can be joined with no change in direction are merged.
+
+        Non-linear geometeries result in an empty GeometryCollection.
+
+        Parameters
+        ----------
+        directed : bool, default False
+            Only combine lines if possible without changing point order.
+            Requires GEOS >= 3.11.0
+
+        Returns
+        -------
+        GeoSeries
+
+        Examples
+        --------
+        >>> from shapely.geometry import MultiLineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         MultiLineString([[(0, 2), (0, 10)], [(0, 10), (5, 10)]]),
+        ...         MultiLineString([[(0, 2), (0, 10)], [(0, 11), (5, 10)]]),
+        ...         MultiLineString(),
+        ...         MultiLineString([[(0, 0), (1, 0)], [(0, 0), (3, 0)]]),
+        ...         Point(0, 0),
+        ...     ]
+        ... )
+        >>> s
+        0    MULTILINESTRING ((0 2, 0 10), (0 10, 5 10))
+        1    MULTILINESTRING ((0 2, 0 10), (0 11, 5 10))
+        2                          MULTILINESTRING EMPTY
+        3       MULTILINESTRING ((0 0, 1 0), (0 0, 3 0))
+        4                                    POINT (0 0)
+        dtype: geometry
+
+        >>> s.line_merge()
+        0                   LINESTRING (0 2, 0 10, 5 10)
+        1    MULTILINESTRING ((0 2, 0 10), (0 11, 5 10))
+        2                       GEOMETRYCOLLECTION EMPTY
+        3                     LINESTRING (1 0, 0 0, 3 0)
+        4                       GEOMETRYCOLLECTION EMPTY
+        dtype: geometry
+
+        With ``directed=True``, you can avoid changing the order of points within lines
+        and merge only lines where no change of direction is required:
+
+        >>> s.line_merge(directed=True)
+        0                   LINESTRING (0 2, 0 10, 5 10)
+        1    MULTILINESTRING ((0 2, 0 10), (0 11, 5 10))
+        2                       GEOMETRYCOLLECTION EMPTY
+        3       MULTILINESTRING ((0 0, 1 0), (0 0, 3 0))
+        4                       GEOMETRYCOLLECTION EMPTY
+        dtype: geometry
+        """
+        return _delegate_geo_method("line_merge", self, directed=directed)
+
     #
     # Reduction operations that return a Shapely geometry
     #
@@ -1670,7 +1748,7 @@ GeometryCollection
     # Binary operations that return a pandas Series
     #
 
-    def contains(self, other, align=True):
+    def contains(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that contains `other`.
 
@@ -1693,9 +1771,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if it
             is contained.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -1786,7 +1864,7 @@ GeometryCollection
         """
         return _binary_op("contains", self, other, align)
 
-    def contains_properly(self, other, align=True):
+    def contains_properly(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is completely inside ``other``, with no common
         boundary points.
@@ -1806,9 +1884,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if it
             is contained.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -1907,7 +1985,7 @@ GeometryCollection
         """
         return _binary_op("contains_properly", self, other, align)
 
-    def dwithin(self, other, distance, align=True):
+    def dwithin(self, other, distance, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is within a set distance from ``other``.
 
@@ -1926,9 +2004,9 @@ GeometryCollection
             applied to all geometries. An array or Series will be applied elementwise.
             If np.array or pd.Series are used then it must have same length as the
             GeoSeries.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2018,7 +2096,7 @@ GeometryCollection
         """
         return _binary_op("dwithin", self, other, distance=distance, align=align)
 
-    def geom_equals(self, other, align=True):
+    def geom_equals(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry equal to `other`.
 
@@ -2036,9 +2114,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test for
             equality.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2128,7 +2206,7 @@ GeometryCollection
         """
         return _binary_op("geom_equals", self, other, align)
 
-    def geom_almost_equals(self, other, decimal=6, align=True):
+    def geom_almost_equals(self, other, decimal=6, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` if
         each aligned geometry is approximately equal to `other`.
 
@@ -2146,9 +2224,9 @@ GeometryCollection
             The GeoSeries (elementwise) or geometric object to compare to.
         decimal : int
             Decimal place precision used when testing for approximate equality.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2205,7 +2283,7 @@ GeometryCollection
             "geom_equals_exact", self, other, tolerance=tolerance, align=align
         )
 
-    def geom_equals_exact(self, other, tolerance, align=True):
+    def geom_equals_exact(self, other, tolerance, align=None):
         """Return True for all geometries that equal aligned *other* to a given
         tolerance, else False.
 
@@ -2220,9 +2298,9 @@ GeometryCollection
             The GeoSeries (elementwise) or geometric object to compare to.
         tolerance : float
             Decimal place precision used when testing for approximate equality.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2270,7 +2348,7 @@ GeometryCollection
             "geom_equals_exact", self, other, tolerance=tolerance, align=align
         )
 
-    def crosses(self, other, align=True):
+    def crosses(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that cross `other`.
 
@@ -2288,9 +2366,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             crossed.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2382,7 +2460,7 @@ GeometryCollection
         """
         return _binary_op("crosses", self, other, align)
 
-    def disjoint(self, other, align=True):
+    def disjoint(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry disjoint to `other`.
 
@@ -2399,9 +2477,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             disjoint.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2483,7 +2561,7 @@ GeometryCollection
         """
         return _binary_op("disjoint", self, other, align)
 
-    def intersects(self, other, align=True):
+    def intersects(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that intersects `other`.
 
@@ -2500,9 +2578,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             intersected.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2594,7 +2672,7 @@ GeometryCollection
         """
         return _binary_op("intersects", self, other, align)
 
-    def overlaps(self, other, align=True):
+    def overlaps(self, other, align=None):
         """Returns True for all aligned geometries that overlap *other*, else False.
 
         Geometries overlaps if they have more than one but not all
@@ -2612,9 +2690,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if
             overlaps.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2705,7 +2783,7 @@ GeometryCollection
         """
         return _binary_op("overlaps", self, other, align)
 
-    def touches(self, other, align=True):
+    def touches(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that touches `other`.
 
@@ -2723,9 +2801,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             touched.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2817,7 +2895,7 @@ GeometryCollection
         """
         return _binary_op("touches", self, other, align)
 
-    def within(self, other, align=True):
+    def within(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is within `other`.
 
@@ -2839,9 +2917,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if each
             geometry is within.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2931,7 +3009,7 @@ GeometryCollection
         """
         return _binary_op("within", self, other, align)
 
-    def covers(self, other, align=True):
+    def covers(self, other, align=None):
         """
         Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is entirely covering `other`.
@@ -2953,9 +3031,9 @@ GeometryCollection
         ----------
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to check is being covered.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3045,7 +3123,7 @@ GeometryCollection
         """
         return _binary_op("covers", self, other, align)
 
-    def covered_by(self, other, align=True):
+    def covered_by(self, other, align=None):
         """
         Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is entirely covered by `other`.
@@ -3066,9 +3144,9 @@ GeometryCollection
         ----------
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to check is being covered.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3159,7 +3237,7 @@ GeometryCollection
         """
         return _binary_op("covered_by", self, other, align)
 
-    def distance(self, other, align=True):
+    def distance(self, other, align=None):
         """Returns a ``Series`` containing the distance to aligned `other`.
 
         The operation works on a 1-to-1 row-wise manner:
@@ -3172,9 +3250,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             distance to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
 
         Returns
@@ -3271,9 +3349,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             distance to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
         densify : float (default None)
             A value between 0 and 1, that splits each subsegment of a line string
             into equal length segments, making the approximation less coarse.
@@ -3393,9 +3471,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             distance to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
         densify : float (default None)
             A value between 0 and 1, that splits each subsegment of a line string
             into equal length segments, making the approximation less coarse.
@@ -3494,7 +3572,7 @@ GeometryCollection
     # Binary operations that return a GeoSeries
     #
 
-    def difference(self, other, align=True):
+    def difference(self, other, align=None):
         """Returns a ``GeoSeries`` of the points in each aligned geometry that
         are not in `other`.
 
@@ -3511,9 +3589,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             difference to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3605,7 +3683,7 @@ GeometryCollection
         """
         return _binary_geo("difference", self, other, align)
 
-    def symmetric_difference(self, other, align=True):
+    def symmetric_difference(self, other, align=None):
         """Returns a ``GeoSeries`` of the symmetric difference of points in
         each aligned geometry with `other`.
 
@@ -3626,9 +3704,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             symmetric difference to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3720,7 +3798,7 @@ GeometryCollection
         """
         return _binary_geo("symmetric_difference", self, other, align)
 
-    def union(self, other, align=True):
+    def union(self, other, align=None):
         """Returns a ``GeoSeries`` of the union of points in each aligned geometry with
         `other`.
 
@@ -3738,9 +3816,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the union
             with.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3834,7 +3912,7 @@ GeometryCollection
         """
         return _binary_geo("union", self, other, align)
 
-    def intersection(self, other, align=True):
+    def intersection(self, other, align=None):
         """Returns a ``GeoSeries`` of the intersection of points in each
         aligned geometry with `other`.
 
@@ -3852,9 +3930,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             intersection with.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4016,7 +4094,7 @@ GeometryCollection
         clipped_geometry = geometry_array.clip_by_rect(xmin, ymin, xmax, ymax)
         return GeoSeries(clipped_geometry, index=self.index, crs=self.crs)
 
-    def shortest_line(self, other, align=True):
+    def shortest_line(self, other, align=None):
         """
         Returns the shortest two-point line between two geometries.
 
@@ -4037,9 +4115,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             shortest line with.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4119,7 +4197,7 @@ GeometryCollection
         """
         return _binary_geo("shortest_line", self, other, align)
 
-    def snap(self, other, tolerance, align=True):
+    def snap(self, other, tolerance, align=None):
         """Snaps an input geometry to reference geometry’s vertices.
 
         Vertices of the first geometry are snapped to vertices of the second. geometry,
@@ -4148,9 +4226,9 @@ GeometryCollection
             The Geoseries (elementwise) or geometric object to snap to.
         tolerance : float or array like
             Maximum distance between vertices that shall be snapped
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4505,7 +4583,7 @@ GeometryCollection
         """
         return _delegate_geo_method("simplify", self, *args, **kwargs)
 
-    def relate(self, other, align=True):
+    def relate(self, other, align=None):
         """
         Returns the DE-9IM intersection matrices for the geometries
 
@@ -4519,9 +4597,9 @@ GeometryCollection
         other : BaseGeometry or GeoSeries
             The other geometry to computed
             the DE-9IM intersection matrices from.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4610,7 +4688,7 @@ GeometryCollection
         """
         return _binary_op("relate", self, other, align)
 
-    def project(self, other, normalized=False, align=True):
+    def project(self, other, normalized=False, align=None):
         """
         Return the distance along each geometry nearest to *other*
 
@@ -4631,9 +4709,9 @@ GeometryCollection
         normalized : boolean
             If normalized is True, return the distance normalized to
             the length of the object.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index d5d361e5..790d371d 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -2117,12 +2117,12 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
 
         >>> groceries_w_communities = groceries.sjoin(chicago)
         >>> groceries_w_communities[["OBJECTID", "community", "geometry"]].head()
-             OBJECTID    community                           geometry
-        0          16       UPTOWN  MULTIPOINT ((-87.65661 41.97321))
-        87        365       UPTOWN  MULTIPOINT ((-87.65465 41.96138))
-        90        373       UPTOWN  MULTIPOINT ((-87.65598 41.96297))
-        140       582       UPTOWN  MULTIPOINT ((-87.67417 41.96977))
-        1          18  MORGAN PARK  MULTIPOINT ((-87.68136 41.69713))
+           OBJECTID       community                           geometry
+        0        16          UPTOWN  MULTIPOINT ((-87.65661 41.97321))
+        1        18     MORGAN PARK  MULTIPOINT ((-87.68136 41.69713))
+        2        22  NEAR WEST SIDE  MULTIPOINT ((-87.63918 41.86847))
+        3        23  NEAR WEST SIDE  MULTIPOINT ((-87.65495 41.87783))
+        4        27         CHATHAM  MULTIPOINT ((-87.62715 41.73623))
 
         Notes
         -----
@@ -2217,9 +2217,9 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
 
         >>> groceries_w_communities = groceries.sjoin_nearest(chicago)
         >>> groceries_w_communities[["Chain", "community", "geometry"]].head(2)
-                     Chain community                               geometry
-        0   VIET HOA PLAZA    UPTOWN  MULTIPOINT ((1168268.672 1933554.35))
-        87      JEWEL OSCO    UPTOWN  MULTIPOINT ((1168837.98 1929246.962))
+                       Chain    community                                geometry
+        0     VIET HOA PLAZA       UPTOWN   MULTIPOINT ((1168268.672 1933554.35))
+        1  COUNTY FAIR FOODS  MORGAN PARK  MULTIPOINT ((1162302.618 1832900.224))
 
 
         To include the distances:
@@ -2227,10 +2227,10 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         >>> groceries_w_communities = groceries.sjoin_nearest(chicago, \
 distance_col="distances")
         >>> groceries_w_communities[["Chain", "community", \
-"distances"]].head(2)  # doctest: +SKIP
-                     Chain community  distances
-        0   VIET HOA PLAZA    UPTOWN        0.0
-        87      JEWEL OSCO    UPTOWN        0.0
+"distances"]].head(2)
+                       Chain    community  distances
+        0     VIET HOA PLAZA       UPTOWN        0.0
+        1  COUNTY FAIR FOODS  MORGAN PARK        0.0
 
         In the following example, we get multiple groceries for Uptown because all
         results are equidistant (in this case zero because they intersect).
@@ -2240,7 +2240,7 @@ distance_col="distances")
 distance_col="distances", how="right")
         >>> uptown_results = \
 chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
-        >>> uptown_results[["Chain", "community"]]  # doctest: +SKIP
+        >>> uptown_results[["Chain", "community"]]
                     Chain community
         30  VIET HOA PLAZA    UPTOWN
         30      JEWEL OSCO    UPTOWN
diff --git a/geopandas/io/file.py b/geopandas/io/file.py
index 5d11dbc7..a3cd6ddc 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -6,6 +6,7 @@ import warnings
 
 import numpy as np
 import pandas as pd
+from geopandas.io.util import vsi_path
 from pandas.api.types import is_integer_dtype
 
 import shapely
@@ -57,6 +58,7 @@ def _import_fiona():
             FIONA_GE_19 = Version(Version(fiona.__version__).base_version) >= Version(
                 "1.9.0"
             )
+
         except ImportError as err:
             fiona = False
             fiona_import_error = str(err)
@@ -175,16 +177,6 @@ def _is_url(url):
         return False
 
 
-def _is_zip(path):
-    """Check if a given path is a zipfile"""
-    parsed = fiona.path.ParsedPath.from_uri(path)
-    return (
-        parsed.archive.endswith(".zip")
-        if parsed.archive
-        else parsed.path.endswith(".zip")
-    )
-
-
 def _read_file(filename, bbox=None, mask=None, rows=None, engine=None, **kwargs):
     """
     Returns a GeoDataFrame from a file or URL.
@@ -320,22 +312,7 @@ def _read_file_fiona(
         # Opening a file via URL or file-like-object above automatically detects a
         # zipped file. In order to match that behavior, attempt to add a zip scheme
         # if missing.
-        if _is_zip(str(path_or_bytes)):
-            parsed = fiona.parse_path(str(path_or_bytes))
-            if isinstance(parsed, fiona.path.ParsedPath):
-                # If fiona is able to parse the path, we can safely look at the scheme
-                # and update it to have a zip scheme if necessary.
-                schemes = (parsed.scheme or "").split("+")
-                if "zip" not in schemes:
-                    parsed.scheme = "+".join(["zip"] + schemes)
-                path_or_bytes = parsed.name
-            elif isinstance(parsed, fiona.path.UnparsedPath) and not str(
-                path_or_bytes
-            ).startswith("/vsi"):
-                # If fiona is unable to parse the path, it might have a Windows drive
-                # scheme. Try adding zip:// to the front. If the path starts with "/vsi"
-                # it is a legacy GDAL path type, so let it pass unmodified.
-                path_or_bytes = "zip://" + parsed.name
+        path_or_bytes = vsi_path(str(path_or_bytes))
 
     if from_bytes:
         reader = fiona.BytesCollection
diff --git a/geopandas/io/util.py b/geopandas/io/util.py
new file mode 100644
index 00000000..2bfbee8d
--- /dev/null
+++ b/geopandas/io/util.py
@@ -0,0 +1,118 @@
+"""Vendored, cut down version of pyogrio/util.py for use with fiona"""
+
+import re
+import sys
+from urllib.parse import urlparse
+
+
+def vsi_path(path: str) -> str:
+    """
+    Ensure path is a local path or a GDAL-compatible vsi path.
+
+    """
+
+    # path is already in GDAL format
+    if path.startswith("/vsi"):
+        return path
+
+    # Windows drive letters (e.g. "C:\") confuse `urlparse` as they look like
+    # URL schemes
+    if sys.platform == "win32" and re.match("^[a-zA-Z]\\:", path):
+        if not path.split("!")[0].endswith(".zip"):
+            return path
+
+        # prefix then allow to proceed with remaining parsing
+        path = f"zip://{path}"
+
+    path, archive, scheme = _parse_uri(path)
+
+    if scheme or archive or path.endswith(".zip"):
+        return _construct_vsi_path(path, archive, scheme)
+
+    return path
+
+
+# Supported URI schemes and their mapping to GDAL's VSI suffix.
+SCHEMES = {
+    "file": "file",
+    "zip": "zip",
+    "tar": "tar",
+    "gzip": "gzip",
+    "http": "curl",
+    "https": "curl",
+    "ftp": "curl",
+    "s3": "s3",
+    "gs": "gs",
+    "az": "az",
+    "adls": "adls",
+    "adl": "adls",  # fsspec uses this
+    "hdfs": "hdfs",
+    "webhdfs": "webhdfs",
+    # GDAL additionally supports oss and swift for remote filesystems, but
+    # those are for now not added as supported URI
+}
+
+CURLSCHEMES = {k for k, v in SCHEMES.items() if v == "curl"}
+
+
+def _parse_uri(path: str):
+    """
+    Parse a URI
+
+    Returns a tuples of (path, archive, scheme)
+
+    path : str
+        Parsed path. Includes the hostname and query string in the case
+        of a URI.
+    archive : str
+        Parsed archive path.
+    scheme : str
+        URI scheme such as "https" or "zip+s3".
+    """
+    parts = urlparse(path)
+
+    # if the scheme is not one of GDAL's supported schemes, return raw path
+    if parts.scheme and not all(p in SCHEMES for p in parts.scheme.split("+")):
+        return path, "", ""
+
+    # we have a URI
+    path = parts.path
+    scheme = parts.scheme or ""
+
+    if parts.query:
+        path += "?" + parts.query
+
+    if parts.scheme and parts.netloc:
+        path = parts.netloc + path
+
+    parts = path.split("!")
+    path = parts.pop() if parts else ""
+    archive = parts.pop() if parts else ""
+    return (path, archive, scheme)
+
+
+def _construct_vsi_path(path, archive, scheme) -> str:
+    """Convert a parsed path to a GDAL VSI path"""
+
+    prefix = ""
+    suffix = ""
+    schemes = scheme.split("+")
+
+    if "zip" not in schemes and (archive.endswith(".zip") or path.endswith(".zip")):
+        schemes.insert(0, "zip")
+
+    if schemes:
+        prefix = "/".join(
+            "vsi{0}".format(SCHEMES[p]) for p in schemes if p and p != "file"
+        )
+
+        if schemes[-1] in CURLSCHEMES:
+            suffix = f"{schemes[-1]}://"
+
+    if prefix:
+        if archive:
+            return "/{}/{}{}/{}".format(prefix, suffix, archive, path.lstrip("/"))
+        else:
+            return "/{}/{}{}".format(prefix, suffix, path)
+
+    return path
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index 97425da3..9ac00d7f 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -7,6 +7,7 @@ import pandas as pd
 
 from geopandas import GeoDataFrame
 from geopandas.array import _check_crs, _crs_mismatch_warn
+from geopandas._compat import PANDAS_GE_30
 
 
 def sjoin(
@@ -80,12 +81,12 @@ def sjoin(
 
     >>> groceries_w_communities = geopandas.sjoin(groceries, chicago)
     >>> groceries_w_communities.head()  # doctest: +SKIP
-            OBJECTID     Ycoord     Xcoord  ... GonorrF GonorrM Tuberc
-    0          16  41.973266 -87.657073  ...   170.8   468.7   13.6
-    87        365  41.961707 -87.654058  ...   170.8   468.7   13.6
-    90        373  41.963131 -87.656352  ...   170.8   468.7   13.6
-    140       582  41.969131 -87.674882  ...   170.8   468.7   13.6
-    1          18  41.696367 -87.681315  ...   800.5   741.1    2.6
+       OBJECTID       community                           geometry
+    0        16          UPTOWN  MULTIPOINT ((-87.65661 41.97321))
+    1        18     MORGAN PARK  MULTIPOINT ((-87.68136 41.69713))
+    2        22  NEAR WEST SIDE  MULTIPOINT ((-87.63918 41.86847))
+    3        23  NEAR WEST SIDE  MULTIPOINT ((-87.65495 41.87783))
+    4        27         CHATHAM  MULTIPOINT ((-87.62715 41.73623))
     [5 rows x 95 columns]
 
     See also
@@ -106,7 +107,9 @@ def sjoin(
 
     indices = _geom_predicate_query(left_df, right_df, predicate, distance)
 
-    joined = _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix)
+    joined, _ = _frame_join(
+        left_df, right_df, indices, None, how, lsuffix, rsuffix, predicate
+    )
 
     return joined
 
@@ -185,19 +188,16 @@ def _geom_predicate_query(left_df, right_df, predicate, distance):
         l_idx, r_idx = sindex.query(
             input_geoms, predicate=predicate, sort=False, distance=distance
         )
-        indices = pd.DataFrame({"_key_left": l_idx, "_key_right": r_idx})
     else:
         # when sindex is empty / has no valid geometries
-        indices = pd.DataFrame(columns=["_key_left", "_key_right"], dtype=float)
+        l_idx, r_idx = np.array([], dtype=np.intp), np.array([], dtype=np.intp)
 
     if original_predicate == "within":
         # within is implemented as the inverse of contains
         # flip back the results
-        indices = indices.rename(
-            columns={"_key_left": "_key_right", "_key_right": "_key_left"}
-        )
+        r_idx, l_idx = l_idx, r_idx
 
-    return indices
+    return l_idx, r_idx
 
 
 def _reset_index_with_suffix(df, suffix, other):
@@ -206,7 +206,12 @@ def _reset_index_with_suffix(df, suffix, other):
     column names.
     """
     index_original = df.index.names
-    df_reset = df.reset_index()
+    if PANDAS_GE_30:
+        df_reset = df.reset_index()
+    else:
+        # we already made a copy of the dataframe in _frame_join before getting here
+        df_reset = df
+        df_reset.reset_index(inplace=True)
     column_names = df_reset.columns.to_numpy(copy=True)
     for i, label in enumerate(index_original):
         # if the original label was None, add suffix to auto-generated name
@@ -251,8 +256,16 @@ def _process_column_names_with_suffix(
             return f"{x}_{suffix}"
         return x
 
-    lrenamer = partial(renamer, suffix=lsuffix, geometry=left_df._geometry_column_name)
-    rrenamer = partial(renamer, suffix=rsuffix, geometry=right_df._geometry_column_name)
+    lrenamer = partial(
+        renamer,
+        suffix=lsuffix,
+        geometry=getattr(left_df, "_geometry_column_name", None),
+    )
+    rrenamer = partial(
+        renamer,
+        suffix=rsuffix,
+        geometry=getattr(right_df, "_geometry_column_name", None),
+    )
 
     # TODO retain index name?
     left_renamed = pd.Index([lrenamer(lab) for lab in left])
@@ -284,7 +297,10 @@ def _restore_index(joined, index_names, index_names_original):
     Set back the the original index columns, and restoring their name as `None`
     if they didn't have a name originally.
     """
-    joined = joined.set_index(list(index_names))
+    if PANDAS_GE_30:
+        joined = joined.set_index(list(index_names))
+    else:
+        joined.set_index(list(index_names), inplace=True)
 
     # restore the fact that the index didn't have a name
     joined_index_names = list(joined.index.names)
@@ -295,42 +311,93 @@ def _restore_index(joined, index_names, index_names_original):
     return joined
 
 
-def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
+def _adjust_indexers(indices, distances, original_length, how, predicate):
+    """
+    The left/right indexers from the query represents an inner join.
+    For a left or right join, we need to adjust them to include the rows
+    that would not be present in an inner join.
+    """
+    # the indices represent an inner join, no adjustment needed
+    if how == "inner":
+        if predicate == "within":
+            # except for the within predicate, where we switched to contains
+            # with swapped left/right -> need to re-sort to have consistent result
+            l_idx, r_idx = indices
+            indexer = np.lexsort((r_idx, l_idx))
+            indices = l_idx[indexer], r_idx[indexer]
+            if distances is not None:
+                distances = distances[indexer]
+        return indices, distances
+
+    l_idx, r_idx = indices
+
+    if how == "right":
+        # re-sort so it is sorted by the right indexer
+        indexer = np.lexsort((l_idx, r_idx))
+        l_idx, r_idx = l_idx[indexer], r_idx[indexer]
+        if distances is not None:
+            distances = distances[indexer]
+
+        # switch order
+        r_idx, l_idx = l_idx, r_idx
+
+    # determine which indices are missing and where they would need to be inserted
+    idx = np.arange(original_length)
+    l_idx_missing = idx[~np.isin(idx, l_idx)]
+    insert_idx = np.searchsorted(l_idx, l_idx_missing)
+    # for the left indexer, insert those missing indices
+    l_idx = np.insert(l_idx, insert_idx, l_idx_missing)
+    # for the right indexer, insert -1 -> to get missing values in pandas' reindexing
+    r_idx = np.insert(r_idx, insert_idx, -1)
+    # for the indices, already insert those missing values manually
+    if distances is not None:
+        distances = np.insert(distances, insert_idx, np.nan)
+
+    if how == "right":
+        # switch back
+        l_idx, r_idx = r_idx, l_idx
+
+    return (l_idx, r_idx), distances
+
+
+def _frame_join(
+    left_df, right_df, indices, distances, how, lsuffix, rsuffix, predicate
+):
     """Join the GeoDataFrames at the DataFrame level.
 
     Parameters
     ----------
-    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
+    indices : tuple of ndarray
+        Indices returned by the geometric join. Tuple with with integer
+        indices representing the matches from `left_df` and `right_df`
+        respectively.
+    distances : ndarray, optional
+        Passed trough and adapted based on the indices, if needed.
+    how : string
+        The type of join to use on the DataFrame level.
     lsuffix : string
         Suffix to apply to overlapping column names (left GeoDataFrame).
     rsuffix : string
         Suffix to apply to overlapping column names (right GeoDataFrame).
-    how : string
-        The type of join to use on the DataFrame level.
 
     Returns
     -------
     GeoDataFrame
         Joined GeoDataFrame.
     """
-    # the spatial index only allows limited (numeric) index types, but an
-    # index in geopandas may be any arbitrary dtype. so reset both indices now
-    # and store references to the original indices, to be reaffixed later.
-    # GH 352
-    left_df = left_df.copy(deep=True)
+    if how in ("inner", "left"):
+        right_df = right_df.drop(right_df.geometry.name, axis=1)
+    else:  # how == 'right':
+        left_df = left_df.drop(left_df.geometry.name, axis=1)
+
+    left_df = left_df.copy(deep=False)
     left_nlevels = left_df.index.nlevels
     left_index_original = left_df.index.names
     left_df, left_column_names = _reset_index_with_suffix(left_df, lsuffix, right_df)
 
-    right_df = right_df.copy(deep=True)
+    right_df = right_df.copy(deep=False)
     right_nlevels = right_df.index.nlevels
     right_index_original = right_df.index.names
     right_df, right_column_names = _reset_index_with_suffix(right_df, rsuffix, left_df)
@@ -350,53 +417,25 @@ def _frame_join(join_df, left_df, right_df, how, lsuffix, rsuffix):
     right_index = right_df.columns[:right_nlevels]
 
     # perform join on the dataframes
-    if how == "inner":
-        join_df = join_df.set_index("_key_left")
-        joined = (
-            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",
-                right_index=True,
-                suffixes=(None, None),
-            )
-            .drop(["_key_right"], axis=1)
-        )
-        joined = _restore_index(joined, left_index, left_index_original)
-
-    elif how == "left":
-        join_df = join_df.set_index("_key_left")
-        joined = (
-            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",
-                left_on="_key_right",
-                right_index=True,
-                suffixes=(None, None),
-            )
-            .drop(["_key_right"], axis=1)
-        )
+    original_length = len(right_df) if how == "right" else len(left_df)
+    (l_idx, r_idx), distances = _adjust_indexers(
+        indices, distances, original_length, how, predicate
+    )
+    # the `take` method doesn't allow introducing NaNs with -1 indices
+    # left = left_df.take(l_idx)
+    # therefore we are using the private _reindex_with_indexers as workaround
+    new_index = pd.RangeIndex(len(l_idx))
+    left = left_df._reindex_with_indexers({0: (new_index, l_idx)})
+    right = right_df._reindex_with_indexers({0: (new_index, r_idx)})
+    joined = pd.concat([left, right], axis=1, copy=False)
+
+    if how in ("inner", "left"):
         joined = _restore_index(joined, left_index, left_index_original)
-
     else:  # how == 'right':
-        joined = (
-            left_df.drop(left_df.geometry.name, axis=1)
-            .merge(
-                join_df.merge(
-                    right_df, left_on="_key_right", right_index=True, how="right"
-                ),
-                left_index=True,
-                right_on="_key_left",
-                how="right",
-                suffixes=(None, None),
-            )
-            .drop(["_key_left", "_key_right"], axis=1)
-            .set_geometry(right_df.geometry.name)
-        )
+        joined = joined.set_geometry(right_df.geometry.name)
         joined = _restore_index(joined, right_index, right_index_original)
 
-    return joined
+    return joined, distances
 
 
 def _nearest_query(
@@ -436,15 +475,14 @@ def _nearest_query(
                 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
+        l_idx, r_idx = np.array([], dtype=np.intp), np.array([], dtype=np.intp)
+        if return_distance:
+            distances = np.array([], dtype=np.float64)
+        else:
+            distances = None
+    return (l_idx, r_idx), distances
 
 
 def sjoin_nearest(
@@ -527,9 +565,9 @@ def sjoin_nearest(
 
     >>> groceries_w_communities = geopandas.sjoin_nearest(groceries, chicago)
     >>> groceries_w_communities[["Chain", "community", "geometry"]].head(2)
-                 Chain community                               geometry
-    0   VIET HOA PLAZA    UPTOWN  MULTIPOINT ((1168268.672 1933554.35))
-    87      JEWEL OSCO    UPTOWN  MULTIPOINT ((1168837.98 1929246.962))
+                   Chain    community                                geometry
+    0     VIET HOA PLAZA       UPTOWN   MULTIPOINT ((1168268.672 1933554.35))
+    1  COUNTY FAIR FOODS  MORGAN PARK  MULTIPOINT ((1162302.618 1832900.224))
 
 
     To include the distances:
@@ -537,10 +575,10 @@ def sjoin_nearest(
     >>> groceries_w_communities = geopandas.sjoin_nearest(groceries, chicago, \
 distance_col="distances")
     >>> groceries_w_communities[["Chain", "community", \
-"distances"]].head(2)  # doctest: +SKIP
-                    Chain community  distances
-    0   VIET HOA PLAZA    UPTOWN        0.0
-    87      JEWEL OSCO    UPTOWN        0.0
+"distances"]].head(2)
+                   Chain    community  distances
+    0     VIET HOA PLAZA       UPTOWN        0.0
+    1  COUNTY FAIR FOODS  MORGAN PARK        0.0
 
     In the following example, we get multiple groceries for Uptown because all
     results are equidistant (in this case zero because they intersect).
@@ -550,7 +588,7 @@ distance_col="distances")
 distance_col="distances", how="right")
     >>> uptown_results = \
 chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
-    >>> uptown_results[["Chain", "community"]]  # doctest: +SKIP
+    >>> uptown_results[["Chain", "community"]]
                 Chain community
     30  VIET HOA PLAZA    UPTOWN
     30      JEWEL OSCO    UPTOWN
@@ -577,19 +615,14 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
 
     return_distance = distance_col is not None
 
-    join_df = _nearest_query(
+    indices, distances = _nearest_query(
         left_df, right_df, max_distance, how, return_distance, exclusive
     )
+    joined, distances = _frame_join(
+        left_df, right_df, indices, distances, how, lsuffix, rsuffix, None
+    )
 
     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]
+        joined[distance_col] = distances
 
     return joined


Test patch:
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index ba934c18..058594ef 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -178,6 +178,17 @@ class TestGeomMethods:
         self.l6 = LineString([(5, 5), (5, 100), (100, 5)])
         self.g12 = GeoSeries([self.l5])
         self.g13 = GeoSeries([self.l6])
+        self.g14 = GeoSeries(
+            [
+                MultiLineString([[(0, 2), (0, 10)], [(0, 10), (5, 10)]]),
+                MultiLineString([[(0, 2), (0, 10)], [(0, 11), (5, 10)]]),
+                MultiLineString(),
+                MultiLineString([[(0, 0), (1, 0)], [(0, 0), (3, 0)]]),
+                Point(0, 0),
+            ],
+            crs=4326,
+            index=range(2, 7),
+        )
 
     def _test_unary_real(self, op, expected, a):
         """Tests for 'area', 'length', 'is_valid', etc."""
@@ -299,15 +310,26 @@ class TestGeomMethods:
     #         self._test_binary_topological('intersection', self.g3,
     #                                       self.g3, no_crs_g3)
 
+    def test_alignment_warning(self):
+        with pytest.warns(
+            UserWarning,
+            match="The indices of the left and right GeoSeries' are not equal",
+        ):
+            self.g0.intersection(self.g9, align=None)
+
+        with warnings.catch_warnings(record=True) as record:
+            self.g0.intersection(self.g9, align=True)
+            self.g0.intersection(self.g9, align=False)
+
+            assert len(record) == 0
+
     def test_intersection(self):
         self._test_binary_topological("intersection", self.t1, self.g1, self.g2)
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            self._test_binary_topological(
-                "intersection", self.all_none, self.g1, self.empty
-            )
+        self._test_binary_topological(
+            "intersection", self.all_none, self.g1, self.empty
+        )
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.intersection(self.g9, align=True) == 8)
+        assert len(self.g0.intersection(self.g9, align=True) == 8)
         assert len(self.g0.intersection(self.g9, align=False) == 7)
 
     def test_clip_by_rect(self):
@@ -322,8 +344,7 @@ class TestGeomMethods:
     def test_union_series(self):
         self._test_binary_topological("union", self.sq, self.g1, self.g2)
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.union(self.g9, align=True) == 8)
+        assert len(self.g0.union(self.g9, align=True) == 8)
         assert len(self.g0.union(self.g9, align=False) == 7)
 
     def test_union_polygon(self):
@@ -332,8 +353,7 @@ class TestGeomMethods:
     def test_symmetric_difference_series(self):
         self._test_binary_topological("symmetric_difference", self.sq, self.g3, self.g4)
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.symmetric_difference(self.g9, align=True) == 8)
+        assert len(self.g0.symmetric_difference(self.g9, align=True) == 8)
         assert len(self.g0.symmetric_difference(self.g9, align=False) == 7)
 
     def test_symmetric_difference_poly(self):
@@ -346,8 +366,7 @@ class TestGeomMethods:
         expected = GeoSeries([GeometryCollection(), self.t2])
         self._test_binary_topological("difference", expected, self.g1, self.g2)
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.difference(self.g9, align=True) == 8)
+        assert len(self.g0.difference(self.g9, align=True) == 8)
         assert len(self.g0.difference(self.g9, align=False) == 7)
 
     def test_difference_poly(self):
@@ -480,8 +499,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.contains(self.t1))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.contains(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.contains(self.g9, align=True))
 
         expected = [False, False, True, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.contains(self.g9, align=False))
@@ -491,10 +509,9 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.contains_properly(Point(0.25, 0.25)))
 
         expected = [False, False, False, False, False, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(
-                expected, self.g0.contains_properly(self.g9, align=True)
-            )
+        assert_array_dtype_equal(
+            expected, self.g0.contains_properly(self.g9, align=True)
+        )
 
         expected = [False, False, True, False, False, False, False]
         assert_array_dtype_equal(
@@ -507,8 +524,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.dwithin(self.p0, 6))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=True))
+        assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=True))
         expected = [True, True, True, True, False, False, False]
         assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=False))
 
@@ -539,8 +555,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.crossed_lines.crosses(self.l3))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.crosses(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.crosses(self.g9, align=True))
 
         expected = [False] * 7
         assert_array_dtype_equal(expected, self.g0.crosses(self.g9, align=False))
@@ -550,8 +565,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.disjoint(self.t1))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.disjoint(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.disjoint(self.g9, align=True))
 
         expected = [False, False, False, False, True, False, False]
         assert_array_dtype_equal(expected, self.g0.disjoint(self.g9, align=False))
@@ -588,8 +602,7 @@ class TestGeomMethods:
             index=range(8),
         )
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.relate(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.relate(self.g9, align=True))
 
         expected = Series(
             [
@@ -615,8 +628,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g6.distance(self.na_none))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.distance(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.distance(self.g9, align=True))
 
         val = self.g0.iloc[4].distance(self.g9.iloc[4])
         expected = Series(np.array([0, 0, 0, 0, val, np.nan, np.nan]), self.g0.index)
@@ -636,10 +648,9 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.na_none.hausdorff_distance(self.p0))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(
-                expected, self.g0.hausdorff_distance(self.g9, align=True)
-            )
+        assert_array_dtype_equal(
+            expected, self.g0.hausdorff_distance(self.g9, align=True)
+        )
 
         val_1 = self.g0.iloc[0].hausdorff_distance(self.g9.iloc[0])
         val_2 = self.g0.iloc[2].hausdorff_distance(self.g9.iloc[2])
@@ -665,10 +676,9 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.na_none.frechet_distance(self.p0))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(
-                expected, self.g0.frechet_distance(self.g9, align=True)
-            )
+        assert_array_dtype_equal(
+            expected, self.g0.frechet_distance(self.g9, align=True)
+        )
 
         # expected returns
         val_1 = 1.0
@@ -704,8 +714,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.intersects(self.empty_poly))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.intersects(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.intersects(self.g9, align=True))
 
         expected = [True, True, True, True, False, False, False]
         assert_array_dtype_equal(expected, self.g0.intersects(self.g9, align=False))
@@ -718,8 +727,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g4.overlaps(self.t1))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.overlaps(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.overlaps(self.g9, align=True))
 
         expected = [False] * 7
         assert_array_dtype_equal(expected, self.g0.overlaps(self.g9, align=False))
@@ -729,8 +737,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.touches(self.t1))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.touches(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.touches(self.g9, align=True))
 
         expected = [True, False, False, True, False, False, False]
         assert_array_dtype_equal(expected, self.g0.touches(self.g9, align=False))
@@ -743,8 +750,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.within(self.sq))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.within(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.within(self.g9, align=True))
 
         expected = [False, True, False, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.within(self.g9, align=False))
@@ -761,8 +767,7 @@ class TestGeomMethods:
         assert_series_equal(res, exp)
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.covers(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.covers(self.g9, align=True))
 
         expected = [False, False, True, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.covers(self.g9, align=False))
@@ -778,8 +783,7 @@ class TestGeomMethods:
         assert_series_equal(res, exp)
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.covered_by(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.covered_by(self.g9, align=True))
 
         expected = [False, True, False, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.covered_by(self.g9, align=False))
@@ -1092,8 +1096,7 @@ class TestGeomMethods:
 
         s = GeoSeries([Point(2, 2), Point(0.5, 0.5)], index=[1, 2])
         expected = Series([np.nan, 2.0, np.nan])
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_series_equal(self.g5.project(s), expected)
+        assert_series_equal(self.g5.project(s), expected)
 
         expected = Series([2.0, 0.5], index=self.g5.index)
         assert_series_equal(self.g5.project(s, align=False), expected)
@@ -1802,6 +1805,35 @@ class TestGeomMethods:
         )
         assert_geoseries_equal(expected, self.g1.force_3d([1, 2]))
 
+    def test_line_merge(self):
+        expected = GeoSeries(
+            [
+                LineString([(0, 2), (0, 10), (5, 10)]),
+                MultiLineString([[(0, 2), (0, 10)], [(0, 11), (5, 10)]]),
+                GeometryCollection(),
+                LineString([(0, 0), (1, 0), (3, 0)]),
+                GeometryCollection(),
+            ],
+            crs=4326,
+            index=range(2, 7),
+        )
+        assert_geoseries_equal(expected, self.g14.line_merge())
+
+    @pytest.mark.skipif(shapely.geos_version < (3, 11, 0), reason="requires GEOS>=3.11")
+    def test_line_merge_directed(self):
+        expected = GeoSeries(
+            [
+                LineString([(0, 2), (0, 10), (5, 10)]),
+                MultiLineString([[(0, 2), (0, 10)], [(0, 11), (5, 10)]]),
+                GeometryCollection(),
+                MultiLineString([[(0, 0), (1, 0)], [(0, 0), (3, 0)]]),
+                GeometryCollection(),
+            ],
+            crs=4326,
+            index=range(2, 7),
+        )
+        assert_geoseries_equal(expected, self.g14.line_merge(directed=True))
+
     @pytest.mark.skipif(
         shapely.geos_version < (3, 9, 5), reason="Empty geom bug in GEOS<3.9.5"
     )
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index b9f81338..99b8336b 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -110,11 +110,11 @@ class TestSeries:
         # Test that warning is issued when operating on non-aligned series
 
         # _series_op
-        with pytest.warns(UserWarning, match="The indices .+ different"):
+        with pytest.warns(UserWarning, match="The indices .+ not equal"):
             self.a1.contains(self.a2)
 
         # _geo_op
-        with pytest.warns(UserWarning, match="The indices .+ different"):
+        with pytest.warns(UserWarning, match="The indices .+ not equal"):
             self.a1.union(self.a2)
 
     def test_no_warning_if_aligned(self):
@@ -136,8 +136,7 @@ class TestSeries:
         assert_array_equal(self.g1.geom_equals(self.sq), [False, True])
 
     def test_geom_equals_align(self):
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            a = self.a1.geom_equals(self.a2, align=True)
+        a = self.a1.geom_equals(self.a2, align=True)
         exp = pd.Series([False, True, False], index=["A", "B", "C"])
         assert_series_equal(a, exp)
 
@@ -153,7 +152,7 @@ class TestSeries:
         with warnings.catch_warnings():
             warnings.filterwarnings(
                 "ignore",
-                "The indices of the two GeoSeries are different",
+                "The indices of the left and right GeoSeries' are not equal",
                 UserWarning,
             )
             assert_array_equal(
@@ -170,7 +169,9 @@ class TestSeries:
         assert_array_equal(self.g1.geom_equals_exact(self.sq, 0.001), [False, True])
         with warnings.catch_warnings():
             warnings.filterwarnings(
-                "ignore", "The indices of the two GeoSeries are different", UserWarning
+                "ignore",
+                "The indices of the left and right GeoSeries' are not equal",
+                UserWarning,
             )
             assert_array_equal(
                 self.a1.geom_equals_exact(self.a2, 0.001, align=True),
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 754d924a..f7542a96 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -704,20 +704,12 @@ def test_groupby_metadata(crs, geometry_name):
         **kwargs,
     )
 
-    if compat.PANDAS_GE_22:
-        # merge sort behaviour changed in pandas #54611
-        take_indices = [0, 0, 2, 2, 1]
-        value_right = [0, 2, 0, 2, 1]
-    else:
-        take_indices = [0, 2, 0, 2, 1]
-        value_right = [0, 0, 2, 2, 1]
-
     expected = (
-        df.take(take_indices)
+        df.take([0, 0, 2, 2, 1])
         .set_index("value2", drop=compat.PANDAS_GE_22, append=True)
         .swaplevel()
         .rename(columns={"value1": "value1_left"})
-        .assign(value1_right=value_right)
+        .assign(value1_right=[0, 2, 0, 2, 1])
     )
     assert_geodataframe_equal(res.drop(columns=["index_right"]), expected)
 
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index a7eaf839..59419fc3 100644
--- a/geopandas/tools/tests/test_sjoin.py
+++ b/geopandas/tools/tests/test_sjoin.py
@@ -574,7 +574,6 @@ class TestIndexNames:
         with pytest.raises(ValueError, match="'index_right' cannot be a column name"):
             sjoin(df1, df2, how=how)
 
-    @pytest.mark.xfail
     @pytest.mark.parametrize("how", ["inner", "left", "right"])
     def test_conflicting_column_with_suffix(self, how):
         # test case where the auto-generated index name conflicts
@@ -608,6 +607,7 @@ class TestIndexNames:
                 "column_right",
                 "geometry",
             ]
+        expected = expected.set_geometry("geometry")
         assert_geodataframe_equal(result, expected)
 
 
@@ -1148,12 +1148,8 @@ class TestNearest:
         result5["index_right"] = result5["index_right"].astype("int64")
         assert_geodataframe_equal(result5, result4, check_like=True)
 
-    expected_index_uncapped = (
-        [1, 3, 3, 1, 2] if compat.PANDAS_GE_22 else [1, 1, 3, 3, 2]
-    )
-
     @pytest.mark.parametrize(
-        "max_distance,expected", [(None, expected_index_uncapped), (1.1, [3, 3, 1, 2])]
+        "max_distance,expected", [(None, [1, 3, 3, 1, 2]), (1.1, [3, 3, 1, 2])]
     )
     def test_sjoin_nearest_exclusive(self, max_distance, expected):
         geoms = shapely.points(np.arange(3), np.arange(3))