309 lines
11 KiB
Python
309 lines
11 KiB
Python
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], ignore_index=True),
|
|
crs=gdf_city1.crs)
|
|
|
|
print("Dataframe for all 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.")
|