“Unire” i poligoni di un layer con grande semplicità: è un lavoro (soltanto?) per mapshaper

Come realizzare l’unione tra geometrie sovrapposte e aggregarne gli attributi in modalità molti-a-uno

Published in

tantotanto

12 min readDec 30, 2019

Quando si lavora con dati spaziali uno degli obiettivi può essere quello di creare nuove geometrie in base a come si sovrappongono fra loro. Queste procedure sono spesso indicate usando il linguaggio degli insiemi: intersezione, unione, differenza, ecc..

Qui sotto un’immagine con alcune delle operazioni “classiche” (questa e il testo introduttivo dalla guida di GeoPandas).

In questo articolo verrà descritta una di queste e in particolare l’unione. È un’operazione realizzabile in diverse modalità e con diversi strumenti, ma diventa abbastanza complessa e/o non realizzabile in modo diretto se l’obiettivo è:

realizzare l’unione tra oggetti di uno stesso layer di input;
produrre in output non soltanto l’unione delle geometrie, ma anche quella degli attributi, sempre in base alla sovrapposizione delle geometrie.

Nell'immagine di sotto un esempio: a sinistra le geometrie di input e a destra l’output. I perimetri delle aree fanno da “lama” di ritaglio e laddove le aree si sovrappongono devono essere trasferiti gli attributi di input (in questo esempio i valori A, B e C); eventuali poligoni sovrapposti e coincidenti sono aggregati in un’unica geometria.

Un’applicazione con cui è molto semplice realizzare quanto descritto è il fantastico mapshaper. È uno strumento scritto in JavaScript, per modificare file Shapefile, GeoJSON, TopoJSON, CSV e altri formati, molto noto per le sue eccellenti doti nella semplificazione di geometrie.
Fa molto molto di più e questa è l’occasione per mostrare un esempio.

A Totò Fiandaca (grazie mille) si deve la redazione del testo che descrive come fare la cosa con QGIS e tramite SQL geografico (in particolare con SpatiaLite).

Unione dei poligoni di un layer con mapshaper

Una volta installato è utilizzabile o a riga di comando o tramite un’interfaccia web (questa la versione pubblica ufficiale).

Per replicare l’operazione descritta a seguire è stato predisposto questo file GeoJSON di input d’esempio, che ha un solo campo (id) ed è composto da tre poligoni (con valori di id pari ad A, B e C).

Questo è il comando tipo da lanciare da shell:

mapshaper -i ./inputLayer.geojson -mosaic calc='output=collect(id).toString()' -o ./output.geojson

Per punti:

con -i si imposta il file di input;
poi il comando da eseguire — -mosaic — che in questo caso è appunto l'unione dei poligoni;
al comando mosaic si associa l'opzione calc, che consente di eseguire calcoli nelle aggregazioni molti-a-uno, utilizzando espressioni JavaScript. In questo caso viene generato un nuovo campo denominato output, in cui vengono aggregati i valori del campo id nei poligoni di intersezione risultanti;
infine si definisce con -o il file di output.

Nell'immagine di sotto è rappresentato il processo. Si veda ad esempio come al poligono centrale, frutto dell’intersezione geometrica di tutti e tre i poligoni di input, venga associato il valore A,B,C.

In termini di geometrie si avranno 7 record di output (corrispondenti alle 7 intersezioni geometriche possibili) con questi attributi:

A,B,C
B,C
A,C
C
A
A,B
B

Tantissimi i campi di applicazione di un processo come questo. Questi tre poligoni potrebbero essere gli areali “colpiti” dalla diffusione di inquinanti da tre punti sorgente, oppure le aree che distano una determinata misura da una fontanella d’acqua potabile, ecc..
L’unione, con intersezione geometrica e aggregazione di attributi, consente di produrre un output in cui è possibile leggere per ogni area il contributo di tutti i poligoni di input.

Altre modalità

Come dicevo in apertura, tutto questo è fattibile in altre modalità, ma non mi risulta ce ne sia una così semplice e diretta in altri ambienti. A seguire qualche esempio.

GeoPandas

GeoPandas è un famoso progetto open source per semplificare il lavoro con i dati geospaziali in Python. È un gran bel modulo!

Consente di eseguire il processo di unione (e altri), ma di base è pensato per essere eseguito tra due layer distinti di input.
A seguito di una issue in merito, aperta sul repository di progetto, è stata fornita questa risposta, in cui si evidenzia come non ci sia un modo diretto per farlo:

