Create Raster Image Mosaic

When working with raster data such as satellite imagery, it is sometimes necessary to combine multiple adjacent, possibly overlapping smaller rasters into a single large raster that covers the entire area of interest. This example demonstrates one way to use R to create a raster image mosaic while also applying a resampling algorithm to align the inputs. The approach described below relies on the raster package, which provides extensive capabilities for working with raster data in R.

Overview

In preparation for creating a raster mosaic, there are two key questions to consider:

  1. Are the input rasters aligned on a common grid? That is, if the rasters were extended indefinitely in all directions, would their respective gridlines all match up perfectly?
  2. Are the input rasters adjacent but non-overlapping?

If the answer to both questions is yes, the mosaic operation simply involves piecing the input rasters together like floor tiles. Indeed, for ease of storage and distribution, data providers often deliberately cut global raster datasets into smaller tiles that users can then mosaic back together.

If the inputs are aligned but have overlaps, one must decide how to choose pixel values in the areas of overlap. The simplest rule is to specify a prioritization order for the input rasters. For example, the raster package's merge function behaves much like the mosaic function demonstrated below, but preferentially uses values from rasters that appear first in the function call. A more interesting rule is to apply some pixel-by-pixel comparison function (e.g., always take the maximum value from the overlapping pixels) or aggregation function (e.g., take the mean of the overlapping pixels). Even more sophisticated methods may blend pixels in a way that smooths out transitions between the input rasters, but these are beyond the scope of this example.

Finally, if the inputs are not aligned, they must first be resampled (and possibly reprojected) into a common grid. Although the exact behavior of raster mosaic operations varies across software applications, in general one would consider the resampling step to be a precursor to the actual mosaic operation itself, which fundamentally involves concatenating pre-aligned rasters into a single larger raster, while applying some rule for choosing between overlapping pixels as described above.




Example: Mosaic nine Landsat TM images spanning Ohio

Here we have nine different overlapping raster images derived from Landsat Thematic Mapper tiles (band 5, near-infrared) that together cover the state of Ohio. As a consequence of prior processing, these rasters each have slightly different pixel resolution, and are slightly misaligned with respect to one another. Before performing the mosaic step itself, we thus need to define a common grid and use this to resample each input raster. The raster package resample function can perform resampling using either a nearest neighbor or bilinear method; of these, the latter is more appropriate for continuous numeric data as we have here.

Input rasters

R code and raster output map

require(raster)
 
# read the individual rasters into a list of RasterLayer objects
input.rasters <- lapply(list.files(pattern="^TmB50.*[.]tif$"), raster)
 
# create an empty output raster that spans the full extent of all input
# rasters, and uses the same coordinate reference system; in this case
# we know that all input rasters share the same CRS, so we can
# arbitrarily extract CRS information from the first one
full.extent <- unionExtent(input.rasters)
bounding.raster <- raster(full.extent,
    crs=projection(input.rasters[[1]]))
 
# set the output resolution to match the center tile (somewhat
# arbitrarily); this can also be specified manually if preferred
res(bounding.raster) <- res(input.rasters[[5]])
 
# for each input raster, extract the corresponding sub-extent of the
# empty output raster, and use this as the basis for resampling
# !! note that this may take several minutes to run !!
resampled.rasters <- lapply(input.rasters, function(input.raster) {
    target.raster <- crop(bounding.raster, input.raster)
    resample(input.raster, target.raster, method="bilinear")
})
 
# mosaic all 9 resampled rasters, taking the maximum pixel value in
# cases where there is overlap
raster.mosaic <- mosaic(resampled.rasters, fun=max)
 
# read in Ohio state border shapefile and generate output map
ohio <- readOGR(".", "Ohio")
plot(raster.mosaic, col=grey((0:256)/256))
plot(ohio, border="yellow", add=TRUE)
 
# write the result out to a GeoTIFF, rounding pixel values to the
# nearest integer to preserve the original 1-byte integer format
writeRaster(round(raster.mosaic), filename="ohio-mosaic.tif",
    datatype="INT1U")

Ohio mosaic

Post-Processing Suggestion:The component images are discernible in the mosaic due to differences in the sunlight illumination level on the acquisition dates for each component image. This effect can be minimized, and the mosaic given a more uniform appearance, by applying histogram matching to the mosaic. Third-party image processing software, such as OSSIM, can be used for this.

Download R script and data

Learning More:

scicomp: