file dump

DistrictHeatingNetworkCreator.py

@@ -3,58 +3,18 @@ import geopandas as gpd
 import matplotlib.pyplot as plt
 from shapely.geometry import Polygon, Point, LineString, MultiPoint
 import networkx as nx
-
-
-def plot_network_graph(network_graph):
-    """
-    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
-    nx.draw_networkx_nodes(network_graph, pos, node_color='blue', node_size=50)
-
-    # Draw edges
-    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.title('District Heating Network Graph')
-    plt.axis('off')
-    plt.show()
+from scipy.spatial import cKDTree
+import hub.helpers.constants as cte
 
 
 class DistrictHeatingNetworkCreator:
     def __init__(self, buildings_file, roads_file, central_plant_longitude, central_plant_latitude):
-        """
-        Initialize the class with paths to the buildings and roads data files, and central plant coordinates.
-
-        :param buildings_file: Path to the GeoJSON file containing building data.
-        :param roads_file: Path to the Shapefile containing road data.
-        :param central_plant_longitude: Longitude of the central plant.
-        :param central_plant_latitude: Latitude of the central plant.
-        """
         self.buildings_file = buildings_file
         self.roads_file = roads_file
         self.central_plant_longitude = central_plant_longitude
         self.central_plant_latitude = central_plant_latitude
 
     def run(self):
-        """
-        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()
@@ -62,19 +22,12 @@ class DistrictHeatingNetworkCreator:
         return network_graph
 
     def _load_and_process_data(self):
-        """
-        Load and process the building and road data, and add central plant node.
-        """
-        # Load road data
         self.gdf_road = gpd.read_file(self.roads_file)
-
-        # Load building data
         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
+        building_ids = []
         for building in city['features']:
             coordinates = building['geometry']['coordinates'][0]
             building_polygon = Polygon(coordinates)
@@ -85,25 +38,17 @@ class DistrictHeatingNetworkCreator:
         centroids.append(Point(self.central_plant_longitude, self.central_plant_latitude))
         building_ids.append(1)
 
-        # 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,
-            'type': ['centroid' for _ in centroids]  # Add type for centroids
+            'type': ['centroid' for _ in centroids]
         }, crs='EPSG:4326')
 
     def _find_nearest_roads(self):
-        """
-        Find the nearest road for each building centroid.
-        """
-        # Ensure centroids are in the same CRS as roads
         self.centroids_gdf = self.centroids_gdf.to_crs(self.gdf_road.crs)
-
-        # Process road geometries
         self.gdf_clean = gpd.GeoDataFrame(
             {'geometry': [LineString([coord for coord in line.coords]) for line in self.gdf_road.geometry]})
 
-        # Find the nearest road line and point for each centroid
         self.closest_linestrings = []
         self.nearest_points = []
         for centroid in self.centroids_gdf.geometry:
@@ -113,18 +58,12 @@ class DistrictHeatingNetworkCreator:
             self.nearest_points.append(nearest_point)
 
     def _process_intersections(self):
-        """
-        Process intersections and create final geometries.
-        """
-        # 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)})
 
-        # Identify intersections and create LineStrings based on intersections
         intersects = []
         for geom in self.gdf_clean.geometry:
             intersecting_points = []
@@ -159,7 +98,6 @@ class DistrictHeatingNetworkCreator:
         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]:
@@ -168,24 +106,22 @@ class DistrictHeatingNetworkCreator:
         self.gdf_cleanest.reset_index(drop=True, inplace=True)
 
     def _create_network_graph(self):
-        """
-        Create a NetworkX graph from the processed geospatial data.
-        :return: A NetworkX graph representing the district heating network.
-        """
         g = nx.Graph()
 
-        # Convert centroids to EPSG:4326 for Google Maps compatibility
+        # Add nodes with geometry attribute
         for idx, row in self.centroids_gdf.iterrows():
             building_name = f"Building_{idx}"
             g.add_node((row.geometry.x, row.geometry.y),
+                       geometry=row.geometry,  # Include geometry attribute
                        type='centroid',
                        name=building_name,
                        building_id=str(row.get('building_id')))
 
         for point in self.nearest_points:
-            g.add_node((point.x, point.y), type='nearest_point')
+            g.add_node((point.x, point.y),
+                       geometry=point,  # Include geometry attribute
+                       type='nearest_point')
 
-        # 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):
@@ -198,9 +134,54 @@ class DistrictHeatingNetworkCreator:
                 start_point, end_point = line.boundary
                 g.add_edge((start_point.x, start_point.y), (end_point.x, end_point.y), weight=length)
 
-        # 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)
 
+        # Check and connect isolated components
+        components = list(nx.connected_components(g))
+        if len(components) > 1:
+            components.sort(key=len)
+            main_component = components[-1]
+            for comp in components[:-1]:
+                self._connect_component_to_main(g, comp, main_component)
+
         return g
