hub/pv_generation.py

216 lines
8.1 KiB
Python
Raw Normal View History

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
input_file = "data/selected_buildings.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)
city = 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).enrich()
WeatherFactory('epw', city).enrich()
ExportsFactory('sra', city, sra_output_path).export()
sra_path = (sra_output_path / f'{city.name}_sra.xml').resolve()
subprocess.run(['sra', str(sra_path)])
ResultFactory('sra', city, sra_output_path).enrich()
tilt_angle = 37
solar_parameters = SolarCalculator(city=city,
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.buildings, random_assignation.residential_systems_percentage)
EnergySystemsFactory('montreal_future', city).enrich()
for building in city.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()
print("simulation done")
# PLOTTING STUFF
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import shape, Polygon
import os
# Paths
layers_file = "data/cmm_limites_avec_mtl.geojson" # Polygon layers
main_dir = os.path.abspath(".")
# Read layer polygons
layers_gdf = gpd.read_file(layers_file)
if layers_gdf.crs is None:
layers_gdf.set_crs(epsg=32188, inplace=True)
buildings_data = []
for b in city.buildings:
yearly_pv = b.pv_generation['year'][0]/1000
if isinstance(yearly_pv, list):
yearly_pv = sum(yearly_pv)
# Extract the ground surface polygon from the building
# Assuming the ground surface is at index 0. If not,
# you may need to find it by checking the surface type.
if not b.surfaces:
continue # no surfaces, skip
ground_surface = b.surfaces[0]
coords_3d = ground_surface.solid_polygon.coordinates
# Convert 3D coordinates (X,Y,Z) into 2D (X,Y)
coords_2d = [(c[0], c[1]) for c in coords_3d]
# Create a Shapely polygon from these coordinates
footprint_poly = Polygon(coords_2d)
buildings_data.append({
"geometry": footprint_poly,
"yearly_pv": yearly_pv
})
# Create the GeoDataFrame from the buildings_data
buildings_gdf = gpd.GeoDataFrame(buildings_data, crs="EPSG:26911") # Change CRS if needed
# If layers_gdf is EPSG:32188 and buildings_gdf is EPSG:4326, reproject buildings:
buildings_gdf = buildings_gdf.to_crs(layers_gdf.crs)
# Loop through selected layers (e.g., layer_6, layer_7, etc.)
# Let's say we have a list of layers we want to process:
layers_of_interest = [f"layer_{i}" for i in range(6, 21)] # example from layer_6 to layer_20
for layer_name in layers_of_interest:
layer_poly = layers_gdf[layers_gdf["region"] == layer_name]
if layer_poly.empty:
continue
# There might be multiple polygons per layer_name; dissolve them into one for simplicity
layer_union = layer_poly.unary_union
if layer_union.is_empty:
continue
# Select buildings within this layer polygon
# Spatial join or mask
buildings_in_layer = buildings_gdf[buildings_gdf.geometry.within(layer_union)]
if buildings_in_layer.empty:
continue
# Plot the layer polygon and buildings
fig, ax = plt.subplots(figsize=(20, 20), dpi=600)
# Plot polygon
layer_poly.plot(ax=ax, facecolor="none", edgecolor="black", linewidth=1)
# Buildings colored by pv_generation
# Let's set a colormap normalized by min/max PV
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.pdf"), dpi=300)
plt.close(fig)
# After generating individual layer plots, create a combined plot of all buildings (all layers):
fig, ax = plt.subplots(figsize=(20, 20), dpi=1200)
# If you want all buildings, possibly also plot the entire layers polygon union as context
layers_gdf.plot(ax=ax, facecolor="none", edgecolor="grey", linewidth=0.5)
vmin = buildings_gdf["yearly_pv"].min()
vmax = buildings_gdf["yearly_pv"].max()
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", fontsize=14)
ax.set_aspect('equal', 'box')
plt.tight_layout()
plt.savefig(os.path.join(main_dir, "all_layers_pv_generation.png"), dpi=1200)
plt.close(fig)
# Finally, sum PV production by layer polygon and create a plot of polygons colored by total PV
layer_sums = []
for idx, row in layers_gdf.iterrows():
region_name = row["region"]
poly_geom = row["geometry"]
# Buildings inside this polygon
b_in_poly = buildings_gdf[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", fontsize=14)
ax.set_aspect('equal', 'box')
plt.tight_layout()
plt.savefig(os.path.join(main_dir, "layers_total_pv.png"), dpi=300)
plt.close(fig)
print("All plots generated successfully.")