4

I have a DEM in QGIS and I want to calculate a ruggedness index that calculates the sum change in elevation between a grid cell and its 8 neighbor grid cells. I have seen several methods to do neighborhood analysis like r-neighbor or LecoS plugin, but from what I've seen they only can do some calculations like average, sum, maximum, etc. I need a custom analysis to build my index. The calculation for the index is this:

Index=[(Σ(el(0,0)-el(i,j))²]^½ 

Where el is the name of the elevation cell (i and j are -1, 0 or 1).

I think in ArcGis it's possible to do this custom calculation with DOCELL, but I don't know what's this, I don't have access to this program and I don't know how to use it.
Is there an option to do it with QGIS, R-statistics or anything else?

2

I think (but I didn't prove it, but the raster calc. says its valid) the syntax you're looking for is square brackets for the relative coordinates, in your case:

sqrt((el-el[-1][-1])^2 + (el-el[0][-1])^2 + ...)

Since there is no 'moving window' funcionality in QGIS Raster Calculator, you have to explicitly add all the surrounding cell values, which might become inappropriate if you have, say, 7x7 windows or so.

You might consider to have a look at r.neighbors GRASS algorithm via the toolbox, cp. How to find buffer radius around raster cell that covers a defined sum of raster cells?

| improve this answer | |
  • The suggested syntax for raster Calculator is valid but does not do what you may think. The square brackets aren't cell references but instead contain operators. At least this is so if my interpretation of "el" as "band@1" is correct – mercergeoinfo Nov 7 '19 at 14:52
3

You can do this in R with a custom function passed to the raster::focal function. In writing a function just keep in mind that each focal window represents a vector (string of values) and is, in effect, a linearized matrix. Because of this, as long as you use the na.rm=FALSE argument, you always know exactly where values are in the vector/matrix. The biggest trick is tracking nodata (NA) in the vector.

To unfold the logic here we create a matrix, that would represent a single focal window, and set the center value to 1.

( x <- matrix(nrow=3,ncol=3, runif(n=9)) )
  x[5] <- 1

Now we coerce to a vector. You will see that, given that we simulated a 3x3 window, the focal value (1) is in the fifth position of the vector. Then, say we want the summed squared deviations of every value from the focal value. Rather than indexing each value separately we rely on R recycling the vector. In this example every values is subtracted from the fifth value with each deviation squared, the resulting vector of eight squared deviations is then summed.

( x <- as.vector(x) )

x[5] - x[-5]
sum( (x[5] - x[-5])^2 )

Here is an example using the simple Riley et al., (1999) Terrain Ruggedness Index. The index is the squared-root of the summed, squared deviations from the center cell thus, requiring one to know what the center cell value to derive the deviations from that explicit value.

First, simulate a raster used in the example

library(raster)
r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
            ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
  r[] <- runif(ncell(r), 1000, 5000) 
  r <- focal(r, w=matrix(1/121,nrow=11,ncol=11)) 
plot(r)

Define function to be passed to focal

tri <- function(x) {
  xc <- x[(length(x)+1) / 2]
  x <- x[-(length(x)+1) / 2]
  x <- x[!is.na(x)]
    if(!is.na(xc) & length(x) > 0) {
      x.dev <- vector()
      for(i in 1:length(x)) { x.dev <- append(x.dev, (xc - x[i])^2) }
      return( sqrt(sum(x.dev)) )
    } else {
    return( NA )  
  }
}

Pass tri function to raster::focal

s = c(3,3) # window size
r.tri <- raster::focal(r, w=matrix(1,s[1],s[2]), fun=tri, na.rm=FALSE) 
  plot(r.tri)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.