Identifying overlap zones in R raster package

Question

Package:

• raster

Data:

• A rasterStack with 10 bands.
• Each of the bands contains an image area surrounded by NAs
• Bands are logical, i.e. "1" for image data and "0"/NA for surrounding area
• The "image areas" of each band do not align completely with each other, though most have partial overlaps

Objective:

• Write a fast function that can return either a rasterLayer or cell numbers for each "zone", for instance a pixel containing data only from bands 1 and 2 falls in zone 1, a pixel containing data only from bands 3 and 4 falls in zone 2, etc. If a rasterLayer is returned, I need to be able to match the zone value with band numbers later.

First attempt:

# Possible band combinations
values = integer(0)
for(i in 1:nlayers(myraster)){
 combs = combn(1:nlayers(myraster), i)
 for(j in 1:ncol(combs)){
  values = c(values, list(combs[,j]))
 }
}

# Define the zone finding function
find_zones = function(bands){

 # The intersection of the bands of interest
 a = subset(myraster, 1)
 values(a) = TRUE
 for(i in bands){
  a = a & myraster[[i]]
 }

 # Union of the remaining bands
 b = subset(myraster, 1)
 values(b) = FALSE
 for(i in seq(1:nlayers(myraster))[-bands]){
  b = b | myraster[[i]]
 }

 #plot(a & !b)
 cells = Which(a & !b, cells=TRUE)
 return(cells)
}

# Applying the function
results = lapply(values, find_zones)

My current function takes a very long time to execute. Can you think of a better way? Note that I don't simply want to know how many bands have data at each pixel, I also need to know which bands. The purpose of this is to process different the areas differently afterwards.

Note also that the real-life scenario is a 3000 x 3000 or more raster with potentially more than 10 bands.

---

EDIT

Some sample data consisting of 10 offset image areas:

# Sample data
library(raster)    
for(i in 1:10) {
  start_line = i*10*1000
  end_line = 1000000 - 800*1000 - start_line
  offset = i * 10
  data = c(rep(0,start_line), rep(c(rep(0,offset), rep(1,800), rep(0,200-offset)), 800), rep(0, end_line))
  current_layer = raster(nrows=1000, ncols=1000)
  values(current_layer) = data
  if(i == 1) {
    myraster = stack(current_layer)
  } else {
    myraster = addLayer(myraster, current_layer)
  }
}
NAvalue(myraster) = 0  # You may not want to do this depending on your solution...

Answer 1

EDIT : Answer updated using Nick's trick and matrix multiplication.

---

You could try the following function, optimized by using Nick's trick and matrix multiplication. The bottleneck now is filling up stack with the seperate layers, but I guess the timings are quite OK now. Memory usage is a bit less, but given your data and the nature of R, I don't know if you can nibble of a bit without hampering the performance big time.

> system.time(T1 <- FindBands(myraster,return.stack=T))
   user  system elapsed 
   6.32    2.17    8.48 
> system.time(T2 <- FindBands(myraster,return.stack=F))
   user  system elapsed 
   1.58    0.02    1.59 
> system.time(results <- lapply(values, find_zones))
  Timing stopped at: 182.27 35.13 217.71

