Unire i poligoni di un layer
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 versalità e la piena compatibilità con la sintassi “tidy” è già divenuto lo standard incontastrato.
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()Genero i poligoni sovrapposti (come nel notebook di geopandas scelgo questa soluzione per la certezza della riproducibilità del codice)
s <- st_read("inputLayer.geojson")## Reading layer `inputLayer' from data source `/home/andria/Dropbox (kode)/tech_zedda/playground/unary/inputLayer.geojson' using driver `GeoJSON'
## Simple feature collection with 3 features and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 12.62027 ymin: 37.73195 xmax: 12.92255 ymax: 37.95624
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defsplot(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: 12.62027 ymin: 37.73195 xmax: 12.92255 ymax: 37.95624
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## id n.overlaps origins geometry
## 1 A 1 1 POLYGON ((12.76425 37.75672...
## 1.1 A 2 1, 2 POLYGON ((12.69781 37.88527...
## 2 B 1 2 MULTIPOLYGON (((12.69734 37...
## 1.2 A 2 1, 3 POLYGON ((12.80191 37.81656...
## 1.3 A 3 1, 2, 3 POLYGON ((12.77859 37.86215...
## 2.1 B 2 2, 3 POLYGON ((12.8758 37.86619,...
## 3 C 1 3 MULTIPOLYGON (((12.74093 37...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: 12.6203 ymin: 37.732 xmax: 12.9226 ymax: 37.9562
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## id n.overlaps origins geometry
## 1 A 1 1 POLYGON ((12.76425 37.75676...
## 1.1 A 2 1, 2 POLYGON ((12.69779 37.88527...
## 2 B 1 2 POLYGON ((12.6973 37.8824, ...
## 1.2 A 2 1, 3 POLYGON ((12.8019 37.8166, ...
## 1.3 A 3 1, 2, 3 POLYGON ((12.77858 37.86215...
## 2.1 B 2 2, 3 POLYGON ((12.87577 37.86622...
## 3 C 1 3 POLYGON ((12.8578 37.873, 1...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’outpout desiderato.
out## Simple feature collection with 7 features and 4 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 12.6203 ymin: 37.732 xmax: 12.9226 ymax: 37.9562
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## id n.overlaps origins geometry output
## 1 A 1 1 POLYGON ((12.76425 37.75676... A
## 2 A 2 1, 2 POLYGON ((12.69779 37.88527... A,B
## 3 B 1 2 POLYGON ((12.6973 37.8824, ... B
## 4 A 2 1, 3 POLYGON ((12.8019 37.8166, ... A,C
## 5 A 3 1, 2, 3 POLYGON ((12.77858 37.86215... A,B,C
## 6 B 2 2, 3 POLYGON ((12.87577 37.86622... B,C
## 7 C 1 3 POLYGON ((12.8578 37.873, 1... Cout["output"] %>% plot