列出 for 循环 returns 的输出为空

Question

我编写了一个代码来获取覆盖栅格的不同区域（由 shapefile 分隔）的栅格堆栈的交叉表结果。但是，我得到一个空列表。

这是函数：

transitions <- function(bound, themat) {    # bound = shapefile # themat = rasterstack
    result = vector("list", nrow(bound))    # empty result list
    names(result) = bound@data$GEOCODIGO

    for (i in 1:nrow(bound)) {    # this is the number of polygons to iterate through
        single <- bound[i,]    # selects a single polygon
        clip <- mask(crop(themat, single), single)    # crops the raster to the polygon boundary

        result[i] <- crosstab(clip, digits = 0, long = FALSE, useNA = FALSE)
        return(result)
    }
}

我已经测试了for循环外shapefile/bound中第一个对象的步骤；它运作良好。但我仍然不明白为什么我得到一个空列表。有什么想法吗？

Answer 1

示例数据：

p <- shapefile(system.file("external/lux.shp", package="raster"))
b <- brick(raster(p), nl=2)
values(b) = sample(2, 200, replace=TRUE)

固定函数：

transitions <- function(poly, rast) {    
    result = vector("list", nrow(poly))
    for (i in 1:nrow(poly)) {  
        clip <- mask(crop(rast, poly[i,]), poly[i,])   
        result[[i]] <- crosstab(clip, digits = 0, long = FALSE, useNA = FALSE)
    }
    return(result)
}

transitions(p, b)

另一种方法是使用 extract

e <- extract(b, p)

在交叉表中制表：

ee <- lapply(e, function(x) aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum))

要理解最后一行，您需要解压它。

class(e)
#[1] "list"
length(e)
#[1] 12
e[[1]]
#     layer.1 layer.2
#[1,]       1       1
#[2,]       1       2
#[3,]       2       2
#[4,]       2       1
#[5,]       2       1
#[6,]       1       2
#[7,]       2       2

e是一个长度与多边形数量相同的列表（见length(p)）

让我们把第一个元素聚合起来，得到一个 table 的案例和计数。

x <- e[[1]]

aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum)
#  layer.1 layer.2 count
#1       1       1     1
#2       2       1     2
#3       1       2     2
#4       2       2     2

通过 table 的类似方法（区别在于您可以获得零的 Freq 值

as.data.frame(table(x[,1], x[,2]))
#  Var1 Var2 Freq
#1    1    1    1
#2    2    1    2
#3    1    2    2
#4    2    2    2

现在把你喜欢的函数包装成一个lapply

z <- lapply(e, function(x) aggregate(data.frame(count=rep(1, nrow(x))), data.frame(x), FUN=sum))

更进一步，绑定 data.frames 并将标识符添加到 link 数据返回到多边形

y <- do.call(rbind, z,)
y$id <- rep(1:length(z), sapply(z, nrow))

head(y)
#  Var1 Var2 Freq id
#1    1    1    1  1
#2    2    1    2  1
#3    1    2    2  1
#4    2    2    2  1
#5    1    1    1  2
#6    2    1    2  2

列出 for 循环 returns 的输出为空

List with output from for loop returns empty

for-loop

r

crosstab

r-raster