hub/pv_generation.py

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.")