在R中计算TWI？

Question

我正在尝试计算 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)

Answer 1

我遇到了同样的问题，偶然发现了其他人 (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 分辨率相同，并且抛出了与您遇到的相同的错误.

希望这对你有用！