Verificare il potenziale fotovoltaico di un tetto
Usando un software GIS è possibile stimare se un tetto è o meno idoneo per un impianto fotovoltaico. E' possibile fare la stessa cosa in Python? Te lo dico qui!
- Prima di iniziare
- 1. Lettura dei dati
- 2. Analisi dell'inclinazione dei tetti
- 3. Analisi dell'orientamento dei tetti
- 4. Produzione energetica stimata
- Conclusione
Un po' di tempo fa ho effettuato un'analisi finalizzata a stimare il potenziale fotovoltaico di un tetto usando QGIS. Con questo articolo voglio ripercorrere lo stesso flusso di lavoro usando però Python!
import pathlib
import geopandas as gpd
import rioxarray as rxr
import contextily as cx
import matplotlib.pyplot as plt
from xrspatial.aspect import aspect
from xrspatial.slope import slope
sample_data = pathlib.Path.cwd().parent.joinpath('sample_data/photovoltaic_potential_analysis')
dsm = sample_data.joinpath('dsm.tif')
buildings = sample_data.joinpath('buildings_footprint.shp')
raster = rxr.open_rasterio(dsm, masked=True).squeeze()
raster
vector = gpd.read_file(buildings)
vector
fig, ax = plt.subplots(figsize=(10, 10))
raster.plot.imshow(
cmap='Blues',
)
plt.title("DSM")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")
aoi = vector.plot(alpha=0.75, color="blue", figsize=(10, 10))
cx.add_basemap(aoi, crs=vector.crs.to_string(), source=cx.providers.OpenStreetMap.Mapnik)
aoi.set_title("Area of Interest")
raster_clipped = raster.rio.clip(
geometries=vector.geometry, crs=raster.rio.crs, all_touched=True, drop=True
).squeeze()
fig, ax = plt.subplots(figsize=(10, 10))
raster_clipped.plot.imshow(
cmap='Blues',
)
plt.title("DEM of buildings footprint")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")
A questo punto posso concentrarmi solo sui tetti calcolando la loro pendenza.
rooftop_slope = slope(raster_clipped)
fig, ax = plt.subplots(figsize=(10, 10))
colormap = plt.cm.get_cmap('RdYlGn')
rooftop_slope.plot.imshow(
cmap=colormap.reversed(),
)
plt.title("Rooftop's slope")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")