feature: add district heating network sizing workflow

Reviewed-on: https://nextgenerations-cities.encs.concordia.ca/gitea/s_ranjbar/energy_system_modelling_workflow/pulls/18

district_heating_network.py

@@ -0,0 +1,111 @@
+from scripts.district_heating_network.directory_manager import DirectoryManager
+import subprocess
+from scripts.ep_run_enrich import energy_plus_workflow
+from hub.imports.geometry_factory import GeometryFactory
+from hub.helpers.dictionaries import Dictionaries
+from hub.imports.construction_factory import ConstructionFactory
+from hub.imports.usage_factory import UsageFactory
+from hub.imports.weather_factory import WeatherFactory
+from hub.imports.results_factory import ResultFactory
+from scripts.energy_system_retrofit_report import EnergySystemRetrofitReport
+from scripts.geojson_creator import process_geojson
+from scripts import random_assignation
+from hub.imports.energy_systems_factory import EnergySystemsFactory
+from scripts.energy_system_sizing import SystemSizing
+from scripts.solar_angles import CitySolarAngles
+from scripts.pv_sizing_and_simulation import PVSizingSimulation
+from scripts.energy_system_retrofit_results import consumption_data, cost_data
+from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
+from scripts.costs.cost import Cost
+from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS
+import hub.helpers.constants as cte
+from hub.exports.exports_factory import ExportsFactory
+from scripts.pv_feasibility import pv_feasibility
+import matplotlib.pyplot as plt
+import numpy as np
+from scripts.district_heating_network.district_heating_network_creator import DistrictHeatingNetworkCreator
+from scripts.district_heating_network.district_heating_factory import DistrictHeatingFactory
+import json
+
+#%% --------------------------------------------------------------------------------------------------------------------
+# Manage File Path
+base_path = "./"
+dir_manager = DirectoryManager(base_path)
+
+# Input files directory
+input_files_path = dir_manager.create_directory('input_files')
+geojson_file_path = input_files_path / 'output_buildings.geojson'
+pipe_data_file = input_files_path / 'pipe_data.json'
+
+# Output files directory
+output_path = dir_manager.create_directory('out_files')
+
+# Subdirectories for output files
+energy_plus_output_path = dir_manager.create_directory('out_files/energy_plus_outputs')
+simulation_results_path = dir_manager.create_directory('out_files/simulation_results')
+sra_output_path = dir_manager.create_directory('out_files/sra_outputs')
+cost_analysis_output_path = dir_manager.create_directory('out_files/cost_analysis')
+
+#%% --------------------------------------------------------------------------------------------------------------------
+# Area Under Study
+location = [45.4934614681437, -73.57982834742518]
+
+#%% --------------------------------------------------------------------------------------------------------------------
+# Create geojson of buildings
+process_geojson(x=location[1], y=location[0], diff=0.001)
+
+#%% --------------------------------------------------------------------------------------------------------------------
+# Create ciry and run energyplus workflow
+city = GeometryFactory(file_type='geojson',
+                       path=geojson_file_path,
+                       height_field='height',
+                       year_of_construction_field='year_of_construction',
+                       function_field='function',
+                       function_to_hub=Dictionaries().montreal_function_to_hub_function).city
+ConstructionFactory('nrcan', city).enrich()
+UsageFactory('nrcan', city).enrich()
+WeatherFactory('epw', city).enrich()
+
+# SRA
+ExportsFactory('sra', city, output_path).export()
+sra_path = (output_path / f'{city.name}_sra.xml').resolve()
+subprocess.run(['sra', str(sra_path)])
+ResultFactory('sra', city, output_path).enrich()
+
+# EP Workflow
+energy_plus_workflow(city, energy_plus_output_path)
+
+#%% --------------------------------------------------------------------------------------------------------------------
+# District Heating Network Creator
+central_plant_locations = [(-73.57812571080625, 45.49499447346277)]  # Add at least one location
+
+roads_file = "./input_files/roads.json"
+
+dhn_creator = DistrictHeatingNetworkCreator(geojson_file_path, roads_file, central_plant_locations)
+
+network_graph = dhn_creator.run()
+
+#%% --------------------------------------------------------------------------------------------------------------------
+# Pipe and pump sizing
+
+with open(pipe_data_file, 'r') as f:
+  pipe_data = json.load(f)
+
+factory = DistrictHeatingFactory(
+  city=city,
+  graph=network_graph,
+  supply_temperature=80 + 273,  # in Kelvin
+  return_temperature=60 + 273,  # in Kelvin
+  simultaneity_factor=0.9
+)
+
+factory.enrich()
+factory.sizing()
+factory.calculate_diameters_and_costs(pipe_data)
+pipe_groups, total_cost = factory.analyze_costs()
+
+# Save the pipe groups with total costs to a CSV file
+factory.save_pipe_groups_to_csv('pipe_groups.csv')
+
+#%% --------------------------------------------------------------------------------------------------------------------
+

