使用 R 中土地覆盖栅格的像素计数方法获得面积估计的任何有效方法？

Question

我想使用像素计数法来估算每个多边形的面积。原始土地覆盖图是栅格数据，每个像素 (landcover map legend) 都分配有独特的 class。然而，通过使用像素计数来估计面积对我来说并不直观。也许，我能做的第一件事是提取每个多边形的所有像素，这些像素提供有关每个多边形内土地覆盖 classes 分布的信息。之后，我需要使用像素计数方法进行面积估计，并汇总每个多边形的城市、农业区的所有 land/soil 覆盖范围。对我来说，使用像素计数来获得面积估计并不简单，也不知道如何在 R 中完成这项工作。

请注意，原始土地覆盖地图相当大（大约 92mb），并且很难为土地覆盖光栅生成可重现的示例，因此请原谅我造成的不便。这是获取地表光栅的 R 脚本：

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")

我想聚合每个多边形的所有 land/soil 覆盖信息（总共 403 个要聚合土地覆盖信息的多边形）；以下是我将用于面积估算的多边形：shapefile with polygons are available on the fly:

我裁剪原始地表光栅如下：

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 )

为了理解像素计数的面积估计，我搜索了相关文章并发现了这篇文章： Using remote sensing for crop area estimation 这篇文章中提出的想法与我的想法非常吻合兴趣，但实现提出的方法纯粹是理论上的，我很难在 R 中做到这一点。如果有任何快速而肮脏的方法可以通过像素计数获得面积估计，我很好，其中 land/soil 覆盖信息（例如城市、森林、农业用地）必须聚合到每个多边形 (shapefile with polygons are available on the fly) .

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

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

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

有什么想法可以在 R 中实现这一点吗？对于给定的土地覆盖栅格，如何使用像素计数方法进行面积估计？是否可以将 land/soil 覆盖信息聚合到每个多边形？提前致谢

Answer 1

这是一个通过多边形获取每个土地覆盖 class 的像素数的工作流程（称为制表面积的操作）。这个想法是使用分辨率和土地覆盖栅格的范围来栅格化 shapefile。然后，使用基于 data.table 的非常有效的函数计算制表面积，该函数计算每个多边形中每个 class 的像素数。

# 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

区域统计功能（改编自here）

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)

使用 R 中土地覆盖栅格的像素计数方法获得面积估计的任何有效方法？

Any efficient way to obtain area estimation with pixel counting approach for land cover raster in R?

r

raster

data-manipulation

区域统计功能（改编自here）

示例数据