DistrictHeatingNetworkCreator.py

@@ -1,31 +1,36 @@
 import json
-import geopandas as gpd
 import matplotlib.pyplot as plt
-from shapely.geometry import Polygon, Point, LineString, MultiPoint
+from shapely.geometry import Polygon, Point, LineString
 import networkx as nx
+from typing import List, Tuple
+from rtree import index
 
 
 class DistrictHeatingNetworkCreator:
-    def __init__(self, buildings_file, roads_file):
+    def __init__(self, buildings_file: str, roads_file: str):
         """
         Initialize the class with paths to the buildings and roads data files.
 
         :param buildings_file: Path to the GeoJSON file containing building data.
-        :param roads_file: Path to the Shapefile containing road data.
+        :param roads_file: Path to the GeoJSON file containing roads data.
         """
         self.buildings_file = buildings_file
         self.roads_file = roads_file
 
-    def run(self):
+    def run(self) -> nx.Graph:
         """
         Main method to execute the district heating network creation process.
         :return: NetworkX graph with nodes and edges representing the network.
         """
         self._load_and_process_data()
         self._find_nearest_roads()
-        self._process_intersections()
-        network_graph = self._create_network_graph()
-        return network_graph
+        self._find_nearest_points()
+        self._break_down_roads()
+        self._create_graph()
+        self._create_mst()
+        self._iteratively_remove_edges()
+        self.add_centroids_to_mst()  # Add centroids to the MST
+        return self.final_mst
 
     def _load_and_process_data(self):
         """
@@ -35,156 +40,211 @@ class DistrictHeatingNetworkCreator:
         with open(self.buildings_file, 'r') as file:
             city = json.load(file)
 
-        # Extract centroids and building IDs from building data
-        centroids = []
-        building_ids = []  # List to store building IDs
+        self.centroids = []
+        self.building_ids = []
         buildings = city['features']
         for building in buildings:
             coordinates = building['geometry']['coordinates'][0]
             building_polygon = Polygon(coordinates)
             centroid = building_polygon.centroid
-            centroids.append(centroid)
-            building_ids.append(building['id'])  # Extract building ID
-
-        # Convert centroids to a GeoDataFrame and include building IDs
-        self.centroids_gdf = gpd.GeoDataFrame({
-            'geometry': [Point(centroid.x, centroid.y) for centroid in centroids],
-            'building_id': building_ids  # Add building IDs as a column
-        }, crs='EPSG:4326')
+            self.centroids.append(centroid)
+            self.building_ids.append(building['id'])
 
         # Load road data
-        self.gdf_road = gpd.read_file(self.roads_file)
+        with open(self.roads_file, 'r') as file:
+            roads = json.load(file)
 
-        # Ensure centroids are in the same CRS as roads
-        self.centroids_gdf = self.centroids_gdf.to_crs(self.gdf_road.crs)
+        line_features = [feature for feature in roads['features'] if feature['geometry']['type'] == 'LineString']
+
+        self.lines = [LineString(feature['geometry']['coordinates']) for feature in line_features]
+        self.cleaned_lines = [LineString([line.coords[0], line.coords[-1]]) for line in self.lines]
 
     def _find_nearest_roads(self):
         """
         Find the nearest road for each building centroid.
         """
-        # Process road geometries
-        self.gdf_clean = gpd.GeoDataFrame(
-            {'geometry': [LineString([coord for coord in line.coords]) for line in self.gdf_road.geometry]})
+        self.closest_roads = []
+        unique_roads_set = set()
+
+        # Create spatial index for roads
+        idx = index.Index()
+        for pos, line in enumerate(self.cleaned_lines):
+            idx.insert(pos, line.bounds)
+
+        for centroid in self.centroids:
+            min_distance = float('inf')
+            closest_road = None
+            for pos in idx.nearest(centroid.bounds, 10):
+                road = self.cleaned_lines[pos]
+                distance = road.distance(centroid)
+                if distance < min_distance:
+                    min_distance = distance
+                    closest_road = road
+
+            if closest_road and closest_road.wkt not in unique_roads_set:
+                unique_roads_set.add(closest_road.wkt)
+                self.closest_roads.append(closest_road)
+
+    def _find_nearest_points(self):
+        """
+        Find the nearest point on each closest road for each centroid.
+        """
+        def find_nearest_point_on_line(point: Point, line: LineString) -> Point:
+            return line.interpolate(line.project(point))
 
-        # Find the nearest road line and point for each centroid
-        self.closest_linestrings = []
         self.nearest_points = []
-        for centroid in self.centroids_gdf.geometry:
-            closest_road = min(self.gdf_clean.geometry, key=lambda x: x.distance(centroid))
-            self.closest_linestrings.append(closest_road)
-            nearest_point = closest_road.interpolate(closest_road.project(centroid))
-            self.nearest_points.append(nearest_point)
+        for centroid in self.centroids:
+            min_distance = float('inf')
+            closest_road = None
+            for road in self.closest_roads:
+                distance = centroid.distance(road)
+                if distance < min_distance:
+                    min_distance = distance
+                    closest_road = road
 
-    def _process_intersections(self):
+            if closest_road:
+                nearest_point = find_nearest_point_on_line(centroid, closest_road)
+                self.nearest_points.append(nearest_point)
+
+    def _break_down_roads(self):
         """