+
+    def _connect_component_to_main(self, graph, component, main_component):
+        main_component_nodes = [graph.nodes[node] for node in main_component if 'geometry' in graph.nodes[node]]
+
+        # Create cKDTree for efficient nearest neighbor search
+        tree = cKDTree([(node['geometry'].x, node['geometry'].y) for node in main_component_nodes])
+
+        # For each node in the component, find the closest street node in the main component
+        for node in component:
+            if 'geometry' in graph.nodes[node]:  # Check for geometry attribute
+                node_geometry = graph.nodes[node]['geometry']
+                distance, idx = tree.query((node_geometry.x, node_geometry.y))
+                closest_node_geometry = main_component_nodes[idx]['geometry']
+
+                # Add edge to the graph
+                graph.add_edge((node_geometry.x, node_geometry.y),
+                               (closest_node_geometry.x, closest_node_geometry.y), weight=distance)
+
+
+def plot_network_graph(network_graph, central_plant_id=1):
+    plt.figure(figsize=(12, 12))
+    pos = {node: (node[0], node[1]) for node in network_graph.nodes()}
+
+    # Node colors based on type
+    node_colors = ['red' if data.get('building_id') == str(central_plant_id) else 'green' if data.get(
+        'type') == 'centroid' else 'blue' for node, data in network_graph.nodes(data=True)]
+
+    # Node sizes, larger for central plant
+    node_sizes = [100 if data.get('building_id') == str(central_plant_id) else 50 for node, data in
+                  network_graph.nodes(data=True)]
+
+    nx.draw_networkx_nodes(network_graph, pos, node_color=node_colors, node_size=node_sizes)
+    nx.draw_networkx_edges(network_graph, pos, edge_color='gray', width=1)
+
+    plt.title('District Heating Network Graph')
+    plt.axis('off')
+    plt.savefig('network_graph_visualization.png', format='png', dpi=300)  # Save as PNG with high dpi for clarity
+    plt.show()

ThermalModeling.py

@@ -0,0 +1,76 @@
+import networkx as nx
+import numpy as np
+
+class ThermalModeling:
+    def __init__(self, graph, T_initial=70, Tg=3, cp=4200, rho=980, U=500, dx=20, delta_t=60):
+        """
+        Initialize the ThermalModeling class with a networkx graph.
+
+        :param graph: A directed networkx graph where edges have 'length', 'diameter', and 'mass_flow_rate_actual'.
+        :param T_initial: Initial temperature at each node in Celsius.
+        :param Tg: Ground temperature in Celsius.
+        :param cp: Isobaric specific heat capacity of water at 60 C in J/(kg*K).
+        :param rho: Density of water in kg/m3.
+        :param U: Heat transfer coefficient for all pipes.
+        :param dx: Number of segments for calculating temperature drops in a pipe.
+        :param delta_t: Time step for the calculation.
+        """
+        self.graph = graph
+        self.T_initial = T_initial
+        self.Tg = Tg
+        self.cp = cp
+        self.rho = rho
+        self.U = U
+        self.dx = dx
+        self.delta_t = delta_t
+
+        # Initialize node temperatures and flow rates
+        for node in self.graph.nodes():
+            self.graph.nodes[node]['temperature'] = T_initial
+            self.graph.nodes[node]['total_mass_flow_in'] = 0
+            self.graph.nodes[node]['weighted_temp_sum'] = 0
+
+    def adjust_flow_directions(self):
+        """
+        Adjust the directions of flow based on the mass flow rate. If the mass flow rate is negative,
+        the direction of the flow is reversed.
+        """
+        to_reverse = [(u, v) for u, v, d in self.graph.edges(data=True) if d['mass_flow_rate_actual'] < 0]
+        for u, v in to_reverse:
+            attrs = self.graph.edges[u, v]
+            self.graph.remove_edge(u, v)
+            self.graph.add_edge(v, u, **attrs)
+            # Update the mass flow rate to be positive after reversing
+            self.graph.edges[v, u]['mass_flow_rate_actual'] = abs(attrs['mass_flow_rate_actual'])
+
+    def calculate_temperatures(self):
+        """
+        Calculate and update temperatures for all nodes based on the network graph, considering weighted averages
+        for nodes with multiple incoming temperatures.
+        """
+        self.adjust_flow_directions()  # Adjust flow directions based on mass flow rates
+
+        # Calculate weighted temperatures for incoming flows
+        for u, v, d in self.graph.edges(data=True):
+            length = d['weight']
+            diameter = d['Diameter']
+            A = np.pi * diameter**2 / 4
+            delta_x = length / self.dx
+            mass_flow_rate = abs(d['mass_flow_rate_actual'])
+
+            C1 = 2 * self.delta_t * self.U / (A * self.rho * self.cp)
+            C2 = 2 * mass_flow_rate * self.delta_t / (self.rho * A * delta_x)
+            C = 1 / (1 + C1 + C2)
+
+            T_in = self.graph.nodes[u]['temperature']
+            T_out = C * (T_in + C1 * self.Tg)  # Simplified model for demonstration
+
+            # Update weighted temperature sum and total mass flow for the target node
+            self.graph.nodes[v]['total_mass_flow_in'] += mass_flow_rate
+            self.graph.nodes[v]['weighted_temp_sum'] += T_out * mass_flow_rate
+
+        # Calculate final temperatures based on weighted averages
+        for node in self.graph.nodes():
+            if self.graph.nodes[node]['total_mass_flow_in'] > 0:  # To avoid division by zero
+                weighted_average_temp = self.graph.nodes[node]['weighted_temp_sum'] / self.graph.nodes[node]['total_mass_flow_in']
+                self.graph.nodes[node]['temperature'] = weighted_average_temp

enrich.py

@@ -0,0 +1,16 @@
+def enrich(graph, city):
+    """
+    Enrich the graph nodes with hourly building demand data.
+
+    :param graph: The networkx graph of the district heating network.
+    :param buildings: A list of building objects, each with a 'name' and 'heating_demand' attribute.
+    """
+    for node in graph.nodes:
+        node_data = graph.nodes[node]
+        # Check if the node has a 'building_id' attribute before comparing
+        if 'building_id' in node_data:
+            for building in city.buildings:
+                if node_data['building_id'] == building.name:
+                    # Assuming `building.heating_demand` is properly structured for direct assignment
+                    graph.nodes[node]["Demand"] = building.heating_demand[cte.HOUR]
+                    graph.nodes[node]["Demand"] = building.heating_peak_load[cte.YEAR]