5 Geometry operations

Prerequisites

• This chapter uses the same packages as Chapter 4 but with the addition of spDataLarge, which was installed in Chapter 2:

library(sf)
library(terra)
library(dplyr)
library(spData)
library(spDataLarge)

5.1 Introduction

The previous three chapters have demonstrated how geographic datasets are structured in R (Chapter 2) and how to manipulate them based on their non-geographic attributes (Chapter 3) and spatial properties (Chapter 4). This chapter extends these skills. After reading it — and attempting the exercises at the end — you should understand and have control over the geometry column in sf objects and the geographic location of pixels represented in rasters.

Section 5.2 covers transforming vector geometries with ‘unary’ and ‘binary’ operations. Unary operations work on a single geometry in isolation. This includes simplification (of lines and polygons), the creation of buffers and centroids, and shifting/scaling/rotating single geometries using ‘affine transformations’ (Sections 5.2.1 to 5.2.4). Binary transformations modify one geometry based on the shape of another. This includes clipping and geometry unions, covered in Sections 5.2.5 and 5.2.6, respectively. Type transformations (from a polygon to a line, for example) are demonstrated in Section 5.2.7.

Section 5.3 covers geometric transformations on raster objects. This involves changing the size and number of the underlying pixels, and assigning them new values. It teaches how to change the resolution (also called raster aggregation and disaggregation), the extent and the origin of a raster. These operations are especially useful if one would like to align raster datasets from diverse sources. Aligned raster objects share a one-to-one correspondence between pixels, allowing them to be processed using map algebra operations, described in Section 4.3.2. The final Section 5.4 connects vector and raster objects. It shows how raster values can be ‘masked’ and ‘extracted’ by vector geometries. Importantly it shows how to ‘polygonize’ rasters and ‘rasterize’ vector datasets, making the two data models more interchangeable.

5.2 Geometric operations on vector data

This section is about operations that in some way change the geometry of vector (sf) objects. It is more advanced than the spatial data operations presented in the previous chapter (in Section 4.2), because here we drill down into the geometry: the functions discussed in this section work on objects of class sfc in addition to objects of class sf.

5.2.1 Simplification

Simplification is a process for generalization of vector objects (lines and polygons) usually for use in smaller scale maps. Another reason for simplifying objects is to reduce the amount of memory, disk space and network bandwidth they consume: it may be wise to simplify complex geometries before publishing them as interactive maps. The sf package provides st_simplify(), which uses the GEOS implementation of the Douglas-Peucker algorithm to reduce the vertex count. st_simplify() uses the dTolerance to control the level of generalization in map units (see Douglas and Peucker 1973 for details). Figure 5.1 illustrates simplification of a LINESTRING geometry representing the river Seine and tributaries. The simplified geometry was created by the following command:

seine_simp = st_simplify(seine, dTolerance = 2000)  # 2000 m

FIGURE 5.1: Comparison of the original and simplified geometry of the seine object.

The resulting seine_simp object is a copy of the original seine but with fewer vertices. This is apparent, with the result being visually simpler (Figure 5.1, right) and consuming less memory than the original object, as verified below:

object.size(seine)
#> 18096 bytes
object.size(seine_simp)
#> 9112 bytes

Simplification is also applicable for polygons. This is illustrated using us_states, representing the contiguous United States. As we show in Chapter 6, GEOS assumes that the data is in a projected CRS and this could lead to unexpected results when using a geographic CRS. Therefore, the first step is to project the data into some adequate projected CRS, such as US National Atlas Equal Area (epsg = 2163) (on the left in Figure 5.2):

us_states2163 = st_transform(us_states, 2163)

st_simplify() works equally well with projected polygons:

us_states_simp1 = st_simplify(us_states2163, dTolerance = 100000)  # 100 km

A limitation with st_simplify() is that it simplifies objects on a per-geometry basis. This means the ‘topology’ is lost, resulting in overlapping and ‘holy’ areal units illustrated in Figure 5.2 (middle panel). ms_simplify() from rmapshaper provides an alternative that overcomes this issue. By default it uses the Visvalingam algorithm, which overcomes some limitations of the Douglas-Peucker algorithm (Visvalingam and Whyatt 1993). The following code chunk uses this function to simplify us_states2163. The result has only 1% of the vertices of the input (set using the argument keep) but its number of objects remains intact because we set keep_shapes = TRUE:²²

# proportion of points to retain (0-1; default 0.05)
us_states2163$AREA = as.numeric(us_states2163$AREA)
us_states_simp2 = rmapshaper::ms_simplify(us_states2163, keep = 0.01,
                                          keep_shapes = TRUE)

Finally, the visual comparison of the original dataset and the two simplified versions shows differences between the Douglas-Peucker (st_simplify) and Visvalingam (ms_simplify) algorithm outputs (Figure 5.2):

FIGURE 5.2: Polygon simplification in action, comparing the original geometry of the contiguous United States with simplified versions, generated with functions from sf (center) and rmapshaper (right) packages.

5.2.2 Centroids

