溶解 rasterToPolygons 的输出
Dissolve output of rasterToPolygons
在 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)
但是,我希望每个具有相邻边或角的多边形成为一个较大多边形的一部分,从而减少多边形总数。有什么办法可以做到这一点?
旧答案:
您可以使用参数 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 个(多)多边形。如果你想分开 non-connected 部分,你可以做
pols <- disaggregate(pol)
pols
#class : SpatialPolygonsDataFrame
#features : 80
请注意,对角线相邻的多边形彼此分开,因为它们不能用于 有效的 单个多边形(它将是 self-intersecting)。
这可以通过使用 spdep
包中的 poly2nb
函数来定义每个多边形的邻居,使用下面创建的函数来创建区域分配向量,使用 [= maptools
包中的 13=] 将 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)
在 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)
但是,我希望每个具有相邻边或角的多边形成为一个较大多边形的一部分,从而减少多边形总数。有什么办法可以做到这一点?
旧答案:
您可以使用参数 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 个(多)多边形。如果你想分开 non-connected 部分,你可以做
pols <- disaggregate(pol)
pols
#class : SpatialPolygonsDataFrame
#features : 80
请注意,对角线相邻的多边形彼此分开,因为它们不能用于 有效的 单个多边形(它将是 self-intersecting)。
这可以通过使用 spdep
包中的 poly2nb
函数来定义每个多边形的邻居,使用下面创建的函数来创建区域分配向量,使用 [= maptools
包中的 13=] 将 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)