Analisi di visibilità
Un’analisi di visibilità consente di stabilire quali sono le porzioni di paesaggio visibili da un osservatore posto in un determinato luogo ed ad una determinata quota.
- Prima di iniziare
- 1. Lettura dei dati
- 2. Analisi su un singolo punto
- 3. Analisi multipla
- Conclusione
L’analisi di visibilità è un argomento molto interessante da trattare con Python. In un mio articolo l’ho affrontata usando QGIS, qui oggi vediamo come raggiungere lo stesso obiettivo con Python.
import pathlib
import geopandas as gpd
import rioxarray as rxr
import rasterio as rio
from rasterio.plot import show
import matplotlib.pyplot as plt
from xrspatial import hillshade, viewshed
from xarray.plot import imshow
Fonti dati
Come DEM ho usato un ritaglio del progetto DEM20 di ISPRA su cui ho posizionato con QGIS undici punti; nell'articolo del blog uso l'undicesimo, qui proverò ad usarli tutti.
sample_data = pathlib.Path.cwd().parent.joinpath('sample_data/viewshed')
dem = sample_data.joinpath('ispra_dtm_20.tif')
points = sample_data.joinpath('poi.shp')
points_data = gpd.read_file(points)
points_data
raster = rxr.open_rasterio(dem).squeeze()
raster
Il vettore dei punti è molto semplice, ha una colonna con le geometrie ed una che identifica il singolo punto; quest'ultima informazione la userò più avanti. Il DEM è rettangolare e monobanda ed in effetti mi aspetto che sia monobanda poichè è quella che contiene le informazioni altimetriche in questo caso, con squeeze() estraggo la banda.
Il DEM ha come valore minimo di quota 0 metri, ma quale è la quota massima?
max_h = int(raster.max().data)
print(f'La quota massima è {max_h} metri')
raster_data = rio.open(dem)
fig, ax = plt.subplots(figsize=(10, 10))
show(
source=raster_data,
cmap='tab20c_r',
contour=True,
ax=ax,
transform=raster_data.transform
)
show(
source=raster_data,
cmap='terrain_r',
ax=ax,
transform=raster_data.transform
)
for _index, row in points_data.iterrows():
coordinates = row.geometry.xy
ax.scatter(
*coordinates,
s=1000,
marker="*",
facecolor='red',
edgecolor='black'
)
plt.text(
x=coordinates[0][0],
y=coordinates[1][0],
s=row.fid,
fontdict=dict(color='black', size=10),
bbox=dict(facecolor='white', alpha=0.75)
)
plt.show()
single_point = points_data[points_data['fid'] == 11]
single_point
Ho selezionato solo il punto con id 11, è quello che ho usato nell'articolo del mio blog. Prima di partire con l'analisi è bene sapere che il risultato di una analisi di visibilità, detta anche viewshed, è un raster con valori che vanno da 0 a 180. Questi sono gradi e fanno riferimento alla visuale verticale dal punto di osservazione. Un valore pari a 0° sta ad indicare che quel pixel è proprio sotto il punto di osservazione, un valore pari a 90° sta a significare che il pixel è sull'orizzonte del punto di osservazione ed un valore pari a 180°, di conseguenza, indica che il pixel è proprio sopra il punto di osservazione.
make_hillshade = hillshade(raster)
fig, ax = plt.subplots(figsize=(20, 10))
imshow(
make_hillshade,
ax=ax,
cmap='gist_gray',
)
imshow(
raster,
ax=ax,
cmap='terrain',
vmin=1,
vmax=max_h,
levels=100,
alpha=0.75
)
for _index, row in single_point.iterrows():
coordinates = row.geometry.xy
ax.scatter(
*coordinates,
s=1000,
marker="*",
facecolor='red',
edgecolor='black'
)
plt.text(
x=coordinates[0][0],
y=coordinates[1][0],
s=row.fid,
fontdict=dict(color='black', size=10),
bbox=dict(facecolor='white', alpha=0.75)
)
plt.show()
Per una migliore resa visiva ho inserito la hillshade sotto al DEM tematizzato. Questa volta anzicchè usare la funzione show di rasterio ho preferito usare imshow di xarray per sperimentarla un po' visto che non l'avevo mai usata. E' interessante il fatto che aggiunga automaticamente la colorbar per ogni raster, devo capire come fare in modo che non venga visualizzata sempre e come disporla in altro modo.
La quota dell'osservatore è un dato fondamentale per portare a termine una analisi di visibilità. Può essere la quota del punto in cui si trova l'osservatore o una quota fittizia compatibile con la topografia dell'area(l'osservatore fittizio non può essere a 100 metri di altezza se stiamo esaminando un'area in cui la quota minima è 500 metri). Quello che manca nei miei dati è proprio la quota!
Ma non è un problema perchè mi viene in soccorso un'esperianza già fatta di recente e che vado a riprodurre:
x_coord = single_point.iloc[0].geometry.xy[0][0]
y_coord = single_point.iloc[0].geometry.xy[1][0]
dataarray_value = raster.sel(x=x_coord, y=y_coord, method="nearest")
pixel_value = dataarray_value.data
pixel_value
Ora che ho tutto quello che mi serve posso effettuare l'analisi di visibilità.
view = viewshed(raster, x=x_coord, y=y_coord, observer_elev=pixel_value)
view
view_sel = view.where(view >= 0)
view_sel
fig, ax = plt.subplots(figsize=(20, 10))
imshow(view_sel, ax=ax)
plt.show()
Analisi effettuata, ora sono curioso di vedere quali aree ricadono nel campo visivo dell'osservatore.
fig, ax = plt.subplots(figsize=(20, 10))
imshow(
make_hillshade,
ax=ax,
cmap='gist_gray',
)
imshow(
raster,
ax=ax,
cmap='terrain',
vmin=1,
vmax=max_h,
levels=100,
alpha=0.75
)
imshow(view_sel, ax=ax, cmap='Wistia')
ax.scatter(
*single_point.iloc[0].geometry.xy,
s=1000,
marker="*",
facecolor='red',
edgecolor='black'
)
plt.text(
x=coordinates[0][0],
y=coordinates[1][0],
s=row.fid,
fontdict=dict(color='black', size=10),
bbox=dict(facecolor='white', alpha=0.75)
)
plt.show()