input_files/pipe_data.json

@@ -0,0 +1,191 @@
+[
+  {
+    "DN": 16,
+    "inner_diameter": 16.1,
+    "outer_diameter": 21.3,
+    "thickness": 2.6,
+    "cost_per_meter": 320
+  },
+  {
+    "DN": 20,
+    "inner_diameter": 21.7,
+    "outer_diameter": 26.9,
+    "thickness": 2.6,
+    "cost_per_meter": 320
+  },
+  {
+    "DN": 25,
+    "inner_diameter": 27.3,
+    "outer_diameter": 33.7,
+    "thickness": 3.2,
+    "cost_per_meter": 320
+  },
+  {
+    "DN": 32,
+    "inner_diameter": 37.2,
+    "outer_diameter": 42.4,
+    "thickness": 2.6,
+    "cost_per_meter": 350
+  },
+  {
+    "DN": 40,
+    "inner_diameter": 43.1,
+    "outer_diameter": 48.3,
+    "thickness": 2.6,
+    "cost_per_meter": 375
+  },
+  {
+    "DN": 50,
+    "inner_diameter": 54.5,
+    "outer_diameter": 60.3,
+    "thickness": 2.9,
+    "cost_per_meter": 400
+  },
+  {
+    "DN": 65,
+    "inner_diameter": 70.3,
+    "outer_diameter": 76.1,
+    "thickness": 2.9,
+    "cost_per_meter": 450
+  },
+  {
+    "DN": 80,
+    "inner_diameter": 82.5,
+    "outer_diameter": 88.9,
+    "thickness": 3.2,
+    "cost_per_meter": 480
+  },
+  {
+    "DN": 90,
+    "inner_diameter": 100.8,
+    "outer_diameter": 108,
+    "thickness": 3.6,
+    "cost_per_meter": 480
+  },
+  {
+    "DN": 100,
+    "inner_diameter": 107.1,
+    "outer_diameter": 114.3,
+    "thickness": 3.6,
+    "cost_per_meter": 550
+  },
+  {
+    "DN": 110,
+    "inner_diameter": 125.8,
+    "outer_diameter": 133,
+    "thickness": 3.6,
+    "cost_per_meter": 550
+  },
+  {
+    "DN": 125,
+    "inner_diameter": 132.5,
+    "outer_diameter": 139.7,
+    "thickness": 3.6,
+    "cost_per_meter": 630
+  },
+  {
+    "DN": 140,
+    "inner_diameter": 151,
+    "outer_diameter": 159,
+    "thickness": 4,
+    "cost_per_meter": 700
+  },
+  {
+    "DN": 150,
+    "inner_diameter": 160.3,
+    "outer_diameter": 168.3,
+    "thickness": 4,
+    "cost_per_meter": 700
+  },
+  {
+    "DN": 180,
+    "inner_diameter": 184.7,
+    "outer_diameter": 193.7,
+    "thickness": 4.5,
+    "cost_per_meter": 700
+  },
+  {
+    "DN": 200,
+    "inner_diameter": 210.1,
+    "outer_diameter": 219.1,
+    "thickness": 4.5,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 250,
+    "inner_diameter": 263,
+    "outer_diameter": 273,
+    "thickness": 5,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 300,
+    "inner_diameter": 312.7,
+    "outer_diameter": 323.9,
+    "thickness": 5.6,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 350,
+    "inner_diameter": 344.4,
+    "outer_diameter": 355.6,
+    "thickness": 5.6,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 400,
+    "inner_diameter": 393.8,
+    "outer_diameter": 406.4,
+    "thickness": 6.3,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 450,
+    "inner_diameter": 444.4,
+    "outer_diameter": 457,
+    "thickness": 6.3,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 500,
+    "inner_diameter": 495.4,
+    "outer_diameter": 508,
+    "thickness": 6.3,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 600,
+    "inner_diameter": 595.8,
+    "outer_diameter": 610,
+    "thickness": 7.1,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 700,
+    "inner_diameter": 696.8,
+    "outer_diameter": 711,
+    "thickness": 7.1,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 800,
+    "inner_diameter": 797,
+    "outer_diameter": 813,
+    "thickness": 8,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 900,
+    "inner_diameter": 894,
+    "outer_diameter": 914,
+    "thickness": 10,
+    "cost_per_meter": 860
+  },
+  {
+    "DN": 1000,
+    "inner_diameter": 996,
+    "outer_diameter": 1016,
+    "thickness": 10,
+    "cost_per_meter": 860
+  }
+]

