diff --git a/DistrictHeatingNetworkCreator.py b/DistrictHeatingNetworkCreator.py index b436d22..9959c97 100644 --- a/DistrictHeatingNetworkCreator.py +++ b/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') diff --git a/Scripts/road_processor.py b/Scripts/road_processor.py index 1bc324a..486ef07 100644 --- a/Scripts/road_processor.py +++ b/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. diff --git a/main.py b/main.py index d9c7327..98c78c0 100644 --- a/main.py +++ b/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()