Centroid operations identify the center of geographic objects. Like statistical measures of central tendency (including mean and median definitions of ‘average’), there are many ways to define the geographic center of an object. All of them create single point representations of more complex vector objects.

The most commonly used centroid operation is the geographic centroid. This type of centroid operation (often referred to as ‘the centroid’) represents the center of mass in a spatial object (think of balancing a plate on your finger). Geographic centroids have many uses, for example to create a simple point representation of complex geometries, or to estimate distances between polygons. They can be calculated with the sf function st_centroid() as demonstrated in the code below, which generates the geographic centroids of regions in New Zealand and tributaries to the River Seine, illustrated with black points in Figure 5.3.

nz_centroid = st_centroid(nz)
seine_centroid = st_centroid(seine)

Sometimes the geographic centroid falls outside the boundaries of their parent objects (think of a doughnut). In such cases point on surface operations can be used to guarantee the point will be in the parent object (e.g., for labeling irregular multipolygon objects such as island states), as illustrated by the red points in Figure 5.3. Notice that these red points always lie on their parent objects. They were created with st_point_on_surface() as follows:²³

nz_pos = st_point_on_surface(nz)
seine_pos = st_point_on_surface(seine)

FIGURE 5.3: Centroids (black points) and ‘points on surface’ (red points) of New Zealand’s regions (left) and the Seine (right) datasets.

Other types of centroids exist, including the Chebyshev center and the visual center. We will not explore these here but it is possible to calculate them using R, as we’ll see in Chapter 10.

5.2.3 Buffers

Buffers are polygons representing the area within a given distance of a geometric feature: regardless of whether the input is a point, line or polygon, the output is a polygon. Unlike simplification (which is often used for visualization and reducing file size) buffering tends to be used for geographic data analysis. How many points are within a given distance of this line? Which demographic groups are within travel distance of this new shop? These kinds of questions can be answered and visualized by creating buffers around the geographic entities of interest.

Figure 5.4 illustrates buffers of different sizes (5 and 50 km) surrounding the river Seine and tributaries. These buffers were created with commands below, which show that the command st_buffer() requires at least two arguments: an input geometry and a distance, provided in the units of the CRS (in this case meters):

seine_buff_5km = st_buffer(seine, dist = 5000)
seine_buff_50km = st_buffer(seine, dist = 50000)

FIGURE 5.4: Buffers around the Seine dataset of 5 km (left) and 50 km (right). Note the colors, which reflect the fact that one buffer is created per geometry feature.

The third and final argument of st_buffer() is nQuadSegs, which means ‘number of segments per quadrant’ and is set by default to 30 (meaning circles created by buffers are composed of \(4 \times 30 = 120\) lines). This argument rarely needs to be set. Unusual cases where it may be useful include when the memory consumed by the output of a buffer operation is a major concern (in which case it should be reduced) or when very high precision is needed (in which case it should be increased).

5.2.4 Affine transformations

Affine transformation is any transformation that preserves lines and parallelism. However, angles or length are not necessarily preserved. Affine transformations include, among others, shifting (translation), scaling and rotation. Additionally, it is possible to use any combination of these. Affine transformations are an essential part of geocomputation. For example, shifting is needed for labels placement, scaling is used in non-contiguous area cartograms (see Section 8.6), and many affine transformations are applied when reprojecting or improving the geometry that was created based on a distorted or wrongly projected map. The sf package implements affine transformation for objects of classes sfg and sfc.

nz_sfc = st_geometry(nz)

Shifting moves every point by the same distance in map units. It could be done by adding a numerical vector to a vector object. For example, the code below shifts all y-coordinates by 100,000 meters to the north, but leaves the x-coordinates untouched (left panel of Figure 5.5).

nz_shift = nz_sfc + c(0, 100000)

Scaling enlarges or shrinks objects by a factor. It can be applied either globally or locally. Global scaling increases or decreases all coordinates values in relation to the origin coordinates, while keeping all geometries topological relations intact. It can be done by subtraction or multiplication of asfg or sfc object.

Local scaling treats geometries independently and requires points around which geometries are going to be scaled, e.g., centroids. In the example below, each geometry is shrunk by a factor of two around the centroids (middle panel in Figure 5.5). To achieve that, each object is firstly shifted in a way that its center has coordinates of 0, 0 ((nz_sfc - nz_centroid_sfc)). Next, the sizes of the geometries are reduced by half (* 0.5). Finally, each object’s centroid is moved back to the input data coordinates (+ nz_centroid_sfc).

nz_centroid_sfc = st_centroid(nz_sfc)
nz_scale = (nz_sfc - nz_centroid_sfc) * 0.5 + nz_centroid_sfc

Rotation of two-dimensional coordinates requires a rotation matrix:

\[ R = \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \\ \end{bmatrix} \]

It rotates points in a clockwise direction. The rotation matrix can be implemented in R as:

rotation = function(a){
  r = a * pi / 180 #degrees to radians
  matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
} 

