Repository: geopandas/geopandas

Issue:
ENH: Line_merge
#### Is your feature request related to a problem?

I would like to merge multilinestrings into linestrings, that is using line_merge from GEOS.

#### Describe the solution you'd like

Line_merge is implemented in pygeos and shapely, but not in GeoPandas

#### API breaking implications

None

#### Describe alternatives you've considered

I'm using GeoPandas.apply(shapely.linemerge)

Thanks

Discussion:
theroggy: resolves #1711

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ab8cbb04..7393cf5a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,7 @@ API changes:
 
 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 +42,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 +68,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/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..4691fe1e 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -1606,6 +1606,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
     #
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/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..7fe0ef62 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."""
@@ -1802,6 +1813,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_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))