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:

In [1]:

from pathlib import Path
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

import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore')

Fonti dati¶

In [2]:

sample_data = Path.cwd().parent.parent.joinpath('sample_data/photovoltaic_potential_analysis')
dsm = sample_data.joinpath('dsm.tif')
buildings = sample_data.joinpath('buildings_footprint.shp')

1. Lettura dei dati¶

In [3]:

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

Out[3]:

<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:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

In [4]:

vector = gpd.read_file(buildings)
vector

Out[4]:

	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¶

In [5]:

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

Out[5]:

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

No description has been provided for this image

In [6]:

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

Out[6]:

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.

In [7]:

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

Out[7]:

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

A questo punto posso concentrarmi solo sui tetti calcolando la loro pendenza.

In [8]:

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

Out[8]:

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

Le pendenze fanno riferimento al piano cartografico. Quindi 0° fa riferimento ad un tetto piano, con 85.8° di pendenza abbiamo un tetto quasi verticale e quindi ortogonale al piano cartografico.

La resa di un pannello fotovoltaico varia in base all'inclinazione rispetto ai raggi solari. Alle latitudini italiane la resa ottimale in fase di produzione di energia si ha quando i pannelli sono installati con una pendenza che varia tra i 30° ed i 40°. La pendenza ottimale è direttamente proporzionale alla latitudine; più ci allontaniamo dall'equatore, più essa aumenta. Quindi nel sud Italia la pendenza ottimale si aggira intorno ai 30° mentre nel nord è intorno ai 40°. Non mi interessa entrare troppo nel dettaglio tecnico legato all'installazione di un impianto fotovoltaico ma dare una stima, una indicazione di produzione energetica, ma il passaggio sull'inclinazione è importante perchè con queste informazioni andrò ad effettuare le analisi successive suddividendo le pendenze dei tetti in tre classi:

0 - 3°
3° - 40°
40° - 90°

Questa suddivisione mi consentirà di individuare sia i tetti piani(quelli che rientrano nella prima classe) che i tetti con la giusta inclinazione per una installazione in pendenza ottimale dell'impianto, consentendomi di escludere quelli con pendenze non ottimali(terza ed ultima classe).

Qui magari, con questa esclusione, storce il naso un installatore di impianti fotovoltaici perchè ci sono soluzioni tecniche che consentono di ovviare il problema della pendenza eccessiva del tetto. Conosco questi aspetti ma a me interessa dare una visione di metodologia di analisi con questo articolo.

2.1 Tetti piani¶

In [9]:

plane_rooftop = rooftop_slope.where(rooftop_slope <= 3)

fig, ax = plt.subplots(figsize=(10, 10))
colormap = plt.cm.get_cmap('RdYlGn')

plane_rooftop.plot.imshow(
    cmap=colormap.reversed(),
)
plt.title("Plane rooftops")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")

Out[9]:

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

2.2 Tetti con inclinazione ottimale¶

In [10]:

optimum_rooftop_inclination = rooftop_slope.where((rooftop_slope > 3) & (rooftop_slope <= 40))

fig, ax = plt.subplots(figsize=(10, 10))
colormap = plt.cm.get_cmap('RdYlGn')

optimum_rooftop_inclination.plot.imshow(
    cmap=colormap.reversed(),
)
plt.title("Rooftops with slope between 3° and 40°")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")

Out[10]:

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

3. Analisi dell'orientamento dei tetti¶

Individuate le pendenze del tetto, passiamo a considerare l'aspetto relativo all'esposizione, o orientamento, del tetto.

Per il discorso lagato alla resa migliore in fase di produzione di energia, anche l'orientamento del tetto è un fattore chiave. Per un tetto piano è insignificante parlare di esposizione, concetto che diventa fondamentale per un tetto inclinato. Anche qui ci sono soluzioni tecniche che possono risolvere problemi di esposizione ma, come scritto poco sopra, non è questo il luogo per approfondire questi aspetti.

In [11]:

rooftop_aspect = aspect(raster_clipped)

fig, ax = plt.subplots(figsize=(10, 10))

rooftop_aspect.plot.imshow(
    cmap='Set2',
)
plt.title("Rooftop's aspect")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")

Out[11]:

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

L'esposizione ottimale è il SUD(tra 135° e 225°), quella sicuramente da scartare è il NORD(tra 0° e 45°, tra 315° e 360°). Rosa dei venti