[…]there is right now no direct way to do this in geopandas. You could do a “cumulative” intersection manually rather easily in a loop, but that would need some more complexity to keep track of the values.[…]

Per l’occasione è stato creato un bel notebook in cui è mostrato come sia possibile eseguire l’intersezione tra poligoni sovrapposti; da notare la maggiore complessità e sopratutto come comunque non sia gestita l’aggregrazione molti-a-uno degli attributi (nel senso che c’è da scrivere questa parte di codice).

QGIS

A partire dalla versione QGIS 3.2 Bonn è disponibile il geoalgoritmo Union che svolge la parte geometrica di ritaglio del processo descritto sopra, ma non raggruppa le geometrie omologhe e di conseguenza gli attributi.
Una singola area condivisa tra n poligoni sarà ritagliata e produrrà in output n elementi geometrici. Nell'immagine di sotto ad esempio, in corrispondenza dell'area selezionata in giallo, tre record e tre geometrie; complessivamente così verranno quindi prodotti 12 record.

Per raggruppare le geometrie e gli attributi occorre un altro passaggio, oppure creare un modello grafico, come quello che trovate qui; in questo si utilizzano anche lo snapping e la riparazione delle geometrie, che servono in casi molto più complessi di quello attuale (vedi issue correlata).

La query utilizzata nel modello è la seguente:

SELECT group_concat ("id") AS ID, st_union (geometry) AS geometry
FROM input1
GROUP BY geometry

Fatto il raggruppamento si avranno in output ancora una volta 7 record.

SQL

Non esiste una funzione di SQL spaziale (SpatiaLite o PostGIS) che risolva in modo diretto la cosa; occorre analizzare step by step il da farsi:

estrarre i perimetri dai poligoni di input;
estrarre i punti di intersezione dei perimetri delle varie feature;
aggiungere i punti alla geometria dei perimetri;
“splittare” la geometria dei perimetri;
“poligonalizzare” a partire dai perimetri “splittati”;
associare gli attributi alle varie geometrie “poligonalizzate”.

Ad esempio in ambiente SpatiaLite si potrebbe partire con l’importazione del file di esempio:

# importa file inputLayer e crea file spatialite
ogr2ogr -append -t_srs EPSG:4326 -f SQLite ./dbOverlayUnion.sqlite ./inputLayer.geojson -nln "inputLayer" -dsco SPATIALITE=YES

E infine lanciare la seguente query:

-- algoritmo proiezioni, the SpatiaLite way
-- di Salvatore Fiandaca
-- e-mail: pigrecoinfinito@gmail.com
-- è stato di aiuto: http://blog.cleverelephant.ca/2019/07/postgis-overlays.html-- crea geotabella polygonize
CREATE TABLE polygonize AS
SELECT St_Polygonize(t.geom) as geom
FROM
(SELECT id, St_Union(st_boundary(geometry)) as geom
FROM inputLayer) t;
SELECT RecoverGeometryColumn('polygonize','geom',4326,'MULTIPOLYGON','XY');-- crea geotabella dalle componenti elementari della geotabella polygonize
SELECT DropGeoTable('elementi');
SELECT ElementaryGeometries( 'elementi' ,'geom' , 'polygonize' ,'out_pk' , 'out_multi_id', 1 ) as num, 'polygon splitted' as label;-- crea poligoni di output con attributi
SELECT DropGeoTable( "OUTPUT");
CREATE TABLE OUTPUT AS
SELECT Group_Concat(id) as id, e.geom
FROM inputLayer p, elementi e
where st_intersects (ST_PointOnSurface(e.geom), p.geometry) = 1
GROUP BY e.geom;
SELECT RecoverGeometryColumn('OUTPUT','geom',4326,'POLYGON','XY');-- aggiorna statistiche e VACUUM
UPDATE geometry_columns_statistics set last_verified = 0;
SELECT UpdateLayerStatistics('geometry_table_name');
VACUUM;

In output ancora una volta 7 record.

Nota bene: Totò ha poi dettagliato tutto in questo post https://pigrecoinfinito.com/2020/01/02/spatialite-overlays/

SAGA

Sezione a cura di Ludovico Frate (grazie).

In SAGA GIS è possibile utilizzare il comando Polygon Self-Intersection per ottenere 7 geometrie (record) derivanti dall'unione del vettore originale.
I passaggi da eseguire sono:

