rr-rastergeotiff

Using mosaic in r to merge multiple geotiff


I have 50 geotiff files in the same folder. All of them represent elevation data in different parts of the world. I would like to merge certain geotiff files, and I found mosaic in R might help us. I have moved those geotiff into the same folder, and I wrote a R script show below:

setwd()
a<-list.files(pattern="*.tiff",file.name=TRUE)
combind<-merge(a,fun=mean)

However, this script returned an error: error in as.data.frame(y)

May I ask how could I improve my script?


Solution

  • Edit 2023-06: use terra::vrt() or terra::mosaic()

    Since gdalUtils seems unmaintained for a while now, you can use a similar approach as in my original answer just using the more up-to-date terra package which is the successor to the raster package:

    library(terra)
    
    vrt(
      x = list.files(path = "folder/to/images", pattern = "*.tif$", full.names = TRUE), 
      filename = "dem.vrt"
    )
    # afterwards read it as if it was a normal raster:
    dem <- rast("dem.vrt")
    

    Another possibility is to use terra::mosaic(). This gives you more control over the aggregation function in overlapping pixels. For that you first create a raster collection using sprc() from a list of rast objects and then mosaic to combine them into one raster.

    library(terra)
    # vector of file names
    fls <- list.files("your/folder", ".tif$", full.names = TRUE)
    # list of rast objects
    r_lst <- lapply(fls, rast)
    # create spatial raster collection
    coll <- sprc(r_list)
    
    # combine all rasters
    mosaic(coll, function = "mean")
    

    Original answer using gdalUtils

    You can make use of the powerful GDAL functions. From my experience these are much faster than pure R code.

    My approach would be with library(gdalUtils):

    First, build a virtual raster file (vrt):

    library(gdalUtils)
    setwd(...)
    gdalbuildvrt(gdalfile = "*.tif", # uses all tiffs in the current folder
                 output.vrt = "dem.vrt")
    

    Then, copy the virtual raster to a actual physical file:

    gdal_translate(src_dataset = "dem.vrt", 
                   dst_dataset = "dem.tif", 
                   output_Raster = TRUE # returns the raster as Raster*Object
                                        # if TRUE, you should consider to assign 
                                        # the whole function to an object like dem <- gddal_tr..
                   options = c("BIGTIFF=YES", "COMPRESSION=LZW"))
    

    Another pure (and probably slower) raster package solution would be:

    f <- list.files(path = "your/path", pattern = ".tif$", full.names = TRUE)
    rl <- lapply(f, raster)
    
    do.call(merge, c(rl, tolerance = 1))
    

    you have to adjust the tolerance since the raster files will probably not have the same origin.