计算R中的TWI？

Question

我正在尝试为R中的DEM计算地形湿度指数（TWI）

对于TWI，我想使用dynatopmodel软件包中的upslope.area函数。运行该函数时，出现以下错误：“栅格的x和y单元格分辨率不同”。但是，当查看栅格的属性时，x和y像元大小均为1。有人遇到过类似的问题并找到了解决方案吗？

这是我的代码：twi

栅格属性：类别：RasterLayer尺寸：22967、30492、700309764（nrow，ncol，ncell）分辨率：1，1（x，y）范围：41539.28、72031.28，-5208329，-5185362（xmin，xmax，ymin，ymax）协调。参考：+ proj = longlat + datum = WGS84 + no_defs + ellps = WGS84 + towgs84 = 0,0,0数据源：D：\ Marion GIS \ DWH \ Marion_DSM.tif名字：Marion_DSM值：-33、1248（最小值，最大值）

谢谢！

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分辨率，并抛出与您遇到的错误相同的错误。

希望这对您有用！