r - 在R中计算TWI?
问题描述
我正在尝试计算 R 中 DEM 的地形湿度指数 (TWI)。
对于 TWI,我想使用 dynatopmodel 包中的 upslope.area 函数。运行该函数时,出现以下错误:“光栅具有不同的 x 和 y 单元格分辨率”。但是当查看我的栅格的属性时,x 和 y 像元大小都是 1。有没有人遇到过类似的问题并找到了解决方案?
这是我的代码:
twi <- upslope.area(marion_dem, log-TRUE, atb=TRUE, deg=0.1, fill.sinks=TRUE)
Raster properties:
class : RasterLayer
dimensions : 22967, 30492, 700309764 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 41539.28, 72031.28, -5208329, -5185362 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : D:\Marion GIS\DWH\Marion_DSM.tif
names : Marion_DSM
values : -33, 1248 (min, max)
解决方案
我遇到了同样的问题,偶然发现了其他人(Ivan Lizarazo)制造的修复程序,并且在某种程度上为我工作。
upslope <- function (dem, log = TRUE, atb = FALSE, deg = 0.12, fill.sinks = TRUE)
{
if (!all.equal(xres(dem), yres(dem))) {
stop("Raster has differing x and y cell resolutions. Check that it is in a projected coordinate system (e.g. UTM) and use raster::projectRaster to reproject to one if not. Otherwise consider using raster::resample")
}
if (fill.sinks) {
capture.output(dem <- invisible(raster::setValues(dem, topmodel::sinkfill(raster::as.matrix(dem), res = xres(dem), degree = deg))))
}
topidx <- topmodel::topidx(raster::as.matrix(dem), res = xres(dem))
a <- raster::setValues(dem, topidx$area)
if (log) {
a <- log(a)
}
if (atb) {
atb <- raster::setValues(dem, topidx$atb)
a <- addLayer(a, atb)
names(a) <- c("a", "atb")
}
return(a)
}
create_layers <- function (dem, fill.sinks = TRUE, deg = 0.1)
{
layers <- stack(dem)
message("Building upslope areas...")
a.atb <- upslope(dem, atb = TRUE, fill.sinks = fill.sinks, deg = deg)
layers <- addLayer(layers, a.atb)
names(layers) <- c("filled.elevations", "upslope.area", "twi")
return(layers)
}
运行它们的“upslope”和“Create_layers”函数,然后对数据使用“Create_layers”函数。输出列表中的第 3 层包含 TWI 值;但是,当它们需要以弧度计算时,它们似乎根据以度为单位的斜率计算不正确。这导致我在倾斜角度较浅的区域有很多 NA 值。这是一个相对简单的修复,使用我们刚刚创建的输出列表 (upslope.area) 中的第二层和“坡度”栅格(如果有的话)。
twi.man <- ln(upslope.area / tan(slope / 180))
我不确定这是否是 dynatopmodel 包中的 upslope.area() 函数的问题,但我的 DEM 还定义了投影 crs 并且 x 和 y 分辨率相同,并且抛出了您遇到的相同错误。
希望这对你有用!
推荐阅读
- mysql - Mysql子查询建议
- java - 如何在 android 应用程序中显示我的动画 svg 文件?
- javascript - 如何在 React Hooks 中应用输入验证
- python - 从 json 输出中删除一个字符
- vaticle-typedb - 使用 Workbase 可视化架构
- javafx - 我在场景构建器中安装 fontawesome 库时遇到问题
- angular - Angular 9 上的 Observable 类型不存在属性
- swift - Swift macOS 如何禁用菜单栏按钮
- python - [(a and b) in c] 和 [a in c and b in c] 的区别
- opengl - 防止 GLSL vec4 输出的 w 乘法?