The rotation function accepts one argument a - a rotation angle in degrees. Rotation could be done around selected points, such as centroids (right panel of Figure 5.5). See vignette("sf3") for more examples.

nz_rotate = (nz_sfc - nz_centroid_sfc) * rotation(30) + nz_centroid_sfc

FIGURE 5.5: Illustrations of affine transformations: shift, scale and rotate.

Finally, the newly created geometries can replace the old ones with the st_set_geometry() function:

nz_scale_sf = st_set_geometry(nz, nz_scale)

5.2.5 Clipping

Spatial clipping is a form of spatial subsetting that involves changes to the geometry columns of at least some of the affected features.

Clipping can only apply to features more complex than points: lines, polygons and their ‘multi’ equivalents. To illustrate the concept we will start with a simple example: two overlapping circles with a center point one unit away from each other and a radius of one (Figure 5.6).

b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points
b = st_buffer(b, dist = 1) # convert points to circles
plot(b)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y")) # add text

FIGURE 5.6: Overlapping circles.

Imagine you want to select not one circle or the other, but the space covered by both x and y. This can be done using the function st_intersection(), illustrated using objects named x and y which represent the left- and right-hand circles (Figure 5.7).

x = b[1]
y = b[2]
x_and_y = st_intersection(x, y)
plot(b)
plot(x_and_y, col = "lightgrey", add = TRUE) # color intersecting area

FIGURE 5.7: Overlapping circles with a gray color indicating intersection between them.

The subsequent code chunk demonstrates how this works for all combinations of the ‘Venn’ diagram representing x and y, inspired by Figure 5.1 of the book R for Data Science (Grolemund and Wickham 2016).

FIGURE 5.8: Spatial equivalents of logical operators.

To illustrate the relationship between subsetting and clipping spatial data, we will subset points that cover the bounding box of the circles x and y in Figure 5.8. Some points will be inside just one circle, some will be inside both and some will be inside neither. st_sample() is used below to generate a simple random distribution of points within the extent of circles x and y, resulting in output illustrated in Figure 5.9.

bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
plot(box)
plot(x, add = TRUE)
plot(y, add = TRUE)
plot(p, add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"))

FIGURE 5.9: Randomly distributed points within the bounding box enclosing circles x and y.

The logical operator way would find the points inside both x and y using a spatial predicate such as st_intersects(), whereas the intersection method simply finds the points inside the intersecting region created above as x_and_y. As demonstrated below the results are identical, but the method that uses the clipped polygon is more concise:

sel_p_xy = st_intersects(p, x, sparse = FALSE)[, 1] &
  st_intersects(p, y, sparse = FALSE)[, 1]
p_xy1 = p[sel_p_xy]
p_xy2 = p[x_and_y]
identical(p_xy1, p_xy2)
#> [1] TRUE

5.2.6 Geometry unions

As we saw in Section 3.2.3, spatial aggregation can silently dissolve the geometries of touching polygons in the same group. This is demonstrated in the code chunk below in which 49 us_states are aggregated into 4 regions using base and tidyverse functions (see results in Figure 5.10):

regions = aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
                    FUN = sum, na.rm = TRUE)
regions2 = us_states %>% group_by(REGION) %>%
  summarize(pop = sum(total_pop_15, na.rm = TRUE))

FIGURE 5.10: Spatial aggregation on contiguous polygons, illustrated by aggregating the population of US states into regions, with population represented by color. Note the operation automatically dissolves boundaries between states.

What is going on in terms of the geometries? Behind the scenes, both aggregate() and summarize() combine the geometries and dissolve the boundaries between them using st_union(). This is demonstrated in the code chunk below which creates a united western US:

us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west)

The function can take two geometries and unite them, as demonstrated in the code chunk below which creates a united western block incorporating Texas (challenge: reproduce and plot the result):

texas = us_states[us_states$NAME == "Texas", ]
texas_union = st_union(us_west_union, texas)

5.2.7 Type transformations

Geometry casting is a powerful operation that enables transformation of the geometry type. It is implemented in the st_cast function from the sf package. Importantly, st_cast behaves differently on single simple feature geometry (sfg) objects, simple feature geometry column (sfc) and simple features objects.

Let’s create a multipoint to illustrate how geometry casting works on simple feature geometry (sfg) objects:

multipoint = st_multipoint(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))

In this case, st_cast can be useful to transform the new object into linestring or polygon (Figure 5.11):

linestring = st_cast(multipoint, "LINESTRING")
polyg = st_cast(multipoint, "POLYGON")

FIGURE 5.11: Examples of linestring and polygon casted from a multipoint geometry.

Conversion from multipoint to linestring is a common operation that creates a line object from ordered point observations, such as GPS measurements or geotagged media. This allows spatial operations such as the length of the path traveled. Conversion from multipoint or linestring to polygon is often used to calculate an area, for example from the set of GPS measurements taken around a lake or from the corners of a building lot.

The transformation process can be also reversed using st_cast:

multipoint_2 = st_cast(linestring, "MULTIPOINT")
multipoint_3 = st_cast(polyg, "MULTIPOINT")
all.equal(multipoint, multipoint_2, multipoint_3)
#> [1] TRUE

