Merge rasters with different origins in R

Question

I have a list of rasters with the same projection and resolution but different origins and I need to merge them in a single raster file.

Here's my code

l <- list.files('path', full.names=TRUE)
lst <- lapply(l, raster)
r_merged <- do.call(raster::merge, lst)

Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  
 : different origin

Here's what I'd like to get. Blue squares are original rasters, I want to get a single raster as the red square.

Answer 1

There was an issue with projectRaster's alignOnly that was fixed in raster 3.4-8 on github on Dec 22nd, 2020. To align and merge two rasters (r1 & r2) this code should work.

template<- projectRaster(from = r2, to= r1, alignOnly=TRUE)
#template is an empty raster that has the projected extent of r2 but is aligned with r1 (i.e. same resolution, origin, and crs of r1)
r2_aligned<- projectRaster(from = r2, to= template)
r_merged<- merge(r1,r2_aligned) 
r_merged2<- mosaic(r1,r2_aligned, fun=mean, na.rm=TRUE)

In merge function if objects overlap, the values get priority in the same order as the arguments (but NA values are ignored), but in mosaic a function is applied to compute cell values in areas where layers overlap.

For combining more than 2 rasters you must first align them all using projectRaster, and then you can put them in a list and use do.call

#For merge
x <- list(r1, r2, r3)
m <- do.call(merge, x)

#For mosaic
x <- list(r1, r2, r3)
names(x)[1:2] <- c('x', 'y')
x$fun <- mean
x$na.rm <- TRUE

m <- do.call(mosaic, x)

Answer 2

I would suggest that you use the stack function from the raster package. For instance, as the example in the documentation shows:

# file with one layer
fn <- system.file("external/test.grd", package="raster")
s <- stack(fn, fn)
r <- raster(fn)
s <- stack(r, fn) 
nlayers(s)

Answer 3

I ran into a similar issue and developed this code to accomplish this task. The first function "reproject_align_raster" is called within "combine_rasters." combine_rasters takes a list of rasters and will align them all to a common grid and then merge them.

#' Reprojects/resamples and aligns a raster
#'
#' Reprojects/resamples and aligns a raster by matching a raster a raster to a specified origin, resolution, and coordinate reference system, or that of a reference raster. Useful for preparing adjacent areas before using raster::merge.
#' @param rast raster to be reprojected or resampled
#' @param ref_rast reference raster with desired properties (Alternatively can supply desired_origin, desired_res, and desired_crs)
#' @param desired_origin desired origin of output raster as a vector with length 2 (x,y)
#' @param desired_res  desired resolution of output raster. Either an integer or a vector of length 2 (x,y)
#' @param desired_crs desired coordinate reference system of output raster (CRS class)
#' @param method resampling method. Either "bilinear" for bilinear interpolation (the default), or "ngb" for using the nearest neighbor
#' @importFrom  raster crs
#' @importFrom  raster extent
#' @importFrom  raster origin
#' @importFrom  raster projectExtent
#' @importFrom raster raster
#' @importFrom raster resample
#' @importFrom raster projectRaster
#' @export

reproject_align_raster<- function(rast, ref_rast=NULL, desired_origin, desired_res, desired_crs, method= "bilinear"){

  if (!is.null(ref_rast)) {
    desired_origin<- origin(ref_rast) #Desired origin
    desired_res<- res(ref_rast) #Desired resolution
    desired_crs<- crs(ref_rast) #Desired crs
  } #Set parameters based on ref rast if it was supplied
  if(length(desired_res)==1){
    desired_res<- rep(desired_res,2)}

  if(identical(crs(rast), desired_crs) & identical(origin(rast), desired_origin) & identical(desired_res, res(rast))){
    message("raster was already aligned")
    return(rast)} #Raster already aligned

  if(identical(crs(rast), desired_crs)){
    rast_orig_extent<- extent(rast)} else{
      rast_orig_extent<- extent(projectExtent(object = rast, crs = desired_crs))} #reproject extent if crs is not the same
  var1<- floor((rast_orig_extent@xmin - desired_origin[1])/desired_res[1])
  new_xmin<-desired_origin[1]+ desired_res[1]*var1 #Calculate new minimum x value for extent
  var2<- floor((rast_orig_extent@ymin - desired_origin[2])/desired_res[2])
  new_ymin<-desired_origin[2]+ desired_res[2]*var2 #Calculate new minimum y value for extent
  n_cols<- ceiling((rast_orig_extent@xmax-new_xmin)/desired_res[1]) #number of cols to be in output raster
  n_rows<- ceiling((rast_orig_extent@ymax-new_ymin)/desired_res[2]) #number of rows to be in output raster
  new_xmax<- new_xmin+(n_cols*desired_res[1]) #Calculate new max x value for extent
  new_ymax<- new_ymin+(n_rows*desired_res[2]) #Calculate new max y value for extent
  rast_new_template<- raster(xmn=new_xmin, xmx =new_xmax,  ymn=new_ymin, ymx= new_ymax, res=desired_res, crs= desired_crs) #Create a blank template raster to fill with desired properties
  if(!identical(desired_origin,origin(rast_new_template))){
    message("desired origin does not match output origin")
    stop()} #Throw error if origin doesn't match
  if(identical(crs(rast),desired_crs)){
    rast_new<- raster::resample(x=rast, y=rast_new_template, method = method)} else{
      rast_new<- projectRaster(from=rast, to=rast_new_template, method = method)} #Use projectRaster if crs doesn't match and resample if they do
  if(!identical(desired_origin,origin(rast_new))){
    message("desired origin does not match output origin")
    stop()} #Throw error if origin doesn't match
  return(rast_new)
}

    #' Combines multiple rasters into one
    #'
    #' Combines multiple rasters into one by using merge after first reprojecting or resampling and aligning rasters by matching them up with a a specified origin, resolution, and coordinate reference system, or that of a reference raster.
    #' @param raster_list a list of rasters
    #' @param rast raster to be reprojected or resampled
    #' @param ref_rast reference raster with desired properties (Alternatively can supply desired_origin, desired_res, and desired_crs)
    #' @param desired_origin desired origin of output raster as a vector with length 2 (x,y)
    #' @param desired_res  desired resolution of output raster. Either an integer or a vector of length 2 (x,y)
    #' @param desired_crs desired coordinate reference system of output raster (CRS class)
    #' @param method resampling method. Either "bilinear" for bilinear interpolation (the default), or "ngb" for using the nearest neighbor
    #' @param display_progress Logical specifying whether or not to indicate progress
    #' @importFrom  raster merge
    #' @export
    
    combine_rasters<- function(raster_list, ref_rast=NULL, desired_origin, desired_res, desired_crs, method= "bilinear", display_progress=TRUE){
      raster_list2<- vector("list", length = length(raster_list))
      for (i in 1:length(raster_list)) {
        if(display_progress){
          print(paste("Reprojecting", as.character(i), "of", as.character(length(raster_list))))}
        raster_list2[[i]]<- reproject_align_raster(raster_list[[i]], ref_rast=ref_rast, desired_origin=desired_origin, desired_res=desired_res, desired_crs=desired_crs, method= method)}
      for (j in 1:length(raster_list2)) {
        if(display_progress){
          print(paste("Combining raster", as.character(j), "of", as.character(length(raster_list2))))}
        if(j==1){
          combined_raster<- raster_list2[[j]]} else{
            combined_raster<- raster::merge(combined_raster, raster_list2[[j]])}}
      return(combined_raster)
      }


  [1]: https://gis.stackexchange.com/questions/333464/re-projecting-a-raster-in-r-to-matching-resolution-crs-and-origin-but-differe/335605