r - 提取R中给定缓冲区附近的像素值和坐标

Question

我想使用raster包x中的提取函数获取随机坐标 ( 我尝试组织结果但没有成功。ystatusbuffer=6ptsdata.frame

library(raster)  

#create some GeoTIFF rasters
r <- raster(ncol=10, nrow=10)
s <- stack(lapply(1:8, function(i) setValues(r, runif(ncell(r)))))
f1 <- file.path(tempdir(), "sl1.tif")
f2 <- file.path(tempdir(), "sl2.tif")
writeRaster(s[[1:4]], f1, overwrite=TRUE)
writeRaster(s[[5:8]], f2, overwrite=TRUE)

# 10 random points in the rasters
set.seed(5)
pts <- sampleRandom(s[[1]], 10, xy=TRUE)[,1:2]
status<-c(rep(c("A","B"),5))
pts<-as.data.frame(cbind(pts,status))
i<-c(1,2)
pts[ , i]<-apply(pts[ , i], 2,            
                    function(x) as.numeric(as.character(x)))

#read all rasters 
f <- c(f1, f2)
ras <- lapply(f, brick)


# extract raster values in 10 random coordinates and 6 meters around and organize the results 
RES<-NULL
for(i in 1:length(ras)){
value <- raster::extract(ras[[i]],pts[,1:2], buffer=6)
RES<-rbind(RES,cbind(pts,coordinates(value),value)) #create a data frame of the results
}


RES

Error in data.frame(..., check.names = FALSE) : 
  arguments imply differing number of rows: 10, 4

我当然有不同的行数！！我想创建一个带有随机坐标（xy in pts）的最终数据帧输出，邻域点的 xy（x2 和 y2 像素坐标在 6m 缓冲区周围），状态（状态的重复pts，我认为邻域具有相同的状态pts父坐标）和每一层的值，如：

      x   y  x2   y2  status      sl1.1      sl1.2      sl1.3      sl1.4 ...
1  -162  45  -165  48   A      0.47991386 0.04220410 0.79925156 0.04536868 0.47991386 ...
...

标签： rrasterr-raster

Answer 1

简化的示例数据（不要在示例中写入磁盘！）

library(raster)  
set.seed(5)
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)))))
s2 <- stack(lapply(1:4, function(i) setValues(r, runif(ncell(r)))))
ras <- list(s1, s2)
pts <- data.frame(pts=sampleRandom(s1, 10, 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])

res <- list()
for (i in 1:length(ras)) {
    v <- raster::extract(ras[[i]], pts[,1:2], buffer=6)
    # 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)

添加原始点的坐标

pts$point <- 1:nrow(pts)
res <- merge(res, pts)
head(res)
#  point layer    x    y    layer.1    layer.2    layer.3    layer.4 pts.x pts.y status
#1     1     1  7.5 37.5 0.72070097 0.98188917 0.44275430 0.77354202   7.5  32.5      A
#2     1     1  2.5 32.5 0.44473056 0.36640641 0.78783480 0.25482562   7.5  32.5      A
#3     1     1  7.5 32.5 0.05936247 0.17737527 0.08365037 0.02629751   7.5  32.5      A
#4     1     1 12.5 32.5 0.27514918 0.01776222 0.05997353 0.48397076   7.5  32.5      A
#5     1     1  7.5 27.5 0.23519875 0.24867338 0.20373370 0.23957656   7.5  32.5      A
#6     1     2  2.5 32.5 0.33440265 0.98600510 0.94576856 0.85867224   7.5  32.5      A