Io andrò ad estrarre solo i tetti con esposizione da EST a SUD-EST, da SUD-EST a SUD-OVEST ad OVEST, quindi:

90° - 135°
135° - 225°
225° - 270°

3.1 Tetti esposti a SUD¶

In [12]:

sud_rooftop_aspect = rooftop_aspect.where((rooftop_aspect >= 135) & (rooftop_aspect <= 225))

fig, ax = plt.subplots(figsize=(10, 10))

sud_rooftop_aspect.plot.imshow(
    cmap='Set2',
)
plt.title("Rooftops facing to SUD")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")

Out[12]:

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

3.2 Tetti esposti ad EST¶

In [13]:

est_rooftop_aspect = rooftop_aspect.where((rooftop_aspect >= 90) & (rooftop_aspect < 135))

fig, ax = plt.subplots(figsize=(10, 10))

est_rooftop_aspect.plot.imshow(
    cmap='Set2',
)
plt.title("Rooftops facing to EST")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")

Out[13]:

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

3.2 Tetti esposti ad OVEST¶

In [14]:

ovest_rooftop_aspect = rooftop_aspect.where((rooftop_aspect > 225) & (rooftop_aspect <= 270))

fig, ax = plt.subplots(figsize=(10, 10))

ovest_rooftop_aspect.plot.imshow(
    cmap='Set2',
)
plt.title("Rooftops facing to OVEST")
plt.ylabel("Y coordinates (meters)")
plt.xlabel("X coordinates (meters)")

Out[14]:

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

4. Produzione energetica stimata¶

Ora che ho una chiara idea delle aree di intervento procederò a stimare la produzione energetica dei tetti presenti nell'area in esame. Prima però è necessario affinare ulteriormente i dati.

4.1 Put them all togheter¶

In [15]:

target_slope_rooftop = rooftop_slope.where(rooftop_slope <= 40)

fig, axs = plt.subplots(1, 2, figsize=(20, 10))
colormap = plt.cm.get_cmap('RdYlGn')

target_slope_rooftop.plot.imshow(
    cmap=colormap.reversed(),
    ax=axs[0],
)
axs[0].set_title("Target rooftops based on slope")

rooftop_slope.plot.imshow(
    cmap=colormap.reversed(),
    ax=axs[1]
)
axs[1].set_title("Rooftop's slope")

Out[15]:

Text(0.5, 1.0, "Rooftop's slope")

In [16]:

target_aspect_rooftop = rooftop_aspect.where((rooftop_aspect >= 90) & (rooftop_aspect <= 270))

fig, axs = plt.subplots(1, 2, figsize=(20, 10))

target_aspect_rooftop.plot.imshow(
    cmap='Set2',
    ax=axs[0]
)
axs[0].set_title("Target rooftops based on aspect")

rooftop_aspect.plot.imshow(
    cmap='Set2',
    ax=axs[1]
)
axs[1].set_title("Rooftop's aspect")

Out[16]:

Text(0.5, 1.0, "Rooftop's aspect")

4.2 Individuazione dei tetti con pendenza ed esposizione ottimali¶

In [17]:

target_slope_rooftop_mask = ((target_slope_rooftop / target_slope_rooftop) >= 1).astype('uint8')
target_aspect_rooftop_mask = ((target_aspect_rooftop / target_aspect_rooftop) >= 1).astype('uint8')

fig, axs = plt.subplots(1, 2, figsize=(20, 10))

target_slope_rooftop_mask.plot.imshow(
    cmap="gray",
    ax=axs[0],

)
axs[0].set_title("Target rooftops based on slope")

target_aspect_rooftop_mask.plot.imshow(
    cmap="gray",
    ax=axs[1]
)
axs[1].set_title("Target rooftops based on aspect")

Out[17]:

Text(0.5, 1.0, 'Target rooftops based on aspect')

In [18]:

target_rooftop = target_slope_rooftop_mask * target_aspect_rooftop_mask

fig, ax = plt.subplots(figsize=(10, 10))

target_rooftop.plot.imshow(
    cmap='gray',
)
plt.title("Target rooftops")

Out[18]:

Text(0.5, 1.0, 'Target rooftops')

4.3 Stima del potenziale produttivo¶

