如何计算R中的土地覆盖面积

Question

基本上，我以 ASCII 的形式计算了一个全局分布概率模型，比如：

gdpm

。

gdpm

的值都在0到1之间。

然后我从形状文件导入了本地地图：

shape <- file.choose()  
map <- readOGR(shape, basename(file_path_sans_ext(shape)))

下一步，我光栅化

gdpm

，并使用本地地图进行裁剪：

ldpm <- mask(gdpm, map)

然后，我将这个连续模型重新分类为离散模型（我将模型分为6个级别）：

recalc <- matrix(c(0, 0.05, 0, 0.05, 0.2, 1, 0.2, 0.4, 2, 0.4, 0.6, 3, 0.6, 0.8, 4, 0.8, 1, 5), ncol = 3, byrow = TRUE) 

ldpmR <- reclassify(ldpm, recalc)

我已经有了一个裁剪并重新分类的栅格，现在我需要总结土地覆盖，也就是说，对于每个级别，我想计算它在本地地图每个区域中所占的面积比例。（我不知道如何用术语来描述它）。我找到并遵循了一个例子（RobertH）：

ext <- raster::extract(ldpmR, map)

tab <- sapply(ext, function(x) tabulate(x, 10))
tab <- tab / colSums(tab)

但我不确定它是否有效，因为

tab

的输出令人困惑。 那么如何正确计算土地覆盖面积呢？如何在每个多边形内应用正确的方法？

我的原始数据太大，我只能提供替代栅格（我认为这个例子应该应用不同的重分类矩阵）：

示例光栅

或者您可以生成测试栅格（RobertH）：

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster")) 
writeRaster(s, file='testtif', format='GTiff', bylayer=T, overwrite=T)
f <- list.files(pattern="testtif_..tif")

我还有一个关于绘制栅格的问题：

r <- as(r, "SpatialPixelsDataFrame")
r <- as.data.frame(r)
colnames(r) <- c("value", "x", "y")

我进行此转换是为了使用 ggplot2 进行栅格绘图，有更简洁的方法吗？

Answer 1

loki的答案是可以的，但这可以通过光栅方式完成，这样更安全。重要的是要考虑坐标是角度（经度/纬度）还是平面（投影）

示例数据

library(raster)
r <- raster(system.file("external/test.grd", package="raster"))
r <- r / 1000
recalc <- matrix(c(0, 0.05, 0, 0.05, 0.2, 1, 0.2, 0.4, 2, 0.4, 0.6, 3, 0.6, 0.8, 4, 0.8, 2, 5), ncol = 3, byrow = TRUE) 
r2 <- reclassify(r, recalc)

方法 1. 仅适用于平面数据

f <- freq(r2, useNA='no')
apc <- prod(res(r))
f <- cbind(f, area=f[,2] * apc)
f

#     value count    area
#[1,]     1    78  124800
#[2,]     2  1750 2800000
#[3,]     3   819 1310400
#[4,]     4   304  486400
#[5,]     5   152  243200

方法 2. 对于角度数据（但也适用于平面数据）

a <- area(r2)
z <- zonal(a, r2, 'sum')
z
#     zone     sum
#[1,]    1  124800
#[2,]    2 2800000
#[3,]    3 1310400
#[4,]    4  486400
#[5,]    5  243200

如果你想通过多边形进行总结，你可以这样做：

# example polygons
a <- rasterToPolygons(aggregate(r, 25))

方法1

# extract values (slow)
ext <- extract(r2, a)

# tabulate values for each polygon
tab <- sapply(ext, function(x) tabulate(x, 5))
# adjust for area (planar data only)
tab <- tab * prod(res(r2))

# check the results, by summing over the regions
rowSums(tab)
#[1]  124800 2800000 1310400  486400  243200    

方法2

x <- rasterize(a, r2)
z <- crosstab(x, r2)
z <- cbind(z, area = z[,3] * prod(res(r2)))

检查结果：

aggregate(z[, 'area', drop=F], z[,'Var2', drop=F], sum)
  Var2    area
#1    1  124800
#2    2 2800000
#3    3 1310400
#4    4  486400
#5    5  243200

请注意，如果您正在处理经/纬度数据，则无法使用 prod(res(r)) 来获取像元大小。在这种情况下，我认为您将需要使用区域函数并循环遍历类。

你还问了关于绘图的问题。绘制 Raster* 对象的方法有多种。基本的有：

 image(r2)
 plot(r2)
 spplot(r2)

 library(rasterVis); 
 levelplot(r2)

更棘手的方法：

 library(ggplot2) # using a rasterVis method
 theme_set(theme_bw())
 gplot(r2) + geom_tile(aes(fill = value)) +
      facet_wrap(~ variable) +
      scale_fill_gradient(low = 'white', high = 'blue') +
      coord_equal()


 library(leaflet)
 leaflet() %>% addTiles() %>%
 addRasterImage(r2, colors = "Spectral", opacity = 0.8)       

Answer 2

好像可以通过像素数得到面积。
让我们从一个可重现的示例开始：

r <- raster(system.file("external/test.grd", package="raster"))
plot(r)

由于该栅格中的值与您的数据处于不同的范围，因此让我们将它们调整为您的值：

r <- r / 1000
r[r>1,] <- 1

之后，我们将应用您的重新分类：

recalc <- matrix(c(0, 0.05, 0, 0.05, 0.2, 1, 0.2, 0.4, 2, 0.4, 0.6, 3, 0.6, 0.8, 4, 0.8, 1, 5), ncol = 3, byrow = TRUE) 
r2 <- reclassify(r, recalc)
plot(r2)

现在，我们如何获得面积？
由于您正在使用投影光栅，因此您可以简单地使用像素数和光栅分辨率。因此，我们首先需要检查投影的分辨率和地图单位：

res(r)
# [1] 40 40
crs(r)
# CRS arguments:
#   +init=epsg:28992
# +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea
# +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000
# +y_0=463000 +ellps=bessel +units=m +no_defs

现在，我们知道我们正在处理 40 x40 米的像素，因为我们有公制 CRS。

让我们使用这些信息来计算每个类别的面积。

app <- res(r)[1] * res(r)[2] # area per pixel

table(r2[]) * app
#      1       2       3       4       5 
# 124800 2800000 1310400  486400  243200 

对于地理参考栅格的绘制，我想向您推荐SO 上的一个较旧问题

Answer 3

如何重新计算<- matrix(c(0, 0.05, 0, 0.05, 0.2, 1, 0.2, 0.4, 2, 0.4, 0.6, 3, 0.6, 0.8, 4, 0.8, 1, 5), ncol = 3, byrow = TRUE) run successfully with ncol = 3? It only worked when ncol was equal to the number of scales (ex: 0.05, 0, 0.05, 0.2 -> ncol =4)

但是只有当 rcl 有 2 或 3 列时才可以运行 reclassify(r, recalc) ..