5

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.

enter image description here

  • 1
    Please explain what you mean by merge. raster::merge is a mosaicking tool – Loïc Dutrieux Jan 17 '17 at 15:00
  • 1
    By running lapply(lst, origin) you'll be able to see if the difference in origin detected by compareRaster() comes from a minor precision issue which you'll be able to ignore or if your rasters really do have different origins – Loïc Dutrieux Jan 18 '17 at 10:39
  • there are a number of R wrappers for gdal. You can use gdalbuildvrt to build a virtual raster, then use wrapped gdal dataset methods to read the data. – Llaves Aug 11 '20 at 22:26
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) 
r_merged2<- mosaic(r1,r2, 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)
1

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)
  • Sorry @dmci, actually the origins are not the same, see edits. – Quechua Jan 18 '17 at 8:24
  • If the extents are different, then I would suggest checking these solutions – dmci Dec 6 '17 at 14:39
1

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

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.