For single simple feature geometries (sfg), st_cast also provides geometry casting from non-multi-types to multi-types (e.g., POINT to MULTIPOINT) and from multi-types to non-multi-types. However, only the first element of the old object would remain in the second group of cases.

Geometry casting of simple features geometry column (sfc) and simple features objects works the same as for single geometries in most of the cases. One important difference is the conversion between multi-types to non-multi-types. As a result of this process, multi-objects are split into many non-multi-objects.

Table 5.1 shows possible geometry type transformations on simple feature objects. Each input simple feature object with only one element (first column) is transformed directly into another geometry type. Several of the transformations are not possible, for example, you cannot convert a single point into a multilinestring or a polygon (so the cells [1, 4:5] in the table are NA). On the other hand, some of the transformations are splitting the single element input object into a multi-element object. You can see that, for example, when you cast a multipoint consisting of five pairs of coordinates into a point.

TABLE 5.1: Geometry casting on simple feature geometries (see Section 2.1) with input type by row and output type by column

	POI	MPOI	LIN	MLIN	POL	MPOL	GC
POI(1)	1	1	1	NA	NA	NA	NA
MPOI(1)	4	1	1	1	1	NA	NA
LIN(1)	5	1	1	1	1	NA	NA
MLIN(1)	7	2	2	1	NA	NA	NA
POL(1)	5	1	1	1	1	1	NA
MPOL(1)	10	1	NA	1	2	1	1
GC(1)	9	1	NA	NA	NA	NA	1
Note: Values like (1) represent the number of features; NA means the operation is not possible. Abbreviations: POI, LIN, POL and GC refer to POINT, LINESTRING, POLYGON and GEOMETRYCOLLECTION. The MULTI version of these geometry types is indicated by a preceding M, e.g., MPOI is the acronym for MULTIPOINT.

Let’s try to apply geometry type transformations on a new object, multilinestring_sf, as an example (on the left in Figure 5.12):

multilinestring_list = list(matrix(c(1, 4, 5, 3), ncol = 2), 
                            matrix(c(4, 4, 4, 1), ncol = 2),
                            matrix(c(2, 4, 2, 2), ncol = 2))
multilinestring = st_multilinestring((multilinestring_list))
multilinestring_sf = st_sf(geom = st_sfc(multilinestring))
multilinestring_sf
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS:           NA
#>                             geom
#> 1 MULTILINESTRING ((1 5, 4 3)...

You can imagine it as a road or river network. The new object has only one row that defines all the lines. This restricts the number of operations that can be done, for example it prevents adding names to each line segment or calculating lengths of single lines. The st_cast function can be used in this situation, as it separates one mutlilinestring into three linestrings:

linestring_sf2 = st_cast(multilinestring_sf, "LINESTRING")
linestring_sf2
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS:           NA
#>                    geom
#> 1 LINESTRING (1 5, 4 3)
#> 2 LINESTRING (4 4, 4 1)
#> 3 LINESTRING (2 2, 4 2)

FIGURE 5.12: Examples of type casting between MULTILINESTRING (left) and LINESTRING (right).

The newly created object allows for attributes creation (see more in Section 3.2.5) and length measurements:

linestring_sf2$name = c("Riddle Rd", "Marshall Ave", "Foulke St")
linestring_sf2$length = st_length(linestring_sf2)
linestring_sf2
#> Simple feature collection with 3 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS:           NA
#>                    geom         name length
#> 1 LINESTRING (1 5, 4 3)    Riddle Rd   3.61
#> 2 LINESTRING (4 4, 4 1) Marshall Ave   3.00
#> 3 LINESTRING (2 2, 4 2)    Foulke St   2.00

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:

• Georectification based on known ground control points;
• Orthorectification, which also accounts for local topography
• Image registration is used to combine images of the same thing but shot from different sensors the process of aligning one image with another (in terms of coordinate system, and resolution)

R is rather 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 argument to FALSE. This will return a raster object containing the cells whose midpoints overlap with clip.

elev = rast(system.file("raster/elev.tif", package = "spData"))
clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]
#> class       : SpatRaster 
#> dimensions  : 2, 1, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : elev 
#> min value   :   18 
#> max value   :   24

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 raster 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).

elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_2 = extend(elev, c(1, 2))

FIGURE 5.13: Original raster (left) and the same raster (right) extended by one row on the top and bottom and two columns on the left and right.

Performing an algebraic operation on two objects with differing extents in R, the terra package returns an error.

elev_3 = elev + elev_2
#> Error: [+] extents do not match

However, we can 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 of a raster is the cell corner closest to the coordinates (0, 0). The origin() function returns the coordinates of the origin. In the below example a cell corner exists with coordinates (0, 0), but that is not necessarily the case.

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)
plot(elev, add = TRUE) # and add the original raster

FIGURE 5.14: Rasters with identical values but different origins.

