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_without_consumption import PvSystemProduction from scripts import random_assignation import subprocess from pathlib import Path import geopandas as gpd import matplotlib.pyplot as plt from shapely.geometry import Polygon import os import pandas as pd 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" 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) 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() 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: PvSystemProduction(building=building, pv_system=None, battery=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).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 simulations done") # PLOTTING STUFF layers_file = "data/cmm_limites_avec_mtl.geojson" main_dir = os.path.abspath(".") layers_gdf = gpd.read_file(layers_file) if layers_gdf.crs is None: layers_gdf.set_crs(epsg=32188, inplace=True) 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) else: yearly_pv = yearly_pv 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 }) gdf = gpd.GeoDataFrame(buildings_data, crs="EPSG:26911") gdf = gdf.to_crs(layers_gdf.crs) 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], ignore_index=True), crs=gdf_city1.crs) 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(os.path.join(main_dir, 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(os.path.join(main_dir, "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(os.path.join(main_dir, "layers_total_pv_all_4_cities.png"), dpi=300) plt.close(fig) print("All plots generated successfully for the 4 cities.")