Un po' di tempo fa ho affrontato in un mio articolo il problema dell'aggiunta della coordinata altimetrica ad un dataset di punti usando QGIS. In questo JupyterNotebook affronterò il problema usando Python.

Prima di iniziare

Librerie

Per raggiungere l'obiettivo userò essenzialmente geopandas, shapely e rasterio. In verità tra gli import troverai anche pathlib e matplotlib ma non sono essenziali per la corretta riuscita del procedimento.

import pathlib

import pandas as pd
import geopandas as gpd
import rasterio as rio
from rasterio.plot import show, show_hist
import matplotlib.pyplot as plt
from shapely.geometry import Point

Fonti dati

Come DEM ho usato un tile del progetto TIN Italy a sua volta ritagliato e su cui ho generato con QGIS quattro punti random.

sample_data = pathlib.Path.cwd().parent.joinpath('sample_data/add_z')
dem = sample_data.joinpath('clip_dem.tif')
points = sample_data.joinpath('points.shp')

1. Quota altimetrica costante

E' il caso più semplice: ho un certo insieme di punti e tutti devono acquisire la stezza Z.

point_data = gpd.read_file(points)
point_data
fid geometry
0 1.0 POINT (405270.221 4512122.298)
1 2.0 POINT (405226.798 4512675.034)
2 3.0 POINT (404962.921 4512294.350)
3 4.0 POINT (404991.893 4512570.669)
point_data.plot(figsize=(10, 10))
<AxesSubplot:>
point_data.has_z
0    False
1    False
2    False
3    False
dtype: bool

Usando has_z ho verificato se i quattro punti sono effettivamente piani e quel False per ognuno di loro mi conferma che lo sono. Non resta che aggiungere la terza quota.

point_z_list = []

for _index, row in point_data.iterrows():
    point_geometry = row.geometry.coords
    # Estraggo la coppia di coordinate
    coordinates_plane = point_geometry[0]
    # Creo una lista con la coppia di coordinate ed aggiungo la z
    coordinates_3d = list(coordinates_plane)
    coordinates_3d.append(10)
    # Converto la lista in PointZ ed aggiungo l'oggetto alla lista di punti
    point_geometry_3d = Point(tuple(coordinates_3d))
    point_z_list.append([row.fid, point_geometry_3d])

# Ricostruisco il GeoDataFrame di punti
df = pd.DataFrame(point_z_list)
df.rename(columns={0:'fid', 1: 'geometry'}, inplace=True)
gdf = gpd.GeoDataFrame(df, geometry=df.geometry, crs=point_data.crs)
gdf
fid geometry
0 1.0 POINT Z (405270.221 4512122.298 10.000)
1 2.0 POINT Z (405226.798 4512675.034 10.000)
2 3.0 POINT Z (404962.921 4512294.350 10.000)
3 4.0 POINT Z (404991.893 4512570.669 10.000)
gdf.has_z
0    True
1    True
2    True
3    True
dtype: bool

Interrogando i dati nuovamente con has_z risultano essere 3D.

2. Quota altimetrica da DEM

Caso più complesso e sicuramente corrispondente alla realtà.

Il primo step è leggere la fonte dati raster. In verità dovrei leggere anche la fonte vettoriale ma siccome l'ho già fatto nell'esempio precedente evito di duplicare il passaggio.

raster_data = rio.open(dem)
raster_data
<open DatasetReader name='/home/max/DEV/Personal/PyGISBlog/sample_data/add_z/clip_dem.tif' mode='r'>
figure, axes = plt.subplots(figsize=(10, 10))
show(
    source=raster_data,
    cmap='terrain'
)
<AxesSubplot:>
show_hist(
    source=raster_data,
    title='Distribuzione delle quote',
    bins=10,
    alpha=0.5,
)

L'istogramma delle quote ci consente di capire che c'è una forte presenza di quote zero che corrispondono all'azzurro del raster precedentemente stampato.

point_z_list = []

for _index, row in point_data.iterrows():
    point_geometry = row.geometry.xy
    # Estraggo la coppia di coordinate
    x_coords = point_geometry[0][0]
    y_coords = point_geometry[1][0]
    # Individuo la riga e la colonna del raster in cui ricade
    # la coppia di coordinate.
    line, column = raster_data.index(x_coords, y_coords)
    # Estraggo il valore della cella dall'intersezione della
    # riga e colonna individuate al passaggio precedente.
    # E' il valore di quota che sto cercando.
    pixel_value = raster_data.read(1)[line, column]
    # Creo il PointZ ed aggiungo l'oggetto alla lista di punti
    point_geometry_3d = Point((x_coords, y_coords, pixel_value))
    point_z_list.append([row.fid, point_geometry_3d])

