提取具有 NA 值的给定缓冲区附近的像素值和坐标

Question

我想获取随机坐标（pts）的邻域（例如buffer=6米）中的值（像素值）、坐标（x和y）和属性（status），使用提取raster 包中的函数。我尝试在没有 NA 值的 data.frame 中组织结果，@Robert Hijmans 在Extracting pixels values and coordinates in neighborhood of given buffer in R 中解决了这个问题。

但是，如果我在一个栅格之外有一些坐标（并且我为此目的创建了s2 栅格），则脚本不起作用。我尝试删除列表中不完整的元素（NA 值、不同数量的元素/列），但最终结果不匹配。

在我的新方法中：

library(raster)  
r <- raster(ncol=10, nrow=10, crs="+proj=utm +zone=1 +datum=WGS84", xmn=0, xmx=50, ymn=0, ymx=50)
s1 <- stack(lapply(1:4, function(i) setValues(r, runif(ncell(r)))))
r2 <- raster(ncol=10, nrow=10, crs="+proj=utm +zone=1 +datum=WGS84", xmn=0, xmx=100, ymn=0, ymx=100) # Large raster for produce NAs
s2 <- stack(lapply(1:4, function(i) setValues(r2, runif(ncell(2)))))
ras <- list(s1, s2)
pts <- data.frame(pts=sampleRandom(s2, 100, xy=TRUE)[,1:2], status=rep(c("A","B"),5))

# get xy from buffer cells
cell <- extract(r, pts[,1:2], buffer=6, cellnumbers=T)
xy <- xyFromCell(r, do.call(rbind, cell)[,1])
xy<-xy[complete.cases(xy),] # Remove NA coordinates


# lopp for extract pixel values and coordinates
res <- list()
for (i in 1:length(ras)) {
    v <- raster::extract(ras[[i]], pts[,1:2], buffer=6)
    delete.NULLs1  <-  function(x.list){   # delele one single column in a list 
    x.list[unlist(lapply(x.list, function(x) length(unique(x))) != 1)]} 
    delete.NULLs2  <-  function(x.list){   # delele different number of elements in a list
    x.list[unlist(lapply(x.list, length)) >= 5]}
    delete.NULLs3  <-  function(x.list){   # delele null/empty entries in a list
    x.list[unlist(lapply(x.list, length) != 0)]}
    v <- delete.NULLs1(v)
    v <- delete.NULLs2(v)
    v <- delete.NULLs3(v)
    # add point id
    for (j in 1:length(v)) {
        v[[j]] <- cbind(point=j, v[[j]])
    }
    #add layer id and xy
    res[[i]] <- cbind(layer=i, xy, do.call(rbind, v))
}
res <- do.call(rbind, res)

我的输出总是：

Error in cbind(layer = i, xy, do.call(rbind, v)) : 
  number of rows of matrices must match (see arg 3)

delete.NULLs 函数后，我丢失了坐标/栅格列表对应关系。请问有什么想法吗？

Answer 1

这是我可能的处理方法

示例数据

library(raster)  
r <- raster(ncol=10, nrow=10, crs="+proj=utm +zone=1 +datum=WGS84", xmn=0, xmx=50, ymn=0, ymx=50)
s1 <- stack(lapply(1:4, function(i) setValues(r, runif(ncell(r)))))
r2 <- raster(ncol=10, nrow=10, crs="+proj=utm +zone=1 +datum=WGS84", xmn=0, xmx=100, ymn=0, ymx=100) # Large raster for produce NAs
s2 <- stack(lapply(1:4, function(i) setValues(r2, runif(ncell(2)))))
ras <- list(s1, s2)
pts <- data.frame(pts=sampleRandom(s2, 100, xy=TRUE)[,1:2], status=rep(c("A","B"),5))

# get xy from buffer cells
cell <- extract(r, pts[,1:2], buffer=6, cellnumbers=T)
xy <- xyFromCell(r, do.call(rbind, cell)[,1])
xy<-xy[complete.cases(xy),] # Remove NA coordinates

更新算法

res <- list()
for (i in 1:length(ras)) {
    v <- raster::extract(ras[[i]], pts[,1:2], buffer=6)
    # find invalid cases (NA or zero rows), a bit tricky
    k <- sapply(sapply(v, nrow), function(i) ifelse(is.null(i), FALSE, i>0))
    # jump out of loop if there is no data
    if (!any(k)) next
    # remove the elements from the list that have no data
    v <- v[k]
    k <- which(k)
    # add point id
    for (j in 1:length(k)) {
        kj <- k[j]
        v[[j]] <- cbind(point=kj, xy[kj,1], xy[kj,2], v[[j]])
    }
    v <- do.call(rbind, v)
    colnames(v)[2:3] <- c("x", "y")
    #add layer id and xy
    res[[i]] <- cbind(layer=i, v)
}
res <- do.call(rbind, res)