【问题标题】:Dissolve output of rasterToPolygons融合 rasterToPolygons 的输出
【发布时间】:2018-04-10 00:34:40
【问题描述】:

当在raster 包中使用rasterToPolygons 时,每个满足公式条件的单元格都会变成自己的多边形:

library(raster)
r <- raster(nrow=18, ncol=36)
r[] <- runif(ncell(r)) * 10
r[r>8] <- NA
pol <- rasterToPolygons(r, fun=function(x){x>6})
plot(pol)

但是,我希望每个具有相邻边或角的多边形成为一个较大多边形的一部分,从而减少总多边形的数量。有没有办法做到这一点?

【问题讨论】:

    标签: geospatial spatial r-raster rgeo


    【解决方案1】:

    旧答案:

    您可以使用参数dissolve=TRUE

    library(raster)
    r <- raster(nrow=18, ncol=36)
    r[] <- sample(2, ncell(r), replace=TRUE)
    pol <- rasterToPolygons(r, dissolve=TRUE)
    plot(pol)
    

    新答案

    如果你不关心值,你可以这样做

    您的示例数据

    library(raster)
    r <- raster(nrow=18, ncol=36)
    r[] <- runif(ncell(r)) * 10
    r[r>8] <- NA
    

    将您想要的所有值单元格设置为一个值,将所有其他值设置为NA

    x <- reclassify(r, rbind(c(-Inf, 6, NA), c(6, Inf, 1)))
    
    pol <- rasterToPolygons(x, dissolve=TRUE)
    

    请注意,pol 现在只有 1 个(多)多边形。如果你想分离非连接部分,你可以这样做

    pols <- disaggregate(pol)
    pols
    #class       : SpatialPolygonsDataFrame 
    #features    : 80 
    

    请注意,对角相邻的多边形彼此分开,因为它们不能用于有效单个多边形(它将是自相交的)。

    【讨论】:

    • 溶解选项仅在此处有效,因为您已强制栅格中的每个单元格具有相同的值。我正在寻找一种将所有相邻输出多边形组合在一起的方法,而不管它们继承的栅格值如何。
    • 我现在更新了答案,因为我认为我更好地理解了你的问题
    【解决方案2】:

    这可以通过使用spdep 包中的poly2nb 函数来定义每个多边形的邻居,使用下面创建的函数来创建区域分配向量,使用来自maptoolsspCbind 来完成包将regions绑定到pol,然后最终使用maptools中的unionSpatialPolygons函数在regions上溶解。创建函数的基本结构是if 至少一个多边形的邻居已分配到组then 分配多边形和邻居到该组else 分配多边形和邻居到新组。

    library(raster)
    library(spdep)
    library(maptools)
    
    r <- raster(nrow=18, ncol=36)
    r[] <- runif(ncell(r)) * 10
    r[r>8] <- NA
    pol <- rasterToPolygons(r, fun=function(x){x>6}, dissolve = T)
    plot(pol)
    
    nb <- poly2nb(pol)
    
    create_regions <- function(data) {
      group <- rep(NA, length(data))
      group_val <- 0
      while(NA %in% group) {
        index <- min(which(is.na(group)))
        nb <- unlist(data[index])
        nb_value <- group[nb]
        is_na <- is.na(nb_value)
        if(sum(!is_na) != 0){
          prev_group <- nb_value[!is_na][1]
          group[index] <- prev_group
          group[nb[is_na]] <- prev_group
        } else {
          group_val <- group_val + 1
          group[index] <- group_val
          group[nb] <- group_val
        }
      }
      group
    }
    
    region <- create_regions(nb)
    pol_rgn <- spCbind(pol, region)
    pol2 <- unionSpatialPolygons(pol_rgn, region)
    plot(pol2)
    

    【讨论】:

    • 这行得通,但是对于手头的任务来说它是不必要的复杂,并且可能不能很好地扩展到大型栅格
    猜你喜欢
    • 1970-01-01
    • 2014-06-10
    • 1970-01-01
    • 2019-03-08
    • 2017-06-16
    • 1970-01-01
    • 1970-01-01
    • 2012-07-04
    • 2014-04-24
    相关资源
    最近更新 更多