计算多个光栅中像素覆盖的面积

Question

我有一个96个分类栅格的列表（每个栅格都有一个相关变量，即rcp、period和month），其值范围为1-16，我想计算每个栅格中每个类别所覆盖的面积，如果栅格中不存在该类别，则返回NA。

这是我现在创建的函数

a <- function(x){
rs <- x #categorical raster
b <- getValues(area(x, weights=FALSE))
b <- aggregate(b, by=list(getValues(x)), sum, na.rm=T)
b <- as.data.frame(b)
names(b) <- c("CLASS","AREA")
return(b)
}

问题是，它返回的数据帧仅包含现有光栅值，而不包含缺少的值。见下文：

  CLASS       AREA
1     1 145084.052
2     5  39425.336
3     6  37912.591
4    10  10089.541
5    11   3150.571
6    15   1451.912
7    16   4289.296

如何返回包含所有类别（1-16）的数据帧？并将所有输出合并为一个？列名应为rcp、period和month。

以下是我目前的代码：

mthLs <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
rcpLs <- c("rcp26", "rcp45", "rcp60", "rcp85")
periodLs <- c("2020_2049", "2040_2069")

for(rcp in rcpLs){
  rcpDir <- paste0(iDir, "/", rcp)
     for(period in periodLs){
         for (mth in 1:12) {
          rs <- raster(paste0(iDir, "/", rcp, "/", period, "/cng_", mth, ".tif", sep=""))
          b <- a(rs) #using function above
         }
       im <- cbind(RCP=rep(rcp,times=nrow(b)), PERIOD=rep(period,times=nrow(b)), MONTH=rep(mth,times=nrow(b)), b)

      } 
  write.csv(im, paste(oDir, "/output.csv", sep=""), quote=T, row.names=F)    
}

Answer 1

示例数据

library(raster)
r <- raster(nc=5, nr=5)
set.seed(20171214)
r1 <- setValues(r, sample(16, ncell(r), replace=TRUE))
r2 <- setValues(r, sample(16, ncell(r), replace=TRUE))
r3 <- setValues(r, sample(16, ncell(r), replace=TRUE))

如果可以制作光栅堆栈（具有相同范围和分辨率的光栅层），则可以执行以下操作：

s <- stack(r1,r2,r3)
x <- freq(s)

m <- data.frame(value=1:16)
for (i in 1:length(x)) {
    y <- x[[i]]
    colnames(y)[2] <- paste0("v", i)
    m <- merge(m, y, all.x=TRUE)
}

m
#   value v1 v2 v3
#1      1  3  2  1
#2      2 NA  3 NA
#3      3  2  2  1
#4      4  1  3  2
#5      5  1 NA  2
#6      6  1  3  3
#7      7 NA  3  4
# ... 

或者，在没有RasterStack的情况下，创建一个RasterLayer对象列表，例如

# z <- lapply(filenames, raster)

但是对于这个例子，我们可以更手动地执行：

z <- list(r1, r2, r3)

m <- data.frame(value=1:16)
for (i in 1:length(z)) {
    y <- freq(z[[i]])  
    colnames(y)[2] <- paste0("v", i)
    m <- merge(m, y, all.x=TRUE)
}

或者——也许是最简单的——使用文件名为ff的向量

m <- data.frame(value=1:16)
for (i in 1:length(ff)) {
    y <- freq(raster(ff[i]))  
    colnames(y)[2] <- paste0("v", i)
    m <- merge(m, y, all.x=TRUE)
}

Answer 2

这个可能会有帮助

txt <- 
"CLASS,REA
1,145084.052
5,39425.336
6,37912.591
10,10089.541
11,3150.571
15,1451.912
16,4289.296"

b <- read.csv(textConnection(txt),sep =  ",")

classes <- 1:16

merge(data.frame(CLASS = classes), b, by='CLASS', all.x=T)

结果：

   CLASS        REA
1      1 145084.052
2      2         NA
3      3         NA
4      4         NA
5      5  39425.336
6      6  37912.591
7      7         NA
8      8         NA
9      9         NA
10    10  10089.541
11    11   3150.571
12    12         NA
13    13         NA
14    14         NA
15    15   1451.912
16    16   4289.296