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:

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

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:
    _FillValue:    -32768.0
    scale_factor:  1.0
    add_offset:    0.0

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

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.

2.1 Quota dell'osservatore

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
array(450, dtype=int16)

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
<xarray.DataArray (y: 1176, x: 1464)>
array([[87.78651279, 87.78502168, 87.78352919, ..., -1.        ,
        -1.        , -1.        ],
       [87.78533465, 87.78384115, 87.78234627, ..., -1.        ,
        -1.        , -1.        ],
       [87.78415626, 87.78266038, 87.7811631 , ..., -1.        ,
        -1.        , -1.        ],
       ...,
       [87.47039509, 87.46816885, 87.46593917, ..., -1.        ,
        -1.        , -1.        ],
       [87.47149388, 87.46927054, 87.46704378, ..., -1.        ,
        -1.        , -1.        ],
       [87.47259367, 87.47037324, 87.46814939, ..., -1.        ,
        -1.        , -1.        ]])
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:
    _FillValue:    -32768.0
    scale_factor:  1.0
    add_offset:    0.0
view_sel = view.where(view >= 0)
view_sel
<xarray.DataArray (y: 1176, x: 1464)>
array([[87.78651279, 87.78502168, 87.78352919, ...,         nan,
                nan,         nan],
       [87.78533465, 87.78384115, 87.78234627, ...,         nan,
                nan,         nan],
       [87.78415626, 87.78266038, 87.7811631 , ...,         nan,
                nan,         nan],
       ...,
       [87.47039509, 87.46816885, 87.46593917, ...,         nan,
                nan,         nan],
       [87.47149388, 87.46927054, 87.46704378, ...,         nan,
                nan,         nan],
       [87.47259367, 87.47037324, 87.46814939, ...,         nan,
                nan,         nan]])
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:
    _FillValue:    -32768.0
    scale_factor:  1.0
    add_offset:    0.0
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()