Calculating TWI in R?

user_12272596 0
Viewed 1485

I am trying to calculate a Topographic Wetness Index (TWI) for a DEM in R.

For the TWI, I want to use the upslope.area function in the dynatopmodel package. When running the function I get the following error: "Raster has differing x and y cell resolutions". But when looking at the properties of my raster, the x and y cell size are both 1. Has anybody had a similar problem and found a solution?

Here is my code:

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)
gis r
1 Answers

I was having the same issue and stumbled upon a fix that someone else (Ivan Lizarazo) manufactured and somewhat worked for me.

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)
}

Run both their 'upslope' and 'Create_layers' functions and then use the 'Create_layers' function on your data. The 3rd layer in the output list contains the TWI values; however, they seem to be incorrectly calculated based on the slope in degrees when they need to be calculated in radians. This caused me to have a lot of NA values in areas with shallow slope angles. This is a relatively simple fix using the 2nd layer in the output list (upslope.area) we just created and a 'slope' raster if you have one. 

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

I'm not sure if this is an issue with the upslope.area() function in the dynatopmodel package, but my DEM also had a projected crs defined and identical x and y resolutions and was throwing the same error you were experiencing.

Hope this works for you!

user_123339470
Related