基本上,我计算了一个ASCII格式的全局分布概率模型,比如:gdpm
. gdpm
的值都在0和1之间。
然后我从shape文件导入了一个本地Map:
shape <- file.choose()
map <- readOGR(shape, basename(file_path_sans_ext(shape)))
字符串
下一步,我光栅化gdpm
,并使用本地Map进行裁剪:
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)
型
的数据
我已经有了一个裁剪和重分类的栅格,现在我需要汇总土地覆盖,也就是说,到每一个级别,我想计算它在本地Map的每个区域中的面积比例。(我不知道如何用术语来描述它)。我找到并遵循了一个例子(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制作一个光栅图,有更简洁的方法吗?
2条答案
按热度按时间qvsjd97n1#
loki的回答是可以的,但这可以通过光栅方式来完成,这更安全。而且考虑坐标是角坐标(经度/纬度)还是平面坐标(投影)是很重要的。
示例数据
字符串
方法1.仅适用于平面数据
型
方法2.用于Angular 数据(但也适用于平面数据)
型
如果你想按多边形汇总,你可以这样做:
型
方法1
型
方法2
检查结果:
注意,如果你正在处理经度/纬度数据,你不能使用prod(res(r))来获取单元格大小。在这种情况下,你需要使用面积函数和循环类,我想。
您还询问了绘图。有许多方法可以绘制Raster* 对象。基本方法如下:
更棘手的方法:
wvt8vs2t2#
看起来你可以通过像素数来得到面积。
让我们从一个可重复的例子开始:
字符串
的数据
由于此栅格中的值位于数据之外的另一个范围内,因此让我们将其调整为您的值:
型
之后,我们会申请您的重新分类:
型
的
现在,我们如何得到该地区?
由于您使用的是投影光栅,因此可以简单地使用像素数和光栅分辨率。因此,我们首先需要检查投影的分辨率和Map单位:
型
现在,我们知道我们正在处理40 x40米的像素,因为我们有一个度量CRS。
让我们使用这些信息来计算每个类的面积。
型
有关地理配准栅格的绘图,请参阅an older question here on SO