main.py

@@ -1,4 +1,5 @@
 from pathlib import Path
+from scripts.district_heating_network.directory_manager import DirectoryManager
 import subprocess
 from scripts.ep_run_enrich import energy_plus_workflow
 from hub.imports.geometry_factory import GeometryFactory
@@ -22,23 +23,31 @@ import hub.helpers.constants as cte
 from hub.exports.exports_factory import ExportsFactory
 from scripts.pv_feasibility import pv_feasibility
 import matplotlib.pyplot as plt
-import numpy as np
-# Specify the GeoJSON file path
-data = {}
-input_files_path = (Path(__file__).parent / 'input_files')
-input_files_path.mkdir(parents=True, exist_ok=True)
-geojson_file = process_geojson(x=-73.58001358793511, y=45.496445294438715, diff=0.0001)
+from scripts.district_heating_network.district_heating_network_creator import DistrictHeatingNetworkCreator
+from scripts.district_heating_network.road_processor import road_processor
+from scripts.district_heating_network.district_heating_factory import DistrictHeatingFactory
+
+base_path = Path(__file__).parent
+dir_manager = DirectoryManager(base_path)
+
+# Input files directory
+input_files_path = dir_manager.create_directory('input_files')
 geojson_file_path = input_files_path / 'output_buildings.geojson'
-output_path = (Path(__file__).parent / 'out_files').resolve()
-output_path.mkdir(parents=True, exist_ok=True)
-energy_plus_output_path = output_path / 'energy_plus_outputs'
-energy_plus_output_path.mkdir(parents=True, exist_ok=True)
-simulation_results_path = (Path(__file__).parent / 'out_files' / 'simulation_results').resolve()
-simulation_results_path.mkdir(parents=True, exist_ok=True)
-sra_output_path = output_path / 'sra_outputs'
-sra_output_path.mkdir(parents=True, exist_ok=True)
-cost_analysis_output_path = output_path / 'cost_analysis'
-cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
+
+# Output files directory
+output_path = dir_manager.create_directory('out_files')
+
+# Subdirectories for output files
+energy_plus_output_path = dir_manager.create_directory('out_files/energy_plus_outputs')
+simulation_results_path = dir_manager.create_directory('out_files/simulation_results')
+sra_output_path = dir_manager.create_directory('out_files/sra_outputs')
+cost_analysis_output_path = dir_manager.create_directory('out_files/cost_analysis')
+
+# Select city area
+location = [45.53067276979674, -73.70234652694087]
+process_geojson(x=location[1], y=location[0], diff=0.001)
+
+# Create city object
 city = GeometryFactory(file_type='geojson',
                        path=geojson_file_path,
                        height_field='height',
@@ -83,4 +92,4 @@ UsageFactory('nrcan', city).enrich()
 # # Save the plot
 # plt.savefig('plot_nrcan.png', bbox_inches='tight')
 # plt.close()
-print('test')
+print('test')

pv_assessment.py

@@ -14,7 +14,7 @@ import hub.helpers.constants as cte
 from hub.exports.exports_factory import ExportsFactory
 from scripts.pv_sizing_and_simulation import PVSizingSimulation
 # Specify the GeoJSON file path
-geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0005)
+geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
 file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson')
 # Specify the output path for the PDF file
 output_path = (Path(__file__).parent / 'out_files').resolve()
@@ -34,7 +34,7 @@ ExportsFactory('sra', city, output_path).export()
 sra_path = (output_path / f'{city.name}_sra.xml').resolve()
 subprocess.run(['sra', str(sra_path)])
 ResultFactory('sra', city, output_path).enrich()
