Mosaic, or merge, rasters in R

Please note: For some reason WordPress is garbling up this post. The full unannotated script for this is presented at the very bottom.

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) ?writeRaster

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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s