Importazione del file geojson utilizzando il comando Import Shapes (Tools -> Import/Export->GDAL/OGR->Import Shapes);

Eseguire l’unione eseguendo il comando Polygon Self-Intersection (Tools -> Shapes->Polygons->Polygons Self-Intersection) ed inserire in Polygons il geojson appena caricato, in identifier la colonna id.

In output i soliti 7 record.

R

Sezione a cura di Andrea Zedda (grazie). Qui il notebook.

Carico la libreria sf, una libreria che manipola geometrie spaziali utilizzando lo standard OGC “Simple feature”.
Sf è una libreria relativamente nuova in ambiente R ma data la semplicità d’uso, la versatilità e la piena compatibilità con la sintassi “tidy” è già divenuto lo standard incontrastato.

library(sf)## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()s <- st_read("inputLayer.geojson")
plot(s)

In in sf il concetto di unione è intersezione sembrano essere stati interpretati diversamente da geopandas e altri software. Per ottenere l’input lancio quindi semplicemente la funziona st_intersection che genera le features derivanti dalle sovrapposizioni e produce di default 2 campi: n.overlap che calcola quante volte la geometria si sovrappone con altre geometrie, e origins, una lista che contiene la combinazione degli indici delle features che overlappano.

i <- s %>% st_intersection()
i## Simple feature collection with 7 features and 3 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##     id n.overlaps origins                       geometry
## 1    A          1       1 POLYGON ((1.5 0.1339746, 1....
## 1.1  A          2    1, 2 POLYGON ((0.1339746 1.5, 0....
## 2    B          1       2 POLYGON ((2 2, 2 2, 2 2, 1....
## 1.2  A          2    1, 3 POLYGON ((2 1, 1.99863 0.94...
## 1.3  A          3 1, 2, 3 GEOMETRYCOLLECTION (LINESTR...
## 2.1  B          2    2, 3 POLYGON ((2 2, 1.99863 1.94...
## 3    C          1       3 MULTIPOLYGON (((3 1, 2.9986...

L’output generato sembra corretto, sono state prodotte 7 features tante quante sono le aree che si sovrappongono. Osservando la colonna delle gemotrie notiamo che non tutte sono dei poligoni ma abbiamo un multypolygon e addirittura una geometrycollection.

i %>% 
  ggplot()+
  geom_sf(aes(fill=row.names(i)), show.legend = F)+
  facet_wrap(.~row.names(i), )+
  theme_bw()

Questo è dovuto ad una delle caratteristiche dello standard simple features che per velocizzare alcune operazione elabora delle coordinate arrotondate. È un caso abbastanza raro ma quando si fanno delle operazioni unarie sulle geometrie è meglio settare la precisione.

i <-  st_intersection(st_set_precision(s, 10000))
i## Simple feature collection with 7 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##     id n.overlaps origins                       geometry
## 1    A          1       1 POLYGON ((1.5 0.134, 1.454 ...
## 1.1  A          2    1, 2 POLYGON ((0.134 1.5, 0.1613...
## 2    B          1       2 POLYGON ((1.5 1.866, 1.454 ...
## 1.2  A          2    1, 3 POLYGON ((2 1, 1.9986 0.947...
## 1.3  A          3 1, 2, 3 POLYGON ((1.5 1.866, 1.5446...
## 2.1  B          2    2, 3 POLYGON ((2 2, 1.9986 1.947...
## 3    C          1       3 POLYGON ((3 1, 2.9986 0.947...

Risolto il problema possiamo infine eseguire una seconda operazione che ci consente di portarci nell’output le elaborazioni sugli attributi. In questo caso inserisco nel pipe una funzione di mutate che crea una colonna output contenente una stringa degli attributi concatenati. Per far ciò si è fatto riferimento al campo origins che contiene come detto la posizione delle geometrie nell’input iniziale. Essendo origins una lista si estraggono i valori utilizzando la funzione sapply.

out <-  i %>% mutate(output=sapply(origins,
                     function(x) paste0(as.character(s$id)[x], collapse = ",")
))

Si ottiene quindi l’output desiderato.