-energy_plus_workflow(city)
+energy_plus_workflow(city, output_path=output_path)
 solar_angles = CitySolarAngles(city.name,
                                city.latitude,
                                city.longitude,

scripts/district_heating_network/directory_manager.py

@@ -0,0 +1,16 @@
+from pathlib import Path
+
+
+class DirectoryManager:
+    def __init__(self, base_path):
+        self.base_path = Path(base_path)
+        self.directories = {}
+
+    def create_directory(self, relative_path):
+        full_path = self.base_path / relative_path
+        full_path.mkdir(parents=True, exist_ok=True)
+        self.directories[relative_path] = full_path
+        return full_path
+
+    def get_directory(self, relative_path):
+        return self.directories.get(relative_path, None)

scripts/district_heating_network/district_heating_factory.py

@@ -0,0 +1,248 @@
+import CoolProp.CoolProp as CP
+import math
+import logging
+import numpy as np
+import csv
+
+
+class DistrictHeatingFactory:
+  """
+  DistrictHeatingFactory class
+
+  This class is responsible for managing the district heating network, including
+  enriching the network graph with building data, calculating flow rates,
+  sizing pipes, and analyzing costs.
+  """
+
+  def __init__(self, city, graph, supply_temperature, return_temperature, simultaneity_factor):
+    """
+    Initialize the DistrictHeatingFactory object.
+
+    :param city: The city object containing buildings and their heating demands.
+    :param graph: The network graph representing the district heating network.
+    :param supply_temperature: The supply temperature of the heating fluid in the network (°C).
+    :param return_temperature: The return temperature of the heating fluid in the network (°C).
+    :param simultaneity_factor: The simultaneity factor used to adjust flow rates for non-building pipes.
+    """
+    self._city = city
+    self._network_graph = graph
+    self._supply_temperature = supply_temperature
+    self._return_temperature = return_temperature
+    self.simultaneity_factor = simultaneity_factor
+    self.fluid = "Water"  # The fluid used in the heating network
+
+  def enrich(self):
+    """
+    Enrich the network graph nodes with the whole building object from the city buildings.
+
+    This method associates each building node in the network graph with its corresponding
+    building object from the city, allowing access to heating demand data during calculations.
+    """
+    for node_id, node_attrs in self._network_graph.nodes(data=True):
+      if node_attrs.get('type') == 'building':
+        building_name = node_attrs.get('name')
+        building_found = False
+        for building in self._city.buildings:
+          if building.name == building_name:
+            self._network_graph.nodes[node_id]['building_obj'] = building
+            building_found = True
+            break
+        if not building_found:
+          logging.error(msg=f"Building with name '{building_name}' not found in city.")
+
+  def calculate_flow_rates(self, A, Gext):
+    """
+    Solve the linear system to find the flow rates in each branch.
+
+    :param A: The incidence matrix representing the network connections.
+    :param Gext: The external flow rates for each node in the network.
+    :return: The calculated flow rates for each edge, or None if an error occurs.
+    """
+    try:
+      G = np.linalg.lstsq(A, Gext, rcond=None)[0]
+      return G
+    except np.linalg.LinAlgError as e:
+      logging.error(f"Error solving the linear system: {e}")
+      return None
+
+  def switch_nodes(self, A, edge_index, node_index, edge):
+    """
+    Switch the in and out nodes for the given edge in the incidence matrix A.
+
+    :param A: The incidence matrix representing the network connections.
+    :param edge_index: The index of edges in the incidence matrix.
+    :param node_index: The index of nodes in the incidence matrix.
+    :param edge: The edge (u, v) to switch.
+    """
+    u, v = edge
+    i = node_index[u]
+    j = node_index[v]
+    k = edge_index[edge]
+    A[i, k], A[j, k] = -A[i, k], -A[j, k]
+
+  def sizing(self):
+    """
+    Calculate the hourly mass flow rates, assign them to the edges, and determine the pipe diameters.
+
+    This method generates the flow rates for each hour, adjusting the incidence matrix as needed to
+    ensure all flow rates are positive. It also applies the simultaneity factor to non-building pipes.
+    """
+    num_nodes = self._network_graph.number_of_nodes()
+    num_edges = self._network_graph.number_of_edges()
+    A = np.zeros((num_nodes, num_edges))  # Initialize incidence matrix
+    node_index = {node: i for i, node in enumerate(self._network_graph.nodes())}
+    edge_index = {edge: i for i, edge in enumerate(self._network_graph.edges())}
+
+    # Initialize mass flow rate attribute for each edge
+    for u, v, data in self._network_graph.edges(data=True):
+      self._network_graph.edges[u, v]['mass_flow_rate'] = {"hour": [], "peak": None}
+
+    # Get the length of the hourly demand for the first building (assuming all buildings have the same length)
+    building = next(iter(self._city.buildings))
+    num_hours = len(building.heating_demand['hour'])
+
+    # Loop through each hour to generate Gext and solve AG = Gext
+    for hour in range(8760):
+      Gext = np.zeros(num_nodes)
+
+      # Calculate the hourly mass flow rates for each edge and fill Gext
+      for edge in self._network_graph.edges(data=True):
+        u, v, data = edge
+        for node in [u, v]:
+          if self._network_graph.nodes[node].get('type') == 'building':
+            building = self._network_graph.nodes[node].get('building_obj')
+            if building and "hour" in building.heating_demand:
+              hourly_demand = building.heating_demand["hour"][hour]  # Get demand for current hour
+              specific_heat_capacity = CP.PropsSI('C', 'T', (self._supply_temperature + self._return_temperature) / 2,
+                                                  'P', 101325, self.fluid)
+              mass_flow_rate = hourly_demand / 3600 / (
+                    specific_heat_capacity * (self._supply_temperature - self._return_temperature))
+              Gext[node_index[node]] += mass_flow_rate
+
+        # Update incidence matrix A
+        i = node_index[u]
+        j = node_index[v]
+        k = edge_index[(u, v)]
+        A[i, k] = 1
+        A[j, k] = -1
+
+      # Solve for G (flow rates)
+      G = self.calculate_flow_rates(A, Gext)
+      if G is None:
+        return
+
+      # Check for negative flow rates and adjust A accordingly
+      iterations = 0
+      max_iterations = num_edges * 2
+      while any(flow_rate < 0 for flow_rate in G) and iterations < max_iterations:
+        for idx, flow_rate in enumerate(G):
+          if flow_rate < 0:
+            G[idx] = -G[idx]  # Invert the sign directly
+        iterations += 1
+
+      # Store the final flow rates in the edges for this hour
+      for idx, (edge, flow_rate) in enumerate(zip(self._network_graph.edges(), G)):
+        u, v = edge
+        if not (self._network_graph.nodes[u].get('type') == 'building' or self._network_graph.nodes[v].get(
+            'type') == 'building'):
+          flow_rate *= self.simultaneity_factor  # Apply simultaneity factor for non-building pipes
+        data = self._network_graph.edges[u, v]
+        data['mass_flow_rate']["hour"].append(flow_rate)  # Append the calculated flow rate
+
+    # Calculate the peak flow rate for each edge
+    for u, v, data in self._network_graph.edges(data=True):
+      data['mass_flow_rate']['peak'] = max(data['mass_flow_rate']['hour'])
+
+  def calculate_diameters_and_costs(self, pipe_data):
+    """
+    Calculate the diameter and costs of the pipes based on the maximum flow rate in each edge.
+
+    :param pipe_data: A list of dictionaries containing pipe specifications, including inner diameters
+                      and costs per meter for different nominal diameters (DN).
+    """
+    for u, v, data in self._network_graph.edges(data=True):
+      flow_rate = data.get('mass_flow_rate', {}).get('peak')
+      if flow_rate is not None:
+        try:
+          # Calculate the density of the fluid
+          density = CP.PropsSI('D', 'T', (self._supply_temperature + self._return_temperature) / 2, 'P', 101325,
+                               self.fluid)
+          velocity = 0.9  # Desired fluid velocity in m/s
+          # Calculate the diameter of the pipe required for the given flow rate
+          diameter = math.sqrt((4 * abs(flow_rate)) / (density * velocity * math.pi)) * 1000  # Convert to mm
+          self._network_graph.edges[u, v]['diameter'] = diameter
+
+          # Match to the closest nominal diameter from the pipe data
+          closest_pipe = self.match_nominal_diameter(diameter, pipe_data)
+          self._network_graph.edges[u, v]['nominal_diameter'] = closest_pipe['DN']
+          self._network_graph.edges[u, v]['cost_per_meter'] = closest_pipe['cost_per_meter']
+        except Exception as e:
+          logging.error(f"Error calculating diameter or matching nominal diameter for edge ({u}, {v}): {e}")
+
+  def match_nominal_diameter(self, diameter, pipe_data):
+    """
+    Match the calculated diameter to the closest nominal diameter.
+
+    :param diameter: The calculated diameter of the pipe (in mm).
+    :param pipe_data: A list of dictionaries containing pipe specifications, including inner diameters
+                      and costs per meter for different nominal diameters (DN).
+    :return: The dictionary representing the pipe with the closest nominal diameter.
+    """
+    closest_pipe = min(pipe_data, key=lambda x: abs(x['inner_diameter'] - diameter))
+    return closest_pipe
+
+  def analyze_costs(self):
+    """
+    Analyze the costs based on the nominal diameters of the pipes.
+
+    This method calculates the total cost of piping for each nominal diameter group
+    and returns a summary of the grouped pipes and the total cost.
+
+    :return: A tuple containing the grouped pipe data and the total cost of piping.
+    """
+    pipe_groups = {}
+    total_cost = 0  # Initialize total cost
+
+    for u, v, data in self._network_graph.edges(data=True):
+      dn = data.get('nominal_diameter')
+      if dn is not None:
+        pipe_length = self._network_graph.edges[u, v].get('length', 1) * 2  # Multiply by 2 for supply and return
+        cost_per_meter = data.get('cost_per_meter', 0)
+
+        if dn not in pipe_groups:
+          pipe_groups[dn] = {
+            'DN': dn,
+            'total_length': 0,
+            'cost_per_meter': cost_per_meter
+          }
+        pipe_groups[dn]['total_length'] += pipe_length
+        group_cost = pipe_length * cost_per_meter
+        total_cost += group_cost  # Add to total cost
+
+    # Calculate total cost for each group
+    for group in pipe_groups.values():
+      group['total_cost'] = group['total_length'] * group['cost_per_meter']
+
+    return pipe_groups, total_cost  # Return both the grouped data and total cost
+
+  def save_pipe_groups_to_csv(self, filename):
+    """
+    Save the pipe groups and their total lengths to a CSV file.
+
+    :param filename: The name of the CSV file to save the data to.
+    """
+    pipe_groups, _ = self.analyze_costs()
+
+    with open(filename, mode='w', newline='') as file:
+      writer = csv.writer(file)
+      # Write the header
+      writer.writerow(["Nominal Diameter (DN)", "Total Length (m)", "Cost per Meter", "Total Cost"])
+
+      # Write the data for each pipe group
+      for group in pipe_groups.values():
+        writer.writerow([
+          group['DN'],
+          group['total_length'],
+          group['cost_per_meter'],
+          group['total_cost']
+        ])

scripts/district_heating_network/district_heating_network_creator.py

@@ -0,0 +1,372 @@
+import json
+import math
+import logging
+import matplotlib.pyplot as plt
+import networkx as nx
+from shapely.geometry import Polygon, Point, LineString
+from typing import List, Tuple
+from rtree import index
+
+# Configure logging
+logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s")
+logging.getLogger("numexpr").setLevel(logging.ERROR)
+
+def haversine(lon1, lat1, lon2, lat2):
+  """
+  Calculate the great-circle distance between two points
+  on the Earth specified by their longitude and latitude.
+  """
+  R = 6371000  # Radius of the Earth in meters
+  phi1 = math.radians(lat1)
+  phi2 = math.radians(lat2)
+  delta_phi = math.radians(lat2 - lat1)
+  delta_lambda = math.radians(lon2 - lon1)
+
+  a = math.sin(delta_phi / 2.0) ** 2 + \
+      math.cos(phi1) * math.cos(phi2) * \
+      math.sin(delta_lambda / 2.0) ** 2
+
+  c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
+  return R * c  # Output distance in meters
+
+class DistrictHeatingNetworkCreator:
+  def __init__(self, buildings_file: str, roads_file: str, central_plant_locations: List[Tuple[float, float]]):
+    """
+    Initialize the class with paths to the buildings and roads data files, and central plant locations.
+
+    :param buildings_file: Path to the GeoJSON file containing building data.
+    :param roads_file: Path to the GeoJSON file containing roads data.
+    :param central_plant_locations: List of tuples containing the coordinates of central plant locations.
+    """
+    if len(central_plant_locations) < 1:
+      raise ValueError("The list of central plant locations must have at least one member.")
+
+    self.buildings_file = buildings_file
+    self.roads_file = roads_file
+    self.central_plant_locations = central_plant_locations
+
+  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.
+    """
+    try:
+      self._load_and_process_data()
+      self._find_nearest_roads()
+      self._find_nearest_points()
+      self._break_down_roads()
+      self._create_graph()
+      self._create_mst()
+      self._iteratively_remove_edges()
+      self._add_centroids_to_mst()
+      self._convert_edge_weights_to_meters()
+      self._create_final_network_graph()
+      return self.network_graph
+    except Exception as e:
+      logging.error(f"Error during network creation: {e}")
+      raise
+
+  def _load_and_process_data(self):
+    """
+    Load and process the building and road data.
+    """
+    try:
+      # Load building data
+      with open(self.buildings_file, 'r') as file:
+        city = json.load(file)
+
+      self.centroids = []
+      self.building_names = []
+      self.building_positions = []
+      buildings = city['features']
+      for building in buildings:
+        coordinates = building['geometry']['coordinates'][0]
+        building_polygon = Polygon(coordinates)
+        centroid = building_polygon.centroid
+        self.centroids.append(centroid)
+        self.building_names.append(str(building['id']))
+        self.building_positions.append((centroid.x, centroid.y))
+
+      # Add central plant locations as centroids
+      for i, loc in enumerate(self.central_plant_locations, start=1):
+        centroid = Point(loc)
+        self.centroids.append(centroid)
+        self.building_names.append(f'central_plant_{i}')
+        self.building_positions.append((centroid.x, centroid.y))
+
+      # Load road data
+      with open(self.roads_file, 'r') as file:
+        roads = json.load(file)
+
+      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]
+    except Exception as e:
+      logging.error(f"Error loading and processing data: {e}")
+      raise
+
+  def _find_nearest_roads(self):
+    """
+    Find the nearest road for each building centroid.
+    """
+    try:
+      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)
+    except Exception as e:
+      logging.error(f"Error finding nearest roads: {e}")
+      raise
+
+  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))
+
+    try:
+      self.nearest_points = []
+      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
+
+        if closest_road:
+          nearest_point = find_nearest_point_on_line(centroid, closest_road)
+          self.nearest_points.append(nearest_point)
+    except Exception as e:
+      logging.error(f"Error finding nearest points: {e}")
+      raise
+
+  def _break_down_roads(self):
+    """
+    Break down roads into segments connecting 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
+
+    try:
+      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)
+    except Exception as e:
+      logging.error(f"Error breaking down roads: {e}")
+      raise
+
+  def _create_graph(self):
+    """
+    Create a NetworkX graph from the cleaned lines.
+    """
+    try:
+      self.G = nx.Graph()
+      for line in self.cleaned_lines:
+        coords = list(line.coords)
+        for i in range(len(coords) - 1):
+          u = coords[i]
+          v = coords[i + 1]
+          self.G.add_edge(u, v, weight=Point(coords[i]).distance(Point(coords[i + 1])))
+    except Exception as e:
+      logging.error(f"Error creating graph: {e}")
+      raise
+
+  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
+
+    try:
+      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)
+    except Exception as e:
+      logging.error(f"Error creating MST: {e}")
+      raise
+
+  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]
+
+    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
+
+    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
+
+    try:
+      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)))
+    except Exception as e:
+      logging.error(f"Error iteratively removing edges: {e}")
+      raise
+
+  def _add_centroids_to_mst(self):
+    """
+    Add centroids to the final MST graph and connect them to their associated node on the graph.
+    """
+    try:
+      for i, centroid in enumerate(self.centroids):
+        building_name = self.building_names[i]
+        pos = (centroid.x, centroid.y)
+        node_type = 'building' if 'central_plant' not in building_name else 'generation'
+        self.final_mst.add_node(pos, type=node_type, name=building_name, pos=pos)
+
+        nearest_point = None
+        min_distance = float('inf')
+        for node in self.final_mst.nodes():
+          if self.final_mst.nodes[node].get('type') != 'building' and self.final_mst.nodes[node].get('type') != 'generation':
+            distance = centroid.distance(Point(node))
+            if distance < min_distance:
+              min_distance = distance
+              nearest_point = node
+
+        if nearest_point:
+          self.final_mst.add_edge(pos, nearest_point, weight=min_distance)
+    except Exception as e:
+      logging.error(f"Error adding centroids to MST: {e}")
+      raise
+
+  def _convert_edge_weights_to_meters(self):
+    """
+    Convert all edge weights in the final MST graph to meters using the Haversine formula.
+    """
+    try:
+      for u, v, data in self.final_mst.edges(data=True):
+        distance = haversine(u[0], u[1], v[0], v[1])
+        self.final_mst[u][v]['weight'] = distance
+    except Exception as e:
+      logging.error(f"Error converting edge weights to meters: {e}")
+      raise
+
+  def _create_final_network_graph(self):
+    """
+    Create the final network graph with the required attributes from the final MST.
+    """
+    self.network_graph = nx.Graph()
+    node_id = 1
+    node_mapping = {}
+    for node in self.final_mst.nodes:
+      pos = node
+      if 'type' in self.final_mst.nodes[node]:
+        if self.final_mst.nodes[node]['type'] == 'building':
+          name = self.final_mst.nodes[node]['name']
+          node_type = 'building'
+        elif self.final_mst.nodes[node]['type'] == 'generation':
+          name = self.final_mst.nodes[node]['name']
+          node_type = 'generation'
+      else:
+        name = f'junction_{node_id}'
+        node_type = 'junction'
+      self.network_graph.add_node(node_id, name=name, type=node_type, pos=pos)
+      node_mapping[node] = node_id
+      node_id += 1
+    for u, v, data in self.final_mst.edges(data=True):
+      u_new = node_mapping[u]
+      v_new = node_mapping[v]
+      length = data['weight']
+      self.network_graph.add_edge(u_new, v_new, length=length)
+
+  def plot_network_graph(self):
+    """
+    Plot the network graph using matplotlib and networkx.
+    """
+    plt.figure(figsize=(15, 10))
+    pos = {node: data['pos'] for node, data in self.network_graph.nodes(data=True)}
+    nx.draw_networkx_nodes(self.network_graph, pos, node_color='blue', node_size=50)
+    nx.draw_networkx_edges(self.network_graph, pos, edge_color='gray')
+    plt.title('District Heating Network Graph')
+    plt.axis('off')
+    plt.show()

