问题描述
我想获取邻域中的值(像素值)、坐标(x
和 y
)和属性(status
)(对于示例在随机坐标 (pts
) 的 buffer=6
米中,使用光栅包中的提取功能.我尝试在 data.frame
中组织结果,但没有成功.
I want to get values (pixels values), coordinates (x
and y
) and attribute (status
) in the neighborhood (for example in a buffer=6
meters) of random coordinates (pts
), using extract function in raster package. I try to organize the results in data.frame
without success.
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
我当然有不同的行数!!我想用随机坐标(pts
中的 xy)、邻域点的 xy(6m 缓冲区周围的 x2 和 y2 像素坐标)、状态( 的重复)创建最终数据帧输出pts
状态,我认为邻域具有与 pts
父坐标相同的状态)并且每个层的值如下:
I have a different number of rows of course!! I've like to create a final data frame output with random coordinates (xy in pts
), xy of neighborhood points (x2 and y2 pixels coordinates around of 6m buffer), status (repetition of pts
status, I consider that neighborhood has the same status of pts
father coordinate) and each layers values like:
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 ...
...
推荐答案
简化的示例数据(不要在示例中写入磁盘!)
Simplified example data (do not write to disk in your examples!)
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
这篇关于提取R中给定缓冲区附近的像素值和坐标的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!