Repository: pysal/momepy

Issue:
remove_false_nodes throws an unexpected error
### Describe the problem

When trying to remove false nodes from a [street network obtained from Ordnance Survey](https://www.ordnancesurvey.co.uk/products/os-open-roads) I get the following error:

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[56], line 1
----> 1 streets_BR = mm.remove_false_nodes(streets_BR)
      2 len(streets_BR)

File [~\AppData\Local\miniconda3\envs\ox\Lib\site-packages\momepy\preprocessing.py:260](http://localhost:8889/lab/tree/OneDrive%20-%20University%20of%20Strathclyde/FB/Birmingham/_scripts/~/AppData/Local/miniconda3/envs/ox/Lib/site-packages/momepy/preprocessing.py#line=259), in remove_false_nodes(gdf)
    256 node_coords = shapely.get_coordinates(target_nodes)
    257 coords = np.array(loop_geom.coords)
    258 new_start = (
    259     node_coords[0]
--> 260     if (node_coords[0] != coords[0]).all()
    261     else node_coords[1]
    262 )
    263 new_start_idx = np.where(coords == new_start)[0][0]
    264 rolled_coords = np.roll(coords[:-1], -new_start_idx, axis=0)

ValueError: operands could not be broadcast together with shapes (2,) (3,)
```

Any solution?

### Steps to reproduce

```
import geopandas as gpd
import momepy as mm
import pandas as pd

sk = gpd.read_file('data/OS Open Roads/SK_RoadLink.shp')
so = gpd.read_file('data/OS Open Roads/SO_RoadLink.shp')
sp = gpd.read_file('data/OS Open Roads/SP_RoadLink.shp')

streets = pd.concat([sk, so, sp], ignore_index=True)

streets = streets.reset_index(drop=True).drop_duplicates('identifier')

bound_BR = gpd.read_file('data/boundary.shp')

streets_BR = streets[streets.within(bound_BR)]

streets_BR = mm.remove_false_nodes(streets_BR)
```

### Versions of your packages

momepy = 0.8.0
geopandas = 1.0.1
libpysal = 4.12.1
networkx = 3.3

### Your operating system

Windows 10

### Additional context

_No response_

Discussion:
@jGaboardi: ~~This seems to be happening in `momepy` actually so I'll transfer this issue over to there.~~

I have no idea why I thought this was opened in `tobler`... LOL@jGaboardi: ~~Seems I don't have the proper permissions for the transfer.~~

Couldn't transfer because it's already here...@martinfleis: Can you try using [`force_2d`](https://geopandas.org/en/latest/docs/reference/api/geopandas.GeoSeries.force_2d.html) prior using the function?@ale-v: Thanks. It solved the issue!@martinfleis: You probably had mixed LineString and LineString Z geometries after the concatenation. 

I am curious if we should do some check internally here to avoid these cryptic errors. Just to provide an informative error rather than forcing 2 or 3D. @jGaboardi thoughts?@jGaboardi: Yes, more checks & informative messages are always great ideas.

@ale-v -- Any chance you can tease out the subset of offending data? We can use that in a test.@martinfleis: That would be a pain with licensing. We can just call force_3d on half of the data to reproduce this. @jGaboardi: Boiled down example this is ***not*** reproducible:

```python
Python 3.12.7 | packaged by conda-forge | (main, Oct  4 2024, 15:57:01) [Clang 17.0.6 ]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.28.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import geopandas, momepy, numpy, shapely
   ...: 
   ...: line_1 = shapely.LineString((shapely.Point(1, 1), shapely.Point(2, 2)))
   ...: line_2 = shapely.LineString((shapely.Point(2, 2), shapely.Point(3, 3)))
   ...: 
   ...: line_2d_input = numpy.array([line_1, line_2])
   ...: line_3d_input = numpy.array([line_1, shapely.force_3d(line_2)])
   ...: 
   ...: gdf_2d_input = geopandas.GeoDataFrame(geometry=line_2d_input)
   ...: gdf_3d_input = geopandas.GeoDataFrame(geometry=line_3d_input)
   ...: 
   ...: 
   ...: (
   ...:     line_1,
   ...:     line_2,
   ...:     line_2d_input,
   ...:     line_3d_input,
   ...:     momepy.remove_false_nodes(line_2d_input),
   ...:     momepy.remove_false_nodes(line_3d_input),
   ...:     momepy.remove_false_nodes(gdf_2d_input),
   ...:     momepy.remove_false_nodes(gdf_3d_input),
   ...: )

Out[1]: 

(<LINESTRING (1 1, 2 2)>,
 <LINESTRING (2 2, 3 3)>,
 array([<LINESTRING (1 1, 2 2)>, <LINESTRING (2 2, 3 3)>], dtype=object),
 array([<LINESTRING (1 1, 2 2)>, <LINESTRING Z (2 2 0, 3 3 0)>],
       dtype=object),
 0    LINESTRING (1 1, 2 2, 3 3)
 dtype: geometry,
 0    LINESTRING (1 1, 2 2, 3 3)
 dtype: geometry,
                      geometry
 0  LINESTRING (1 1, 2 2, 3 3),
                      geometry
 0  LINESTRING (1 1, 2 2, 3 3))
```@jGaboardi: Ideas about what to do different here to reproduce?@martinfleis: It is in the code handling and rotating loops, so we need a loop that is formed by 3D geometry. Because `shapely.get_coordinates` ignores Z by default, while `geom.coords` does not. So it needs a loop of some sort.@jGaboardi: Good call. The following reproduces the error:

```python
In [1]: import momepy, numpy, shapely
   ...: 
   ...: line_1 = shapely.LineString((shapely.Point(1, 1), shapely.Point(2, 2)))
   ...: line_2 = shapely.LineString((shapely.Point(2, 2), shapely.Point(3, 3)))
   ...: line_3 = shapely.LineString((shapely.Point(2, 2), shapely.Point(1, 2)))
   ...: line_4 = shapely.LineString((shapely.Point(1, 2), shapely.Point(1, 3)))
   ...: line_5 = shapely.LineString((shapely.Point(1, 3), shapely.Point(2, 2)))
   ...: 
   ...: line_2d_input = numpy.array([line_1, line_2, line_3, line_4, line_5])
   ...: momepy.remove_false_nodes(line_2d_input)
   ...: 
   ...: line_3d_input = numpy.array([line_1, line_2, shapely.force_3d(line_3), line_4, line_5])
   ...: momepy.remove_false_nodes(line_3d_input)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[1], line 14
     11 momepy.remove_false_nodes(line_2d_input)
     13 line_3d_input = numpy.array([line_1, line_2, shapely.force_3d(line_3), line_4, line_5])
---> 14 momepy.remove_false_nodes(line_3d_input)

File ~/github_repos/momepy/momepy/preprocessing.py:260, in remove_false_nodes(gdf)
    256 node_coords = shapely.get_coordinates(target_nodes)
    257 coords = np.array(loop_geom.coords)
    258 new_start = (
    259     node_coords[0]
--> 260     if (node_coords[0] != coords[0]).all()
    261     else node_coords[1]
    262 )
    263 new_start_idx = np.where(coords == new_start)[0][0]
    264 rolled_coords = np.roll(coords[:-1], -new_start_idx, axis=0)

ValueError: operands could not be broadcast together with shapes (2,) (3,) 
```@jGaboardi: > You probably had mixed LineString and LineString Z geometries after the concatenation.
> 
> I am curious if we should do some check internally here to avoid these cryptic errors. Just to provide an informative error rather than forcing 2 or 3D. @jGaboardi thoughts?

Pushing up a proposed fix shortly.

Solution patch:
diff --git a/momepy/preprocessing.py b/momepy/preprocessing.py
index 196331a4..839cdca6 100644
--- a/momepy/preprocessing.py
+++ b/momepy/preprocessing.py
@@ -166,8 +166,8 @@ def preprocess(
 def remove_false_nodes(gdf):
     """
     Clean topology of existing LineString geometry by removal of nodes of degree 2.
-
-    Returns the original gdf if there's no node of degree 2.
+    Returns the original gdf if there's no node of degree 2. Some geometries may
+    be forced to 2D where a Z coordinate is present.
 
     Parameters
     ----------
@@ -178,7 +178,7 @@ def remove_false_nodes(gdf):
     -------
     gdf : GeoDataFrame, GeoSeries
 
-    See also
+    See Also
     --------
     momepy.extend_lines
     momepy.close_gaps
@@ -247,14 +247,16 @@ def remove_false_nodes(gdf):
             ),
             lines=False,
         )
+
         loops = combined[combined.is_ring]
+
         node_ix, loop_ix = loops.sindex.query(nodes.geometry, predicate="intersects")
         for ix in np.unique(loop_ix):
             loop_geom = loops.geometry.iloc[ix]
             target_nodes = nodes.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)
+                coords = shapely.get_coordinates(loop_geom)
                 new_start = (
                     node_coords[0]
                     if (node_coords[0] != coords[0]).all()


Test patch:
diff --git a/momepy/tests/test_preprocess.py b/momepy/tests/test_preprocess.py
index 98f8ee44..01803374 100644
--- a/momepy/tests/test_preprocess.py
+++ b/momepy/tests/test_preprocess.py
@@ -2,7 +2,7 @@
 import numpy as np
 import pytest
 from geopandas.testing import assert_geodataframe_equal
-from shapely import affinity
+from shapely import affinity, force_3d
 from shapely.geometry import LineString, MultiPoint, Point, Polygon
 from shapely.ops import polygonize
 
@@ -67,6 +67,16 @@ def test_remove_false_nodes(self):
         df = self.df_streets.drop([4, 7, 17, 22])
         assert_geodataframe_equal(df, mm.remove_false_nodes(df))
 
+        # check 3d coords in loop
+        line_1 = LineString((Point(1, 1), Point(2, 2)))
+        line_2 = LineString((Point(2, 2), Point(3, 3)))
+        line_3 = force_3d(LineString((Point(2, 2), Point(1, 2))))
+        line_4 = LineString((Point(1, 2), Point(1, 3)))
+        line_5 = LineString((Point(1, 3), Point(2, 2)))
+        with_3d = np.array([line_1, line_2, line_3, line_4, line_5])
+        no_3d = mm.remove_false_nodes(with_3d)
+        assert len(no_3d) == 3
+
     def test_CheckTessellationInput(self):
         df = self.df_buildings
         df.loc[144, "geometry"] = Polygon([(0, 0), (0, 1), (1, 0)])