R中土地覆被面积的计算
基本上，我以ASCII的形式计算了一个全局分布概率模型，例如：gdpm。gdpm的值都在0到1之间。
然后，我从shape文件中导入了一个本地映射：
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制作一个栅格图，有更简洁的方法吗？
回答 2
Stack Overflow用户
发布于 2017-12-03 12:40:08
洛基的回答是好的，但这可以用光栅的方式，这是比较安全的。考虑坐标是角(经纬度)还是平面(投影)是很重要的。
示例数据
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请注意，如果您正在处理lon/lat数据，则不能使用prod(res(r))获取单元格大小。我认为，在这种情况下，您需要使用area函数并对类进行循环。
你还问过密谋的事。有许多方法来绘制一个光栅*对象。基本原则是：
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) Stack Overflow用户
发布于 2017-12-03 10:21:42
看上去你可以用像素的数量来得到面积。
让我们从一个可重复的例子开始：
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 关于地理坐标光栅的绘制，我想请您参考这里有个老问题
https://stackoverflow.com/questions/47616990
复制相似问题
领取专属 10元无门槛券
AI混元助手 在线答疑
洞察 腾讯核心技术
剖析业界实践案例
腾讯云开发者
