Avete mai avuto la necessità di eseguire un’operazione tanto semplice quando indispensabile come il dover estrarre l’area non condivisa fra due geometrie sovrapposte? In sostanza, l’opposto di un’operazione di Intersect?

symmetrical_difference
Differenza simmetrica

Tale operazione può essere eseguita con lo strumento Symmetrical Difference. L’unico problema è che si tratta di un tool che necessita di una licenza di livello Advanced, non sempre a disposizione di tutti gli utenti. Nessun problema! In questo post vi faremo vedere come sfruttare gli oggetti Geometria per rendere possibili questa e molte altre operazioni utilissime, superando il limite della licenza. Il tutto sfruttando Python e il modulo arcpy.Ma cominciamo subito con un esempio. Il seguente codice crea, a partire da uno shapefile poligonale contenente una sola geometria rettangolare, un oggetto Geometria:

import arcpy
geom_1 = arcpy.CopyFeatures_management("geom_1.shp", arcpy.Geometry())

Abbiamo sfruttato lo strumento Copy Features per copiare le sole geometrie contenute nello shapefile geom_1.shp in una variabile che abbiamo denominato geom_1. Il tutto specificando come parametro di output della funzione la chiamata a arcpy.Geometry(), che crea l’oggetto Geometria. Quest’oggetto ha una serie di metodi (intesi in senso pythonico), fra cui uno chiamato symmetricDifference(), che andremo a utilizzare per il nostro scopo.

Andiamo oltre e procediamo con la creazione del secondo oggetto Geometria che sarà sovrapposto in parte al primo precedentemente creato (proprio come nell’immagine a inizio post).

geom_2 = arcpy.CopyFeatures_management("geom_2.shp", arcpy.Geometry())

Nulla di nuovo: abbiamo replicato il codice per la creazione della prima geometria, semplicemente sostituendo il nome della variabile (geom_2 anziché geom_1), e lo shapefile da copiare. Ora abbiamo due oggetti Geometria da poter usare come input per ottenere una terza Geometria che sia l’unione delle due geometrie meno l’intersezione fra di esse.

Un’illustrazione renderà più chiaro il funzionamento del metodo symmetricDifference.

symmetrical_difference
simmetricDifference

Prima di procedere col nostro codice, vanno precisate alcune cose da tenere a mente quando si usano gli oggetti Geometria:

  1. ogni poligono, inteso come entità geometrica individuale, che compone lo shapefile in input, corrisponderà ad un oggetto Geometria di tipo poligonale. Questo, insieme agli altri (eventuali) oggetti, vengono raccolti in una lista di Python. In sostanza, geom_1 e geom_2 sono due liste, ciascuna avente al proprio interno un solo oggetto Geometria (poligonale in questo caso come gli shapefile da cui sono derivati). In pratica se stampassimo (con la funzione print) il valore di una delle due variabili, otterremmo qualcosa di simile a:

[Polygon object at 0x4d34bb0[0x26e2040]]

  1. gli oggetti Geometria rappresentano una posizione nello spazio e la forma ad essa associata, pertanto non contengono attributi.

Chiariti questi punti, finalizziamo ora il codice.

geom_diff = geom_1[0].symmetricDifference(geom_2[0])
arcpy.CopyFeatures_management(geom_diff,"geom_diff.shp")

Fatto! Abbiamo usato il metodo symmetricDifference() per ottenere un nuovo oggetto Geometria che rappresenta l’unione delle due geometrie meno l’intersezione fra di esse (proprio quello che fa lo strumento Symmetrical Difference dell’Analysis Toolbox in ArcMap), e l’abbiamo salvato in uno shapefile geom_diff.shp usando ancora lo strumento Copy Features (in maniera opposta rispetto alle righe 2 e 3 nel codice). Notare che alla riga 4 abbiamo dovuto specificare un indice per la variabile geom_1 e geom_2, dato che in realtà si tratta di due liste contenenti ciascuna un solo elemento, l’oggetto Geometria. Ricordiamo che [0] rappresenta un indice che punta al primo (e in questo caso unico) elemento all’interno di una lista in Python.

Riportiamo di seguito il codice completo:

import arcpy
geom_1 = arcpy.CopyFeatures_management("geom_1.shp", arcpy.Geometry())
geom_2 = arcpy.CopyFeatures_management("geom_2.shp", arcpy.Geometry())
geom_diff = geom_1[0].difference(geom_2[0])
arcpy.CopyFeatures_management(geom_diff,"geom_diff.shp")

Abbiamo visto in questo post come sfruttare gli oggetti Geometria e il loro metodo symmetricDifference() per replicare il funzionamento dello strumento Symmetrical Difference presente in ArcMap, ma utilizzabile solo con una licenza di livello Advanced. Questo è utile in svariate casistiche (clausole di contenimento, individuazione delle interferenze, ecc.). Ovviamente non è il solo metodo degli oggetti Geometria: equals() consente di confrontare due geometrie e restituisce true se sono identiche, ed è utile per verificare che non vi siano duplicati geometrici all’interno di un proprio dataset; generalize() crea una geometria che è la versione generalizzata (semplificata) di quella in input, utile quando si hanno troppi vertici e si vuole “snellire” il dato; l’opposto è il metodo densify(), che invece aggiunge vertici. Per un elenco completo, si faccia riferimento a questa pagina dell’help online di ArcGIS.

Gli oggetti Geometria possono essere utilizzati come input o output al posto di creare nuove Feature Class temporanee, applicarvi dei cursori, e alla fine cancellarle durante una fase intermedia di un processo di geoprocessing. Possono anche essere creati da zero chiamando la classe Geometry(), perfino “al volo”, se necessario, inserendoli come parametri durante una chiamata a una funzione di geoprocessing.

Nei prossimi post vedremo come mantenere le informazioni (gli attributi) degli shapefile di partenza agganciati agli oggetti Geometria tramite l’uso dei cursori.

Leave a Comment