import subprocess from pathlib import Path import geopandas as gpd import pandas as pd import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import Normalize from shapely.geometry import Polygon from hub.imports.energy_systems_factory import EnergySystemsFactory from hub.imports.geometry_factory import GeometryFactory from hub.helpers.dictionaries import Dictionaries from hub.imports.construction_factory import ConstructionFactory from hub.imports.results_factory import ResultFactory from hub.exports.exports_factory import ExportsFactory from hub.imports.weather_factory import WeatherFactory from pv_assessment.solar_calculator import SolarCalculator from pv_assessment.pv_system_assessment_with_lcoe import PvSystemAssessment from scripts import random_assignation input_file_1 = "./data/updated_geojson_layers/output_layer_6.geojson" input_file_2 = "./data/updated_geojson_layers/output_layer_20.geojson" input_file_3 = "./data/updated_geojson_layers/output_layer_38.geojson" input_file_4 = "./data/updated_geojson_layers/output_layer_53.geojson" demand_file = "data/energy_demand_data.csv" output_path = (Path(__file__).parent.parent / 'hub/out_files').resolve() output_path.mkdir(parents=True, exist_ok=True) sra_output_path = output_path / 'sra_outputs' sra_output_path.mkdir(parents=True, exist_ok=True) pv_assessment_path = output_path / 'pv_outputs' pv_assessment_path.mkdir(parents=True, exist_ok=True) figures_path = (Path(__file__).parent.parent / 'hub/figures').resolve() figures_path.mkdir(parents=True, exist_ok=True) def process_city(input_file, name_suffix): city_obj = GeometryFactory( "geojson", input_file, height_field="height", year_of_construction_field="contr_year", function_field="function_c", adjacency_field="adjacency", lot_area_field='lot_area', build_area_field='build_area', function_to_hub=Dictionaries().montreal_function_to_hub_function ).city ConstructionFactory('nrcan', city_obj).enrich() WeatherFactory('epw', city_obj).enrich() ResultFactory('archetypes', city_obj, demand_file).enrich() ExportsFactory('sra', city_obj, sra_output_path).export() sra_path = (sra_output_path / f'{city_obj.name}_sra.xml').resolve() subprocess.run(['sra', str(sra_path)]) ResultFactory('sra', city_obj, sra_output_path).enrich() tilt_angle = 37 solar_parameters = SolarCalculator(city=city_obj, surface_azimuth_angle=180, tilt_angle=tilt_angle, standard_meridian=-75) solar_angles = solar_parameters.solar_angles solar_parameters.tilted_irradiance_calculator() random_assignation.call_random(city_obj.buildings, random_assignation.residential_systems_percentage) EnergySystemsFactory('montreal_future', city_obj).enrich() for building in city_obj.buildings: PvSystemAssessment(building=building, pv_system=None, battery=None, electricity_demand=None, tilt_angle=tilt_angle, solar_angles=solar_angles, pv_installation_type='rooftop', simulation_model_type='explicit', module_model_name=None, inverter_efficiency=0.95, system_catalogue_handler=None, roof_percentage_coverage=0.75, facade_coverage_percentage=0, csv_output=False, output_path=pv_assessment_path, run_lcoe=True).enrich() return city_obj city1 = process_city(input_file_1, "_city1") city2 = process_city(input_file_2, "_city2") city3 = process_city(input_file_3, "_city3") city4 = process_city(input_file_4, "_city4") print("All cities processed.") def city_to_gdf(city_obj): buildings_data = [] for b in city_obj.buildings: yearly_pv = b.pv_generation['year'] if isinstance(yearly_pv, list): yearly_pv = sum(yearly_pv) lcoe_val = b.pv_generation['LCOE_PV'] yearly_self_sufficiency = b.self_sufficiency['year']/1000.0 if 'year' in b.self_sufficiency else np.nan if not b.surfaces: continue ground_surface = b.surfaces[0] coords_3d = ground_surface.solid_polygon.coordinates coords_2d = [(c[0], c[1]) for c in coords_3d] footprint_poly = Polygon(coords_2d) buildings_data.append({ "geometry": footprint_poly, "yearly_pv": yearly_pv, "LCOE": lcoe_val, "self_sufficiency": yearly_self_sufficiency }) gdf = gpd.GeoDataFrame(buildings_data, crs="EPSG:26911") return gdf gdf_city1 = city_to_gdf(city1) gdf_city2 = city_to_gdf(city2) gdf_city3 = city_to_gdf(city3) gdf_city4 = city_to_gdf(city4) all_buildings_gdf = gpd.GeoDataFrame(pd.concat([gdf_city1, gdf_city2, gdf_city3, gdf_city4], ignore_index=True), crs=gdf_city1.crs) print("Dataframe for all 4 cities created, now plotting and saving figures") layers_file = "data/cmm_limites_avec_mtl.geojson" layers_gdf = gpd.read_file(layers_file) if layers_gdf.crs is None: layers_gdf.set_crs(epsg=32188, inplace=True) all_buildings_gdf = all_buildings_gdf.to_crs(layers_gdf.crs) if 'self_sufficiency' in all_buildings_gdf.columns and not all_buildings_gdf['self_sufficiency'].isna().all(): vmin = min(0, all_buildings_gdf['self_sufficiency'].min()) vmax = max(0, all_buildings_gdf['self_sufficiency'].max()) fig, ax = plt.subplots(1, 1, figsize=(14, 10)) cmap = plt.cm.viridis norm = Normalize(vmin=vmin, vmax=vmax) all_buildings_gdf.plot(column='self_sufficiency', cmap=cmap, linewidth=0.8, edgecolor='grey', legend=False, ax=ax) sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) sm._A = [] cbar = fig.colorbar(sm, ax=ax, fraction=0.03, pad=0.04) cbar.set_label('Self-Sufficiency (kWh/year)', fontsize=12) ax.grid(color='lightgrey', linestyle='--', linewidth=0.5, alpha=0.7) ax.set_title('Building Self-Sufficiency Levels (All 4 Cities)', fontsize=16, fontweight='bold', pad=20) ax.set_xlabel('Longitude', fontsize=12) ax.set_ylabel('Latitude', fontsize=12) plt.tight_layout() plt.savefig(figures_path / "self_sufficiency.png", dpi=300) plt.close(fig) else: print("No self_sufficiency data available.") if 'yearly_pv' in all_buildings_gdf.columns and not all_buildings_gdf['yearly_pv'].isna().all(): fig2, ax2 = plt.subplots(1, 1, figsize=(14, 10)) vmin_pv = all_buildings_gdf['yearly_pv'].min() vmax_pv = all_buildings_gdf['yearly_pv'].max() all_buildings_gdf.plot(column='yearly_pv', cmap="plasma", edgecolor="white", linewidth=0.5, legend=True, legend_kwds={'label': "Yearly PV Generation (kWh)", 'orientation': "vertical"}, vmin=vmin_pv, vmax=vmax_pv, ax=ax2) ax2.set_title("Yearly PV Generation (All 4 Cities)", fontsize=16, fontweight='bold') ax2.set_aspect('equal', 'box') plt.tight_layout() plt.savefig(figures_path / "yearly_pv_generation.png", dpi=300) plt.close(fig2) else: print("No yearly_pv data available.") lcoe_gdf = all_buildings_gdf.dropna(subset=['LCOE']) if not lcoe_gdf.empty: vmin_lcoe = lcoe_gdf['LCOE'].min() vmax_lcoe = lcoe_gdf['LCOE'].max() cmap_lcoe = plt.cm.viridis norm_lcoe = Normalize(vmin=vmin_lcoe, vmax=vmax_lcoe) fig3, ax3 = plt.subplots(1, 1, figsize=(14, 10)) lcoe_gdf.plot(column='LCOE', cmap=cmap_lcoe, linewidth=0.8, edgecolor='grey', legend=False, ax=ax3) sm_lcoe = plt.cm.ScalarMappable(cmap=cmap_lcoe, norm=norm_lcoe) sm_lcoe._A = [] cbar_lcoe = fig3.colorbar(sm_lcoe, ax=ax3, fraction=0.03, pad=0.04) cbar_lcoe.set_label('LCOE (currency/kWh)', fontsize=12) ax3.grid(color='lightgrey', linestyle='--', linewidth=0.5, alpha=0.7) ax3.set_title('LCOE Values for Buildings (All 4 Cities)', fontsize=16, fontweight='bold', pad=20) ax3.set_xlabel('Longitude', fontsize=12) ax3.set_ylabel('Latitude', fontsize=12) plt.tight_layout() plt.savefig(figures_path / "lcoe_values.png", dpi=300) plt.close(fig3) else: print("No buildings have LCOE values computed.") layers_of_interest = [f"layer_{i}" for i in range(6, 21)] for layer_name in layers_of_interest: layer_poly = layers_gdf[layers_gdf["region"] == layer_name] if layer_poly.empty: continue layer_union = layer_poly.unary_union if layer_union.is_empty: continue buildings_in_layer = all_buildings_gdf[all_buildings_gdf.geometry.within(layer_union)] if buildings_in_layer.empty: continue fig, ax = plt.subplots(figsize=(20, 20), dpi=600) layer_poly.plot(ax=ax, facecolor="none", edgecolor="black", linewidth=1) vmin = buildings_in_layer["yearly_pv"].min() vmax = buildings_in_layer["yearly_pv"].max() buildings_in_layer.plot(column="yearly_pv", ax=ax, cmap="viridis", edgecolor="white", linewidth=0.5, legend=True, legend_kwds={'label': "Yearly PV Generation (kWh)", 'orientation': "vertical"}, vmin=vmin, vmax=vmax) ax.set_title(f"Yearly PV Generation - {layer_name}", fontsize=14) ax.set_aspect('equal', 'box') plt.tight_layout() plt.savefig(figures_path / f"{layer_name}_pv_generation_multi_city.pdf", dpi=300) plt.close(fig) fig, ax = plt.subplots(figsize=(20, 20), dpi=1200) layers_gdf.plot(ax=ax, facecolor="none", edgecolor="grey", linewidth=0.5) vmin = all_buildings_gdf["yearly_pv"].min() vmax = all_buildings_gdf["yearly_pv"].max() all_buildings_gdf.plot(column="yearly_pv", ax=ax, cmap="plasma", edgecolor="none", legend=True, legend_kwds={'label': "Yearly PV Generation (kWh)", 'orientation': "vertical"}, vmin=vmin, vmax=vmax) ax.set_title("Yearly PV Generation - All Layers (All 4 Cities)", fontsize=14) ax.set_aspect('equal', 'box') plt.tight_layout() plt.savefig(figures_path / "all_layers_pv_generation_all_4_cities.png", dpi=1200) plt.close(fig) layer_sums = [] for idx, row in layers_gdf.iterrows(): region_name = row["region"] poly_geom = row["geometry"] b_in_poly = all_buildings_gdf[all_buildings_gdf.geometry.within(poly_geom)] total_pv = b_in_poly["yearly_pv"].sum() layer_sums.append({ "region": region_name, "total_pv": total_pv, "geometry": poly_geom }) layers_pv_gdf = gpd.GeoDataFrame(layer_sums, crs=layers_gdf.crs) fig, ax = plt.subplots(figsize=(10, 10), dpi=300) vmin = layers_pv_gdf["total_pv"].min() vmax = layers_pv_gdf["total_pv"].max() layers_pv_gdf.plot(column="total_pv", ax=ax, cmap="YlOrRd", edgecolor="black", legend=True, legend_kwds={'label': "Total PV Generation per Layer (kWh)", 'orientation': "vertical"}, vmin=vmin, vmax=vmax) ax.set_title("Total PV Generation by Layer (All 4 Cities)", fontsize=14) ax.set_aspect('equal', 'box') plt.tight_layout() plt.savefig(figures_path / "layers_total_pv_all_4_cities.png", dpi=300) plt.close(fig) print("All plots generated successfully for the 4 cities and saved in the 'figures' folder.")