4

I am storing ~ 30 raster files in a folder and want to create a raster stack. I know that not all of these rasters have completely the same extent. They do have the same resolution and CRS.

Thus I try to set the same extent for all of my rasters in this folder with the following loop in R.

#reference raster to get the extent    
a <- raster("F:/SDM/DATA/env_tiff/e_veg_d10.tif")

# where the new rasters will be stored
    outpath <- "F:/SDM/DATA/temp/"
    dir.create(outpath)

    files <- list.files(path="F:/SDM/DATA/env_tiff", pattern=".tif$")   

    # add output directory
    outfiles <- paste0(outpath, files)
    # change extensions
    extension(outfiles) <- 'tif'

    for(i in 1:length(files)) {
      e <- extent(a)
      r <-raster(files[i])
      rc <- crop(r, e)
      rc<- mask(rc, e)
      rw <- writeRaster(rc, outfiles[i], overwrite=TRUE)
      print(outfiles[i])
    }

    #creates a raster stack with all environmental rasters in this file
    env_data <- list.files(path="F:/SDM/DATA/temp", pattern = 'tif$')
    env_data<- stack(env_data)

For my loop I get the following Error message:

Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘mask’ for signature ‘"RasterLayer", "Extent"’

When I comment the mask part, the loop works, but the extents are still different.

I am also not sure if I choose the right extent, since I don't know the minimum extent. I spent a lot of time searching a way to handle this multiple-extent-problem within my data and found the crop() and mask() functions.

  • 1
    There is a funciton "align_rasters" in the gdalUtils package that makes this a bit easier. – Jeffrey Evans Nov 8 '16 at 21:36
4

You can manage multi-extent-problem resampling your data before mask() function. This work for aligned and non-aligned pixels (for non-aligned, choose wisely method argument). Also, you can use extendto align boundaries of your data. I'd made an reproducible example for recreate our problem:

library(raster)

a <- raster(xmn=-100, xmx=100, ymn=-90, ymx=90,res=10) 

# for reproducible example
r1 <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90,res=10)
r2 <- raster(xmn=-190, xmx=180, ymn=-80, ymx=10,res=10)
r3 <- raster(xmn=-50, xmx=130, ymn=-80, ymx=90,res=10)

a <- setValues(a, 1:ncell(a))
r1 <- setValues(r1, 1:ncell(r1))
r2 <- setValues(r2, 1:ncell(r2))
r3 <- setValues(r3, 1:ncell(r3))

files <- list()
files[[1]] <- r1
files[[2]] <- r2
files[[3]] <- r3

results <- list()

for(i in 1:length(files)) {
  e <- extent(a)
  r <-files[[i]] # raster(files[i])
  rc <- crop(r, e)
  if(sum(as.matrix(extent(rc))!=as.matrix(e)) == 0){ # edited
    rc <- mask(rc, a) # You can't mask with extent, only with a Raster layer, RStack or RBrick
  }else{
    rc <- extend(rc,a)
    rc<- mask(rc, a)
  }

  # commented for reproducible example      
  results[[i]] <- rc # rw <- writeRaster(rc, outfiles[i], overwrite=TRUE)
  # print(outfiles[i])

}

env_data<- stack(results)

The result will depends of the nature of you data, in this case I exaggerated the problem. After this code, we got this:

plot(env_data)

enter image description here


Edit:

Test your layers

library(raster)

a <- raster(xmn=-100, xmx=100, ymn=-90, ymx=90,res=10) 

# for reproducible example
r1 <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90,res=10)
r2 <- raster(xmn=-190, xmx=180, ymn=-80, ymx=10,res=10)
r3 <- raster(xmn=-50, xmx=130, ymn=-80, ymx=90,res=10)
# Other posibilities
r4 <- raster(xmn=-200, xmx=-150, ymn=-80, ymx=90,res=10)
r5 <- raster(xmn=-100, xmx=100, ymn=-90, ymx=90,res=5)
r6 <- a

files <- list()
files[[1]] <- r1
files[[2]] <- r2
files[[3]] <- r3
files[[4]] <- r4
files[[5]] <- r5
files[[6]] <- r6

library(testthat)

test <- list()

for(i in 1:length(files)){
  test[[i]] <- capture_warnings(compareRaster(a,files[[i]], res=T, orig=T, stopiffalse=F, showwarning=T))
}

test

[[1]]
[1] "different extent"            "different number or columns"

[[2]]
[1] "different extent"            "different number or columns"
[3] "different number or rows"   

[[3]]
[1] "different extent"            "different number or columns"
[3] "different number or rows"   

[[4]]
[1] "different extent"            "different number or columns"
[3] "different number or rows"   

[[5]]
[1] "different number or columns" "different number or rows"   
[3] "different resolution"       

[[6]]
character(0)
  • Raster stacks not working with different extents is not entirely true, it depends on what you want to do. You can create a stack with mismatched extents by using the "quick=TRUE" argument. – Jeffrey Evans Nov 8 '16 at 21:32
  • Thank you for ansering! The raster cells are aligned but don't match at the borders. What I plan to do is to extract values to points. – parallax Nov 8 '16 at 21:43
  • Actually, this code apply to your purposes. I use it in Landsat 8 time series, especially in subsets near to borders. – aldo_tapia Nov 9 '16 at 1:38
  • I'll try this code and give a feedback soon, Thank you aldo_tapia! – parallax Nov 9 '16 at 19:16
  • I tried the code but I am not sure how to integrate it in my code above. I am a newbie in R and coding in general, so could you make some more comments? – parallax Nov 10 '16 at 22:03
0

Crop can write straight to file; I don't think the mask step is necessary? you should make sure your extent aligns with pixel edges though - see here Improving Clip Raster workflow in QGIS?

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.