scripts/district_heating_network/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

simulation_result_test.py

@@ -1,54 +0,0 @@
-from pathlib import Path
-import subprocess
-from scripts.ep_run_enrich import energy_plus_workflow
-from hub.imports.geometry_factory import GeometryFactory
-from hub.helpers.dictionaries import Dictionaries
-from hub.imports.construction_factory import ConstructionFactory
-from hub.imports.usage_factory import UsageFactory
-from hub.imports.weather_factory import WeatherFactory
-from hub.imports.results_factory import ResultFactory
-from scripts.energy_system_retrofit_report import EnergySystemRetrofitReport
-from scripts.geojson_creator import process_geojson
-from scripts import random_assignation
-from hub.imports.energy_systems_factory import EnergySystemsFactory
-from scripts.energy_system_sizing import SystemSizing
-from scripts.solar_angles import CitySolarAngles
-from scripts.pv_sizing_and_simulation import PVSizingSimulation
-from scripts.energy_system_retrofit_results import consumption_data, cost_data
-from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
-from scripts.costs.cost import Cost
-from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, CURRENT_STATUS
-import hub.helpers.constants as cte
-from hub.exports.exports_factory import ExportsFactory
-from scripts.pv_feasibility import pv_feasibility
-
-# Specify the GeoJSON file path
-input_files_path = (Path(__file__).parent / 'input_files')
-input_files_path.mkdir(parents=True, exist_ok=True)
-geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
-geojson_file_path = input_files_path / 'output_buildings.geojson'
-output_path = (Path(__file__).parent / 'out_files').resolve()
-output_path.mkdir(parents=True, exist_ok=True)
-energy_plus_output_path = output_path / 'energy_plus_outputs'
-energy_plus_output_path.mkdir(parents=True, exist_ok=True)
-simulation_results_path = (Path(__file__).parent / 'out_files' / 'simulation_results').resolve()
-simulation_results_path.mkdir(parents=True, exist_ok=True)
-sra_output_path = output_path / 'sra_outputs'
-sra_output_path.mkdir(parents=True, exist_ok=True)
-cost_analysis_output_path = output_path / 'cost_analysis'
-cost_analysis_output_path.mkdir(parents=True, exist_ok=True)
-city = GeometryFactory(file_type='geojson',
-                       path=geojson_file_path,
-                       height_field='height',
-                       year_of_construction_field='year_of_construction',
-                       function_field='function',
-                       function_to_hub=Dictionaries().montreal_function_to_hub_function).city
-ConstructionFactory('nrcan', city).enrich()
-UsageFactory('nrcan', city).enrich()
-WeatherFactory('epw', city).enrich()
-energy_plus_workflow(city, energy_plus_output_path)
-random_assignation.call_random(city.buildings, random_assignation.residential_new_systems_percentage)
-EnergySystemsFactory('montreal_future', city).enrich()
-for building in city.buildings:
-  EnergySystemsSimulationFactory('archetype1', building=building, output_path=simulation_results_path).enrich()
-