[GIS] Statistiche zonali condizionate

[pigreco]

1. ho un raster delle pendenze riclassificato con valori 1,2,3,4;
2. un vettore poligonale che copre il raster;
3. che tool/procedura usare per calcolare una sorta di statistiche zonali? ovvero: al poligonale devo aggiungere 4 attributi (pend_1, pend_2, pend_3 e pend_4) uno per ogni valore riclassificato e popolarli con il conteggio dei relativi valori dei pixel che vi ricadono dentro (pend_1 = conteggio dei soli pixel con valore 1; pend_2 = conteggio dei soli pixel con valore 2 e cosi via);

sotto uno screenshot:

io lo ho risolto velocemente (ma per piccole aree) convertendo il raster in vettore (da pixel a punti) e usando (4 volte) le espressioni di QGIS:

aggregate(
layer:='punti',
aggregate:='count',
expression:="value",
filter:= "value"=1 and intersects($geometry, geometry(@parent)))

come potrei ovviare per raster molto più grandi?

allego dati di esempio zonal.zip

[aborruso]

Caro @pigreco , proverei a usare gdal e WhiteboxTools (o grass).

Ma ne riparliamo con più calma

[aborruso]

Python mi sembra la via più comoda.

Genera output.zip in 0.4 secondi.

import rasterio
import geopandas as gpd
from rasterstats import zonal_stats
import os

input_raster = 'testOneR.tif'
input_shapefile = 'zone.shp'
output_shapefile = 'output_shapefile.shp'

# Carica lo shapefile
gdf = gpd.read_file(input_shapefile)

# Loop per i valori dei pixel da 1 a 4
for value in range(1, 5):
    temp_raster_file = f'temp_raster_{value}.tif'

    # Crea un raster temporaneo per ogni valore
    with rasterio.open(input_raster) as src:
        data = src.read(1)
        data[data != value] = src.nodata
        profile = src.profile
        with rasterio.open(temp_raster_file, 'w', **profile) as dst:
            dst.write(data, 1)

    # Calcola le statistiche zonali
    stats = zonal_stats(input_shapefile, temp_raster_file, stats="count")
    gdf[f'count_{value}'] = [stat['count'] if stat else None for stat in stats]

    # Rimuovi il file temporaneo
    os.remove(temp_raster_file)

# Salva il nuovo shapefile con le colonne aggiunte
gdf.to_file(output_shapefile)

[pigreco]

Per chi usasse solo QGIS:

tra gli strumenti di processing c'è l'algoritmo Istogramma zonale:

che risolve brillantemente il quesito:

[aborruso]

Non ti resta che scrivere il post

[aborruso]

È molto carino anche così con rio e zonalstats. Il file vettoriale deve essere in formato geojson

rio zonalstats --categorical zone.geojson -r testOneR.tif >out.geojson

https://pythonhosted.org/rasterstats/cli.html?highlight=command#command-line-interface

[pigreco]

Ho aggiunto prefix TEXT e ottengo :

rio zonalstats --categorical zone.geojson prefix Andrea_ -r testOneR.tif >out.geojson

[pigreco]

ricetta fatta e pubblicata

https://tansignari.opendatasicilia.it/ricette/mappe/istogramma_zonale/