7

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?

AndrewGB
  • 16,126
  • 5
  • 18
  • 49
Bing-Hong Huang
  • 125
  • 1
  • 9

1 Answers1

14

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.

loki
  • 9,816
  • 7
  • 56
  • 82
  • Thank you for kindly reply. The gdalUtils finally return warning message ""C:\Program Files (x86)\QGIS 2.18\bin\gdal_translate.exe" -of "GTiff" "dem.vrt" "elevation.tif"' had status 1. Therefore, I am trying the second code. – Bing-Hong Huang May 08 '18 at 16:04
  • Perhaps you have to set `options = c("BIGTIFF=YES")` – loki May 08 '18 at 19:06
  • I had the time to check. In my case `gdalbuildvrt(gdalfile = ".tif$" ...)` failed. I edited my answer. Perhapy you might check. Just substitute `".tif$"` with `".tif"`. I also added some options. – loki May 09 '18 at 07:14