Note that changing the resolution (next section) frequently 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 (disagg()) the resolution of one raster.²⁵ As an example, we here change the spatial resolution of dem (found in the spDataLarge 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.):

dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
dem_agg = aggregate(dem, fact = 5, fun = mean)

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

By contrast, the disagg() function increases the resolution. However, we have to specify a method on how to fill the new cells. The disagg() function provides two methods. The default one (method = "near") simply gives all output cells the value of the 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 = disagg(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 compareGeom() 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.

5.3.4 Resampling

The above methods of aggregation and disaggregation are only suitable when we want to change the resolution of our raster by the aggregation/disaggregation factor. However, what to do when we have two or more rasters with different resolutions and origins? This is the role of resampling – a process of computing values for new pixel locations. In short, this process takes the values of our original raster and recalculates new values for a target raster with custom resolution and origin.

Several methods for recalculating (estimating) values for a raster with different resolutions/origins exist (Figure 5.17). It includes:

• Nearest neighbor - assigns the value of the nearest cell of the original raster to the cell of the target one. It is fast and usually suitable for categorical rasters.
• Bilinear interpolation - assigns a weighted average of the four nearest cells from the original raster to the cell of the target one (Figure 5.16). The fastest method for continuous rasters.
• Cubic interpolation - uses values of 16 nearest cells of the original raster to determine the output cell value, applying third-order polynomial functions. Used for continuous rasters. It results in a more smoothed surface than the bilinear interpolation, but is also more computationally demanding.
• Cubic spline interpolation - also uses values of 16 nearest cells of the original raster to determine the output cell value, but applies cubic splines (piecewise third-order polynomial functions) to derive the results. Used for continuous rasters.
• Lanczos windowed sinc resampling - uses values of 36 nearest cells of the original raster to determine the output cell value. Used for continuous rasters.²⁶

As you can find in the above explanation, only nearest neighbor is suitable for categorical rasters, while all the methods can be used (with different outcomes) for the continuous rasters. Additionally, each successive method requires more processing time.

To apply resampling, the terra package provides a resample() function. It accepts an input raster (x), a raster with target spatial properties (y), and a resampling method (method).

We need a raster with target spatial properties to see how the resample() function works. For this example, we create target_rast, but you would often use an already existing raster object.

target_rast = rast(xmin = 794600, xmax = 798200, 
                   ymin = 8931800, ymax = 8935400,
                   resolution = 150, crs = "EPSG:32717")

Next, we need to provide our two raster objects as the first two arguments and one of the resampling methods described above.

dem_resampl = resample(dem, y = target_rast, method = "bilinear")

Figure 5.17 shows a comparison of different resampling methods on the dem object.

FIGURE 5.17: Visual comparison of the original raster and five different resampling methods.

As you will see in section 6.6, raster reprojection is a special case of resampling when our target raster has a different CRS than the original raster.

Most geometry operations in terra are user-friendly, rather fast, and work on large raster objects. However, there could be some cases, when terra is not the most performant either for extensive rasters or many raster files, and some alternatives should be considered.

The most established alternatives come with the GDAL library. It contains several utility functions, including:

• gdalinfo - lists various information about a raster file, including its resolution, CRS, bounding box, and more
• gdal_translate - converts raster data between different file formats
• gdal_rasterize - converts vector data into raster files
• gdalwarp - allows for raster mosaicing, resampling, cropping, and reprojecting

All of the above functions are written in C++, but can be called in R using the gdalUtilities package. Importantly, all of these functions expect a raster file path as an input and often return their output as a raster file (for example, gdalUtilities::gdal_translate("my_file.tif", "new_file.tif", t_srs = "EPSG:4326")) This is very different from the usual terra approach, which expects SpatRaster objects as inputs.

5.4 Raster-vector interactions

This section focuses on interactions between raster and vector geographic data models, introduced in Chapter 2. It includes four main techniques: raster cropping and masking using vector objects (Section 5.4.1); extracting raster values using different types of vector data (Section 5.4.2); and raster-vector conversion (Sections 5.4.3 and 5.4.4). The above concepts are demonstrated using data used in previous chapters to understand their potential real-world applications.

5.4.1 Raster cropping

Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case raster cropping and masking are useful for unifying the spatial extent of input data. Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.

We will use two objects to illustrate raster cropping:

• A SpatRaster object srtm representing elevation (meters above sea level) in south-western Utah.
• A vector (sf) object zion representing Zion National Park.

Both target and cropping objects must have the same projection. The following code chunk therefore not only reads the datasets from the spDataLarge package9 installed in Chapter 2), it also reprojects zion (see Section 6 for more on reprojection):

srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))

We will use crop() from the terra package to crop the srtm raster. It reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument, as demonstrated in the command below (which generates Figure 5.18(B) — note the smaller extent of the raster background):

srtm_cropped = crop(srtm, vect(zion))

Related to crop() is the terra function mask(), which sets values outside of the bounds of the object passed to its second argument to NA. The following command therefore masks every cell outside of the Zion National Park boundaries (Figure 5.18(C)):

