I am looking for a method in a R-Package to come up with an index value that represents the similarity and dissimilarity between multiple, different scaled, input rasters. These multiple input rasters could be: elevation, slope, rainfall, temperature, etc.

The multiple input rasters should be normalized/standardized. Then, using a certain method(?), generate one raster/map with values ranging between 1 (high similarity) and 0 (low similarity) between the multiple input rasters. I have been reading into methods, considering e.g. Principal Component Analysis (PCA) and Multivariate Environmental Similarity Surfaces (MESS), but could not come up with a solution to get what I want so far. Any tips on which method and R-package could perform the above?

EDIT: Let me provide some more context: For my research, I need to place camera traps to monitor wildlife. But in order to make my research statistically robust, I would like to place camera's at locations that have similar conditions, such as similar yearly rainfall, temperature, elevation, slope, etc. (which probably means there will be some Multi-collinearity).

the output should then return an index value (between 0 and 1), where 0 = pixels with dissimilar conditions and 1 = most similar conditions (thus suitable for camera traps).

library(raster) # copied from RobertH r <- raster(nrow=5, ncol=5) set.seed(20181901) s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))

scalefun <- function(r) {
(r - minValue(r)) / (maxValue(r) - minValue(r))}

ss <- scalefun(s)
...some method on ss

If my request is not feasible, are other options then to perform a:

  • Unsupervised classification (on output ss), but then I would receive classes instead of an index value. (see example code below)?
  • Simply go criteria based, by selecting thresholds (e.g. Elevation > 100 & Elevation < 200)? But I prefer to not use that option..

Unsupervised classification (multi-collinearity not considered yet)

test <- unsuperClass(ss, nSamples = 10000, nClasses = 5, nStarts = 30,
         nIter = 100, norm = F, clusterMap = T,
         algorithm = "Hartigan-Wong")
  • I don't get what you want. If you want an output in 0 - 1 range, you must make a comparison with some control data
    – aldo_tapia
    Jan 19 '18 at 18:19
  • 1
    How are you going to measure pixel-wise similarity (a local property) across single pixels of different measurement units? Is 10mm rain similar to 23C temperature and 1001mB pressure? That seems to be what you are asking.
    – Spacedman
    Jan 19 '18 at 18:29
  • If you normalise within raster layers so each layer is N(0,1) you could use the standard deviation of the values in each pixel across layers. That would be low if all values are equally above average or equally below average relative to that variables global mean and sd... Does that sound useful?
    – Spacedman
    Jan 19 '18 at 18:32
  • @Spacedman, thank you for trying to make sense of my question. I added a short explanation (see EDIT), hopefully my request makes more sense now. Any help is much appreciated.
    – JSD
    Jan 19 '18 at 20:21

Always provide some example data, please:

r <- raster(nrow=5, ncol=5)
s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))

I have edited my answer based on the edited question. I think you should start with (at least) one location you want to include, and then compute similarity to the given location

#scale raster data
s <- scale(s) 

#arbitrary location you like:
xy <- cbind(0,0)
v <- extract(s, xy) 
d <- abs(s - as.vector(v))
dd <- sum(d)

# similar sites
plot(dd < 2)
plot(dd < 3)
  • Thank you for your quick reply! I added some more context in my question above EDIT,to explain, in short, my research idea. Does your answer then still provides the solution? Your help is much appreciated.
    – JSD
    Jan 19 '18 at 20:20

Here is my approach:





temp <- srtm

temp[] <- srtm[]*-.1 + 30 + rnorm(ncell(temp),sd = 1.5)

slope <- terrain(srtm, unit = 'tangent')

rain <- srtm

df <- as.data.frame(srtm, xy = T)

rain[] <- df$x/max(df$x) * 5 + (df$y/max(df$y))^(5/3) * 5 + rnorm(ncell(temp),sd = 0.01)

s <- stack(srtm, slope, temp, rain)

names(s) <- c('dem','slope','temp','pp')


enter image description here

Using the approach of mister @RobertH

# assuming only positive values

scalefun <- function(r, precision = 2) {
  v <- round((r - minValue(r)) / (maxValue(r) - minValue(r)),precision)
  v[is.na(v)] <- 0 # avoid NA in flat terrains (slope layer)

ss <- scalefun(s)


enter image description here

precision argument round values. It's useful because you can find similar places by the number of occurrences:

# search similar places


cond <- as.data.frame(ss) %>% group_by(dem,slope,temp,pp) %>%
  summarise(n = n()) %>% arrange(desc(n)) %>% head(1)


## # A tibble: 1 x 5
## # Groups: dem, slope, temp [1]
##      dem slope  temp    pp     n
##    <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.0600     0 0.750 0.630    25

You have 25 places with similar conditions to choose. You can either explore beyond the first element.


cond <- as.vector(unlist(cond))

selected <- ss[[1]] == cond[1] & ss[[2]] == cond[2] & ss[[3]] == cond[3] & ss[[4]] == cond[4]

selected <- as.data.frame(selected,xy = T)

head(selected[selected[,3] == T,])
##            x       y layer
## 46746 626940 -415080  TRUE
## 47318 626880 -415140  TRUE
## 47560 625530 -415170  TRUE
## 48190 627210 -415230  TRUE
## 50170 626340 -415440  TRUE
## 51631 627120 -415590  TRUE

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.