5.3 Geometric operations on raster data

Geometric raster operations include the shift, flipping, mirroring, scaling, rotation or warping of images.These operations are necessary for a variety of applications including georeferencing, used to allow images to be overlaid on an accurate map with a known CRS (Liu and Mason 2009).A variety of georeferencing techniques exist, including:

• Georeferencing based on known ground control points.
• Orthorectification also georeferences an image, but additionally takes into account local topography.
• Image (co-)registration is the process of aligning one image with another (in terms of coordinate reference system, origin and resolution).Registration becomes necessary for images from the same scene but shot from different sensors or from different angles or at different points in time.
  R is unsuitable for the first two points since these often require manual intervention which is why they are usually done with the help of dedicated GIS software (see also Chapter 9).On the other hand, aligning several images is possible in R and this section shows among others how to do so.This often includes changing the extent, the resolution and the origin of an image.A matching projection is of course also required but is already covered in Section 6.6.In any case, there are other reasons to perform a geometric operation on a single raster image.For instance, in Chapter 13 we define metropolitan areas in Germany as 20 km² pixels with more than 500,000 inhabitants.The original inhabitant raster, however, has a resolution of 1 km² which is why we will decrease (aggregate) the resolution by a factor of 20 (see Section 13.5).Another reason for aggregating a raster is simply to decrease run-time or save disk space.Of course, this is only possible if the task at hand allows a coarser resolution.Sometimes a coarser resolution is sufficient for the task at hand.

5.3.1 Geometric intersections

In Section 4.3.1 we have shown how to extract values from a raster overlaid by other spatial objects.To retrieve a spatial output, we can use almost the same subsetting syntax.The only difference is that we have to make clear that we would like to keep the matrix structure by setting the drop-parameter to FALSE.This will return a raster object containing the cells whose midpoints overlap with clip.

data("elev", package = "spData")
clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
              res = 0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]
#> class      : RasterLayer 
#> dimensions : 2, 1, 2  (nrow, ncol, ncell)
#> resolution : 0.5, 0.5  (x, y)
#> extent     : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#> source     : memory
#> names      : layer 
#> values     : 18, 24  (min, max)

For the same operation we can also use the intersect() and crop() command.

5.3.2 Extent and origin

When merging or performing map algebra on rasters, their resolution, projection, origin and/or extent have to match. Otherwise, how should we add the values of one raster with a resolution of 0.2 decimal degrees to a second with a resolution of 1 decimal degree?The same problem arises when we would like to merge satellite imagery from different sensors with different projections and resolutions.We can deal with such mismatches by aligning the rasters.

In the simplest case, two images only differ with regard to their extent.Following code adds one row and two columns to each side of the raster while setting all new values to an elevation of 1000 meters (Figure 5.13).

data(elev, package = "spData")
elev_2 = extend(elev, c(1, 2), value = 1000)
plot(elev_2)

Figure 5.13: Original raster extended by one row on each side (top, bottom) and two columns on each side (right, left).

Performing an algebraic operation on two objects with differing extents in R, the raster package returns the result for the intersection, and says so in a warning.

elev_3 = elev + elev_2
#> Warning in elev + elev_2: Raster objects have different extents. Result for
#> their intersection is returned

However, we can also align the extent of two rasters with extend().Instead of telling the function how many rows or columns should be added (as done before), we allow it to figure it out by using another raster object.Here, we extend the elev object to the extent of elev_2.The newly added rows and column receive the default value of the value parameter, i.e., NA.

elev_4 = extend(elev, elev_2)

The origin is the point closest to (0, 0) if you moved towards it (starting from any given cell corner of the raster) in steps of x and y resolution.

origin(elev_4)
#> [1] 0 0

If two rasters have different origins, their cells do not overlap completely which would make map algebra impossible.To change the origin , use origin().²¹Looking at Figure 5.14 reveals the effect of changing the origin.

# change the origin
origin(elev_4) = c(0.25, 0.25)
plot(elev_4)
# and add the original raster
plot(elev, add = TRUE)

Figure 5.14: Rasters with identical values but different origins.

Note that changing the resolution frequently (next section) also changes the origin.

5.3.3 Aggregation and disaggregation

Raster datasets can also differ with regard to their resolution.To match resolutions, one can either decrease (aggregate()) or increase (disaggregate()) the resolution of one raster.²²As an example, we here change the spatial resolution of dem (found in the RQGIS package) by a factor of 5 (Figure 5.15).Additionally, the output cell value should correspond to the mean of the input cells (note that one could use other functions as well, such as median(), sum(), etc.):

data("dem", package = "RQGIS")
dem_agg = aggregate(dem, fact = 5, fun = mean)

Figure 5.15: Original raster (left). Aggregated raster (right).

By contrast, thedisaggregate() function increases the resolution.However, we have to specify a method on how to fill the new cells.The disaggregate() function provides two methods.The first (nearest neighbor, method = "") simply gives all output cells the value of the nearest input cell, and hence duplicates values which leads to a blocky output image.

The bilinear method, in turn, is an interpolation technique that uses the four nearest pixel centers of the input image (salmon colored points in Figure 5.16) to compute an average weighted by distance (arrows in Figure 5.16 as the value of the output cell - square in the upper left corner in Figure 5.16).

dem_disagg = disaggregate(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)
#> [1] FALSE

Figure 5.16: Bilinear disaggregation in action.

Comparing the values of dem and dem_disagg tells us that they are not identical (you can also use compareRaster() or all.equal()).However, this was hardly to be expected, since disaggregating is a simple interpolation technique.It is important to keep in mind that disaggregating results in a finer resolution; the corresponding values, however, are only as accurate as their lower resolution source.

The process of computing values for new pixel locations is also called resampling.In fact, the raster package provides a resample() function.It lets you align several raster properties in one go, namely origin, extent and resolution.By default, it uses the bilinear-interpolation.

# add 2 rows and columns, i.e. change the extent
dem_agg = extend(dem_agg, 2)
dem_disagg_2 = resample(dem_agg, dem)

Finally, in order to align many (possibly hundreds or thousands of) images stored on disk, you could use the gdalUtils::align_rasters() function.However, you may also use raster with very large datasets.This is because raster:

• Lets you work with raster datasets that are too large to fit into the main memory (RAM) by only processing chunks of it.
• Tries to facilitate parallel processing.For more information, see the help pages of beginCluster() and clusteR().Additionally, check out the Multi-core functions section in vignette("functions", package = "raster").