-        Process intersections and create final geometries.
+        Break down roads into segments connecting nearest points.
         """
-        # Create additional GeoDataFrames for points and nearest points
-        self.gdf_pts = gpd.GeoDataFrame(
-            {'geometry': [Point(coord) for line in self.gdf_clean.geometry for coord in line.coords]})
-        self.gdf_pts2 = gpd.GeoDataFrame({'geometry': self.nearest_points})
+        def break_down_roads(closest_roads: List[LineString], nearest_points_list: List[Point]) -> List[LineString]:
+            new_segments = []
+            for road in closest_roads:
+                coords = list(road.coords)
+                points_on_road = [point for point in nearest_points_list if road.distance(point) < 0.000000001]
+                sorted_points = sorted(points_on_road, key=lambda point: road.project(point))
+                sorted_points.insert(0, Point(coords[0]))
+                sorted_points.append(Point(coords[-1]))
+                for i in range(len(sorted_points) - 1):
+                    segment = LineString([sorted_points[i], sorted_points[i + 1]])
+                    new_segments.append(segment)
+            return new_segments
 
-        # Combine nearest points and road points into one GeoDataFrame
-        self.gdf_pts3 = gpd.GeoDataFrame({'geometry': self.nearest_points + list(self.gdf_pts.geometry)})
+        self.new_segments = break_down_roads(self.closest_roads, self.nearest_points)
+        self.cleaned_lines = [line for line in self.cleaned_lines if line not in self.closest_roads]
+        self.cleaned_lines.extend(self.new_segments)
 
-        # Identify intersections and create LineStrings based on intersections
-        self.gdf_clean["intersect"] = [
-            [y for y in range(len(self.gdf_pts2)) if self.gdf_pts2.geometry[y].distance(geom) <= 1.0] for geom in
-            self.gdf_clean.geometry]
-        self.gdf_cleaner = self.gdf_clean[self.gdf_clean["intersect"].apply(len).gt(0)].reset_index(drop=True)
-
-        self.test_list = []
-        for idx, row in self.gdf_cleaner.iterrows():
-            for i in range(len(row["intersect"]) + 1):
-                if i == 0:
-                    self.test_list.append(
-                        LineString([row['geometry'].coords[0], self.gdf_pts3.geometry[row['intersect'][i]]]))
-                elif i < len(row['intersect']):
-                    self.test_list.append(LineString(
-                        [self.gdf_pts3.geometry[row['intersect'][i - 1]], self.gdf_pts3.geometry[row['intersect'][i]]]))
-                else:
-                    self.test_list.append(
-                        LineString([self.gdf_pts3.geometry[row['intersect'][i - 1]], row['geometry'].coords[-1]]))
-
-        self.gdf_cleanest = gpd.GeoDataFrame({'geometry': self.test_list})
-
-        points = [coord for geom in self.gdf_cleanest.geometry for coord in geom.coords]
-        gdf_pts_cnt = self.gdf_pts.copy()
-        gdf_pts_cnt["count"] = gdf_pts_cnt.geometry.apply(lambda x: points.count(x.coords[0]))
-        gdf_pts_reset = gdf_pts_cnt[gdf_pts_cnt["count"] > 1].reset_index(drop=True)
-        gdf_pts_drop = gdf_pts_cnt[gdf_pts_cnt["count"] <= 1].reset_index(drop=True)
-
-        # Remove unnecessary geometries from gdf_cleanest
-        for idx, geom in self.gdf_cleanest.iterrows():
-            for coord in geom.geometry.coords:
-                if coord in [pt.coords[0] for pt in gdf_pts_drop.geometry]:
-                    self.gdf_cleanest = self.gdf_cleanest.drop(idx)
-
-        self.gdf_cleanest.reset_index(drop=True, inplace=True)
-
-    def _create_network_graph(self):
+    def _create_graph(self):
         """
