Extracción eficiente de todos los sub-polígonos generados por las características de auto-intersección en un MultiPolygon
intersection sp (3)
Entrada
Modifico un poco los datos de la maqueta para ilustrar la capacidad de tratar con múltiples atributos.
library(tibble)
library(dplyr)
library(sf)
ncircles <- 9
rmax <- 120
x_limits <- c(-70,70)
y_limits <- c(-30,30)
set.seed(100)
xy <- data.frame(
id = paste0("id_", 1:ncircles),
val = paste0("val_", 1:ncircles),
x = runif(ncircles, min(x_limits), max(x_limits)),
y = runif(ncircles, min(y_limits), max(y_limits)),
stringsAsFactors = FALSE) %>%
as_tibble()
polys <- st_as_sf(xy, coords = c(3,4)) %>%
st_buffer(runif(ncircles, min = 1, max = 20))
plot(polys[1])
Operación básica
Luego define las siguientes dos funciones.
-
cur: el índice actual del polígono base -
x: el índice de polígonos, que se interseca concur -
input_polys: la característica simple de los polígonos -
keep_columns: el vector de nombres de atributos necesarios para mantener después del cálculo geométrico
get_difference_region() obtiene la diferencia entre el polígono base y otros polígonos intersectados; get_intersection_region() obtiene las intersecciones entre los polígonos intersectados.
library(stringr)
get_difference_region <- function(cur, x, input_polys, keep_columns=c("id")){
x <- x[!x==cur] # remove self
len <- length(x)
input_poly_sfc <- st_geometry(input_polys)
input_poly_attr <- as.data.frame(as.data.frame(input_polys)[, keep_columns])
# base poly
res_poly <- input_poly_sfc[[cur]]
res_attr <- input_poly_attr[cur, ]
# substract the intersection parts from base poly
if(len > 0){
for(i in 1:len){
res_poly <- st_difference(res_poly, input_poly_sfc[[x[i]]])
}
}
return(cbind(res_attr, data.frame(geom=st_as_text(res_poly))))
}
get_intersection_region <- function(cur, x, input_polys, keep_columns=c("id"), sep="&"){
x <- x[!x<=cur] # remove self and remove duplicated obj
len <- length(x)
input_poly_sfc <- st_geometry(input_polys)
input_poly_attr <- as.data.frame(as.data.frame(input_polys)[, keep_columns])
res_df <- data.frame()
if(len > 0){
for(i in 1:len){
res_poly <- st_intersection(input_poly_sfc[[cur]], input_poly_sfc[[x[i]]])
res_attr <- list()
for(j in 1:length(keep_columns)){
pred_attr <- str_split(input_poly_attr[cur, j], sep, simplify = TRUE)
next_attr <- str_split(input_poly_attr[x[i], j], sep, simplify = TRUE)
res_attr[[j]] <- paste(sort(unique(c(pred_attr, next_attr))), collapse=sep)
}
res_attr <- as.data.frame(res_attr)
colnames(res_attr) <- keep_columns
res_df <- rbind(res_df, cbind(res_attr, data.frame(geom=st_as_text(res_poly))))
}
}
return(res_df)
}
Primer nivel
Diferencia
Veamos el efecto de la función de diferencia en los datos de la maqueta.
flag <- st_intersects(polys, polys)
first_diff <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_difference_region(i, flag[[i]], polys, keep_column = c("id", "val"))
first_diff <- rbind(first_diff, cur_df)
}
first_diff_sf <- st_as_sf(first_diff, wkt="geom")
first_diff_sf
plot(first_diff_sf[1])
Intersección
first_inter <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_intersection_region(i, flag[[i]], polys, keep_column=c("id", "val"))
first_inter <- rbind(first_inter, cur_df)
}
first_inter <- first_inter[row.names(first_inter %>% select(-geom) %>% distinct()),]
first_inter_sf <- st_as_sf(first_inter, wkt="geom")
first_inter_sf
plot(first_inter_sf[1])
Segundo nivel
Use la intersección del primer nivel como entrada y repita el mismo proceso.
Diferencia
flag <- st_intersects(first_inter_sf, first_inter_sf)
# Second level difference region
second_diff <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_difference_region(i, flag[[i]], first_inter_sf, keep_column = c("id", "val"))
second_diff <- rbind(second_diff, cur_df)
}
second_diff_sf <- st_as_sf(second_diff, wkt="geom")
second_diff_sf
plot(second_diff_sf[1])
Intersección
second_inter <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_intersection_region(i, flag[[i]], first_inter_sf, keep_column=c("id", "val"))
second_inter <- rbind(second_inter, cur_df)
}
second_inter <- second_inter[row.names(second_inter %>% select(-geom) %>% distinct()),] # remove duplicated shape
second_inter_sf <- st_as_sf(second_inter, wkt="geom")
second_inter_sf
plot(second_inter_sf[1])
Obtenga las distintas intersecciones del segundo nivel y utilícelas como entrada del tercer nivel. Podríamos obtener que los resultados de la intersección del tercer nivel sean NULL , entonces el proceso debería terminar.
Resumen
Ponemos todos los resultados de la diferencia en una lista cerrada, y ponemos todos los resultados de la intersección en una lista abierta. Entonces nosotros tenemos:
- Cuando la lista abierta está vacía, detenemos el proceso.
- Los resultados están cerca de la lista.
Por lo tanto, obtenemos el código final aquí (las dos funciones básicas deben ser declaradas):
# init
close_df <- data.frame()
open_sf <- polys
# main loop
while(!is.null(open_sf)) {
flag <- st_intersects(open_sf, open_sf)
for(i in 1:length(flag)) {
cur_df <- get_difference_region(i, flag[[i]], open_sf, keep_column = c("id", "val"))
close_df <- rbind(close_df, cur_df)
}
cur_open <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_intersection_region(i, flag[[i]], open_sf, keep_column = c("id", "val"))
cur_open <- rbind(cur_open, cur_df)
}
if(nrow(cur_open) != 0) {
cur_open <- cur_open[row.names(cur_open %>% select(-geom) %>% distinct()),]
open_sf <- st_as_sf(cur_open, wkt="geom")
}
else{
open_sf <- NULL
}
}
close_sf <- st_as_sf(close_df, wkt="geom")
close_sf
plot(close_sf[1])
Partiendo de un shapefile que contiene un número bastante grande (alrededor de 20000) de polígonos potencialmente superpuestos parcialmente, necesitaría extraer todos los subpolígonos originados por la intersección de sus diferentes "límites".
En la práctica, a partir de algunos datos de maqueta:
library(tibble)
library(dplyr)
library(sf)
ncircles <- 9
rmax <- 120
x_limits <- c(-70,70)
y_limits <- c(-30,30)
set.seed(100)
xy <- data.frame(
id = paste0("id_", 1:ncircles),
x = runif(ncircles, min(x_limits), max(x_limits)),
y = runif(ncircles, min(y_limits), max(y_limits))) %>%
as_tibble()
polys <- st_as_sf(xy, coords = c(2,3)) %>%
st_buffer(runif(ncircles, min = 1, max = 20))
plot(polys[1])
Necesitaría derivar un multipolígono de sf o sp contenga TODOS y SOLAMENTE los polígonos generados por las intersecciones, algo como:
(tenga en cuenta que los colores están ahí solo para ejemplificar el resultado esperado, en el que cada área de "color diferente" es un polígono separado que no se superpone a ningún otro polígono)
Sé que podría abrirme camino analizando un polígono a la vez, identificando y guardando todas sus intersecciones y luego "borrando" esas áreas del multipolígono completo y proceder en un ciclo, pero eso es bastante lento.
Creo que debería haber una solución más eficiente para esto, pero no puedo resolverlo, ¡así que cualquier ayuda sería apreciada! (Ambas soluciones basadas en sf y sp son bienvenidas)
ACTUALIZACIÓN :
Al final, me di cuenta de que incluso al ir "un polígono a la vez", la tarea no es nada simple. Estoy realmente luchando en este problema aparentemente "fácil"! ¿Alguna pista? ¡Incluso se agradecería una solución lenta o sugerencias para iniciar un camino adecuado!
ACTUALIZACIÓN 2 :
Quizás esto aclare las cosas: la funcionalidad deseada sería similar a la que se describe aquí:
ACTUALIZACIÓN 3 :
Otorgué la recompensa a @ shuiping-chen (¡gracias!), Cuya respuesta resolvió correctamente el problema en el conjunto de datos de ejemplo proporcionado. Sin embargo, el "método" tiene que ser generalizado a las situaciones donde son posibles las intersecciones "cuádruple" o "n-uple". ¡Intentaré trabajar en eso en los próximos días y publicar una solución más general si me las arreglo!
Esto ahora se ha implementado en R package sf como el resultado predeterminado cuando se llama st_intersection con un solo argumento (sf o sfc), vea https://r-spatial.github.io/sf/reference/geos_binary_ops.html para los ejemplos . (No estoy seguro de que el campo de origins contenga índices útiles; idealmente, deberían apuntar a índices solo en x , en este momento son una especie de autorreferencia).
No estoy seguro de si te ayuda ya que no está en R, pero creo que hay una buena manera de resolver este problema usando Python. Hay una biblioteca llamada GeoPandas ( http://geopandas.org/index.html ) que le permite realizar operaciones de geo fácilmente. En pasos lo que deberías hacer es lo siguiente:
- Cargar todos los polígonos en geopandas GeoDataFrames
- Bucle todos los GeoDataFrames ejecutando una operación de superposición de unión ( http://geopandas.org/set_operations.html )
El ejemplo exacto se muestra en la documentación.
Antes de la operación - 2 polígonos
Después de la operación - 9 polígonos
Si hay algo que no está claro, no dude en hacérmelo saber! ¡Espero eso ayude!