The function returns either a rasterStack with the different level combinations present in the plot (that's not all possible level combinations, so you have some gain there already), or a matrix with the level number and level names. This allows you to do something like :

levelnames <- attr(T2,"levels")[T2]

to get the level names for each cell point. As shown below, you can easily put that matrix inside a rasterLayer object.

The function :

 FindBands <- function(x,return.stack=F){
    dims <- dim(x)
    Values <- getValues(x)
    nn <- colnames(Values)

    vec <- 2^((1:dims[3])-1)
    #Get all combinations and the names
    id <- unlist(
                lapply(1:10,function(x) combn(1:10,x,simplify=F))
              ,recursive=F)

    nameid <- sapply(id,function(i){
      x <- sum(vec[i])
      names(x) <- paste(i,collapse="-")
      x
    })
    # Nicks approach
    layers <- Values %*% vec
    # Find out which levels we need
    LayerLevels <- unique(sort(layers))
    LayerNames <- c("No Layer",names(nameid[nameid %in% LayerLevels]))

    if(return.stack){
        myStack <- lapply(LayerLevels,function(i){
          r <- raster(nr=dims[1],nc=dims[2])
          r[] <- as.numeric(layers == i)
          r
          } )
        myStack <- stack(myStack)
        layerNames(myStack) <- LayerNames
        return(myStack)

    } else {

      LayerNumber <- match(layers,LayerLevels)
      LayerNumber <- matrix(LayerNumber,ncol=dims[2],byrow=T)
      attr(LayerNumber,"levels") <- LayerNames
      return(LayerNumber)
    }    
}

Proof of concept, using the data of RobertH :

r <- raster(nr=10, nc=10)
r[]=0
r[c(20:60,90:93)] <- 1
s <- list(r)
r[]=0
r[c(40:70,93:98)] <- 1
s <- c(s, r)
r[]=0
r[50:95] <- 1
s <- (c(s, r))
aRaster <- stack(s)


> X <- FindBands(aRaster,return.stack=T)
> plot(X)

> X <- FindBands(aRaster,return.stack=F)
> X
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    1    1    1    1    1    1    1    1     1
 [2,]    1    1    1    1    1    1    1    1    1     2
 [3,]    2    2    2    2    2    2    2    2    2     2
 [4,]    2    2    2    2    2    2    2    2    2     4
 [5,]    4    4    4    4    4    4    4    4    4     8
 [6,]    8    8    8    8    8    8    8    8    8     8
 [7,]    7    7    7    7    7    7    7    7    7     7
 [8,]    5    5    5    5    5    5    5    5    5     5
 [9,]    5    5    5    5    5    5    5    5    5     6
[10,]    6    6    8    7    7    3    3    3    1     1
attr(,"levels")
[1] "No Layer" "1"        "2"        "3"        "1-2"      "1-3"
       "2-3"      "1-2-3"   

> XX <- raster(ncol=10,nrow=10)
> XX[] <- X
> plot(XX)

Answer 2

How about this?

library(raster)
#setting up some data

r <- raster(nr=10, nc=10)
r[]=0
r[c(20:60,90:93)] <- 1
s <- list(r)
r[]=0
r[c(40:70,93:98)] <- 1
s <- c(s, r)
r[]=0
r[50:95] <- 1
s <- (c(s, r))
plot(stack(s))

# write a vectorized function that classifies the data
# 
fun=function(x,y,z)cbind(x+y+z==0, x==1&y+z==0, y==1&x+z==0, z==1&x+y==0, x==0&y+z==2, y==0&x+z==2, z==0&x+y==2,x+y+z==3)

z <- overlay(s[[1]], s[[2]], s[[3]], fun=fun)
# equivalent to
#s <- stack(s)
#z <- overlay(s[[1]], s[[2]], s[[3]], fun=fun)

ln <- c("x+y+z==0", "x==1&y+z==0", "y==1&x+z==0", "z==1&x+y==0", "x==0&y+z==2", "y==0&x+z==2", "z==0&x+y==2", "x+y+z==3")
layerNames(z) <- ln
x11()
plot(z)

more generic:

s <- stack(s)
fun=function(x)as.numeric(paste(which(x==1), collapse=""))
x <- calc(s,fun)

this is not good when nlayers(s) has double digits ("1", "2" is the same as "12", and in those cases you could use the function below (fun2) instead:

fun2=function(x)as.numeric(paste(c(9, x), collapse=""))
x2 <- calc(s,fun2)

unique(x)
# [1]   1   2   3  12  13  23 123

unique(x2)
# [1] 9000 9001 9010 9011 9100 9101 9110 9111

for the toy example only:

plot(x)
text(x)
p=rasterToPolygons(x)
plot(p, add=T)

Answer 3

I'm not familiar with raster, but from what I grasp from the above, you essentially have a 10*3000*3000 array, right?

If so, for each position in the raster (second and third indices, currow and curcol), you can calculate a unique identifier for its 'zone' by using binary: run i over the 'bands' (first index) and sum r[i,currow, curcol]*2^(i-1). Depending on the internal workings of raster, it should be possible to have a rather quick implementation of this.

This results in a new 'raster' of size 3000*3000 holding the unique identifiers of each position. Finding the unique values in there gives you back the zones that actually occur in your data, and reversing the binary logic should give you the bands that belong to a given zone.

Pardon me if my interpretation of raster is incorrect: then please ignore my musings. Either way not a complete solution.

Answer 4

I've written code for @Nick Sabbe's suggestion, which I think is very concise and relatively fast. This assumes that the input rasterStack already has logical 1 or 0 data:

# Set the channels to 2^i instead of 1
bands = nlayers(myraster)
a = stack()
for (i in 1:bands) {
  a = addLayer(a, myraster[[i]] * 2^i)
}
coded = sum(a)
#plot(coded)
values = unique(coded)[-1]
remove(a, myraster)

# Function to retrieve which coded value means which channels
which_bands = function(value) {
  single = numeric()
  for (i in bands:1) {
    if ((0 < value) & (value >= 2^i)) {
     value = value - 2^i
      single = c(single, i)
    }
  }
  return(single)
}