Bring rasters to same extent by filling with NA - in R

Question

I have several cropped rasters with different geometry/outlines. Specifically spatial yield maps from several years of the same field, but the extent varies - the measurements were not always overall the whole field, but in some years only part of it. I want to calculate a mean value of those maps and combine them into one mean-value-raster. That does mean however, that not for every pixel in let's say 5 layers/rasters there is a value. I could accept these missing values to be NA, so the final mean value would only be calculated by let's say 3 rasters for parts of the field, where the maps are not overlapping.

I thought of extending the raster with 'extend{raster}', filling the non-overlapping parts with NA values:

y <- extend(y, shape, value=NA) #Shape is a rectangular shape that enframes all yield map rasters

That works fine, for all rasters. But they still don't have the same extent. Even if I adjust the extent by setExtent() or extent() <- extent() to the extent of the rectangular shapefile or even to one of the other extended rasters, I still get:

Error in compareRaster(x) : different number or columns

..when I want to stack them and use calc(y, fun=mean,...). The original raster extents are too different to resample. But they do have the same resolution and CRS.

Has anyone an idea how to solve this?

Answer 1

If you want to achieve an operation such as calc(y, fun=mean,...), you could get the minimal common extent of your rasters, and crop them all to that extent before stacking them together and applying the operation.

Suppose you have three rasters:

# Generate 3 dummy rasters with different extents
r1 <- raster( crs="+proj=utm +zone=31")
extent(r1) <- extent(0, 100, 0, 500)
res(r1) <- c(5, 5)
values(r1) <- sample(10, ncell(r1), replace=TRUE)

r2 <- raster( crs="+proj=utm +zone=31")
extent(r2) <- extent(10, 120, -10, 400) 
res(r2) <- c(5, 5)
values(r2) <- runif(ncell(r2), 1, 10)

r3 <- raster( crs="+proj=utm +zone=31")
extent(r3) <- extent(50, 150, 30, 200) 
res(r3) <- c(5, 5)
values(r3) <- runif(ncell(r3), 1, 10)

A first way to do that:

# Summing your rasters will only work where they are not NA
r123 = r1+r2+r3 # r123 has the minimal common extent
r1 = crop(r1, r123) # crop to that minimal extent
r2 = crop(r2, r123)
r3 = crop(r3, r123)

s123 = stack(r1, r2, r3)
s123.mean = calc(s123, fun=mean)

Another one:

# Manually compute the minimal extent
xmin <- max(bbox(r1)[1,1], bbox(r2)[1,1], bbox(r3)[1,1])
xmax <- min(bbox(r1)[1,2], bbox(r2)[1,2], bbox(r3)[1,2])  
ymin <- max(bbox(r1)[2,1], bbox(r2)[2,1], bbox(r3)[2,1])  
ymax <- min(bbox(r1)[2,2], bbox(r2)[2,2], bbox(r3)[2,2])  
newextent=c(xmin, xmax, ymin, ymax)
r1 = crop(r1, newextent)
r2 = crop(r2, newextent)
r3 = crop(r3, newextent)

s123 = stack(r1, r2, r3)
s123.mean = calc(s123, fun=mean)

I am not too sure of what you exactly want, but if you really want to keep all your rasters complete (no crop operation), you could also compute the largest extent (same as the second method above, just invert min and max) and extend all you rasters to that one before stacking them (same thing, you just end up with larger rasters, filled with NA's... not sure this is really desired):

xmin <- min(bbox(r1)[1,1], bbox(r2)[1,1], bbox(r3)[1,1])
xmax <- max(bbox(r1)[1,2], bbox(r2)[1,2], bbox(r3)[1,2])  
ymin <- min(bbox(r1)[2,1], bbox(r2)[2,1], bbox(r3)[2,1])  
ymax <- max(bbox(r1)[2,2], bbox(r2)[2,2], bbox(r3)[2,2])  
newextent=c(xmin, xmax, ymin, ymax)
r1 = extend(r1, newextent)
r2 = extend(r2, newextent)
r3 = extend(r3, newextent)
s123 = stack(r1, r2, r3)
s123.mean = calc(s123, fun=mean)

Does that help?

NOTE

You might not want to play with setExtent() or extent() <- extent(), as you could end with wrong geographic coordinates of your rasters (i.e., the extent would be modified, but not the content, so you'd actually translate your rasters likely without aiming at that as the resulting overlap between them would be physically meaningless).