srtm_masked = mask(srtm, vect(zion))

Importantly, we want to use both crop() and mask() together in most cases. This combination of functions would (a) limit the raster’s extent to our area of interest and then (b) replace all of the values outside of the area to NA.

srtm_cropped = crop(srtm, vect(zion))
srtm_final = mask(srtm_cropped, vect(zion))

Changing the settings of mask() yields different results. Setting updatevalue = 0, for example, will set all pixels outside the national park to 0. Setting inverse = TRUE will mask everything inside the bounds of the park (see ?mask for details) (Figure 5.18(D)).

srtm_inv_masked = mask(srtm, vect(zion), inverse = TRUE)

FIGURE 5.18: Illustration of raster cropping and raster masking.

5.4.2 Raster extraction

Raster extraction is the process of identifying and returning the values associated with a ‘target’ raster at specific locations, based on a (typically vector) geographic ‘selector’ object. The results depend on the type of selector used (points, lines or polygons) and arguments passed to the terra::extract() function, which we use to demonstrate raster extraction. The reverse of raster extraction — assigning raster cell values based on vector objects — is rasterization, described in Section 5.4.3.

The basic example is of extracting the value of a raster cell at specific points. For this purpose, we will use zion_points, which contain a sample of 30 locations within the Zion National Park (Figure 5.19). The following command extracts elevation values from srtm and creates a data frame with points’ IDs (one value per vector’s row) and related srtm values for each point. Now, we can add the resulting object to our zion_points dataset with the cbind() function:

data("zion_points", package = "spDataLarge")
elevation = terra::extract(srtm, vect(zion_points))
zion_points = cbind(zion_points, elevation)

FIGURE 5.19: Locations of points used for raster extraction (left) and their 1km buffers (right).

Raster extraction also works with line selectors. Then, it extracts one value for each raster cell touched by a line. However, the line extraction approach is not recommended to obtain values along the transects as it is hard to get the correct distance between each pair of extracted raster values.

In this case, a better approach is to split the line into many points and then extract the values for these points. To demonstrate this, the code below creates zion_transect, a straight line going from northwest to southeast of the Zion National Park, illustrated in Figure 5.20(A) (see Section 2.2 for a recap on the vector data model):

zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
  st_linestring() %>% 
  st_sfc(crs = crs(srtm)) %>% 
  st_sf()

The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an ‘elevation profile’ of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.

The first step is to add a unique id for each transect. Next, with the st_segmentize() function we can add points along our line(s) with a provided density (dfMaxLength) and convert them into points with st_cast().

zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")

Now, we have a large set of points, and we want to derive a distance between the first point in our transects and each of the subsequent points. In this case, we only have one transect, but the code, in principle, should work on any number of transects:

zion_transect = zion_transect %>% 
  group_by(id) %>% 
  mutate(dist = st_distance(geometry)[, 1]) 

Finally, we can extract elevation values for each point in our transects and combine this information with our main object.

zion_elev = terra::extract(srtm, vect(zion_transect))
zion_transect = cbind(zion_transect, zion_elev)

The resulting zion_transect can be used to create elevation profiles, as illustrated in Figure 5.20(B).

FIGURE 5.20: Location of a line used for raster extraction (left) and the elevation along this line (right).

The final type of geographic vector object for raster extraction is polygons. Like lines, polygons tend to return many raster values per polygon. This is demonstrated in the command below, which results in a data frame with column names ID (the row number of the polygon) and srtm (associated elevation values):

zion_srtm_values = terra::extract(x = srtm, y = vect(zion))

Such results can be used to generate summary statistics for raster values per polygon, for example to characterize a single region or to compare many regions. The generation of summary statistics is demonstrated in the code below, which creates the object zion_srtm_df containing summary statistics for elevation values in Zion National Park (see Figure 5.21(A)):

group_by(zion_srtm_values, ID) %>% 
  summarize(across(srtm, list(min = min, mean = mean, max = max)))
#> # A tibble: 1 × 4
#>      ID srtm_min srtm_mean srtm_max
#>   <dbl>    <dbl>     <dbl>    <dbl>
#> 1     1     1122     1818.     2661

The preceding code chunk used the tidyverse to provide summary statistics for cell values per polygon ID, as described in Chapter 3. The results provide useful summaries, for example that the maximum height in the park is around 2,661 meters above see level (other summary statistics, such as standard deviation, can also be calculated in this way). Because there is only one polygon in the example a data frame with a single row is returned; however, the method works when multiple selector polygons are used.

The similar approach works for counting occurrences of categorical raster values within polygons. This is illustrated with a land cover dataset (nlcd) from the spDataLarge package in Figure 5.21(B), and demonstrated in the code below:

nlcd = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
zion2 = st_transform(zion, st_crs(nlcd))
zion_nlcd = terra::extract(nlcd, vect(zion2))
zion_nlcd %>% 
  group_by(ID, levels) %>%
  count()
