0

I am trying to stack multiple rasters together in one RasterStack object but am having problems with aligning the extents and dimensions/resolution so that it can work. All rasters have the same CRS and projections.

The four raster files I am using can be downloaded from this dropbox link: https://www.dropbox.com/sh/1fafynlmoveqr2z/AACdQlHMzwRgfdnKKxP8mcm_a?dl=0

I have tried using the following, but I am not sure if this is the best/correct method:

    library(raster)
    library(sp)
    library(rgeos)
    library(rgdal)
    library(dismo)

    # Load country administrative boundaries:
    data("wrld_simpl")

    r1 <- raster("./NDVI_metu_saudi_arabia.2021001_NDVI.MOD44CQ.250m.tif")
    r2 <- raster("./slope_1KMmd_SRTM.tif")
    r3 <- raster("./wc2.1_30s_bio_1.tif")

    # Mask raster with country border shapefile
    mskd_rstr <- mask(
      r1, 
      wrld_simpl[177,]  # country border
    )
    
    # Trim raster to country boundary. Make all pixels outside border == NA.
    trmmd_rstr <- trim(mskd_rstr)
    
    # Check result.
    plot(trmmd_rstr)
    
    # Write to disk.
    writeRaster(
      trmmd_rstr, 
      file = paste0(
        "r1", 
        "_SAU_trmmd_rstr.tif"
      )
    )

# Repeat for r2, r3, r4 so that all rasters are the same and can be stacked.

This method takes extremely long (especially when using trim()). Is there a better alternative?

Ultimately, I want to be able to stack them together without getting errors about the extent or dimensions. For example by using:

predictors <- stack(r1, r2, r3, r3)

1 Answer 1

0

you can align the global raster to the extent and resolution of the NDVI raster and stack them together. Once you have the stack you could mask and trim in one go all the layers. To align the extent and resolution I recommend gdalwarp from gdalUtils package. I've written some code that can help you, with the data you provided:

library(raster)
library(rgdal)
library(gdalUtils)

setwd("your/path/to/raster/file")
#path of the raster to be aligned
input <- list.files(full.names = T,pattern = 'tif')[-1]#left out the NDVI raster
#base raster that give extent and resolution
base <- raster("./NDVI_metu_saudi_arabia.2021001_NDVI.MOD44CQ.250m.tif")
ext <- extent(base)
ext <- c(ext@xmin, ext@ymin, ext@xmax, ext@ymax)
r <-c( 0.002248181, 0.002248301)#the resolution of the NDVI raster
#name of the output. They will be saved in the getwd() directory
output <- paste0(basename(tools::file_path_sans_ext(input)),"_allign.tif")
#allign

myallign <- function(input,output,res,ext){
mask_out <- gdalwarp(
srcfile = input,
dstfile = output,
t_srs=proj4string(base),
tr = res,
te = ext,
output_Raster = T,
co = "compress=DEFLATE",
overwrite = TRUE,
verbose = TRUE
)
}

for (i in seq_along(input)) {
alligned <- myallign(input[i],output[i],res=r,ext = ext)
}#40 sec on my PC
#now you can easily stack the file
s <- stack(output)
s <- stack(base,s)

Now you can mask and trim the s stack to obtain the final output. I warn you that probably these last operations will be the slowest one, but at least you have to perform them one time. There are probably a lot of ways to speed up computing in such a task. One could be to implement your code in parallel with foreach or ParLapply. My favourite function for parallel processing is sfLapply form snowfall package, it is very effective and fast. The gdal solution is still quite fast and somehow elegant

2
  • Thanks Elia, this worked well except for that I was encountering an error with the gdalUtils package: ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db Googling this error, I found other people having the same issue with gdalUtils. In order to get around it, I used the gdal software in terminal instead. Code in next message for anyone who experience the same issue.
    – rhyncogale
    Apr 21, 2021 at 13:08
  • "/usr/local/Cellar/gdal/3.2.2_3/bin/gdalwarp" -overwrite -te 33.44490261 12.211649492 62.747693178 32.444114516 -tr 0.00224818095504066 0.00224830148060896 -t_srs "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" -of "GTiff" -co "compress=DEFLATE" "./elevation_1KMmd_SRTM.tif" "elevation_1KMmd_SRTM_allign.tif"
    – rhyncogale
    Apr 21, 2021 at 13:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.