DistrictHeatingNetworkCreator.py

@@ -1,7 +1,6 @@
 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
 
 
@@ -11,7 +10,7 @@ class DistrictHeatingNetworkCreator:
         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
@@ -23,9 +22,12 @@ class DistrictHeatingNetworkCreator:
         """
         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()
+        return self.final_mst
 
     def _load_and_process_data(self):
         """
@@ -36,154 +38,199 @@ class DistrictHeatingNetworkCreator:
             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 = []  # List to store 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'])  # Extract 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']
+
+        # Create a list of LineString objects and their properties
+        self.lines = []
+        for feature in line_features:
+            # Create a LineString from coordinates
+            linestring = LineString(feature['geometry']['coordinates'])
+            self.lines.append(linestring)
+
+        self.cleaned_lines = []
+        for line in self.lines:
+            coords = list(line.coords)
+            cleaned_line = LineString([coords[0], coords[-1]])
+            self.cleaned_lines.append(cleaned_line)
 
     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()
+
+        # Loop through each centroid
+        for centroid in self.centroids:
+            min_distance = float('inf')  # Start with a large number to ensure any real distance is smaller
+            closest_road = None
+
+            # Loop through each road and calculate the distance to the current centroid
+            for line in self.cleaned_lines:
+                distance = line.distance(centroid)
+                # Check if the current road is closer than the ones previously checked
+                if distance < min_distance:
+                    min_distance = distance
+                    closest_road = line
+
+            # Add the closest road to the list if it's not already added
+            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, line):
+            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)
 
-    def _process_intersections(self):
+        # Find the nearest point on each closest road for each centroid
+        for centroid in self.centroids:
+            # Find the closest road for this centroid
+            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
+
+            # Find the nearest point on the closest road
+            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})
 
-        # Combine nearest points and road points into one GeoDataFrame
-        self.gdf_pts3 = gpd.GeoDataFrame({'geometry': self.nearest_points + list(self.gdf_pts.geometry)})
+        def break_down_roads(closest_roads, nearest_points_list):
+            new_segments = []
+            for road in closest_roads:
+                # Get coordinates of the road
+                coords = list(road.coords)
+                # Find all nearest points for this road
+                points_on_road = [point for point in nearest_points_list if road.distance(point) < 0.000000001]
+                # Sort nearest points along the road
+                sorted_points = sorted(points_on_road, key=lambda point: road.project(point))
+                # Add the start node to the sorted points
+                sorted_points.insert(0, Point(coords[0]))
+                # Add the end node to the sorted points
+                sorted_points.append(Point(coords[-1]))
+                # Create new segments
+                for i in range(len(sorted_points) - 1):
+                    segment = LineString([sorted_points[i], sorted_points[i + 1]])
+                    new_segments.append(segment)
+            return 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)
+        # Create new segments
+        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)
 
-        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()
 
-        # 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
+        # Add edges to the graph from the cleaned lines
+        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])))
 
-        for point in self.nearest_points:
-            G.add_node((point.x, point.y), type='nearest_point')
+    def _create_mst(self):
+        """
+        Create a Minimum Spanning Tree (MST) from the graph.
+        """
 
-        # 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 find_paths_between_nearest_points(g, nearest_points):
+            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
 
-        # 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)
+        # Find the edges used to connect the nearest points
+        edges = find_paths_between_nearest_points(self.G, self.nearest_points)
 
-        return G
+        # Create a graph from these edges
+        h = nx.Graph()
+        h.add_weighted_edges_from(edges)
 
-    def plot_network_graph(self, network_graph):
+        # Compute the Minimum Spanning Tree (MST) using Kruskal's algorithm
+        mst = nx.minimum_spanning_tree(h, weight='weight')
+
+        # Perform pathfinding again on the MST to ensure shortest paths within the MST
+        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)
+
+        # Create the final MST graph with these edges
+        self.final_mst = nx.Graph()
+        self.final_mst.add_weighted_edges_from(final_edges)
+
+    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.
+        """
+        nearest_points_tuples = [(point.x, point.y) for point in self.nearest_points]
+
+        def find_edges_to_remove(graph):
+            edges_to_remove = []
+            for u, v in graph.edges():
+                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))
+            return edges_to_remove
+
+        edges_to_remove = find_edges_to_remove(self.final_mst)
+
+        while edges_to_remove:
+            self.final_mst.remove_edges_from(edges_to_remove)
+            # Find and remove nodes with no connections
+            nodes_to_remove = [node for node in self.final_mst.nodes() if self.final_mst.degree(node) == 0]
+            self.final_mst.remove_nodes_from(nodes_to_remove)
+            edges_to_remove = find_edges_to_remove(self.final_mst)
+
+    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()}
+        plt.figure(figsize=(15, 10))
+        pos = {node: (node[0], node[1]) for node in self.final_mst.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')
+        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')

Scripts/road_processor.py

@@ -5,7 +5,7 @@ import json
 
 def road_processor(x, y, diff):
     """
-    Processes a GeoJSON file to find roads that have at least one node within a specified polygon.
+    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.

main.py

@@ -1,13 +1,12 @@
 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.json'
-)
+from Scripts.road_processor import road_processor
+from pathlib import Path
 
-# Create the network graph
-network_graph = network_creator.run()
 
-network_creator.plot_network_graph(network_graph)
+roads_file = road_processor(-73.61038745537758, 45.484399882086215, 0.001)
+
+buildings_file = Path('./input_files/buildings.geojson').resolve()
+
+dhn_creator = DistrictHeatingNetworkCreator(buildings_file, roads_file)
+network_graph = dhn_creator.run()
+dhn_creator.plot_network_graph()