-        Create a NetworkX graph from the processed geospatial data.
-        :return: A NetworkX graph representing the district heating network.
+        Create a NetworkX graph from the cleaned lines.
         """
-        G = nx.Graph()
+        self.G = nx.Graph()
+        for line in self.cleaned_lines:
+            coords = list(line.coords)
+            for i in range(len(coords) - 1):
+                self.G.add_edge(coords[i], coords[i + 1], weight=Point(coords[i]).distance(Point(coords[i + 1])))
 
-        # Convert centroids to EPSG:4326 for Google Maps compatibility
-        for idx, row in self.centroids_gdf.iterrows():
-            building_name = f"Building_{idx + 1}"
-            G.add_node((row.geometry.x, row.geometry.y),
-                       type='centroid',
-                       name=building_name,
-                       building_id=row['building_id'])  # Add building ID as an attribute
+    def _create_mst(self):
+        """
+        Create a Minimum Spanning Tree (MST) from the graph.
+        """
+        def find_paths_between_nearest_points(g: nx.Graph, nearest_points: List[Point]) -> List[Tuple]:
+            edges = []
+            for i, start_point in enumerate(nearest_points):
+                start = (start_point.x, start_point.y)
+                for end_point in nearest_points[i + 1:]:
+                    end = (end_point.x, end_point.y)
+                    if nx.has_path(g, start, end):
+                        path = nx.shortest_path(g, source=start, target=end, weight='weight')
+                        path_edges = list(zip(path[:-1], path[1:]))
+                        edges.extend((u, v, g[u][v]['weight']) for u, v in path_edges)
+            return edges
 
-        for point in self.nearest_points:
-            G.add_node((point.x, point.y), type='nearest_point')
+        edges = find_paths_between_nearest_points(self.G, self.nearest_points)
+        h = nx.Graph()
+        h.add_weighted_edges_from(edges)
+        mst = nx.minimum_spanning_tree(h, weight='weight')
+        final_edges = []
+        for u, v in mst.edges():
+            if nx.has_path(self.G, u, v):
+                path = nx.shortest_path(self.G, source=u, target=v, weight='weight')
+                path_edges = list(zip(path[:-1], path[1:]))
+                final_edges.extend((x, y, self.G[x][y]['weight']) for x, y in path_edges)
+        self.final_mst = nx.Graph()
+        self.final_mst.add_weighted_edges_from(final_edges)
 
-        # Add edges with lengths as weights for the road network
-        for line in self.gdf_cleanest.geometry:
-            length = line.length
-            if isinstance(line.boundary, MultiPoint):
-                # Handle MultiPoint boundary
-                points = list(line.boundary.geoms)
-                for i in range(len(points) - 1):
-                    start_point = points[i]
-                    end_point = points[i + 1]
-                    G.add_edge((start_point.x, start_point.y), (end_point.x, end_point.y), weight=length)
-            else:
-                # Handle typical case with two endpoints
-                start_point, end_point = line.boundary
-                G.add_edge((start_point.x, start_point.y), (end_point.x, end_point.y), weight=length)
+    def _iteratively_remove_edges(self):
+        """
+        Iteratively remove edges that do not have any nearest points and have one end with only one connection.
+        Also remove nodes that don't have any connections and street nodes with only one connection.
+        """
+        nearest_points_tuples = [(point.x, point.y) for point in self.nearest_points]
 
-        # Add edges connecting nearest points to their centroids
-        for point, centroid in zip(self.nearest_points, self.centroids_gdf.geometry):
-            distance = point.distance(centroid)
-            G.add_edge((point.x, point.y), (centroid.x, centroid.y), weight=distance)
+        def find_edges_to_remove(graph: nx.Graph) -> List[Tuple]:
+            edges_to_remove = []
+            for u, v, d in graph.edges(data=True):
+                if u not in nearest_points_tuples and v not in nearest_points_tuples:
+                    if graph.degree(u) == 1 or graph.degree(v) == 1:
+                        edges_to_remove.append((u, v, d))
+            return edges_to_remove
 
-        return G
+        def find_nodes_to_remove(graph: nx.Graph) -> List[Tuple]:
+            nodes_to_remove = []
+            for node in graph.nodes():
+                if graph.degree(node) == 0:
+                    nodes_to_remove.append(node)
+            return nodes_to_remove
 
-    def plot_network_graph(self, network_graph):
+        edges_to_remove = find_edges_to_remove(self.final_mst)
+        self.final_mst_steps = [list(self.final_mst.edges(data=True))]
+
+        while edges_to_remove or find_nodes_to_remove(self.final_mst):
+            self.final_mst.remove_edges_from(edges_to_remove)
+            nodes_to_remove = find_nodes_to_remove(self.final_mst)
+            self.final_mst.remove_nodes_from(nodes_to_remove)
+            edges_to_remove = find_edges_to_remove(self.final_mst)
+            self.final_mst_steps.append(list(self.final_mst.edges(data=True)))
+
+        def find_single_connection_street_nodes(graph: nx.Graph) -> List[Tuple]:
+            single_connection_street_nodes = []
+            for node in graph.nodes():
+                if node not in nearest_points_tuples and graph.degree(node) == 1:
+                    single_connection_street_nodes.append(node)
+            return single_connection_street_nodes
+
+        single_connection_street_nodes = find_single_connection_street_nodes(self.final_mst)
+
+        while single_connection_street_nodes:
+            for node in single_connection_street_nodes:
+                neighbors = list(self.final_mst.neighbors(node))
+                self.final_mst.remove_node(node)
+                for neighbor in neighbors:
+                    if self.final_mst.degree(neighbor) == 0:
+                        self.final_mst.remove_node(neighbor)
+            single_connection_street_nodes = find_single_connection_street_nodes(self.final_mst)
+            self.final_mst_steps.append(list(self.final_mst.edges(data=True)))
+
+    def add_centroids_to_mst(self):
+        """
+        Add centroids to the final MST graph and connect them to their associated node on the graph.
+        """
+        for centroid in self.centroids:
+            centroid_tuple = (centroid.x, centroid.y)
+            self.final_mst.add_node(centroid_tuple, type='centroid')
+            nearest_point = None
+            min_distance = float('inf')
+            for node in self.final_mst.nodes():
+                if self.final_mst.nodes[node].get('type') != 'centroid':
+                    node_point = Point(node)
+                    distance = centroid.distance(node_point)
+                    if distance < min_distance:
+                        min_distance = distance
+                        nearest_point = node
+
+            if nearest_point:
+                self.final_mst.add_edge(centroid_tuple, nearest_point, weight=min_distance)
+
+    def plot_network_graph(self):
         """
         Plot the network graph using matplotlib and networkx.
