<jupyter_start><jupyter_text>Da Point a PointZ> Con questo articolo ti spiegherò come convertire facilmente un vettore di punti 2D in uno 3D usando Python.- toc: true- badges: true- comments: false- author: <NAME>- categories: [point]- image: images/logo.png Un po' di tempo fa ho affrontato in un [mio articolo](https://massimilianomoraca.it/blog/gis/da-point-a-pointz/) 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  LibreriePer 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.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Fonti datiCome [DEM](https://massimilianomoraca.it/blog/gis/3d-open-source/) ho usato un tile del progetto [TIN Italy](http://tinitaly.pi.ingv.it/Download_Area2.html) a sua volta ritagliato e su cui ho generato con QGIS quattro punti random.<jupyter_code>sample_data = pathlib.Path.cwd().parent.joinpath('sample_data/add_z')
dem = sample_data.joinpath('clip_dem.tif')
points = sample_data.joinpath('points.shp')<jupyter_output><empty_output><jupyter_text>1. Quota altimetrica costanteE' il caso più semplice: ho un certo insieme di punti e tutti devono acquisire la stezza Z.<jupyter_code>point_data = gpd.read_file(points)
point_data
point_data.plot(figsize=(10, 10))
point_data.has_z<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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
gdf.has_z<jupyter_output><empty_output><jupyter_text>Interrogando i dati nuovamente con `has_z` risultano essere 3D.  2. Quota altimetrica da DEMCaso 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.<jupyter_code>raster_data = rio.open(dem)
raster_data
figure, axes = plt.subplots(figsize=(10, 10))
show(
    source=raster_data,
    cmap='terrain'
)
show_hist(
    source=raster_data,
    title='Distribuzione delle quote',
    bins=10,
    alpha=0.5,
)<jupyter_output><empty_output><jupyter_text>L'istogramma delle quote ci consente di capire che c'è una forte presenza di quote zero che corrispondono all'azzurro del raster precedentemente stampato.<jupyter_code>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
gdf.has_z<jupyter_output><empty_output><jupyter_text>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](https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html) estensione di [xarray](https://docs.xarray.dev/en/stable/).<jupyter_code>import rioxarray as rxr

new_raster_data = rxr.open_rasterio(dem).squeeze()
new_raster_data
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<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08625204582651391, 'fraction_numerical': 0.00900163666121113, 'mean_word_length': 5.287037037037037, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '654846', 'license': 'Apache-2.0', 'n_tokens_mistral': 2242, 'n_tokens_neox': 2217, 'n_words': 597, 'path': '_notebooks/2022-03-31-da-point-a-pointz.ipynb', 'repo_name': 'MaxDragonheart/PyGIS-Blog', 'size': '6123'}