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.
Prima di iniziare¶
Librerie¶
Per raggiungere l'obiettivo userò le librerie che seguono:
from pathlib import Path
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 = Path.cwd().parent.parent.joinpath('sample_data/viewshed')
dem = sample_data.joinpath('ispra_dtm_20.tif')
points = sample_data.joinpath('poi.shp')
output_folder = Path.cwd().parent.parent.joinpath('output_folder')
Path.mkdir(output_folder, exist_ok=True)
1. Lettura dei dati¶
points_data = gpd.read_file(points)
points_data
| fid | geometry | |
|---|---|---|
| 0 | 1 | POINT (458757.920 4498413.572) |
| 1 | 2 | POINT (462111.858 4502990.699) |
| 2 | 3 | POINT (459359.399 4506106.798) |
| 3 | 4 | POINT (457577.922 4495996.256) |
| 4 | 5 | POINT (447422.726 4500303.775) |
| 5 | 6 | POINT (461931.392 4502754.454) |
| 6 | 7 | POINT (464694.875 4495428.963) |
| 7 | 8 | POINT (457380.393 4500985.601) |
| 8 | 9 | POINT (452798.980 4493684.317) |
| 9 | 10 | POINT (455753.530 4496944.529) |
| 10 | 11 | POINT (456914.431 4498265.587) |
raster = rxr.open_rasterio(dem).squeeze()
raster
<xarray.DataArray (y: 1176, x: 1464)>
[1721664 values with dtype=int16]
Coordinates:
band int64 1
* x (x) float64 4.386e+05 4.387e+05 ... 4.679e+05 4.679e+05
* y (y) float64 4.513e+06 4.513e+06 ... 4.489e+06 4.489e+06
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
_FillValue: -32768
scale_factor: 1.0
add_offset: 0.0Il 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')
La quota massima è 1429 metri
1.1 Visualizzazione dei dati grezzi¶
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()
2. Analisi su un singolo punto¶
single_point = points_data[points_data['fid'] == 11]
single_point
| fid | geometry | |
|---|---|---|
| 10 | 11 | POINT (456914.431 4498265.587) |
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()