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!

Prima di iniziare

Librerie

Per raggiungere l'obiettivo userò le librerie che seguono:

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

Fonti dati 

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')

1. Lettura dei dati 

raster = rxr.open_rasterio(dsm, masked=True).squeeze()
raster

<xarray.DataArray (y: 500, x: 500)>
[250000 values with dtype=float32]
Coordinates:
    band         int64 1
  * x            (x) float64 4.445e+05 4.445e+05 4.445e+05 ... 4.45e+05 4.45e+05
  * y            (y) float64 4.529e+06 4.529e+06 ... 4.529e+06 4.529e+06
    spatial_ref  int64 0
Attributes:
    STATISTICS_MAXIMUM:        62.470001220703
    STATISTICS_MEAN:           38.810858080246
    STATISTICS_MINIMUM:        26.229999542236
    STATISTICS_STDDEV:         7.799208310175
    STATISTICS_VALID_PERCENT:  100
    scale_factor:              1.0
    add_offset:                0.0
vector = gpd.read_file(buildings)
vector
id layer tipoedific uso sup_m2 altezza altez_dsm geometry
0 0 0201N-Edificio Generico Nuovo Edificio generico Abitazione 54 4.67 36.77 POLYGON ((444654.507 4529354.495, 444651.839 4...
1 1 0201N-Edificio Generico Nuovo Edificio generico Abitazione 49 7.60 40.79 POLYGON ((444639.838 4529162.951, 444638.794 4...
2 2 0208N-Baracca Nuovo Baracca Altro 23 3.30 38.11 POLYGON ((444726.330 4529069.704, 444730.929 4...
3 3 0208I-Baracca Invariato Baracca Altro 22 3.19 37.95 POLYGON ((444523.992 4529241.197, 444523.992 4...
4 4 0208I-Baracca Invariato Baracca Altro 55 3.96 38.68 POLYGON ((444518.076 4529233.077, 444518.192 4...
... ... ... ... ... ... ... ... ...
312 312 0201I-Edificio Generico Invariato Edificio generico Abitazione 187 10.57 42.86 POLYGON ((444676.418 4529355.727, 444676.008 4...
313 313 0201M-Edificio Generico Modificato Edificio generico Abitazione 323 7.41 42.35 POLYGON ((444545.683 4529221.436, 444524.605 4...
314 314 0201I-Edificio Generico Invariato Edificio generico Abitazione 370 2.88 35.15 POLYGON ((444723.780 4529165.567, 444719.890 4...
315 315 0201M-Edificio Generico Modificato Edificio generico Abitazione 2389 11.78 45.07 POLYGON ((444848.335 4529111.537, 444850.324 4...
316 316 0201M-Edificio Generico Modificato Edificio generico Abitazione 2389 11.78 45.07 POLYGON ((444860.572 4529113.825, 444863.015 4...

317 rows × 8 columns

1.1 Visualizzazione dei dati grezzi 

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

Text(0.5, 0, '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")

Text(0.5, 1.0, 'Area of Interest')

2. Analisi dell'inclinazione dei tetti

Come primo passo in questa analisi estrarrò dal DSM i soli edifici.

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

Text(0.5, 0, '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)")

Text(0.5, 0, 'X coordinates (meters)')