Enhance OSM getter script

requirements.txt

@@ -9,3 +9,4 @@ Pillow
 seaborn
 tqdm
 opencv-python
+shapely

src/dsm/dsm.hpp

@@ -6,7 +6,7 @@
 
 static constexpr uint8_t DSM_VERSION_MAJOR = 2;
 static constexpr uint8_t DSM_VERSION_MINOR = 2;
-static constexpr uint8_t DSM_VERSION_PATCH = 8;
+static constexpr uint8_t DSM_VERSION_PATCH = 9;
 
 static auto const DSM_VERSION =
     std::format("{}.{}.{}", DSM_VERSION_MAJOR, DSM_VERSION_MINOR, DSM_VERSION_PATCH);

src/dsm/headers/Graph.cpp

@@ -305,8 +305,8 @@ namespace dsm {
         std::string id, lat, lon, highway;
         // osmid;x;y;highway
         std::getline(iss, id, ';');
-        std::getline(iss, lat, ';');
         std::getline(iss, lon, ';');
+        std::getline(iss, lat, ';');
         std::getline(iss, highway, ';');
         auto const nodeId{static_cast<Id>(std::stoul(id))};
         if (highway.find("traffic_signals") != std::string::npos) {
@@ -470,7 +470,10 @@ namespace dsm {
   void Graph::addStreet(std::unique_ptr<Street> street) {
     if (m_streets.contains(street->id())) {
       throw std::invalid_argument(
-          buildLog(std::format("Street with id {} already exists.", street->id())));
+          buildLog(std::format("Street with id {} from {} to {} already exists.",
+                               street->id(),
+                               street->nodePair().first,
+                               street->nodePair().second)));
     }
     // emplace nodes
     const auto srcId{street->nodePair().first};
@@ -488,7 +491,10 @@ namespace dsm {
   void Graph::addStreet(const Street& street) {
     if (m_streets.contains(street.id())) {
       throw std::invalid_argument(
-          buildLog(std::format("Street with id {} already exists.", street.id())));
+          buildLog(std::format("Street with id {} from {} to {} already exists.",
+                               street.id(),
+                               street.nodePair().first,
+                               street.nodePair().second)));
     }
     // emplace nodes
     const auto srcId{street.nodePair().first};

utils/get_osm_data.py

@@ -3,7 +3,7 @@
 The place is passed as a command line argument.
 
 Example:
-python get_osm_data.py --place "postua, vercelli, italy"
+python get_osm_data.py --place "Bologna, Emilia-Romagna, Italy" --exclude-residential True
 
 The output files are:
 - nodes.csv
@@ -13,12 +13,16 @@
 """
 
 from argparse import ArgumentParser
+import logging
 import osmnx as ox
 import networkx as nx
 import matplotlib.patches as mpatches
 import matplotlib.pyplot as plt
+import pandas as pd
+from shapely.geometry import MultiLineString, LineString
+from shapely.ops import linemerge
 
-RGBA_RED = (1, 0, 0, 0.1)
+RGBA_RED = (1, 0, 0, 0.3)
 RGBA_WHITE = (1, 1, 1, 1)
 
 FLAGS_MOTORWAY = ["motorway", "motorway_link"]
@@ -36,8 +40,176 @@
     "residential",
     "living_street",
     "unclassified",
+    "service",
+    "pedestrian",
+    "busway",
 ]
 
+
+def merge_edges(
+    graph: nx.DiGraph, previous_node: int, successive_node: int, node: int
+) -> dict:
+    """
+    Merge two edges into a single edge.
+    The function merges the edges into a single edge if the following conditions are met:
+    - the name of the two edges is the same
+    - the number of lanes is the same
+    - the geometry of the two edges is contiguous
+    - the coordinates of the previous_node and successive_node are in the geometry
+
+    Parameters:
+    ----------------
+    graph (networkx.DiGraph): the graph
+    previous_node (int): the previous node
+    successive_node (int): the successive node
+    node (int): the id of the node which will be removed
+
+    Returns:
+    ----------------
+    dict: the new edge
+    """
+    try:
+        data_u = graph.get_edge_data(previous_node, node)[0]
+        data_v = graph.get_edge_data(node, successive_node)[0]
+        data_u.setdefault("lanes", 1)
+        data_v.setdefault("lanes", 1)
+        if (
+            not (data_u["name"] in data_v["name"] or data_v["name"] in data_u["name"])
+            or data_u["lanes"] != data_v["lanes"]
+        ):
+            return None
+        edge_uv = data_u.copy()
+        # set length as the sum
+        edge_uv["length"] = data_u["length"] + data_v["length"]
+        edge_uv["lanes"] = int(data_u["lanes"])
+        # merge also linestrings
+        edge_uv["geometry"] = data_u["geometry"].union(data_v["geometry"])
+        if isinstance(edge_uv["geometry"], MultiLineString):
+            edge_uv["geometry"] = linemerge(edge_uv["geometry"])
+        else:
+            edge_uv["geometry"] = edge_uv["geometry"]
+        if isinstance(edge_uv["geometry"], LineString):
+            coords = list(edge_uv["geometry"].coords)
+        else:
+            # If it's still a MultiLineString,
+            # handle it by iterating through its individual LineStrings
+            coords = []
+            for line in edge_uv["geometry"]:
+                coords.extend(
+                    list(line.coords)
+                )  # Add coords from each individual LineString
+        # take the list from coordinates of previous_node to coordinates of successive_node
+        u_coords = (graph.nodes[previous_node]["x"], graph.nodes[previous_node]["y"])
+        v_coords = (
+            graph.nodes[successive_node]["x"],
+            graph.nodes[successive_node]["y"],
+        )
+        if u_coords not in coords or v_coords not in coords:
+            return None
+        # cut coords from u_index to v_index
+        edge_uv["geometry"] = LineString(
+            coords[coords.index(u_coords) : coords.index(v_coords)]
+        )
+    except TypeError:
+        # type error means that data_u or data_v cannot be created,
+        # which means that the road is a one-way road
+        # thus, skip the type error
+        return None
+
+    return edge_uv
+
+
+def simplify_graph(graph_original: nx.DiGraph) -> nx.DiGraph:
+    """
+    Simplify the graph by removing nodes that have only two neighborsand are actually the same
+    street.
+    The function merges the edges into a single edge.
+
+    Parameters:
+    ----------------
+    graph_original (networkx.DiGraph): the graph to simplify
+
+    Returns:
+    ----------------
+    networkx.DiGraph: the simplified graph
+    """
+    graph = graph_original.copy()
+    previous_nodes = 0
+    while previous_nodes != len(graph.nodes):
+        logging.info("New cycle: current_nodes=%d", len(graph.nodes))
+        previous_nodes = len(graph.nodes)
+        for node in graph.copy().nodes:
+            # define neighborus as list of predecessors and successors
+            neighbours = list(graph.predecessors(node)) + list(graph.successors(node))
+            if (
+                len(neighbours) != 2
+                or graph.in_degree(node) != graph.out_degree(node)
+                or graph.in_degree(node) > 2
+            ):
+                continue
+            u, v = neighbours
+            if graph.has_edge(u, v):
+                continue
+            edge_uv = merge_edges(graph, u, v, node)
+            edge_vu = merge_edges(graph, v, u, node)
+
+            if not (edge_uv is None and edge_vu is None):
+                if edge_uv is not None:
+                    # print(f"Edges {u} -> {node} and {node} -> {v} can be merged.")
+                    graph.add_edge(
+                        u,
+                        v,
+                        length=edge_uv["length"],
+                        name=edge_uv["name"],
+                        geometry=edge_uv["geometry"],
+                    )
+                    # print(f"Added edge {graph.get_edge_data(u, v)}")
+                if edge_vu is not None:
+                    # print(f"Edges {v} -> {node} and {node} -> {u} can be merged.")
+                    graph.add_edge(
+                        v,
+                        u,
+                        length=edge_vu["length"],
+                        name=edge_vu["name"],
+                        geometry=edge_vu["geometry"],
+                    )
+                    # print(f"Added edge {graph.get_edge_data(v, u)}")
+                graph.remove_node(node)
+                # print(f"Removed node {node}")
+
+    # Remove all nodes that are not in the giant component
+    graph.remove_nodes_from(
+        set(graph.nodes) - max(list(nx.weakly_connected_components(graph)), key=len)
+    )
+    # remove all self-loops
+    graph.remove_edges_from(list(nx.selfloop_edges(graph)))
+    # check if there are edges with same u and v. If true, keep only the one with the bigger lanes
+    edges_to_remove = []
+    seen_edges = {}
+
+    for u, v, data in graph.edges(data=True):
+        lanes = data.get("lanes", 0)
+
+        if (u, v) not in seen_edges:
+            seen_edges[(u, v)] = (lanes, None)  # Store first edge and its lanes count
+        else:
+            existing_lanes, existing_edge = seen_edges[(u, v)]
+
+            if lanes > existing_lanes:
+                edges_to_remove.append(
+                    existing_edge
+                )  # Remove the previous edge if the current one has more lanes
+                seen_edges[(u, v)] = (lanes, (u, v))  # Update to keep current edge
+            else:
+                edges_to_remove.append(
+                    (u, v)
+                )  # Remove the current edge if it has fewer 'lanes'
+
+    graph.remove_edges_from(edges_to_remove)
+
+    return graph
+
+
 if __name__ == "__main__":
     parser = ArgumentParser("Script to get the OSM data of a place.")
     parser.add_argument(
@@ -58,15 +230,36 @@
         help="Exclude residential roads from the data",
     )
     parser = parser.parse_args()
+    logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
+    # set up colored logging
+    logging.addLevelName(
+        logging.INFO, f"\033[1;32m{logging.getLevelName(logging.INFO)}\033[1;0m"
+    )
+    logging.addLevelName(
+        logging.WARNING, f"\033[1;33m{logging.getLevelName(logging.WARNING)}\033[1;0m"
+    )
+    logging.addLevelName(
+        logging.ERROR, f"\033[1;31m{logging.getLevelName(logging.ERROR)}\033[1;0m"
+    )
 
     # get the street network for San Cesario sul Panaro
     G_ALL = ox.graph_from_place(parser.place, network_type="drive")
-    print(f"Graph created with {len(G_ALL.nodes)} nodes and {len(G_ALL.edges)} edges.")
+    logging.info(
+        "Graph created with %d nodes and %d edges.", len(G_ALL.nodes), len(G_ALL.edges)
+    )
 
     gdf_nodes, gdf_edges = ox.graph_to_gdfs(G_ALL)
     gdf_edges["highway"] = gdf_edges["highway"].apply(
         lambda x: x[-1] if isinstance(x, list) else x
     )
+    if "lanes" not in gdf_edges.columns:
+        gdf_edges["lanes"] = 1
+    gdf_edges["lanes"] = gdf_edges["lanes"].apply(
+        lambda x: max(x) if isinstance(x, list) else 1 if pd.isna(x) else x
+    )
+    gdf_edges["name"] = gdf_edges["name"].apply(
+        lambda x: " ".join(x) if isinstance(x, list) else " " if pd.isna(x) else x
+    )
     # gdf_edges = gdf_edges[~gdf_edges["access"].isin(["no", "private"])]
 
     # Make a plot to visualize the removed links
@@ -80,7 +273,19 @@
     # rebuild the graph
     G = ox.graph_from_gdfs(gdf_nodes, gdf_edges)
     G.remove_nodes_from(list(nx.isolates(G)))
-    print(f"Graph filterded: now it has {len(G.nodes)} nodes and {len(G.edges)} edges.")
+    logging.info(
+        "Graph filtered: now it has %d nodes and %d edges.", len(G.nodes), len(G.edges)
+    )
+    G = simplify_graph(G)
+    logging.info(
+        "Graph simplified: now it has %d nodes and %d edges.",
+        len(G.nodes),
+        len(G.edges),
+    )
+    # assert that graph has not isolated nodes
+    assert not list(nx.isolates(G))
+    # assert that graph has not self-loops
+    assert not list(nx.selfloop_edges(G))
 
     fig, ax = ox.plot_graph(
         G_ALL,
@@ -94,23 +299,26 @@
         close=False,
     )
     ax.legend(handles=[removed_patch])
-    plt.show()
+    fig.set_size_inches(16, 9)
+    plt.savefig("removed_nodes_and_edges.png")
 
-    gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)
+    # Plot resulting graph
+    fig, ax = ox.plot_graph(G, show=False, close=False)
+    fig.set_size_inches(16, 9)
+    plt.savefig("final_graph.png")
 
+    gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)
     # notice that osmnid is the index of the gdf_nodes DataFrame, so take it as a column
     gdf_nodes.reset_index(inplace=True)
     gdf_edges.reset_index(inplace=True)
+
+    # assert that there are no edges with the same u and v
+    assert not gdf_edges.duplicated(subset=["u", "v"]).any()
     # Prepare node dataframe
     gdf_nodes = gdf_nodes[["osmid", "x", "y", "highway"]]
     # Prepare edge dataframe
-    if "lanes" not in gdf_edges.columns:
-        gdf_edges["lanes"] = None
-
-    gdf_edges["lanes"] = gdf_edges["lanes"].apply(
-        lambda x: max(x) if isinstance(x, list) else x
-    )
 
+    gdf_edges.to_csv("edges_ALL.csv", sep=";", index=False)
     gdf_edges = gdf_edges[
         ["u", "v", "length", "oneway", "lanes", "highway", "maxspeed", "name"]
     ]