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
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)
| improve this answer | |
  • 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 cols 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(ref_rast),crs(rast))){
    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
| improve this answer | |

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.