Ora che è chiaro quali sono i tetti potenzialmente utilizzabili è possibile calcolarne la produzione.

Il primo passo da seguire è convertire il DataArray creato in precedenza in un DataFrame.

In [19]:

target_rooftop_df = target_rooftop.to_dataframe(name="value").reset_index()
target_rooftop_df = target_rooftop_df[target_rooftop_df['value'] == 1]
target_rooftop_df = target_rooftop_df[['x', 'y', 'value']]
target_rooftop_df

Out[19]:

	x	y	value
1951	444946.808	4529404.911	1
1952	444947.808	4529404.911	1
1953	444948.808	4529404.911	1
1954	444949.808	4529404.911	1
2455	444950.808	4529403.911	1
...	...	...	...
248266	444761.808	4528911.911	1
248267	444762.808	4528911.911	1
248268	444763.808	4528911.911	1
248269	444764.808	4528911.911	1
248767	444762.808	4528910.911	1

25141 rows × 3 columns

In questo modo si è ottenuto un oggetto che contiene i centroidi dei pixel del raster precedente, con i soli valori di interesse. Il passo successivo è convertire in GeoDataFrame il dataset ottenuto in precedenza.

In [20]:

gdf_point = gpd.GeoDataFrame(
    target_rooftop_df.value,
    geometry=gpd.points_from_xy(target_rooftop_df.x,target_rooftop_df.y),
    crs=raster.rio.crs.to_string()
)
gdf_point

Out[20]:

	value	geometry
1951	1	POINT (444946.808 4529404.911)
1952	1	POINT (444947.808 4529404.911)
1953	1	POINT (444948.808 4529404.911)
1954	1	POINT (444949.808 4529404.911)
2455	1	POINT (444950.808 4529403.911)
...	...	...
248266	1	POINT (444761.808 4528911.911)
248267	1	POINT (444762.808 4528911.911)
248268	1	POINT (444763.808 4528911.911)
248269	1	POINT (444764.808 4528911.911)
248767	1	POINT (444762.808 4528910.911)

25141 rows × 2 columns

Si è ottenuto un GeoDataFrame di punti, ma a me interessano i poligoni, quindi creo un buffer intorno ai punti in modo da riprodurre i pixel del raster con cui ho iniziato l'analisi. Quindi devo andare ad individuare la risoluzione del mio raster di partenza:

In [21]:

x, y = raster.rio.resolution()
print(f"Il raster di partenza ha una risoluzione di {x} x {y*-1} metri")

Il raster di partenza ha una risoluzione di 1.0 x 1.0 metri

Ora posso creare i miei poligoni.

In [22]:

gdf_polygon = gdf_point.buffer(distance=0.5, cap_style=3)
polygons = gdf_polygon.unary_union

gdf = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polygons), crs=gdf_polygon.crs).explode(index_parts=True)
gdf.insert(loc=1, column="superficie_mq", value=gdf.geometry.area)
gdf

Out[22]:

		geometry	superficie_mq
