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

Question

基本上，我用ASCII的形式计算了一个全局分布概率模型，例如： gdpm。 gdpm的值全部为0和1之间

然后我输入从形状文件的本地地图：

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

下一步，我光栅化gdpm，并裁剪使用本地地图：

ldpm <- mask(gdpm, map) 

然后，我重新分类该连续模型转换为离散的模型（I划分模型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输出是混乱的。 那么如何正确计算土地覆盖面积呢？我如何在每个多边形内应用正确的方法？

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

Example raster

或者你可以生成测试光栅（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

好像你可以通过像素数获得面积。
让我们先从一个可重复的例子：

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×40米的像素，因为我们有一个指标CRS。

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

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

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

对于地理参考栅格的绘制，我想请您看看an older question here on SO

Answer 2

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 

如果要总结e。通过多边形，你可以做这样的事情：

# 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 

请注意，如果您正在处理您可以输入的lon/lat数据不使用prod（res（r））来获取单元大小。在这种情况下，你会需要使用区域函数并遍历类，我想。

您还问过关于绘图。绘制栅格*对象的方法很多。基本的是：

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)