# Ricostruisco il GeoDataFrame di punti
df = pd.DataFrame(point_z_list)
df.rename(columns={0:'fid', 1: 'geometry'}, inplace=True)
gdf = gpd.GeoDataFrame(df, geometry=df.geometry, crs=point_data.crs)
gdf
fid geometry
0 1.0 POINT Z (405270.221 4512122.298 53.347)
1 2.0 POINT Z (405226.798 4512675.034 0.000)
2 3.0 POINT Z (404962.921 4512294.350 65.504)
3 4.0 POINT Z (404991.893 4512570.669 18.646)
gdf.has_z
0    True
1    True
2    True
3    True
dtype: bool

Anche in questo caso la verifica sulla tridimensionalità del GeoDataFrame ha dato esito positivo.

2bis. Quota altimetrica da DEM usando rioxarray

Aggiornamento del 07-04-2022

Procedura simile alla precedente con la sola differenza legata al tipo di libreria con cui vado a leggere il DEM; andrò ad usare infatti rioxarray estensione di xarray.

import rioxarray as rxr

new_raster_data = rxr.open_rasterio(dem).squeeze()
new_raster_data
<xarray.DataArray (y: 72, x: 55)>
array([[ 0.      ,  0.      ,  0.      , ...,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      , ...,  0.      ,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      , ...,  0.      ,  0.      ,  0.      ],
       ...,
       [69.77669 , 69.81791 , 69.8414  , ..., 15.89126 , 13.74339 ,  8.158463],
       [70.76064 , 69.98181 , 69.8414  , ..., 16.82376 , 13.96572 , 10.95792 ],
       [70.76064 , 69.98181 , 69.95621 , ..., 19.76735 , 16.77215 , 13.75869 ]],
      dtype=float32)
Coordinates:
    band         int64 1
  * x            (x) float64 4.049e+05 4.05e+05 4.05e+05 ... 4.055e+05 4.055e+05
  * y            (y) float64 4.513e+06 4.513e+06 ... 4.512e+06 4.512e+06
    spatial_ref  int64 0
Attributes:
    STATISTICS_MAXIMUM:        87.830078125
    STATISTICS_MEAN:           25.813233406838
    STATISTICS_MINIMUM:        0
    STATISTICS_STDDEV:         26.3760348731
    STATISTICS_VALID_PERCENT:  100
    _FillValue:                -9999.0
    scale_factor:              1.0
    add_offset:                0.0
new_point_z_list = []

for _index, row in point_data.iterrows():
    point_geometry = row.geometry.xy
    # Estraggo la coppia di coordinate
    x_coords = point_geometry[0][0]
    y_coords = point_geometry[1][0]

    # Individuo la riga e la colonna del DataArray in cui ricade
    # la coppia di coordinate.
    value = new_raster_data.sel(x=x_coords, y=y_coords, method="nearest")
    # Estraggo il valore della cella dall'intersezione della
    # riga e colonna individuate al passaggio precedente.
    # E' il valore di quota che sto cercando.
    pixel_value = value.data
    # Creo il PointZ ed aggiungo l'oggetto alla lista di punti
    point_geometry_3d = Point((x_coords, y_coords, pixel_value))
    new_point_z_list.append([row.fid, point_geometry_3d])

# Ricostruisco il GeoDataFrame di punti
new_df = pd.DataFrame(new_point_z_list)
new_df.rename(columns={0:'fid', 1: 'geometry'}, inplace=True)
new_gdf = gpd.GeoDataFrame(new_df, geometry=new_df.geometry, crs=point_data.crs)
new_gdf
fid geometry
0 1.0 POINT Z (405270.221 4512122.298 53.347)
1 2.0 POINT Z (405226.798 4512675.034 0.000)
2 3.0 POINT Z (404962.921 4512294.350 65.504)
3 4.0 POINT Z (404991.893 4512570.669 18.646)

Conclusione

Quelli riportati sono solo due dei metodi possibili per raggiungere l'obiettivo, sono quelli che ho individuato io. Se ne conosci altri, o se hai indicazioni da dare per migliorare il post, fallo scrivendo in Discussioni. Nota bene che per 2 è indispensabile che raster e punti abbiano lo stesso sistema di riferimento!

Newsletter

Se vuoi ricevere aggiornamenti sui prossimi articoli iscriviti alla newsletter!

Sponsorizza

Se vuoi sponsorizza il progetto andando su Sponsor this project nella home di questo repository.

Alla prossima!