0	0	POLYGON ((444503.308 4528932.411, 444503.308 4...	1.0
	1	POLYGON ((444505.308 4528934.411, 444506.308 4...	1.0
	2	POLYGON ((444501.308 4528936.411, 444501.308 4...	2.0
	3	POLYGON ((444504.308 4528937.411, 444504.308 4...	5.0
	4	POLYGON ((444500.308 4528947.411, 444500.308 4...	8.0
	...	...	...
	2284	POLYGON ((444989.308 4529288.411, 444989.308 4...	1.0
	2285	POLYGON ((444980.308 4529322.411, 444980.308 4...	1.0
	2286	POLYGON ((444980.308 4529325.411, 444981.308 4...	1.0
	2287	POLYGON ((444992.308 4529354.411, 444992.308 4...	1.0
	2288	POLYGON ((444985.308 4529360.411, 444984.308 4...	18.0

2289 rows × 2 columns

Ho ottenuto un GeoDataFrame con i dati di mio interesse. Quando scrissi l'articolo del blog, oramai 6 annifa, 1 kwp occupava 7 metri quadri in media; non conosco i risvolti del mercato attuale ma ricordo che con la miniaturizzazione si andava verso la riduzione della superficie utile per produrre 1 kWp. Anche per questo articolo mi atterrò al dato del 2018 ma, prima di calcolare la produzione potenziale di energia dell'are in esame, andrò a filtrare tutto ciò che è inferiore a 7 metri quadri.

In [23]:

final_analysis = gdf[gdf.superficie_mq >= 7].reset_index()
final_analysis = final_analysis[['geometry', 'superficie_mq']]
final_analysis.insert(loc=2, column="potential_production_kw", value=final_analysis.superficie_mq*7)
final_analysis

Out[23]:

	geometry	superficie_mq	potential_production_kw
0	POLYGON ((444500.308 4528947.411, 444500.308 4...	8.0	56.0
1	POLYGON ((444506.308 4528961.411, 444507.308 4...	8.0	56.0
2	POLYGON ((444518.308 4528970.411, 444519.308 4...	35.0	245.0
3	POLYGON ((444518.308 4528967.411, 444519.308 4...	9.0	63.0
4	POLYGON ((444524.308 4529036.411, 444525.308 4...	24.0	168.0
...	...	...	...
641	POLYGON ((444988.308 4529225.411, 444989.308 4...	12.0	84.0
642	POLYGON ((444985.308 4529260.411, 444986.308 4...	16.0	112.0
643	POLYGON ((444985.308 4529279.411, 444986.308 4...	9.0	63.0
644	POLYGON ((444987.308 4529283.411, 444987.308 4...	22.0	154.0
645	POLYGON ((444985.308 4529360.411, 444984.308 4...	18.0	126.0

646 rows × 3 columns

In [24]:

fig, ax = plt.subplots(figsize=(10, 10))

final_view = final_analysis.plot(column="potential_production_kw", cmap="gnuplot", ax=ax, legend=True)
cx.add_basemap(final_view, crs=final_analysis.crs.to_string(), source=cx.providers.OpenStreetMap.Mapnik)
final_view.set_title("Energy production map kWp")

Out[24]:

Text(0.5, 1.0, 'Energy production map kWp')

Conclusione¶

La stima del potenziale produttivo di un tetto è sicuramente un passo fondamentale per capire se il tetto in esame è da considerare per una analisi più approfondita o meno. Nel processo di analisi che ho esposto ho volutamente evitato di scendere troppo nel dettaglio di analisi proprio perchè:

avevo interesse ad esporre un flusso di analisi;
avevo interesse a stimare dei valori.

Ad esempio un aspetto su cui ho volutamente sorvolato è il calcolo delle superfici inclinate che mi avrebbe restituito un valore di produttività per tetto differente.

Verificare il potenziale fotovoltaico di un tetto

Prima di iniziare¶

Librerie¶

Fonti dati¶

Contenuti¶

1. Lettura dei dati¶

1.1 Visualizzazione dei dati grezzi¶

2. Analisi dell'inclinazione dei tetti¶

2.1 Tetti piani¶

2.2 Tetti con inclinazione ottimale¶

3. Analisi dell'orientamento dei tetti¶

3.1 Tetti esposti a SUD¶

3.2 Tetti esposti ad EST¶

3.2 Tetti esposti ad OVEST¶

4. Produzione energetica stimata¶

4.1 Put them all togheter¶

4.2 Individuazione dei tetti con pendenza ed esposizione ottimali¶

4.3 Stima del potenziale produttivo¶

Conclusione¶

Offrimi una birra!

Prima di iniziare¶

Librerie¶

Fonti dati¶

Contenuti¶

1. Lettura dei dati¶

1.1 Visualizzazione dei dati grezzi¶

2. Analisi dell'inclinazione dei tetti¶

2.1 Tetti piani¶

2.2 Tetti con inclinazione ottimale¶

3. Analisi dell'orientamento dei tetti¶

3.1 Tetti esposti a SUD¶

3.2 Tetti esposti ad EST¶

3.2 Tetti esposti ad OVEST¶

4. Produzione energetica stimata¶

4.1 Put them all togheter¶

4.2 Individuazione dei tetti con pendenza ed esposizione ottimali¶

4.3 Stima del potenziale produttivo¶

Conclusione¶

Altri articoli della stessa categoria

Mappatura dell'andamento del PM10 nel Nord Italia

Analisi dell'andamento del PM10 nel Nord Italia

Gestione di un file NetCDF

Siti archeologici italiani

Linea meno lunga in una lista

Offrimi una birra!