3

in R, I have two rasters with different extents and origins (and possibly even resolution). How to efficiently take the subset of each of them only where there are data in both of them?

Here is a code that attempts to do that (through the function getOverlap):

library(raster)

# Generate 3 dummy rasters
r1 <- raster()
extent(r1) <- extent(0, 100, 0, 500)
res(r1) <- c(5, 5)
values(r1) <- sample(10, ncell(r1), replace=TRUE)

r2 <- raster()
extent(r2) <- extent(10, 120, -10, 400) 
res(r2) <- c(5, 5)
values(r2) <- runif(ncell(r2), 1, 50)

r3 <- raster()
extent(r3) <- extent(2.5, 120, -10, 400)  # xmin is the problem
res(r3) <- c(5, 5)
values(r3) <- runif(ncell(r3), 1, 20)


# Function to return a brick of two rasters only where they superpose
getOverlap <- function(a, b) {
  ab <- intersect(a, b)
  rbrick <- brick(resample(a, ab),
                  resample(b, ab))
  rmask <- trim(sum(rbrick))
  mask(crop(rbrick, rmask), rmask)
}

# Function to check if a raster and its subset are identical (conserved)
checkConservation <- function(zbrick, a, b) {
  checka <- 
    (as(crop(a, zbrick[[1]]), 'BasicRaster') == as(zbrick[[1]], 'BasicRaster'))
  print(paste("1st raster is conserved: ",checka, sep=""))
  checkb <- 
    (as(crop(b, zbrick[[2]]), 'BasicRaster') == as(zbrick[[2]], 'BasicRaster'))
  print(paste("2d raster is conserved: ",checkb, sep=""))
}

There is a trouble when the extent happens to be not properly placed in comparison to the resolution (see that r3 is actually not conserved):

z12=getOverlap(r1,r2)
checkConservation(z12, r1, r2)
  [1] "1st raster is conserved: TRUE"
  [1] "2d raster is conserved: TRUE"

z13=getOverlap(r1,r3)
checkConservation(z13, r1, r3)
  [1] "1st raster is conserved: TRUE"
  [1] "2d raster is conserved: FALSE"

Here is an attempt to correct for this by artificially adapting the extent of r3 as a function of the resolution (it works):

r3bis <- raster()
extent0 <- c(2.5, 120, -10, 400) 
resol <- 5
xmin <- round_any(extent0[1], resol, floor) 
xmax <- round_any(extent0[2], resol, ceiling)
ymin <- round_any(extent0[3], resol, floor) 
ymax <- round_any(extent0[4], resol, ceiling)
extent(r3bis) <- c(xmin, xmax, ymin, ymax)
res(r3bis) <- c(resol, resol)
values(r3bis) <- runif(ncell(r3bis), 1, 20)

z13bis=getOverlap(r1,r3bis)
checkConservation(z13bis, r1, r3bis)
  [1] "1st raster is conserved: TRUE"
  [1] "2d raster is conserved: TRUE"

Any better way to do that generically and more concisely than above?

Is there really no built-in function to do that (common?) operation?

(i.e., replacing getOverlap --> superposing two or more raster layers with different extent/origin and automatically selecting (subsetting) areas where there are data in all of them and homogenize the extent/origin or the resulting subsets)?

Ideally, the solution would even generically do that for rasters with different resolution (i.e., take the subset of one raster where there are data in another, coarser one).

Maybe in another R package (code snippet very welcome)?

This question is related to Calculating area of overlap between 2 raster layers with different origins and extents in R and How to create a raster brick with rasters of different extents?, though I haven't really found my answer there.

closed as too broad by PolyGeo Jun 29 at 12:41

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • Welcome to GIS SE! As a new user I recommend taking the Tour and reviewing How to ask a good question?. From the Tour you will see that an important protocol is to ask only one question per question. As a result yours may be closed as too broad. However, there is an edit button beneath it that you can use to focus it down. – PolyGeo Jan 8 '15 at 8:31
  • Thanks @PolyGeo, for the welcome and for the links! What I marked as an "additional question" was actually a variant/additional request (about rasters with different resolutions), but still about what seems to me a single question (= getting efficiently the overlap between two rasters with different properties). Hence I rephrased. I tried to be clear, produced an easy example code and some attempts, to make my problem concrete... I thus hope the reader has the info, and that the question won't be closed. – ztl Jan 8 '15 at 8:52

Browse other questions tagged or ask your own question.