问用像素计数法获得R区土地覆盖栅格面积估计的有效方法吗？
EN

Stack Overflow用户

提问于 2018-05-15 02:49:23

回答 1查看 645关注 0票数 2

我想用像素计数法得到每个多边形的面积估计。最初的陆地覆盖图是栅格数据，并为每个像素(地物图图例)分配了一个独特的类。然而，用像素计数来获得区域估计对我来说是不直观的。也许，我能做的第一件事就是提取每个多边形的所有像素，这些像素表示每个多边形内土地覆盖等级的分布信息。在此基础上，采用像素计数法对城市、农业各多边形的土地/土壤覆盖面积进行估算和汇总。对我来说，用像素计数来获得区域估计并不简单，也不知道如何在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")

我想聚合每个多边形的所有土地/土壤覆盖信息(总共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::extract，如下所示(下面的代码是试用的)：

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

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

有什么办法在R里实现吗？对于给定的土地覆盖栅格，如何利用像素计数方法进行面积估计？是否有可能将土地/土壤覆盖信息聚合到每个多边形？提前感谢

raster

data-manipulation

TDSQL-CMySQL免费体验

云原生数据库TDSQL-C MySQL 2核4G 免费体验15天

回答 1

Stack Overflow用户

回答已采纳

发布于 2018-05-15 04:44:25

下面是一个按多边形(一种称为表格区域的操作)获取每个土地覆盖类像素数的工作流。其想法是利用分辨率和土地覆盖栅格的范围，对shapefile进行栅格化。然后，使用一个基于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)