5

I generated, from point data, a kernel density map using GRASS, and I would like to identify the 95% volume contour in it.

Is it possible in GRASS?

Improve this question

2 Answers 2

Reset to default
4

Here is an approach in R. It is computationally expensive and a bit slow when applied to large rasters. Because of this, I added a point sub-sampling approach that seems tractable. I should note that, because the raster needs to be coerced into a vector, this function is not memory safe. This function returns the percent volume data and not a contour, but that is easy to derive from the output.

I just can't see a way around reading and indexing the rasters values. You have to sum the vector, multiply by p, sort descending and then start summing the sorted vector until the threshold is met. The index is then used to resort the vector to its original order to assign back to the raster array correctly. I added this function raster.vol into the spatialEco package (on CRAN).

Examples 

library(raster)
library(spatialEco)

  r <- raster(ncols=100, nrows=100)
    r[] <- runif(ncell(r), 0, 1)
      r <- focal(r, w=focalWeight(r, 6, "Gauss"))
        r[sample(1000, 1:ncell(r))] <- NA

Raster percent volume

p30 <- raster.vol(r, p=0.30)
p50 <- raster.vol(r, p=0.50)
p80 <- raster.vol(r, p=0.80)

par(mfrow=c(2,2))
  plot(r, col=cm.colors(10), main="original raster")
    plot(p30, breaks=c(0,0.1,1), col=c("cyan","red"),
         legend=FALSE, main="30% volume")
      plot(p50, breaks=c(0,0.1,1), col=c("cyan","red"),
          legend=FALSE, main="50% volume")
     plot(p80, breaks=c(0,0.1,1), col=c("cyan","red"),
         legend=FALSE, main="80% volume")
Improve this answer
3
  • +1 It's good finally to see a valid answer to this question. Replacing the non-sampling block of code with the following will be more than three orders of magnitude faster (16 seconds for a 10K by 10K raster). Its output compares with yours, but I haven't double-checked it for robustness at edge cases, etc. It does preserve NA cells, as one would expect. if( sample == FALSE ) { den <- getValues(x); z <- sort(den[!is.na(den)], decreasing=TRUE); y <- cumsum(as.numeric(z)); i <- sum(y <= p*y[length(y)]); return(setValues(x, den >= z[i])) }
    – whuber
    Feb 28, 2015 at 1:22
  • @whuber Ah, I see it. I have the flu and am not at my sharpest, thanks.
    – Jeffrey Evans
    Feb 28, 2015 at 1:28
  • Get well soon! I'm pleased you're alert enough to be entertaining yourself here, too. :-)
    – whuber
    Feb 28, 2015 at 1:30
-2

Use r.quantile to find the 95% percentile, then with that value fed into r.mapcalc create a new raster with values <= the 95% percentile. i.e.

GRASS 7.0.svn (ITM):~/geodata > r.quantile dem perc=95.0
Computing histogram
 100%
Computing bins
Binning data
 100%
Sorting bins
 100%
Computing quantiles
0:95.000000:735.000000

GRASS 7.0.svn (ITM):~/geodata > r.mapcalc "dem95 = if(dem<=735.0, dem, null())"
 100%

You might choose, in the r.mapcalc expression to set the new raster to a fixed value instead of the same values as the original raster. So

GRASS 7.0.svn (ITM):~/geodata > r.mapcalc "dem95 = if(dem<=735.0, 1, null())"
Improve this answer
1
  • This calculation has little or no meaning in the context of the question. A "95% volume contour" is going to be at a low level (not a high one). Its defining property is that the sum of all values in the raster exceeding this level will be 0.95 times the sum of all values in the raster. A correct answer will use the cumulative sum of the ordered raster values to find the threshold.
    – whuber
    Mar 4, 2014 at 21:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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