Update Feb 5, 2022: The post showing you how to layer rasters on top of each other so that you can extract values is live. You can find it here.
Update Oct 13, 2020: This post has gained a lot of traction in recent weeks. I do want to point out a misconception about what this post is about. This post shows you how to merge rasters. This is not the same as layering rasters on top of each other so that you can extract values. I will how to do that at another time, since there seems to be a lot of interest around doing that.
Please note: People tend to have trouble on this particular post, and have asked me for help either through the comments or the contact me form. As with most things coding, you’re going to need to provide the data and code you use. As you can see in the comments below, plenty of people don’t. Something to keep in mind before contacting me for help.
To do a mosaic of rasters in GIS software, whether Arc or QGIS, is simple enough command-wise. However, I was running into issues about a month ago getting my rasters to mosaic properly. Something always seemed to come up wrong, so I thought I would do it in R to see how my results would differ. I’m glad I did, and so I thought I’d put this up in case anyone else out there is having the same trouble I did. Just as a note, I think the likelihood of an error increases depending on the size of your raster. Since I was using 1/3 arcsecond, i.e. 10 m spatial resolution, the rasters were fairly large.
All of my DEM data was downloaded from Earthexplorer. This mosaic would eventually become the basis for my watershed delineation, as seen in my portfolio.
Now onto the code.
I used the libraries rgdal, gdalutils, and raster. If they’re not present, use the install.packages() function. By the way, I’ve just discovered the package, easypackages. This package will allow you to load multiple packages with just one line. I recommend it.
library(easypackages)
libraries("rgdal", "gdalUtils", "raster")
So in the next line of code, I’m loading my rasters into R. All of my rasters, in this case I’m working with geotiffs, are all in the same folder for ease of loading. I’ve set my working directory to this folder.
setwd("D:/Rasters")
a <- c('1.tif', '2.tif', '3.tif', '4.tif')
The following is perhaps the most annoying part of mosaicing rasters in R. You need to create a template raster. This is basically the blank canvas that R will use to piece all of your rasters together. Think of it like a blank piece of paper. You’ll need the coordinates for the four corners of the blank piece of paper. The coordinates will be entered in the following format: x-minimum, x-maximum, y-minimum, y-maximum. You’ll use the extent() function to set the extent of your template raster.
e <- extent(-85, -83, 39, 41)
template <- raster(e)
As with all things geospatial, you need to set your projection. The projection for this particular raster dataset is NAD83, or espg 4269. If you’ve read my other tutorials, you know that you’ll be using the proj4string() and CRS() functions and that epsg information can be obtained here.
proj4string(template) <- CRS('+init=epsg:4269')
writeRaster(template, file="MiamiWatershed.tif", format="GTiff")
mosaic_rasters(gdalfile=a,dst_dataset="MiamiWatershed.tif",of="GTiff")
gdalinfo("MiamiWatershed.tif")
And you’re all set. Your mosaic raster should now be available in your working directory. Honestly, this will probably be my go to for mosaicing rasters. It seems to be faster than using Arc or Q.
setwd("D:/Raster")
a <- c('1.tif', '2.tif','3.tif','4.tif')
e <- extent(-85, -83, 39, 41)
template <- raster(e)
proj4string(template) <- CRS("+init=epsg:4269")
writeRaster(template, file="MiamiWatershed.tif", format="GTiff")
mosaic_rasters(gdalfile=a,dst_dataset="MiamiWatershed.tif",of="GTiff")
gdalinfo("MiamiWatershed.tif")
The featured image is Mosaic by Jenny Brown. It is used under a Creative Commons 2.0 License.
Teng Keng Vang 
Hi there, I jist Tried your method of mosaicing. However I keep getting an error that reads: “in local(x, filename) all cell values are NA” what does that mean.
Could you provide the code that you used?
What does error on “function X invalif file argument mean”?
This is a generic error that could be applicable to a wide range of problems. If you could provide the code that you used and where the error occurred, I would have a better idea of what went wrong.
Hi I get this error. My guess is that list of rasters is wrong? the “a”
mosaic_rasters(gdalfile=a,dst_dataset=”BigRaster.tif”,of=”GTiff”)
Error in FUN(X[[i]], …) : invalid ‘file’ argument
This is my code. naip_ndvi and naip_evi are my rasters. Img is a polygon I used to crop the original raster file before creating naip_ndvi and naip_evi.
a <- c(naip_ndvi, naip_evi)
img
extent(img)
e <- extent(img)
template <- raster(e)
proj4string(template) <- CRS('+proj=utm +zone=54 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
writeRaster(template, file="BigRaster.tif", format="GTiff", overwrite=TRUE)
mosaic_rasters(gdalfile=a,dst_dataset="BigRaster.tif",of="GTiff")
I hope you can help me!
cheers and thanks for your code.
The extent needs to be set to the rasters you’re trying to mosaic. You can mask the raster after you’ve mosaic’d them.
Thanks for this, very helpful. When running the following code I get the below error.
mosaic_rasters(gdalfile=Viet_NDVI,dst_dataset= “Vietnam_NDVI.tif”, of=”GTiff”)
“Error in mosaic_rasters(gdalfile = Viet_NDVI, dst_dataset = “Vietnam_NDVI.tif”, : Some of the input files are missing: ndvi32-00.tif Some of the input files are missing: ndvi00-00.tif”
I’m pretty sure I followed everything exactly, please let me know how to fix this.
Its hard to tell whats wrong without seeing all the code. You’d have to provide all of the code. Its likely your mistake occurred before this line. Its just that it didn’t become apparent until you tried to run this function.
Hi tK,
Thanks for this post! I run the following script to mosaic 15 .tif files into a new raster. The new raster is created, but only contains data from the very first filer ‘OF.tif’ from object ‘a’ (I can see when plotted and with GDALinfo. Another problem I see is the message saying CRS was discarded.
Below my code:
a <- c('OF.tif', 'OG.tif', 'OH.tif', 'OI.tif', 'OJ.tif', 'PF.tif', 'PG.tif', 'PH.tif', 'PI.tif',
'QF.tif', 'QG.tif', 'QH.tif', 'QI.tif', 'RF.tif', 'RG.tif')
# a large raster was clipped to (42.448°S 170.695°E)-(44.379°S 173.496°E) to export the 15 rasters
e <- extent(170.85, 173.34, -44.40, -42.55)
template <- raster(e)
proj4string(template) <- CRS("+init=epsg:4959") #.tif files are in NZGD 2000/ epgs: 2193
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
Discarded datum New_Zealand_Geodetic_Datum_2000 in CRS definition
writeRaster(template, file=ʺraster_4_farms.tifʺ, format=ʺGTiffʺ, overwrite=TRUE)
Warning messages:
1: In .local(x, filename, …) : all cell values are NA
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
mosaic_rasters(gdalfile=a,dst_dataset=ʺraster_4_farms.tifʺ,of=ʺGTiffʺ)
# list from [1] to [15] with warning:
[15] "C:\\Users\\morenogc\\AppData\\Local\\Temp\\RtmpigJ8o5\\file2a81e2cf25_warped.vrt"
NULL
There were 15 warnings (use warnings() to see them)
gdalinfo(ʺraster_4_farms.tifʺ)
[1] "Driver: GTiff/GeoTIFF"
[2] "Files: raster_4_farms.tif"
[3] "Size is 3922, 7188"
[4] "Origin = (1410416.000000000000000,5300384.000000000000000)"
[5] "Pixel Size = (8.000000000000000,-8.000000000000000)"
[6] "Image Structure Metadata:"
[7] " INTERLEAVE=BAND"
[8] "Corner Coordinates:"
[9] "Upper Left ( 1410416.000, 5300384.000) "
[10] "Lower Left ( 1410416.000, 5242880.000) "
[11] "Upper Right ( 1441792.000, 5300384.000) "
[12] "Lower Right ( 1441792.000, 5242880.000) "
[13] "Center ( 1426104.000, 5271632.000) "
[14] "Band 1 Block=3922×1 Type=Float32, ColorInterp=Gray"
[15] " NoData Value=-32767"
Off the top of my head, it looks like a projection issue. It looks like your raster extent was set to WGS84 (EPSG:4326) whereas your rasters are in a different projection. You need to set the extent to the projection that your rasters are in. You can then do a transformation if you wanted it in a different projection.
Hey tK,
Thanks. I uninstalled and reinstalled the three packages. Thenk, I set the extend with the projection of the rasters:
e gdalinfo(ʺraster_4_farms.tifʺ)
[1] “Driver: GTiff/GeoTIFF” “Files: raster_4_farms.tif”
[3] “Size is 3922, 7188” “Origin = (1410416.000000000000000,5300384.000000000000000)”
[5] “Pixel Size = (8.000000000000000,-8.000000000000000)” “Image Structure Metadata:”
[7] ” INTERLEAVE=BAND” “Corner Coordinates:”
[9] “Upper Left ( 1410416.000, 5300384.000) ” “Lower Left ( 1410416.000, 5242880.000) ”
[11] “Upper Right ( 1441792.000, 5300384.000) ” “Lower Right ( 1441792.000, 5242880.000) ”
[13] “Center ( 1426104.000, 5271632.000) ” “Band 1 Block=3922×1 Type=Float32, ColorInterp=Gray”
[15] ” NoData Value=-32767″
Hmm, I’m not sure why WordPress didn’t notify me about your next comment. Apologies about that. I’m not sure why you reinstalled the packages. You’ll also need to provide the whole script that you used. I’m not really sure what’s going on with the snippets you provided.
Hello, thank you for your code. I was trying to use it on rasters from jp2 file format (those are sentinel 2 data) and get the following error.
Error in FUN(X[[i]], …) : invalid ‘file’ of argument
UVV <- list.files(path ="C:", pattern="*.jp2")
s_UVV <- stack(UVV)
#UUV
UUV <- list.files(path = "C" pattern="*.jp2")
s_UUV <- stack(UUV)
#merge tiles T33UUv und T33UVV
a <- c(s_UVV, s_UUV)
e <- extent(3e+05, 509760, 5890200, 6e+06)
template <- raster(e)
proj4string(template) <- '+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs'
setwd("C:/Users/bovn_551/Documents/Masterarbeit/Sentinel2/2020_05_31")
writeRaster(template, file="MiamiWatershed.tif", format="GTiff")
mosaic_rasters(gdalfile=a,dst_dataset="MiamiWatershed.tif",of="GTiff")
gdalinfo("MiamiWatershed.tif")
The two rasters (s_UUV,s_UVV) look the same as a stacked .tif. So, I think there is no error. But now I was wondering if maybe the mosaic_raster doesn’t support jp2 file formats. Or is the problem that the input is a jp2 and the output is a tif.
Do you by any chance know a solution to this problem. Thanks in advance.
Ignore the other comment. I meant to reply to another comment. Out of the box, gdal doesn’t recognize sentinel-2 images. You need to redirect gdal to the correct drivers. Try the following line of code after importing your libraries:
>gdal_chooseInstallation(“JP2OpenJPEG”)
tried the gdal_chooseInstallation. But when I run the code, I still get the same error as bevor, invalid ‘file’ of argument. I also looked into gdalDrivers() and for the JP2OpenJPEG I get for create=FALSE, copy=TRUE, isRASTER=TRUE. This is the same as prior to gdal_chooseInstallation. So, it seems this is not working. Do you have any other suggestions. I am using Rversion 4.0.2 and the rgdal version is 1.1-17 (also doesn’t change after the gdal_chooseInstallation).
A quick look shows that this seems to be a common error. Unfortunately, since every install and computer configuration is a little bit different, you’ll need to browse stackexchange for possible solutions. I had to do the same when working with MrSID files. I wish you the best in your search.
Hi! Thank you very much for this very helpful post. I am trying to mosaic 4 raster stacks with 500 layers each. I tried to use your method but I am encountering errors. My intuition tells me that there are some extra steps that I am missing due to the fact that I am using stacks. Would really appreciate any insights you can give me.
> a b c d
> abcd
> e template
> ref_crs
> proj4string(template) writeRaster(template, file=”MdC_NDMI_1990_2019″, format = “GTiff”)
Warning message:
In .local(x, filename, …) : all cell values are NA
> mosaic_rasters(gdalfile = abcd, dst_dataset = “MdC_NDMI_1990_2019.tif”, of = “GTiff”)
Error in FUN(X[[i]], …) : invalid ‘file’ argument
Try a raster brick instead of stack. The documentation specifies that the mosaic function only works with a raster layer or brick.