`consolidate_intersections` crashes with dictionary for tolerances

EwoutH commented 4 months ago

Contributing guidelines

[X] I understand the contributing guidelines

Documentation

[X] My problem is not addressed by the documentation or examples

Existing issues

[X] My problem does not appear in an existing issue

What operating system and Python version are you using?

Windows 11, Python 3.12

What OSMnx version are you using?

2.0.0b0

Environment packages and versions

OSMnx version: 2.0.0b0
Shapely version: 2.0.4
GeoPandas version: 1.0.0

How did you install OSMnx?

Pip

Problem description

When using a dictionary for tolerances the consolidate_intersections() function crashes, while it works with a single tolerance value (replace tolerance=tolerance_dict with tolerance=1 in the example below).

See for context:

1150
1160

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[16], line 41
     38 tolerance_dict = {node: tolerance for node, tolerance in tolerance_dict.items() if node in G.nodes}
     40 # Consolidate intersections
---> 41 merged_network = ox.consolidate_intersections(G, tolerance=tolerance_dict, rebuild_graph=True)
     43 # Print the results
     44 print("The simplified combined network has", len(merged_network.nodes), "nodes and", len(merged_network.edges), "edges")

File ~\.virtualenvs\Py312\Lib\site-packages\osmnx\simplification.py:537, in consolidate_intersections(G, tolerance, rebuild_graph, dead_ends, reconnect_edges, node_attr_aggs)
    534         return G
    536     # otherwise
--> 537     return _consolidate_intersections_rebuild_graph(
    538         G,
    539         tolerance,
    540         reconnect_edges,
    541         node_attr_aggs,
    542     )
    544 # otherwise, if we're not rebuilding the graph
    545 if len(G) == 0:
    546     # if graph has no nodes, just return empty GeoSeries

File ~\.virtualenvs\Py312\Lib\site-packages\osmnx\simplification.py:640, in _consolidate_intersections_rebuild_graph(G, tolerance, reconnect_edges, node_attr_aggs)
    635     node_attr_aggs = {"elevation": np.mean}
    637 # STEP 1
    638 # buffer nodes to passed-in distance and merge overlaps. turn merged nodes
    639 # into gdf and get centroids of each cluster as x, y
--> 640 node_clusters = gpd.GeoDataFrame(geometry=_merge_nodes_geometric(G, tolerance))
    641 centroids = node_clusters.centroid
    642 node_clusters["x"] = centroids.x

File ~\.virtualenvs\Py312\Lib\site-packages\osmnx\simplification.py:583, in _merge_nodes_geometric(G, tolerance)
    578 if isinstance(tolerance, dict):
    579     # create series of tolerances reindexed like nodes, then buffer, then
    580     # fill nulls (resulting from missing tolerances) with original points,
    581     # then merge overlapping geometries
    582     tols = pd.Series(tolerance).reindex(gdf_nodes.index)
--> 583     merged = gdf_nodes.buffer(tols).fillna(gdf_nodes["geometry"]).unary_union
    584 else:
    585     # buffer nodes then merge overlapping geometries
    586     merged = gdf_nodes.buffer(tolerance).unary_union

File ~\.virtualenvs\Py312\Lib\site-packages\geopandas\base.py:5046, in GeoPandasBase.buffer(self, distance, resolution, cap_style, join_style, mitre_limit, single_sided, **kwargs)
   4972 def buffer(
   4973     self,
   4974     distance,
   (...)
   4980     **kwargs,
   4981 ):
   4982     """Returns a ``GeoSeries`` of geometries representing all points within
   4983     a given ``distance`` of each geometric object.
   4984 
   (...)
   5044 
   5045     """
-> 5046     return _delegate_geo_method(
   5047         "buffer",
   5048         self,
   5049         distance=distance,
   5050         resolution=resolution,
   5051         cap_style=cap_style,
   5052         join_style=join_style,
   5053         mitre_limit=mitre_limit,
   5054         single_sided=single_sided,
   5055         **kwargs,
   5056     )

