Как переназначить маленькие полигоны ближайшему соседу в R (похоже на функцию gdal_sieve.py)?

Question

Я создал файл пространственного многоугольника из набора данных с сеткой, используя функцию grid2polygons в R. Теперь я хотел бы взять этот файл многоугольника и объединить все многоугольники, содержащие менее 6 пикселей, с ближайшим наибольшим соседом. В ссылке на изображение вы можете увидеть пример маленьких многоугольников, обведенных синим цветом. Я хотел бы, чтобы те слились с большим окружающим многоугольником.

Изображение многоугольника

В Gdal есть функция: gdal_sieve.py, которая может сделать именно это ( http://www.gdal.org/gdal_sieve.html), если вы хотите получить представление о том, что именно я хотел бы сделать. Кто-нибудь знает способ сделать это в R?

Answer 1

... довольно поздно, но это все еще может быть кому-то полезно;)

Мне потребовалось некоторое время, чтобы изучить и протестировать... и еще есть возможности для совершенствования. Тем не менее, это решило мою проблему!

      sieve_slivers = function(x, area.thr=NULL, neigh.case="rook" ){

  if( is.null(area.thr) | area.thr <= 0 ){ stop("Area threshold parameter is wrong or missing") }

  if ( !is(x, "sf") ){ x = st_as_sf(x) }

  if( tolower(neigh.case) == "rook"){ neigh.case = "F***1****" } else { neigh.case = "F***T****" }

  x = x %>%
    mutate( ID = 1:n(), # assign ID first
            area_mq = as.numeric(st_area(.)) )

  newID = do.call(c,
                  lapply(which(x$area_mq <= area.thr), function(g){
                    tmp = st_relate(x[g,], x, pattern = neigh.case )# check: https://github.com/r-spatial/sf/issues/234
                    tmp = x[tmp[[1]],] %>%
                      st_drop_geometry()

                    if( nrow(tmp)>0 ){
                      out = tmp[which.max(tmp$area_mq),"ID"] # in case of 2 or more polygons with the same area, it takes the first
                    } else {
                      out = NA_real_ # for isolated polygons #x$ID[g]
                    }
                    return(out)
                  }) )

  x$ID[which(x$area_mq <= area.thr)] = newID

  # final dissolve
  x = x %>%
    filter( !is.na(ID) ) %>% # remove tiny isolated polygons that could not be merged to anything
    group_by( ID ) %>%
    summarise( area_mq = sum(area_mq, na.rm = TRUE)) %>%
    mutate( ID = NULL)

  return(x)
}