Extracción eficiente de todos los sub-polígonos generados por las características de auto-intersección en un MultiPolygon

intersection sp

Question

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í:

https://it.mathworks.com/matlabcentral/fileexchange/18173-polygon-intersection?requestedDomain=www.mathworks.com

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!

Answer 1

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 con cur
• 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])

Answer 2

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).

Answer 3

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:

1. Cargar todos los polígonos en geopandas GeoDataFrames
2. 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!