我需要从栅格中与我的焦点单元格相邻的单元格中提取值，并且我需要能够将这些值与焦点单元格链接起来。 到目前为止，我可以从相邻单元格中提取值，但结果不是有组织的形式。
library(raster)
r <- raster(matrix(runif(100), 10))
cells <- c(34,22,50,10)
cells_ad <- adjacent(r,cells, directions = 8, pairs=FALSE)
extract (r,cells_ad)
我需要一个数据框，其中包含为每个相邻单元格提取的值，每列是一个相邻单元格（理想情况下，我会知道顺序 - 南、西南、西、西北......），每一行焦点细胞。按照我上面的示例，数据框将有 8 列和 4 行。
需要为大量的光栅文件和数百万个要提取的点执行此操作。所以这里应该考虑计算时间。
感谢您的帮助。
【问题讨论】：
这是一个简单的 Queen-case 相邻函数，包括栅格之外的单元格（如 NA）。
adj <- function(x, rr, ngb, global) {
a <- x + ngb
# the below line is not strictly needed, so it should be faster without it
# a[a < 1 | a > ncell(rr)] <- NA
if (!global) {
col <- colFromCell(rr, x)
if (col == 1) {
a[c(1,4,6)] <- NA
} else if (col==nc) {
a[c(3,5,8)] <- NA
}
}
a
}
您的示例数据
library(raster)
r <- raster(matrix(runif(100), 10))
cells <- c(34,22,50,10)
应用函数
nc <- ncol(r)
nbrs <- c(-nc-1, -nc, -nc+1, -1, +1, +nc-1, +nc, +nc+1)
ad <- t(sapply(cells, adj, rr=raster(r), ngb=nbrs, global=FALSE))
colnames(ad) <- c("NW", "N", "NE", "W", "E", "SW", "S", "SE")
提取值
ad[] <- extract(r, as.vector(ad))
【讨论】：
请参阅下面我的（草率）解决方案。我注意到相邻（）函数省略了 NA 单元格，因此我对其进行了修改以将其删除。提取单元格的方向遵循以下顺序："nw","w","sw","ne","e","se","n","s"。
library(raster)
r <- raster(matrix(runif(100), 10))
cells <- c(34,22,50,10)
# modified adjacent function to exclud the na.omit() part, therefore it returns all adjacent cells around the cell
adjacent_mod = function (x, cells, directions = 4, pairs = TRUE, target = NULL,
sorted = FALSE, include = FALSE, id = FALSE)
{
.isGlobalLonLat <- function(x) {
res <- FALSE
tolerance <- 0.1
scale <- xres(x)
if (isTRUE(all.equal(xmin(x), -180, tolerance=tolerance, scale=scale)) &
isTRUE(all.equal(xmax(x), 180, tolerance=tolerance, scale=scale))) {
if (couldBeLonLat(x, warnings=FALSE)) {
res <- TRUE
}
}
res
}
if (is.character(directions)) {
directions <- tolower(directions)
}
x <- raster(x)
r <- res(x)
xy <- xyFromCell(x, cells)
mat <- FALSE
if (is.matrix(directions)) {
stopifnot(length(which(directions == 0)) == 1)
stopifnot(length(which(directions == 1)) > 0)
d <- .adjacentUD(x, cells, directions, include)
directions <- sum(directions == 1, na.rm = TRUE)
mat <- TRUE
}
else if (directions == 4) {
if (include) {
d <- c(xy[, 1], xy[, 1] - r[1], xy[, 1] + r[1],
xy[, 1], xy[, 1], xy[, 2], xy[, 2], xy[, 2],
xy[, 2] + r[2], xy[, 2] - r[2])
}
else {
d <- c(xy[, 1] - r[1], xy[, 1] + r[1], xy[, 1],
xy[, 1], xy[, 2], xy[, 2], xy[, 2] + r[2], xy[,
2] - r[2])
}
}
else if (directions == 8) {
if (include) {
d <- c(xy[, 1], rep(xy[, 1] - r[1], 3), rep(xy[,
1] + r[1], 3), xy[, 1], xy[, 1], xy[, 2], rep(c(xy[,
2] + r[2], xy[, 2], xy[, 2] - r[2]), 2), xy[,
2] + r[2], xy[, 2] - r[2])
}
else {
d <- c(rep(xy[, 1] - r[1], 3), rep(xy[, 1] + r[1],
3), xy[, 1], xy[, 1], rep(c(xy[, 2] + r[2],
xy[, 2], xy[, 2] - r[2]), 2), xy[, 2] + r[2],
xy[, 2] - r[2])
}
}
else if (directions == 16) {
r2 <- r * 2
if (include) {
d <- c(xy[, 1], rep(xy[, 1] - r2[1], 2), rep(xy[,
1] + r2[1], 2), rep(xy[, 1] - r[1], 5), rep(xy[,
1] + r[1], 5), xy[, 1], xy[, 1], xy[, 2], rep(c(xy[,
2] + r[2], xy[, 2] - r[2]), 2), rep(c(xy[, 2] +
r2[2], xy[, 2] + r[2], xy[, 2], xy[, 2] - r[2],
xy[, 2] - r2[2]), 2), xy[, 2] + r[2], xy[, 2] -
r[2])
}
else {
d <- c(rep(xy[, 1] - r2[1], 2), rep(xy[, 1] + r2[1],
2), rep(xy[, 1] - r[1], 5), rep(xy[, 1] + r[1],
5), xy[, 1], xy[, 1], rep(c(xy[, 2] + r[2],
xy[, 2] - r[2]), 2), rep(c(xy[, 2] + r2[2],
xy[, 2] + r[2], xy[, 2], xy[, 2] - r[2], xy[,
2] - r2[2]), 2), xy[, 2] + r[2], xy[, 2] -
r[2])
}
}
else if (directions == "bishop") {
if (include) {
d <- c(xy[, 1], rep(xy[, 1] - r[1], 2), rep(xy[,
1] + r[1], 2), xy[, 2], rep(c(xy[, 2] + r[2],
xy[, 2] - r[2]), 2))
}
else {
d <- c(rep(xy[, 1] - r[1], 2), rep(xy[, 1] + r[1],
2), rep(c(xy[, 2] + r[2], xy[, 2] - r[2]), 2))
}
directions <- 4
}
else {
stop("directions should be one of: 4, 8, 16, \"bishop\", or a matrix")
}
if (include)
directions <- directions + 1
d <- matrix(d, ncol = 2)
if (.isGlobalLonLat(x)) {
d[, 1] <- (d[, 1] + 180)%%360 - 180
}
if (pairs) {
if (mat) {
cell <- rep(cells, each = directions)
}
else {
cell <- rep(cells, directions)
}
if (id) {
if (mat) {
ID <- rep(1:length(cells), each = directions)
}
else {
ID <- rep(1:length(cells), directions)
}
d <- cbind(ID, cell, cellFromXY(x,
d))
attr(d, "na.action") <- NULL
colnames(d) <- c("id", "from", "to")
if (!is.null(target)) {
d <- d[d[, 3] %in% target, ]
}
}
else {
d <- cbind(cell, cellFromXY(x, d))
attr(d, "na.action") <- NULL
colnames(d) <- c("from", "to")
if (!is.null(target)) {
d <- d[d[, 2] %in% target, ]
}
}
if (sorted) {
d <- d[order(d[, 1], d[, 2]), ]
}
}
else {
d <- as.vector(unique(cellFromXY(x, d)))
if (!is.null(target)) {
d <- intersect(d, target)
}
if (sorted) {
d <- sort(d)
}
}
d
}
# extract adjacent cell numbers
cells_ad <- adjacent_mod(r,cells, directions = 8, pairs=T, sorted=T, id=F)
# table them into a data.frame with the orientation
i=1
for (i in 1:length(unique(cells))) {
if (i == 1) {
df = data.frame(matrix(c(cells[i], cells_ad[which(cells_ad[,1]==cells[i]),2]), ncol=9))
} else {
df = rbind(df, data.frame(matrix(c(cells[i], cells_ad[which(cells_ad[,1]==cells[i]),2]), ncol=9)))
}
}
names(df) = c("cell","nw","w","sw","ne","e","se","n","s")
# extract values from raster r
i=1
df_values = df
for (i in 1:dim(df)[1]) {
df_values[i,] = extract (r, as.numeric(df[i,]))
}
【讨论】：