I have a raster stack that contains terrain bends extracted with different neighborhood windows "crosc_nxn" (e.g. 3x3, 5x5, etc ...).

crosc_nxn <- raster::stack(path.files.list)
> crosc_nxn
class       : RasterStack 
dimensions  : 170, 172, 29240, 5  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 794418.4, 799578.4, 7400299, 7405399  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=22 +south +ellps=WGS84 +units=m +no_defs 
names       : crosc_3, crosc_5, crosc_9, crosc_11, crosc_15 

I want to create a new stack containing only the extreme values of each layer using two quantile ranges ( <= 0.01 or 0.99 >=).

out_prob <- function(x, min= quantile(x, 0.01,na.rm=TRUE),  max= quantile(x, 0.99, na.rm=TRUE)){ifelse(x <= min | x>= max , x, NA)}

range_calc <- calc(crosc_nxn, fun= out_prob)
range_over <- overlay(crosc_nxn, fun= out_prob)

It appears that calc applies the function to a certain pixel and traverses all layers to return the desired value. However I want my function to sweep all the pixels of the layer and then move on to the next one. I'm surprised that overlay has returned the same result.

On the other hand when specifying the quantile for a certain layer the result is consistent with my query, in that layer where I expressed the interval.

    # calc quantile
    > quantile(crosc_nxn, c(0.01, 0.99), na.rm=TRUE)
                       1%         99%
    crosc_3  -0.001324438 0.001358304
    crosc_5  -0.000933790 0.000937456
    crosc_9  -0.000481476 0.000486876
    crosc_11 -0.000356400 0.000362014
    crosc_15 -0.000223900 0.000228338

    # explicit range in function    
    out_prob_ex <- function(x) {ifelse(x <= -0.001324438 | x >= 0.001358304, x, NA)}

    # generic 
    out_prob <- function(x, min= quantile(x, 0.01,na.rm=TRUE),  max= quantile(x, 0.99, na.rm=TRUE)){ifelse(x <= min | x>= max , x, NA)}

    quantile_generic <- overlay(crosc_nxn, fun= out_prob)
    quantile_explicit <- overlay(crosc_nxn, fun= out_prob_ex)

    # comparison using only the first stack layer
    par(mfrow = c(2,2))
    plot(outlier_generic[[1]], main= "generic formula")
    plot(outlier_explicit[[1]], main= "explicit quantile")

enter image description here

I need a generic function for each layer, without having to express the quantiles one by one. I thought overlay would do this and I can not figure it out. Someone has an explanation to help me see the way out.

  • If I understand you correctly, you want an output stack of 10 layers. Each layer is a raster that only has a single value. out[[1]] is the 0.01 quantile of the values in crosc_nxn[[1]], out[[2]] is the .99 quantile of crosc_nxn[[2]], out[[3]] = .01 quantile of crosc_nxn[[3]] and so on?
    – Spacedman
    Oct 25, 2018 at 21:30
  • Dear @Spacedman. I need return a stack of 5 layers. One for each input of crosc_nxn (3x3, 5x5, ... 15x15). The output layer will receive two intervals of each input layer. out[[1]] <- x [[1]] [ x [[1]] <= quantile (x [[1]], 0.01) | > = quantile (x, 0.99) ]
    – viniciovcl
    Oct 26, 2018 at 1:53

3 Answers 3


Lets make a test raster to make it easy to see what's going on:

> r = raster(matrix(1:100,10,10))
> plot(r)

enter image description here

The highest and lowest values are in the left and right columns, so we can see if we are selecting them properly.

This function takes a single raster and returns the parts of the raster that are outside the given quantiles, and NA elsewhere:

> extremes = function(r,q=c(0.01,0.99)){
        qua = quantile(r[],q,na.rm=TRUE)
        Q = r<qua[1] | r > qua[2]

It works by computing the TRUE/FALSE raster that satisfies the condition, replacing FALSE (zero) with NA, and then multiplying by the raster values. That gives NA where FALSE and r values where TRUE. Test:

> plot(extremes(r, c(.2,.8)))

enter image description here

I'm using narrower quantile bands so I can see what's going on. With only 100 pixels I wouldnt see much with a 0.01 and 0.99 level. This looks good so far. Now let's construct a test stack:

> s = stack(list(r,r*2,r*3,r*4))
> plot(s)

enter image description here

To apply a function layer-wise, loop over layers with lapply and stack the returned layer list:

> exs = stack(lapply(1:nlayers(s), function(i){extremes(s[[i]], q=c(.2,.8))}))
> plot(exs)

enter image description here

Which I think is your final answer, once you've run it with your desired quantile levels.

  • thanks @Spacedman. very correct. Many thanks for your generous resolution.
    – viniciovcl
    Oct 29, 2018 at 15:54

Below is a more memory safe version of @Spacedman's extremes function

extremes <- function(r,q=c(0.01,0.99)){
    qua <- quantile(r, q, na.rm=TRUE)
    reclassify(r, cbind(qua[1], qua[2], NA))
  • Is the last mask step needed? The values in Q will be the values from r and the mask will do nothing...
    – Spacedman
    Oct 26, 2018 at 14:31

The question and the answers gave solution if anyone wanted to remove values that is within a set range in a raster or raster stack. So, I thought it is useful to modify the function that will do the opposite (i.e., keep values that falls in between a set range of values). This can be particularly useful for beginners trying to find a solution.

The solution below is based on answer given by @Spacedman

# make a new raster
r = raster(matrix(1:100,10,10))

# function to remove values outside of desired range
extremesOut = function(r,q=c(0.01,0.99)){
  qua = quantile(r[],q,na.rm=TRUE)
  Q = r>qua[1] & r < qua[2]

plot(extremesOut(r, c(.2,.8))) # keep only values that fall between 0.2 and 0.8

s = stack(r,r*2,r*3,r*4)
exs = stack(lapply(1:nlayers(s), function(i){extremesOut(s[[i]], q=c(.2,.8))}))

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.