-
-        :param network_graph: The NetworkX graph to be plotted.
         """
-        plt.figure(figsize=(12, 12))
-        pos = {node: (node[0], node[1]) for node in network_graph.nodes()}
-
-        # Draw nodes and edges
-        nx.draw_networkx_nodes(network_graph, pos, node_color='blue', node_size=50)
-        nx.draw_networkx_edges(network_graph, pos, edge_color='gray')
-
-        # Create a dictionary for node labels for centroids only
-        node_labels = {node: data['name'] for node, data in network_graph.nodes(data=True) if
-                       data.get('type') == 'centroid'}
-
-        # Adjust node label positions to reduce overlap
-        label_pos = {node: (coords[0], coords[1] + 0.03) for node, coords in pos.items()}  # Shift labels up
-
-        # Draw node labels for centroids
-        nx.draw_networkx_labels(network_graph, label_pos, labels=node_labels, font_size=8, verticalalignment='bottom')
-
+        plt.figure(figsize=(15, 10))
+        pos = {node: (node[0], node[1]) for node in self.final_mst.nodes()}
+        nx.draw_networkx_nodes(self.final_mst, pos, node_color='blue', node_size=50)
+        nx.draw_networkx_edges(self.final_mst, pos, edge_color='gray')
         plt.title('District Heating Network Graph')
         plt.axis('off')
         plt.show()

Scripts/road_processor.py

@@ -0,0 +1,56 @@
+from pathlib import Path
+from shapely.geometry import Polygon, Point, shape
+import json
+
+
+def road_processor(x, y, diff):
+    """
+    Processes a .JSON file to find roads that have at least one node within a specified polygon.
+
+    Parameters:
+    x (float): The x-coordinate of the center of the selection box.
+    y (float): The y-coordinate of the center of the selection box.
+    diff (float): The half-width of the selection box.
+
+    Returns:
+    str: The file path of the output GeoJSON file containing the selected roads.
+    """
+    diff += 2 * diff
+    # Define the selection polygon
+    selection_box = Polygon([
+        [x + diff, y - diff],
+        [x - diff, y - diff],
+        [x - diff, y + diff],
+        [x + diff, y + diff]
+    ])
+
+    # Define input and output file paths
+    geojson_file = Path("./input_files/roads.json").resolve()
+    output_file = Path('./input_files/output_roads.geojson').resolve()
+
+    # Initialize a list to store the roads in the region
+    roads_in_region = []
+
+    # Read the GeoJSON file
+    with open(geojson_file, 'r') as file:
+        roads = json.load(file)
+        line_features = [feature for feature in roads['features'] if feature['geometry']['type'] == 'LineString']
+
+        # Check each road feature
+        for feature in line_features:
+            road_shape = shape(feature['geometry'])
+            # Check if any node of the road is inside the selection box
+            if any(selection_box.contains(Point(coord)) for coord in road_shape.coords):
+                roads_in_region.append(feature)
+
+    # Create a new GeoJSON structure with the selected roads
+    output_geojson = {
+        "type": "FeatureCollection",
+        "features": roads_in_region
+    }
+
+    # Write the selected roads to the output file
+    with open(output_file, 'w') as outfile:
+        json.dump(output_geojson, outfile)
+
+    return output_file

input_files/roads/geobase_mtl.cpg

@@ -1 +0,0 @@
-UTF-8

input_files/roads/geobase_mtl.prj

@@ -1 +0,0 @@
-PROJCS["NAD_1983_MTM_8",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",304800.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-73.5],PARAMETER["Scale_Factor",0.9999],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

main.py

@@ -1,13 +1,19 @@
 from DistrictHeatingNetworkCreator import DistrictHeatingNetworkCreator
-# building_file = "./input_files/buildings.geojson"
-# road_file = "./input_files/roads/geobase_mtl.shp"
-# Initialize the class
-network_creator = DistrictHeatingNetworkCreator(
-    './input_files/buildings.geojson',
-    './input_files/roads/geobase_mtl.shp'
-)
+from Scripts.road_processor import road_processor
+from pathlib import Path
+import time
+location = [45.51663850312751, -73.59854314961274]
+start_time = time.perf_counter()
+roads_file = road_processor(location[1], location[0], 0.001)
 
-# Create the network graph
-network_graph = network_creator.run()
+buildings_file = Path('./input_files/buildings.geojson').resolve()
 
-network_creator.plot_network_graph(network_graph)
+dhn_creator = DistrictHeatingNetworkCreator(buildings_file, roads_file)
+network_graph = dhn_creator.run()
+end_time = time.perf_counter()
+
+elapsed_time = end_time - start_time
+
+print(f"The simulation took {elapsed_time:.4f} seconds to run.")
+
+dhn_creator.plot_network_graph()