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?
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this communityHere 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")
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])) }
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())"