#> # A tibble: 7 × 3
#> # Groups:   ID, levels [7]
#>      ID levels      n
#>   <dbl>  <dbl>  <int>
#> 1     1      2   4205
#> 2     1      3  98285
#> 3     1      4 298299
#> 4     1      5 203701
#> # … with 3 more rows

FIGURE 5.21: Area used for continuous (left) and categorical (right) raster extraction.

Polygons usually have irregular shapes, and therefore, a polygon can overlap only some parts of a raster’s cells. To get more detailed results, the extract() function has an argument called exact. With exact = TRUE, we get one more column fraction in the output data frame, which contains a fraction of each cell that is covered by the polygon. This could be useful to calculate a weighted mean for continuous rasters or more precise coverage for categorical rasters. By default, it is FALSE as this operation requires more computations.

5.4.3 Rasterization

Rasterization is the conversion of vector objects into their representation in raster objects. Usually, the output raster is used for quantitative analysis (e.g., analysis of terrain) or modeling. As we saw in Chapter 2 the raster data model has some characteristics that make it conducive to certain methods. Furthermore, the process of rasterization can help simplify datasets because the resulting values all have the same spatial resolution: rasterization can be seen as a special type of geographic data aggregation.

The terra package contains the function rasterize() for doing this work. Its first two arguments are, x, vector object to be rasterized and, y, a ‘template raster’ object defining the extent, resolution and CRS of the output. The geographic resolution of the input raster has a major impact on the results: if it is too low (cell size is too large), the result may miss the full geographic variability of the vector data; if it is too high, computational times may be excessive. There are no simple rules to follow when deciding an appropriate geographic resolution, which is heavily dependent on the intended use of the results. Often the target resolution is imposed on the user, for example when the output of rasterization needs to be aligned to the existing raster.

To demonstrate rasterization in action, we will use a template raster that has the same extent and CRS as the input vector data cycle_hire_osm_projected (a dataset on cycle hire points in London is illustrated in Figure 5.22(A)) and spatial resolution of 1000 meters:

cycle_hire_osm_projected = st_transform(cycle_hire_osm, "EPSG:27700")
raster_template = rast(ext(cycle_hire_osm_projected), resolution = 1000,
                       crs = st_crs(cycle_hire_osm_projected)$wkt)

Rasterization is a very flexible operation: the results depend not only on the nature of the template raster, but also on the type of input vector (e.g., points, polygons) and a variety of arguments taken by the rasterize() function.

To illustrate this flexibility we will try three different approaches to rasterization. First, we create a raster representing the presence or absence of cycle hire points (known as presence/absence rasters). In this case rasterize() requires only one argument in addition to x and y (the aforementioned vector and raster objects): a value to be transferred to all non-empty cells specified by field (results illustrated Figure 5.22(B)).

ch_raster1 = rasterize(vect(cycle_hire_osm_projected), raster_template,
                       field = 1)

The fun argument specifies summary statistics used to convert multiple observations in close proximity into associate cells in the raster object. By default fun = "last" is used but other options such as fun = "length" can be used, in this case to count the number of cycle hire points in each grid cell (the results of this operation are illustrated in Figure 5.22(C)).

ch_raster2 = rasterize(vect(cycle_hire_osm_projected), raster_template, 
                       fun = "length")

The new output, ch_raster2, shows the number of cycle hire points in each grid cell. The cycle hire locations have different numbers of bicycles described by the capacity variable, raising the question, what’s the capacity in each grid cell? To calculate that we must sum the field ("capacity"), resulting in output illustrated in Figure 5.22(D), calculated with the following command (other summary functions such as mean could be used):

ch_raster3 = rasterize(vect(cycle_hire_osm_projected), raster_template, 
                       field = "capacity", fun = sum)

FIGURE 5.22: Examples of point rasterization.

Another dataset based on California’s polygons and borders (created below) illustrates rasterization of lines. After casting the polygon objects into a multilinestring, a template raster is created with a resolution of a 0.5 degree:

california = dplyr::filter(us_states, NAME == "California")
california_borders = st_cast(california, "MULTILINESTRING")
raster_template2 = rast(ext(california), resolution = 0.5,
                        crs = st_crs(california)$wkt)

When considering line or polygon rasterization, one useful additional argument is touches. By default it is FALSE, but when changed to TRUE – all cells that are touched by a line or polygon border get a value. Line rasterization with touches = TRUE is demonstrated in the code below (Figure 5.23(A)).

california_raster1 = rasterize(vect(california_borders), raster_template2,
                               touches = TRUE)

Compare it to a polygon rasterization, with touches = FALSE by default, which selects only cells whose centroids are inside the selector polygon, as illustrated in Figure 5.23(B).

california_raster2 = rasterize(vect(california), raster_template2) 

FIGURE 5.23: Examples of line and polygon rasterizations.

5.4.4 Spatial vectorization

Spatial vectorization is the counterpart of rasterization (Section 5.4.3), but in the opposite direction. It involves converting spatially continuous raster data into spatially discrete vector data such as points, lines or polygons.

Be careful with the wording! In R, vectorization refers to the possibility of replacing for-loops and alike by doing things like 1:10 / 2 (see also Wickham (2019)).

