Объединение небольших полигонов с крупнейшим соседом в R

У меня есть несколько полигонов, представляющих классы растительного покрова, однако мне не нужны такие подробности. Я хотел бы объединить небольшие многоугольники (т. е. <150 м2) с их самым большим соседом. Это будет похоже на «Удалить» в ArcGIS Pro или «Удалить выбранные полигоны» в QGIS. Этот процесс будет повторяться, поэтому я хотел бы написать его на R и сэкономить немного времени.

Единственный способ, которым я мог это сделать, — это добавить столбец, указывающий, какие полигоны были меньше 150 м2, а затем каким-то образом объединить их с ближайшим соседом.

      #Packages
library(terra)
library(sf)
library(dplyr)

#Adding polygons from Terra and converting to sf
v <- vect(system.file("ex/lux.shp", package="terra"))
v <- st_as_sf(v[c(1:12)])


#Adding a column indicating which polygons are under 150m2
mutate(v, area_thresh = case_when(AREA < 150 ~ 1,
                               AREA > 150 ~ 0))

2 ответа

Можешь попробоватьterra::combineGeoms.

Пример данных

      library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))

Решение

      th <- v$AREA < 150
x <- combineGeoms(v[!th], v[th])

Иллюстрация

      v$th <- th
plot(v, "th")
lines(x, lwd=3, col="blue")

Это интересная проблема, так как вам нужно перебирать строки изменяющегося шейп-файла (количество строк будет уменьшаться по мере объединения меньших объектов).

В качестве возможного подхода рассмотрим этот фрагмент кода, основанный на хорошо известном и любимом шейп-файле NC, который поставляется с{sf}.

Обратите внимание, что я использую county_id в качестве уникального идентификатора; объединенный полигон сохранит идентификатор большего округа.

      library(sf)
library(dplyr)

# the one & only NC shapefile... as geometry only
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  select(CNTY_ID)

# limit for merging - in this case median area
limit <- median(st_area(shape))

# initialize iterator
i <- 1

# iterate over rows of shapefile
while (i < nrow(shape)) {
  
  # check if area over or under limit
  if(st_area(shape[i, ]) < limit) {
    
    # find id of largest neighbor
    largest_neighbor <- shape[unlist(st_touches(shape[i, ], shape)), ] %>% 
      mutate(area = st_area(.)) %>% 
      top_n(1, area) %>% 
      pull(CNTY_ID)
    
    # merge offending shape with largest neighbor
    wrk_shape <- st_union(st_geometry(shape[i, ]), 
                          st_geometry(shape[shape$CNTY_ID == largest_neighbor, ])) %>% 
      st_make_valid() %>% 
      st_as_sf() %>% 
      mutate(CNTY_ID = largest_neighbor)
    
    # rename geometry, column to enable binding
    st_geometry(wrk_shape) = attr(shape, "sf_column")
    
    # remove offending shape and unmerged neighbor
    shape <- shape[-i, ]
    shape <- filter(shape, !CNTY_ID == largest_neighbor) 
    
    # add merged shape back in
    shape <- bind_rows(shape, wrk_shape)
    
  }
  
  # don't forget to update iterator :)
  i <- i + 1
    
}

plot(st_geometry(shape))

Другие вопросы по тегам