BUG: fix ordering of closed loops in remove_false_nodes

momepy/preprocessing.py

@@ -19,8 +19,9 @@
 from tqdm.auto import tqdm
 
 from .coins import COINS
+from .graph import node_degree
 from .shape import CircularCompactness
-from .utils import nx_to_gdf
+from .utils import gdf_to_nx, nx_to_gdf
 
 __all__ = [
     "preprocess",
@@ -221,7 +222,7 @@
         df = df.drop(merge_res)
         final = gpd.GeoSeries(new_geoms, crs=df.crs).explode(ignore_index=True)
         if isinstance(gdf, gpd.GeoDataFrame):
-            return pd.concat(
+            combined = pd.concat(
                 [
                     df,
                     gpd.GeoDataFrame(
@@ -231,7 +232,59 @@
                 ignore_index=True,
             )
         else:
-            return pd.concat([df, final], ignore_index=True)
+            combined = pd.concat([df, final], ignore_index=True)
+
+        # re-order closed loops
+        fixed_loops = []
+        fixed_index = []
+        nodes = nx_to_gdf(
+            node_degree(
+                gdf_to_nx(
+                    combined
+                    if isinstance(combined, gpd.GeoDataFrame)
+                    else combined.to_frame("geometry")
+                )
+            ),
+            lines=False,
+        )
+        degree2 = nodes[nodes.degree == 2]
+        loops = combined[combined.is_ring]
+        node_ix, loop_ix = loops.sindex.query(degree2.geometry, predicate="intersects")
+        for ix in np.unique(loop_ix):
+            loop_geom = loops.geometry.iloc[ix]
+            target_nodes = degree2.geometry.iloc[node_ix[loop_ix == ix]]
+            if len(target_nodes) == 2:
+                node_coords = shapely.get_coordinates(target_nodes)
+                coords = np.array(loop_geom.coords)
+                new_start = (
+                    node_coords[0]
+                    if (node_coords[0] != coords[0]).all()
+                    else node_coords[1]
+                )
+                new_start_idx = np.where(coords == new_start)[0][0]
+                rolled_coords = np.roll(coords[:-1], -new_start_idx, axis=0)
+                new_sequence = np.append(rolled_coords, rolled_coords[[0]], axis=0)
+                fixed_loops.append(shapely.LineString(new_sequence))
+                fixed_index.append(ix)
+        fixed_loops = gpd.GeoSeries(fixed_loops, crs=df.crs).explode(ignore_index=True)
+
+        if isinstance(gdf, gpd.GeoDataFrame):
+            return pd.concat(
+                [
+                    combined.drop(loops.iloc[fixed_index].index),
+                    gpd.GeoDataFrame(
+                        {df.geometry.name: fixed_loops},
+                        geometry=df.geometry.name,
+                        crs=df.crs,
+                    ),
+                ],
+                ignore_index=True,
+            )
+        else:
+            return pd.concat(
+                [combined.drop(loops.iloc[fixed_index].index), fixed_loops],
+                ignore_index=True,
+            )
 
     # if there's nothing to fix, return the original dataframe
     return gdf

momepy/tests/test_preprocess.py

@@ -39,6 +39,20 @@ def test_remove_false_nodes(self):
         assert isinstance(fixed, gpd.GeoDataFrame)
         assert self.false_network.crs.equals(fixed.crs)
         assert sorted(self.false_network.columns) == sorted(fixed.columns)
+
+        # check loop order
+        expected = np.array(
+            [
+                [-727238.49292668, -1052817.28071986],
+                [-727253.1752498, -1052827.47329062],
+                [-727223.93217677, -1052829.47624082],
+                [-727238.49292668, -1052817.28071986],
+            ]
+        )
+        np.testing.assert_almost_equal(
+            np.array(fixed.loc[53].geometry.coords), expected
+        )
+
         fixed_series = mm.remove_false_nodes(self.false_network.geometry)
         assert len(fixed_series) == 56
         assert isinstance(fixed_series, gpd.GeoSeries)