File ~\.virtualenvs\Py312\Lib\site-packages\geopandas\base.py:117, in _delegate_geo_method(op, this, **kwargs)
    107         raise ValueError(
    108             f"Index of the Series passed as '{key}' does not match index of the "
    109             f"{klass}. If you want both Series to be aligned, align them before "
   (...)
    113             f"using `{key}.values`."
    114         )
    116 a_this = GeometryArray(this.geometry.values)
--> 117 data = getattr(a_this, op)(**kwargs)
    118 return GeoSeries(data, index=this.index, crs=this.crs)

File ~\.virtualenvs\Py312\Lib\site-packages\geopandas\array.py:811, in GeometryArray.buffer(self, distance, resolution, **kwargs)
    809 if not (isinstance(distance, (int, float)) and distance == 0):
    810     self.check_geographic_crs(stacklevel=5)
--> 811 return GeometryArray(
    812     shapely.buffer(self._data, distance, quad_segs=resolution, **kwargs),
    813     crs=self.crs,
    814 )

File ~\.virtualenvs\Py312\Lib\site-packages\geopandas\array.py:323, in GeometryArray.__init__(self, data, crs)
    321     data = data._data
    322 elif not isinstance(data, np.ndarray):
--> 323     raise TypeError(
    324         "'data' should be array of geometry objects. Use from_shapely, "
    325         "from_wkb, from_wkt functions to construct a GeometryArray."
    326     )
    327 elif not data.ndim == 1:
    328     raise ValueError(
    329         "'data' should be a 1-dimensional array of geometry objects."
    330     )

TypeError: 'data' should be array of geometry objects. Use from_shapely, from_wkb, from_wkt functions to construct a GeometryArray.

Complete minimal reproducible example

import networkx as nx
import osmnx as ox
import geopandas as gpd
from shapely.geometry import Point

# Create a sample graph
G = nx.MultiDiGraph()

# Add nodes with attributes including network type and coordinates
nodes = [
    (1, {"network": "CityNetwork", "x": 0, "y": 0, "geometry": Point(0, 0)}),
    (2, {"network": "CityNetwork", "x": 1, "y": 1, "geometry": Point(1, 1)}),
    (3, {"network": "SurroundingNetwork", "x": 2, "y": 2, "geometry": Point(2, 2)}),
    (4, {"network": "SurroundingNetwork", "x": 3, "y": 3, "geometry": Point(3, 3)})
]

G.add_nodes_from(nodes)

# Add edges
edges = [(1, 2), (2, 3), (3, 4)]
G.add_edges_from(edges)

# Add CRS to the graph (yes this doesn't make sense for a graph, but it's just for the example)
G.graph['crs'] = 3857

# Create the tolerance dictionary
city_name = "CityNetwork"
tolerance_dict = {node: 10 if data["network"] == city_name else 50 for node, data in G.nodes(data=True)}
print(f"The tolerance dictionary is: {tolerance_dict}")

# Ensure the tolerance_dict only includes nodes present in the graph
tolerance_dict = {node: tolerance for node, tolerance in tolerance_dict.items() if node in G.nodes}

# Consolidate intersections
merged_network = ox.consolidate_intersections(G, tolerance=tolerance_dict, rebuild_graph=True)

# Print the results
print("The simplified combined network has", len(merged_network.nodes), "nodes and", len(merged_network.edges), "edges")

EwoutH commented 4 months ago

What's the weirdest thing is that I have this notebook tracked in git, and it was running perfectly just a few weeks ago (May 6): https://github.com/EwoutH/urban-self-driving-effects/blob/main/network/create_network.ipynb

I even explicitly tested it with this exact osmnx 2.0.0b0 version. So it almost has to have been an underlying dependency that got updated. Shapely hasn't had a release since April 17th, and the only other option I can maybe think of is geopandas, which had their 1.0 update last week.

EwoutH commented 4 months ago

Wow it looks like it is geopandas! With geopandas 0.14.4 it actually runs perfectly fine.

martinfleis commented 4 months ago

Yeah, sorry! This is a regression in geopandas 1.0 that'll be fixed in 1.0.1 in a few days.

EwoutH commented 4 months ago

I can confirm this is issue was fixed by https://github.com/geopandas/geopandas/issues/3362, thanks a lot @martinfleis!

gboeing / osmnx