Расчет TWI в R?

Question

Я пытаюсь вычислить индекс топографической влажности (TWI) для матрицы высот в R.

Для TWI я хочу использовать функцию upslope.area в пакете dynatopmodel. При запуске функции я получаю следующую ошибку: «Растр имеет различное разрешение ячеек x и y». Но если посмотреть на свойства моего растра, то размеры ячеек x и y равны 1. Кто-нибудь сталкивался с подобной проблемой и нашел решение?

Вот мой код: twi <- upslope.area (marion_dem, log-TRUE, atb = TRUE, deg = 0.1, fill.sinks = TRUE) </p>

Свойства растра: класс: размеры 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

У меня возникла та же проблема, и я наткнулся на исправление, которое кто-то другой (Иван Лизаразо) создал и несколько сработало для меня.

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 в областях с малыми углами наклона. Это относительно простое исправление, использующее 2-й слой в списке вывода (upslope.area), который мы только что создали, и растр «наклон», если он у вас есть.

twi.man <- ln(upslope.area / tan(slope / 180))

Я не уверен, что это проблема с функцией upslope.area () в пакете dynatopmodel, но в моей DEM также были определены проецируемые crs и идентичные разрешения x и y, и он выбрасывалта же самая ошибка, с которой вы столкнулись.

Надеюсь, это работает для вас!