在 R 或 Grass GIS 中填充栅格孔

Question

示例数据

x <- raster(x=matrix(data=1:36, nrow=6), xmn=-1000, xmx=1000, ymn=-100, ymx=900)
x[c(8, 15, 16, 17, 22, 25, 26, 30, 31)] <- NA
plot(x)

问题

我如何（从算法上）区分栅格中的孔洞，即单元格 c(15:17, 22) 与其他非孔洞间隙（即其余空单元格）界定的区域？

这样就可以只对光栅的孔/岛区域进行操作，用自定义值填充孔等。

实际光栅有大约 30000 个孔，因此速度很重要。我对 R 和 Grass GIS 解决方案都感兴趣。非常感谢您的帮助，非常感谢！

Answer 1

多边形化怎么样？对于速度，我不知道它的价值，但你可以：

x[!is.na(values(x))]<-1
plot(x)
x[is.na(values(x))]<-0

hole <- rasterToPolygons(x, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=T)

现在您必须将单部分多边形变为多部分：

hole2 <- SpatialPolygons(lapply(hole@polygons[[1]]@Polygons, function(xx) Polygons(list(xx),ID=round(runif(1,1,100000000)))))

plot(hole2, add=T)

现在你找到 "real" 个孔，这些孔没有接触到边界

around <- as(extent(x), "SpatialLines")
touch_border <- gIntersects(around, hole2, byid=T) 
extract(x, hole2[!touch_border,],cellnumbers=T)

这为您提供了每个孔的单元格编号。它找到了单元格“8”，你不是说它是一个洞，所以我不确定它是否正确，但它一定非常接近！

如果在 R 中速度慢，请在 GRASS 中执行相同的算法（主要是 rasterToPolygons 调用）

Answer 2

这是一个没有多边形化的解决方案：（它不优雅，但它有效）。但是您必须将 holes/island 重新分类为值（即 999），并将所有其他非岛屿重新分类为 NA。像这样：

x <- raster(x=matrix(rep(NA,36), nrow=6), xmn=-1000, xmx=1000, ymn=-100, ymx=900)
x[c(8, 15, 16, 17, 22, 25, 26, 30, 31)] <- 999

plot(x)

现在我使用 clump() 函数来检查是否有岛屿，该函数最酷的地方在于，它还 returns 这些岛屿的 ID：

#Get Islands with IDs
cl <- clump(x,directions=8)
plot(cl)

然后我根据岛屿的频率创建一个数据框（这只是为了获取每个岛屿的 ID）

freqCl <- as.data.frame(freq(cl))

#remove the (row) which corresponds to the NA values (this is important for the last step)
freqCl <- freqCl[-which(is.na(freqCl$value)),]

检查是否有任何岛屿触及边界：

#Check if the island touches any border and therefore isn't a "real island" (first and last column or row)

noIslandID <- c()
#First row
if(any(rownames(freqCl) %in% cl[1,])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[1,]]
  noIslandID  <- append(noIslandID, eliminate)
}
#Last row
if(any(rownames(freqCl) %in% cl[nrow(cl),])){
  eliminate  <- rownames(freqCl)[rownames(freqCl) %in% cl[nrow(cl),]]
  noIslandID <- append(noIslandID, eliminate)
}
#First col
if(any(rownames(freqCl) %in% cl[,1])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[,1]]
  noIslandID  <- append(noIslandID, eliminate)
}
#Last col
if(any(rownames(freqCl) %in% cl[,ncol(cl)])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[,ncol(cl)]]
  noIslandID  <- append(noIslandID, eliminate)
}

消除那些接触边界的岛屿：

noIslandID <- unique(noIslandID)
IslandID   <- setdiff(rownames(freqCl), noIslandID)

在最后一步中，为初始栅格中的每个 "real island" 分配一个 1：

for(i in 1:length(IslandID)) {
  x[cl[]==as.numeric(IslandID[i])] <- 1
}

plot(x)

Answer 3

在这种情况下，如果一个单元格在四个主要方向中的每一个方向上都至少有一个非缺失值作为边界，则该单元格位于一个洞中。您可以使用下面的代码对此进行测试，该代码还将缺失值替换为 99.

x2 <- x                        ## Create a copy of the raster
notNA <- !is.na(as.matrix(x))  ## Matrix identifying cells that are not NA

nr <- nrow(x) ## Number of rows
nc <- ncol(x) ## Number of columns

for(i in 2:(nr-1)) {       ## Loop over rows, ignoring first and last
    for(j in 2:(nc-1)) {   ## Loop over columns, ignoring first and last
        if(notNA[i,j]) ## Skip cell if not NA
            next
        ## For NA cells, determine if non-NA cell exists in each cardinal direction
        any.north <- any(notNA[1:(i-1),j])
        any.south <- any(notNA[(i+1):nr,j])
        any.west <- any(notNA[i,1:(j-1)])
        any.east <- any(notNA[i,(j+1):nc])
        if(any.north & any.south & any.west & any.east)
            x2[i,j] <- 99 ## If cell is in hole, repalce with new value
    }
}

这是新栅格：

New raster with holes filled by the value 99

在 R 或 Grass GIS 中填充栅格孔

Fill raster holes in R or Grass GIS

r

r-raster

sp

grass

示例数据

问题