out## Simple feature collection with 7 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##   id n.overlaps origins                       geometry output
## 1  A          1       1 POLYGON ((1.5 0.134, 1.454 ...      A
## 2  A          2    1, 2 POLYGON ((0.134 1.5, 0.1613...    A,B
## 3  B          1       2 POLYGON ((1.5 1.866, 1.454 ...      B
## 4  A          2    1, 3 POLYGON ((2 1, 1.9986 0.947...    A,C
## 5  A          3 1, 2, 3 POLYGON ((1.5 1.866, 1.5446...  A,B,C
## 6  B          2    2, 3 POLYGON ((2 2, 1.9986 1.947...    B,C
## 7  C          1       3 POLYGON ((3 1, 2.9986 0.947...      Cout["output"] %>% plot

out %>% 
  ggplot()+
  geom_sf(aes(fill=output), show.legend = F)+
  facet_wrap(.~output)+
  theme_bw()

QGIS (seconda modalità)

Sezione a cura di Ivano Giuliano (grazie).

L’esempio grafico da risolvere mostrato in quest’articolo, è una rappresentazione di quello che viene definito nella teoria degli insiemi un diagramma di Venn, ovvero un diagramma che mostra tutte le possibili relazioni logiche tra un numero finito di insiemi (le tre circonferenze), che vengono a costituire le regioni di interesse. Il metodo si basa proprio su questo concetto, individuazione delle parti di sovrapposizione di regioni e quindi l’insieme di punti in esse ricadenti, nonchè le regioni di provenienza, che andranno a costituire 12 elementi nella prima fase e 7 nella seconda.

Algoritmo Union —> A-B-B,A-C-C,A-C,B-C,B,A (12 elementi)

Algoritmo Aggregate —> {A} - {B} - {B,A} - {C} - {C,A} - {C,B} - {C,B,A} ( 7 insiemi)

Il metodo che si andrà ad esporre qui di seguito, è eseguibile completamente dall’interno di QGIS che prevede due fasi sequenziali con i seguenti geoprocessi.

Geoalgoritmo Union

Riprendendo la definizione : Questo algoritmo controlla le sovrapposizioni tra le features all'interno del livello di Input (nell'esempio : polyover) e crea features separate per quelle parti sovrapposte e non sovrapposte. Ovvero per le parti di aree sovrapposte, verranno create tante features quante sono le effettive parti di sovrapposizione. La tabella degli attributi del livello Union che verrà generato, conterrà tutti i valori di attributo delle features originali per le parti non sovrapposte e valori di attributo per le parti sovrapposte, in totale 12 elementi come indicato qui sopra. L’ algoritmo prevede inoltre la possibilità di utilizzare un livello di sovrapposizione e quindi tutte le sue derivazioni, ma che in questo caso non vengono prese in considerazione.

Geoalgoritmo Aggregate

Riprendendo la definizione : Questo algoritmo a partire da un livello vettoriale in ingresso o tabella, permette l’ aggregazione di features tramite una espressione di raggruppamento Group By Expression . Le features che secondo l’espressione restituiranno lo stesso valore, saranno aggregate. Questo rappresenta la chiave di questa risoluzione al quesito proposto, tramite QGIS. L’algoritmo è molto ben strutturato e permette di effettuare tutta una serie di raggruppamenti mediante modalità parametriche diverse, in questo esempio, viene utilizzata la seguente espressione di raggruppamento:

geom_to_wkt($geometry)

classica funzione utilizzabile nel field calculator che permette in modalità Well-Known-Text (WKT), di rappresentare la geometria . La risultante del geoprocesso, sarà che le geometrie riconosciute, saranno combinate in una geometria multipart per ogni gruppo creato. I relativi attributi di output, verranno cosí calcolati sulla base di ciascuna aggregazione definita.

NOTA BENE: Salvatore Fiandaca ha realizzato (grazie) questo modello grafico QGIS, per questa modalità

Per concludere

L’autore di mapshaper è Matthew Bloch. È uno sviluppatore con doti al di fuori del comune, sia in termini professionali, che in termini di attitudine all'ascolto e allo scambio con gli altri.
Questo comando di unione è stato introdotto da poco (il 18 novembre 2019), in risposta a una richiesta ricevuta nel repository del progetto (da sottolineare che anche precedentemente fosse semplice realizzare la cosa).

È evidente come mapshaper sia (tra quelle di questo articolo) la soluzione più diretta e semplice. Lo è anche perché è stata costruita una funzione specializzata per svolgere questo compito; vista l’utilità dovrebbe/potrebbe essere disponibile anche in altri ambienti.

Non vi resta che installarlo e apprezzarne anche tutte le altre ottime caratteristiche.

Un grosso grazie a Matthew per quello che realizza e mette a disposizione e per il modo in cui lo fa!