Merging multiple rasters in R

Question

I've been trying to find a time-efficient way to merge multiple raster images in R. These are adjacent ASTER scenes from the southern Kilimanjaro region, and my target is to put them together to obtain one large image.

This is what I got so far (object 'ast14dmo' representing a list of RasterLayer objects):

# Loop through single ASTER scenes
for (i in seq(ast14dmo.sd)) {
  if (i == 1) {
    # Merge current with subsequent scene
    ast14dmo.sd.mrg <- merge(ast14dmo.sd[[i]], ast14dmo.sd[[i+1]], tolerance = 1)
  } else if (i > 1 && i < length(ast14dmo.sd)) {
    tmp.mrg <- merge(ast14dmo.sd[[i]], ast14dmo.sd[[i+1]], tolerance = 1)
    ast14dmo.sd.mrg <- merge(ast14dmo.sd.mrg, tmp.mrg, tolerance = 1)
  } else {
    # Save merged image
    writeRaster(ast14dmo.sd.mrg, paste(path.mrg, "/AST14DMO_sd_", z, "m_mrg", sep = ""), format = "GTiff", overwrite = TRUE)
  }
}

As you surely guess, the code works. However, merging takes quite long considering that each single raster object is some 70 mb large. I also tried Reduce and do.call, but that failed since I couldn't pass the argument 'tolerance' which circumvents the different origins of the raster files.

Anybody got an idea of how to speed things up?

Komsa

Answer 1

With "terra" you first create a SpatRasterCollection

library(terra)
r1 <- rast(xmax=-150, ymin=60, ncols=30, nrows=30)
values(r1) <- 1:ncell(r1)
r2 <- rast(xmin=-100, xmax=-50, ymax=50, ymin=30)
res(r2) <- c(xres(r1), yres(r1))
values(r2) <- 1:ncell(r2)

s <- sprc(r1, r2)
m <- merge(s)

---

With the deprecated "raster" package.

You can use do.call

ast14dmo.sd$tolerance <- 1
ast14dmo.sd$filename <- paste(path.mrg, "/AST14DMO_sd_", z, "m_mrg.tif", sep = "")
ast14dmo.sd$overwrite <- TRUE
mm <- do.call(merge, ast14dmo.sd)

Here with some data, from the example in raster::merge

r1 <- raster(xmx=-150, ymn=60, ncols=30, nrows=30)
r1[] <- 1:ncell(r1)
r2 <- raster(xmn=-100, xmx=-50, ymx=50, ymn=30)
res(r2) <- c(xres(r1), yres(r1))
r2[] <- 1:ncell(r2)

x <- list(r1, r2)
names(x) <- c("x", "y")
x$filename <- 'test.tif'
x$overwrite <- TRUE
m <- do.call(merge, x)

Boathouse

Great solution as well, thank you! I just had a quick look at the computation time and it turned out that your approach via do.call works almost twice as fast as Reduce. – Komsa

For some reason this doesn't work with a named list wich can be troublesome – Ionia

In that case, you need to add this line, I think, names(x) <- c("x", "y") – Boathouse

Error while trying to “do.call”: Error in as.data.frame(x) : argument "x" is missing, with no default. I’ve tried raster::merge instead of merge, but got the same error. – Jura

With raster being obsolete and the newer terra package out, would it make sense to use terra::merge in replacement of raster::merge? or is there a better solution as of today? – Bureaucrat

Answer 2

The 'merge' function from the Raster package is a little slow. For large projects a faster option is to work with gdal commands in R.

library(gdalUtils)
library(rgdal)

Build list of all raster files you want to join (in your current working directory).

all_my_rasts <- c('r1.tif', 'r2.tif', 'r3.tif')

Make a template raster file to build onto. Think of this a big blank canvas to add tiles to.

e <- extent(-131, -124, 49, 53)
template <- raster(e)
projection(template) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
writeRaster(template, file="MyBigNastyRasty.tif", format="GTiff")

Merge all raster tiles into one big raster.

mosaic_rasters(gdalfile=all_my_rasts,dst_dataset="MyBigNastyRasty.tif",of="GTiff")
gdalinfo("MyBigNastyRasty.tif")

This should work pretty well for speed (faster than merge in the raster package), but if you have thousands of tiles you might even want to look into building a vrt first.

Pushcart

Isn't there a more efficient way to set the extent without guessing lat/long? – Kokand

@HermanToothrot: without looking at each individual file, to know the final extent of the merged files you either need to make a (generous) guess (educated e.g. by opening border files in GIS). To get the exact extents of the final merged file, you can do a "fast" merge of the extents before doing the real merge. FINAL_EXTENT = Reduce(function(d1, d2) merge(extent(raster(d1)), extent(raster(d2))), RASTER_LIST) – Checkerwork

one downside of mosaic_rasters seems to be that no overlap behaviour can be defined. the behaviour seems to be to keep ALL values of the first raster, even if they are NA values. – Checkerwork

Answer 3

You can use Reduce like this for example :

Reduce(function(...)merge(...,tolerance=1),ast14dmo.sd)

Samoyed

Answer 4

SAGA GIS mosaicking tool (http://www.saga-gis.org/saga_tool_doc/7.3.0/grid_tools_3.html) gives you maximum flexibility for merging numeric layers, and it runs in parallel by default! You only have to translate all rasters/images to SAGA .sgrd format first, then run the command line saga_cmd.

Immediacy

Answer 5

I have tested the solution using gdalUtils as proposed by Matthew Bayly. It works quite well and fast (I have about 1000 images to merge). However, after checking with document of mosaic_raster function here, I found that it works without making a template raster before mosaic the images. I pasted the example codes from the document below:

outdir <- tempdir()
gdal_setInstallation()
valid_install <- !is.null(getOption("gdalUtils_gdalPath"))
if(require(raster) && require(rgdal) && valid_install)
{
layer1 <- system.file("external/tahoe_lidar_bareearth.tif", package="gdalUtils")
layer2 <- system.file("external/tahoe_lidar_highesthit.tif", package="gdalUtils")
mosaic_rasters(gdalfile=c(layer1,layer2),dst_dataset=file.path(outdir,"test_mosaic.envi"),
    separate=TRUE,of="ENVI",verbose=TRUE)
gdalinfo("test_mosaic.envi")

}

Cursorial

Answer 6

I was faced with this same problem and I used

#Read desired files into R
data_name1<-'file_name1.tif' 

r1=raster(data_name1)

data_name2<-'file_name2.tif'

r2=raster(data_name2)

#Merge files
new_data <- raster::merge(r1, r2)

Although it did not produce a new merged raster file, it stored in the data environment and produced a merged map when plotted.

Irreplaceable

Answer 7

I ran into the following problem when trying to mosaic several rasters on top of each other

In vv[is.na(vv)] <- getValues(x[[i]])[is.na(vv)] :
  number of items to replace is not a multiple of replacement length 

As @Robert Hijmans pointed out, it was likely because of misaligned rasters. To work around this, I had to resample the rasters first

library(raster)

x  <- raster("Base_raster.tif")
r1 <- raster("Top1_raster.tif")
r2 <- raster("Top2_raster.tif")

# Resample
x1 <- resample(r1, crop(x, r1))
x2 <- resample(r2, crop(x, r2))

# Merge rasters. Make sure to use the right order
m <- merge(merge(x1, x2), x)

# Write output
writeRaster(m,
            filename = file.path("Mosaic_raster.tif"),
            format = "GTiff",
            overwrite = TRUE)

Acinaciform