Эффективное извлечение всех подполигонов, создаваемых самопересекающимися функциями в MultiPolygon

Ответ 1

Теперь это реализовано в R-пакете sf как результат по умолчанию, когда st_intersection вызывается с единственным аргументом (sf или sfc), см. https://r-spatial.github.io/sf/reference/geos_binary_ops.html для примеров. (Я не уверен, что в поле origins содержатся полезные индексы, в идеале они должны указывать только на индексы только в x, прямо сейчас они вроде себя относятся).

Ответ 2

Ввод

Я немного модифицирую данные макета, чтобы проиллюстрировать способность обрабатывать несколько атрибутов.

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

Основные операции

Затем определите следующие две функции.

• cur: текущий индекс базового многоугольника
• x: индекс многоугольников, который пересекается с cur
• input_polys: простая особенность многоугольников
• keep_columns: вектор имен атрибутов, необходимых для сохранения после геометрического расчета

get_difference_region() получить разницу между базовым многоугольником и другими пересекающимися многоугольниками; get_intersection_region() получить пересечения между пересекающимися многоугольниками.

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

Первый уровень

Разница

Посмотрите эффект разностной функции на данные макета.

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

Пересечение

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

Второй уровень

используйте пересечение первого уровня в качестве входа и повторите тот же процесс.

Разница

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

Пересечение

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

Получите четкие пересечения второго уровня и используйте их как вход третьего уровня. Мы могли бы получить, что результаты пересечения третьего уровня NULL, тогда процесс должен завершиться.

Резюме

Мы помещаем все результаты разности в закрытый список и помещаем все результаты пересечения в открытый список. Тогда имеем:

• Когда открытый список пуст, мы останавливаем процесс
• Результаты закрытого списка

Поэтому мы получаем окончательный код здесь (должны быть объявлены основные две функции):

# 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])

Ответ 3

Не уверен, что это поможет вам, так как это не в R, но я думаю, что есть хороший способ решить эту проблему с помощью Python. Существует библиотека под названием GeoPandas (http://geopandas.org/index.html), которая позволяет вам легко выполнять геоинформацию. В шагах вам нужно будет сделать следующее:

• Загрузка всех полигонов в геоданные GeoDataFrames
• Завершить все GeoDataFrames, выполняющие операцию наложения объединения (http://geopandas.org/set_operations.html)

Точный пример показан в документации.

До операции - 2 многоугольника

После операции - 9 многоугольников

Если есть что-то неясное, не стесняйтесь, дайте мне знать! Надеюсь, это поможет!