对于R中的土地覆盖栅格，是否有任何有效的方法可以通过像素计数方法获得面积估计？

Question

我想通过使用像素计数方法得到每个多边形的面积估计。原始土地覆盖图是光栅数据，为每个像素指定了一个独特的类别（土地覆盖图图例）。然而，使用像素计数进行面积估计对我来说并不直观。也许，我能做的第一件事就是提取每个多边形的所有像素，这些像素表示每个多边形内土地覆盖类别分布的信息。之后，我需要使用像素计数方法来获得面积估计，并聚合每个多边形的城市、农业区域的所有土地/土壤覆盖率。对我来说，使用像素计数来获得面积估计并不简单，也不知道如何在R中实现这一点。

请注意，原始土地覆盖图相当大（约92 mb），很难为土地覆盖光栅生成可复制的示例，请原谅我的不方便。以下是获取土地覆盖物光栅的脚本：

library(raster)
library(R.utils)

url = "https://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/LUISA/PrimaryOutput/Europe/REF-2014/JRC_LUISA_Input_Corine_land_cover_2006_r_ref_2014.zip"
download.file(url, basename(url))
gunzip(basename(url))
rname <- list.files(getwd(), "tif$")
land_cover <- raster::raster("~/LUISA_CLC_land_coverage/clc06_r.tif")

我想聚合每个多边形的所有土地/土壤覆盖信息（总共403个要聚合土地覆盖信息的多边形）；以下是我将用于面积估计的多边形：带有多边形的shapefile可以随时使用：

我裁剪了原始土地覆盖光栅，如下所示：

deu_shp = maptools::readShapePoly("~/adm2.shp",
                                  proj4string=crswgs84,verbose=TRUE)
deu_proj <- spTransform(deu_shp, CRSobj = land_cover@crs)
land_cover_deu <- crop(land_cover, deu_proj )

为了了解使用像素计数进行面积估计，我在谷歌上搜索了相关文章，发现：使用遥感进行作物面积估计，本文提出的想法与我的兴趣非常吻合，但要实现所提出的方法纯粹是理论上的，对我来说很难做到这一点在R中。如果有任何快速而肮脏的方法可以通过像素计数进行面积估计，其中土地/土壤覆盖信息（例如城市、森林、农田）必须聚合到每个多边形（带有多边形的shapefile可在飞行中获得）。

对于像素提取，我可以简单地使用raster::extrac如下（下面的代码是试用版）：

lapply(deu_proj, function(x) {
    pixel_extract = raster::extract(land_cover_deu, deu_proj[x,])
    pixel_extract= as.data.frame(pixel_extract)
})

而上述简单的像元提取对于原始的土地覆盖物栅格来说并不是很有效。我不知道如何对每个多边形中提取的像素进行像素计数，并得到相应的面积估计（如城市、森林、农业面积等的土地覆盖率）。

有没有办法让这一切发生在R？对于给定的土地覆盖栅格，如何使用像素计数方法进行面积估计？是否可以将土地/土壤覆盖信息聚合到每个多边形？提前感谢

Answer 1

这是一个通过多边形（称为表格区域的操作）获取每个土地覆盖类别的像素数量的工作流程。其想法是使用土地覆盖光栅的分辨率和范围对形状文件进行光栅化。然后，使用基于data.table的非常有效的函数计算表格区域，该函数计算每个多边形中每个类别的像素数量。

# add ID field to the shapefile
deu_proj@data$ID <- 1:nrow(deu_proj@data)

# extract extent and resolution of land cover raster
ext <- extent(land_cover_deu)
ext <- paste(ext[1], ext[3], ext[2], ext[4])
res <- paste(res(land_cover_deu)[1], res(land_cover_deu)[2])

# rasterize shapefile using gdal (more efficent than rasterize from raster package)
# you can change GDAL_CACHEMAX according to your RAM
system(paste0("gdal_rasterize --config GDAL_CACHEMAX 8000 -a ID -te ",
               ext," -tr ", res,
               " -ot Int16 /home/deu_proj.shp /home/zones.tif"))

# load raster
zones <- raster("/home/zones.tif")

# zonal stats using myZonal function
Zstat <- myZonal(land_cover_deu, zones)

# reshape output    
results <- data.table::dcast(Zstat, z ~ vals)
colnames(results)[1] <- "ID"

# merge data
deu_proj@data <- plyr::join(deu_proj@data, results, by="ID")

# show results
head(deu_proj@data[, c(8, 17:ncol(deu_proj@data))])

           NAME_2  0   1    2    3    4  5   6   7  8   9 10  11    12  15   16    18    20   21 22    23    24    25  26
1 Alb-Donau-Kreis NA 241 4504  781 1317 NA  NA 357 NA  NA NA 133 57460  NA   NA  7212 19899 1627 NA 16403  6628 15248 135
2       Böblingen NA 369 4947 1599 1140 NA  NA 149 87 116 98 438 17845  NA 1712  1717  4742 3079 NA  3988  2286 14557  NA
3     Baden-Baden NA  45  791  150  306 NA  83  33 NA  NA 38  41   695 338  602   721   468   92 NA  1090  2272  4401  NA
4        Biberach NA 201 4681  462 1264 NA 152 312 NA  48 NA 131 46163  NA   29 23780 19126  920 NA  2291 28679  5140  NA
5   Bodenseekreis NA 268 3219  561  814  7 169  81 NA  NA NA  76  9482 418 5403  4750 19105  804 NA    96  9151  8797  NA
6        Bodensee NA  NA   13    1   NA  9  NA  NA NA  NA NA   3    NA  NA   NA     1     2   NA NA    NA    NA     4  NA
  27   29 30 31 32 34   35  36 37 39 40  41 42 43   45 46   128
1 NA  581 NA NA NA NA  104  NA NA NA NA 397 NA NA 2643 40    NA
2 NA  801 NA NA NA NA   NA  NA NA NA NA  NA NA NA 2001 58    NA
3 NA 1117 NA NA NA NA  157  NA NA NA NA  79 NA NA  486  1    NA
4 NA 2063 NA NA NA NA 1105 273 NA NA NA 384 NA NA 3760 25    NA
5 NA  240 NA NA NA NA  299  NA NA NA NA 193 NA NA 2476 66    45
6 NA   NA NA NA NA NA    4  NA NA NA NA 415 NA NA   20  1 27563

myZonal <- function (x, z, digits = 0, na.rm = TRUE) { 
  vals <- raster::getValues(x) 
  zones <- round(raster::getValues(z), digits = digits) 
  rDT <- data.table::data.table(vals, z = zones) 
  plyr::count(rDT) }

# I have loaded the shapefile using getData but you can use your own shp.
# The only difference will be in the column numbers of "show results" step.
deu_shp <- getData("GADM", country="DEU", level=2)