The simplest form of vectorization is to convert the centroids of raster cells into points. as.points() does exactly this for all non-NA raster grid cells (Figure 5.24). Note, here we also used st_as_sf() to convert the resulting object to the sf class.

elev_point = as.points(elev) %>% 
  st_as_sf()

FIGURE 5.24: Raster and point representation of the elev object.

Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures (isotherms) for example. We will use a real-world digital elevation model (DEM) because the artificial raster elev produces parallel lines (task: verify this and explain why this happens). Contour lines can be created with the terra function as.contour(), which is itself a wrapper around filled.contour(), as demonstrated below (not shown):

dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
cl = as.contour(dem)
plot(dem, axes = FALSE)
plot(cl, add = TRUE)

Contours can also be added to existing plots with functions such as contour(), rasterVis::contourplot() or tmap::tm_iso(). As illustrated in Figure 5.25, isolines can be labelled.

# create hillshade
hs = shade(slope = terrain(dem, "slope", unit = "radians"),
           aspect = terrain(dem, "aspect", unit = "radians"))
plot(hs, col = gray(0:100 / 100), legend = FALSE)
# overlay with DEM
plot(dem, col = terrain.colors(25), alpha = 0.5, legend = FALSE, add = TRUE)
# add contour lines
contour(dem, col = "white", add = TRUE)

FIGURE 5.25: DEM hillshade of the southern flank of Mt. Mongón overlaid by contour lines.

The final type of vectorization involves conversion of rasters to polygons. This can be done with terra::as.polygons(), which converts each raster cell into a polygon consisting of five coordinates, all of which are stored in memory (explaining why rasters are often fast compared with vectors!).

This is illustrated below by converting the grain object into polygons and subsequently dissolving borders between polygons with the same attribute values (also see the dissolve argument in as.polygons()).

grain = rast(system.file("raster/grain.tif", package = "spData"))
grain_poly = as.polygons(grain) %>% 
  st_as_sf()

FIGURE 5.26: Illustration of vectorization of raster (left) into polygon (center) and polygon aggregation (right).

5.5 Exercises

Some of the exercises use a vector (zion_points) and raster dataset (srtm) from the spDataLarge package. They also use a polygonal ‘convex hull’ derived from the vector dataset (ch) to represent the area of interest:

library(sf)
library(terra)
zion_points = read_sf(system.file("vector/zion_points.gpkg", package = "spDataLarge"))
srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
ch = st_combine(zion_points) %>%
  st_convex_hull() %>% 
  st_as_sf()

E1. Generate and plot simplified versions of the nz dataset. Experiment with different values of keep (ranging from 0.5 to 0.00005) for ms_simplify() and dTolerance (from 100 to 100,000) st_simplify().

• At what value does the form of the result start to break down for each method, making New Zealand unrecognizable?
• Advanced: What is different about the geometry type of the results from st_simplify() compared with the geometry type of ms_simplify()? What problems does this create and how can this be resolved?

E2. In the first exercise in Chapter Spatial data operations it was established that Canterbury region had 70 of the 101 highest points in New Zealand. Using st_buffer(), how many points in nz_height are within 100 km of Canterbury?

E3. Find the geographic centroid of New Zealand. How far is it from the geographic centroid of Canterbury?

E4. Most world maps have a north-up orientation. A world map with a south-up orientation could be created by a reflection (one of the affine transformations not mentioned in this chapter) of the world object’s geometry. Write code to do so. Hint: you need to use a two-element vector for this transformation. Bonus: create an upside-down map of your country.

E5. Subset the point in p that is contained within x and y.

• Using base subsetting operators.
• Using an intermediary object created with st_intersection().

E6. Calculate the length of the boundary lines of US states in meters. Which state has the longest border and which has the shortest? Hint: The st_length function computes the length of a LINESTRING or MULTILINESTRING geometry.

E7. Crop the srtm raster using (1) the zion_points dataset and (2) the ch dataset. Are there any differences in the output maps? Next, mask srtm using these two datasets. Can you see any difference now? How can you explain that?

E8. Firstly, extract values from srtm at the points represented in zion_points. Next, extract average values of srtm using a 90 buffer around each point from zion_points and compare these two sets of values. When would extracting values by buffers be more suitable than by points alone?

E9. Subset points higher than 3100 meters in New Zealand (the nz_height object) and create a template raster with a resolution of 3 km. Using these objects:

• Count numbers of the highest points in each grid cell.
• Find the maximum elevation in each grid cell.

E10. Aggregate the raster counting high points in New Zealand (created in the previous exercise), reduce its geographic resolution by half (so cells are 6 by 6 km) and plot the result.

• Resample the lower resolution raster back to a resolution of 3 km. How have the results changed?
• Name two advantages and disadvantages of reducing raster resolution.

E11. Polygonize the grain dataset and filter all squares representing clay.

• Name two advantages and disadvantages of vector data over raster data.
• At which points would it be useful to convert rasters to vectors in your work?