Chapter 10 Combining rasters and vector layers

Last updated: 2020-08-12 00:39:38

Aims

Our aims in this chapter are:

  • Crop and mask a raster according to a vector layer
  • Switch from vector to raster representation, and vice versa
  • Extract raster values from locations defined by a vector layer

We will use the following R packages:

  • sf
  • stars
  • units

10.1 Masking and cropping rasters

10.1.1 Introduction

Masking a raster means turning pixels values outside of a an area of interest—defined using a polygonal layer—into NA (Figure 10.1). Cropping a raster means deleting whole rows and/or columns, so that raster extent is reduced to a new (smaller) rectangular shape, also according to the extent of a vector layer. The [ operator can be used for masking or masking and cropping (the default).

Masking a raster^[http://rpubs.com/etiennebr/visualraster]

Figure 10.1: Masking a raster48

10.1.2 Masking and cropping

For an example of masking and cropping, we will prepare a raster of average NDVI over the period 2000-2019, based on MOD13A3_2000_2019.tif. The following code section combines what we learned in Sections 6.3.2, 6.6.1.3 and 9.3:

We will also read a Shapefile with a polygon of Israel borders, named israel_borders.shp:

Figure 10.2 shows the average NDVI raster and the borders polygon:

Raster and crop/mask geometry

Figure 10.2: Raster and crop/mask geometry

Masking and cropping a raster based on an sf layer can be done with the [ operator, as follows:

As we can see in Figure 10.3, all pixels that were outside of the polygon are now NA (this is called masking). It’s difficult to see, but the raster extent is also slightly reduced according to the extent of borders (this is called cropping):

Raster after cropping and masking

Figure 10.3: Raster after cropping and masking

Zooming in on a small portion of the masked raster (Figure 10.4) demonstrates that pixels whose center is not inside polygon were converted to NA:

Raster after cropping and masking

Figure 10.4: Raster after cropping and masking

10.1.3 Masking-only

Masking-only is done with the same operator ([), using the (non-default) argument crop=FALSE. To demonstrate the difference between masking+cropping (which we just did) and masking only, consider the following example where r_avg is masked and cropped (r_avg1), or just masked (r_avg2), according to the "Negev" administrative area. First we prepare the polygon used for masking and cropping, obtained from a Shapefile named nafot.shp:

Then we produce the masked+cropped (r_avg1) and masked (r_avg2) rasters:

As we can see in Figure 10.5, masking transforms pixel values to NA (left panel) while cropping reduces raster extent by deleting whole rows and columns (right panel).

Cropping and masking (left) vs. masking (right), raster extent is shown in greyCropping and masking (left) vs. masking (right), raster extent is shown in grey

Figure 10.5: Cropping and masking (left) vs. masking (right), raster extent is shown in grey

Note that for plotting the raster extents (grey boxes in Figure 10.5), the above code section uses a combination of st_bbox and st_as_sfc. The combination returns the bounding box of a stars (or sf) layer as an sfc (geometry) object. For example, the following expression returns a geometry column with a single polygon, which is the bounding box of r_avg1:

10.2 Vector layer → raster

10.2.1 Introduction

The st_rasterize function converts a vector layer to a raster, given two arguments:

  • sf—The vector layer to convert
  • template—A raster “template” (if missing it can be generated automatically)

The resulting raster retains the original values (from the template) in pixels that do not overlap with the vector layer (Figure 10.6). In those pixels that do overlap with the vector layer, the value of the (first) vector layer attribute is “burned” into those pixels. The values are burned sequentially, according to the order of features. Therefore, when there are more than one vector features coinciding with the same pixel, the last feature “wins”.

Conversion of points, lines and polygons to raster^[http://rpubs.com/etiennebr/visualraster]

Figure 10.6: Conversion of points, lines and polygons to raster49

10.2.2 Rasterizing polygon attributes

For example, let’s rasterize the nafot vector layer to a raster, using r_avg as a template. The attribute which will be burned into the raster is the Id of the administrative area:

The resulting raster s has the same dimensions of r_avg but different values, as shown on Figure 10.7. We can see that those areas not covered by the nafot layer retained their original value (average NDVI). All pixels that were covered by nafot got the attribute value (Id).

Rasterizing the `Id` of the `nafot` administrative areas layer into the average NDVI raster

Figure 10.7: Rasterizing the Id of the nafot administrative areas layer into the average NDVI raster

Often we want to start with an empty template, so that the only values in the resulting raster are those “burned” from the vector layer. For example, we can use a copy of the r_avg raster, where all of the pixel values are replaced with NA, as the template:

The result is shown in Figure 10.8:

Rasterizing into an empty template

Figure 10.8: Rasterizing into an empty template

10.2.3 Rasterizing point counts

Sometimes we want to add up attribute values from overlapping features coinciding with the same pixel, rather that take the last value. This is useful, for instance, when calculating a point density raster. For example, let’s read the plants.shp and reserve.shp Shapefiles, containing a point layer of rare plants observations and a polygon with the borders of a nature reserve in the Negev:

We will also reproject both layers to UTM, and subset only those plants that are inside the reserve:

Figure 10.9 shows both layers—the polygonal layer reserve and the point layer plants:

Nature reserve and rare plants observations

Figure 10.9: Nature reserve and rare plants observations

How can we calculate a raster expressing the density of plants points across the nature reserve? First, we set up a empty template raster, but this time the area of interest—the nature reserve—is given an initial value of zero:

Then, we create a new attrubute named count, with a value of 1 for each plant observation:

Finally, we rasterize the "count" attribute into the template raster, with an additional option options="MERGE_ALG=ADD". The "MERGE_ALG=ADD" option acts as a flag instructing st_rasterize to add up any overlapping attributes “burned” into the same pixel:

The resulting raster gives the number of rare plant observations per pixel (Figure 10.10):

Density (observations per pixel) of rare plants in the nature reserve

Figure 10.10: Density (observations per pixel) of rare plants in the nature reserve

The raster coloring is almost uniform, because the distribution is highly skewed (try running hist(s[[1]])). We can use a logarithmic scale to visualize the density pattern more clearly (Figure 10.11):

Density (observations per pixel) of rare plants in the nature reserve, with a logarithmic scale

Figure 10.11: Density (observations per pixel) of rare plants in the nature reserve, with a logarithmic scale

10.2.4 Standardizing density units

The units of measurement of the density raster are, at the moment, plants per pixel, which happens to have an area of 0.86 \(km^2\). We can find out the area of a pixel by applying the st_area function (Section 8.3.2.2) on the raster, which returns a raster of area values per pixel:

It is more convenient to work with counts per standard area unit, such as plants per 1 \(km^2\). To do that, we can calculate a raster with pixel areas a:

and then divide the plants count raster s by a:

The result is a raster with counts per unit area, in this case plant observations per \(km^2\) (Figure 10.12):

Density (observations per $km^2$) of rare plants in the nature reserve, with a logarithmic scale

Figure 10.12: Density (observations per \(km^2\)) of rare plants in the nature reserve, with a logarithmic scale

10.3 Raster → Polygons

10.3.1 Raster to polygons conversion

The st_as_sf function makes the Raster→Polygons, using the (default) as_points=FALSE. The function ignores pixels that have a NA value in all layers. A useful option merge=TRUE dissolves all adjacent polygons that have the same raster value (in the first layer) into a single feature. The dissolve algorithm can be specified with the connect8 parameter: 4d, which is the default (connect8=FALSE) or 8d (connect8=TRUE). The attribute table of the polygon layer contains the raster values—with a separate attribute for each layer.

For example, let’s take a subset of the r:

  • Rows 100-102
  • Columns 200-202
  • Layers 1-2

We will replace some of the pixel values with NA:

and also round the values to one digit:

The resulting small raster u is visualized in Figure 10.13. Note that the raster has one pixel with NA in both layers and another pixel with NA in the first layer only.

Sample raster

Figure 10.13: Sample raster

Now, let’s transform the raster u to a polygon layer p, using st_as_sf:

The resulting polygonal layer p has 8 polygons, even though the raster u had 9 pixels. That is because one of the pixels had NA in all layers and therefore was not converted to a polygon:

The p layer is shown in Figure 10.14:

Polygon layer created from a raster

Figure 10.14: Polygon layer created from a raster

Let’s try the merge=TRUE option. Since merging is only affected by the first layer, it makes sense transform each of the layers separately:

As a result, we have two separate polygonal layers p1 and p2:

Figure 10.15 displays both p1 and p2:

Polygon layer created from a rasterPolygon layer created from a raster

Figure 10.15: Polygon layer created from a raster

Do you think there will there be a difference if we use connect8=TRUE?

10.3.2 Example: segmentation

One example of a use case of raster to polygon conversion is the delineation of distinct inter-connected groups of pixels, sharing the same (or a similar) value. This operationis also known as segmentation. In it’s simplest form—detecting segments with exactly the same value—segmentation can be done using st_as_sf with merge=TRUE.

For example, we can derive segments of NDVI>0.2 in the reclassified NDVI raster l_rec_focal from Section 9.4.3. First, let’s recreate the l_rec_focal raster:

Segments in the l_rec_focal raster represent distinct continuous areas with \(NDVI>0.2\). To detect them, we can convert the raster to polygons using st_as_sf with the merge=TRUE:

Then, we filter only those segments where raster value was 1:

The result is a polygon layer where each feature represents a single continuous area where \(NDVI>0.2\):

Figure 10.16) shows the segments on top of the NDVI raster:

Segments with NDVI>0.2

Figure 10.16: Segments with NDVI>0.2

What is the exact number of segments in pol? If we ran the st_as_sf function on l_rec instead of l_rec_focal, do you think the number of segments would be higher or lower?

10.4 Raster → Points

The Raster→Points transformation is done using function st_as_sf with the as_points=TRUE option. Pixel centers—except for pixels with NA in all layers—become points. The attribute table of the point layer contains the raster values.

For example, here is how we can transform the small raster u to a point layer:

The resulting point layer p has 8 points, even though the raster has 9 pixels, because one of the pixels had NA in all layers and therefore was not converted to a point:

The point layer p is shown in Figure 10.17, with the first and second attribute values in black and in red, respectively:

Point layer created from a raster

Figure 10.17: Point layer created from a raster

10.5 Distance to nearest point

Another example of a spatial operator involving a raster and a vector layer is the calculation of a raster of distances to nearest point. For example, we may be interested in mapping the distances to nearest meteorological stations in Israel, perhaps to evaluate where coverage is too sparse and reliable meteorological data are missing.

Since the st_distance function (Section 8.3.2.3) expects vector layers as input, to calculate a raster of distances from nearest point we first need to convert the “template” raster to a point later. Our “template” for the distances raster, hereby grid is the \(926\times926\) \(m^2\) MODIS NDVI raster r_avg, but we could similarly use any other raster:

Given the grid point layer, we can use the st_distance to calculate pairwise distances between every grid point (grid) and meteorological station (rainfall):

The result is a distance matrix named distance. Its rows correspond to grid points and its columns correspond to rainfall points:

In this example, we are interested in the minimal distance, per grid point, among the distances to the various meteorological stations. Therefore we apply the min function on the distance matrix rows, which refer to grid points. The result is a numeric vector of minimal distances, which we assign as an attribute in grid:

The grid point layer now has a new distance attribute:

The distance attribute values reflect distance to nearest station, as is shown in Figure 10.18:

Distance to nearest meteorological station (point grid)

Figure 10.18: Distance to nearest meteorological station (point grid)

To have the results back as a raster, we can rasterize the points (Section 10.2). We use the same raster which we started with as template—this guarantees that every point corresponds exactly to one pixel:

It is also useful to set the raster value measurement units, for example to easily convert the distances from \(m\) to \(km\):

The final raster with distances to the nearest stations is shown in Figure 10.19. The locations of the meteorological stations and contour lines are shown on top of the raster, to emphasize the distance gradient.

Distance to nearest meteorological station (raster)

Figure 10.19: Distance to nearest meteorological station (raster)

10.6 Raster → Lines (contours)

We already saw how a raster can be converted to points or polygons, which simply results in a layer of cell centroids or outlines, respectively. Another common transformation is to calculate contour lines of equal raster values.

To illustrate contour lines calculation, we will reconstruct the Haifa DEM from Chapter 9:

Raster contours can be created using the st_contour function. The function can accept a breaks argument, which determines the break points where contours are generated. Another parameter named contour_lines determines whether contours are returned as a line layer (TRUE) or polygon layer (FALSE, the default).

To decide on breaks values we are interested in, it is useful to examine the range of values in the raster:

According to the above result, we will use breaks from -100 to 600 meters, in steps of 50 meters:

The result, dem_contour, is a LINESTRING layer with one feature per contour line:

The contours are shown in Figure 10.20:

Contour lines

Figure 10.20: Contour lines

10.7 Extracting raster values

10.7.1 Introduction

It is often necessary to “extract” raster values according to overlapping vector features, such as points, lines or polygons (Figure 10.21). For example, given an NDVI raster, we may be interested to calculate the NDVI value observed in particular point locations, or the average NDVI observed over an administrative area polygon.

Extracting raster values^[http://rpubs.com/etiennebr/visualraster]

Figure 10.21: Extracting raster values50

In the next few sections, we will see examples of extracting raster values:

  • From single-band raster to points (Section 10.7.2)
  • From multi-band raster to points (Section 10.7.3)
  • From single-band raster to polygons (Section 10.7.4)
  • From multi-band raster to polygons (Section 10.7.5)

10.7.2 Extracting to points: single-band

10.7.2.1 NDVI in meteorological stations

Raster values can be extracted to points by:

  • Converting the raster to a polygon layer (Section 10.4); and then
  • spatially joining the polygon attributes to the points (Section 8.1).

The spatial join is one-to-one, because each point falls in exactly one raster pixel (or none, if it is beyond raster extent)51.

For example, we can determine the average NDVI (r_avg) in the pixel where each meteorological station (rainfall) falls in as follows:

As a result, the rainfall layer now has an additional attribuite named NDVI:

The attribute values are shown as text labels in Figure 10.22:

Raster values extracted to points

Figure 10.22: Raster values extracted to points

Now that we know the average NDVI for each meteorological station location, we can examine, for instance, the association between NDVI and rainfall in December (Figure 10.23):

Average NDVI as function of average rainfall in December

Figure 10.23: Average NDVI as function of average rainfall in December

What is the number and proportion of rainfall points that have an NA value in the NDVI attribute? Why did those stations get NA?

10.7.2.2 Elevation along GPS track

The following example uses a GPX (GPS Exchange Format) file. A GPX file contains one or more layers, of the following types:

  • waypoints
  • tracks (recorder track as line)
  • routes
  • track_points (recorder track as points)
  • route_points

When reading a file format that contains more than one layer, such as GPX, we may first need to find out which layers the file contains. This can be done using st_layers:

which returns the list of layers contained in the GPX file:

## Driver: GPX 
## Available layers:
##     layer_name     geometry_type features fields
## 1    waypoints             Point        0     23
## 2       routes       Line String        0     12
## 3       tracks Multi Line String        1     13
## 4 route_points             Point        0     25
## 5 track_points             Point     5091     26

Once we know which layer we want to read, we can specify at as the second argument in st_read:

Next, we will prepare a \(1\times1\) \(km^2\) elevation raster named dem1km. The methods we use should be familiar from Chapter 9:

We will also reproject the track point layer to UTM and crop the raster according to its extent, plus a 2 \(km\) buffer:

The elevation raster and the GPS track are shown in Figure 10.24:

GPS track and elevation

Figure 10.24: GPS track and elevation

To extract the elevation values per GPS point, we convert the raster to a polygon layer, then join the polygon attributes (elevation) with the track point:

The track layer contains a field named ele, which contains elevation recorded by the GPS device. Figure 10.25 shows the elevation profile recorded by the GPS (ele, in red) and the elevation profile extracted from the DEM (elevation, in black):

Elevation profile, based on GPS measurements (in red) and a DEM raster (in black)

Figure 10.25: Elevation profile, based on GPS measurements (in red) and a DEM raster (in black)

What is the reason for the stairs-like shape of the black line in Figure 10.25?

10.7.3 Extracting to points: multi-band

As another example, we can extract the NDVI values for different dates from the multi-band raster r. For simplicity, let’s create a subset rainfall1, with just three meteorological stations from the rainfall layer:

Also, we will remove all non-spatial attributes, by creating a new sf layer with nothing but the geometry from rainfall1:

To extract the NDVI values from the first five bands of r, we transform the raster to a polygon layer and perform a spatial join:

As a result, five new attributes were joined:

To work with the data, we can rearrange them to a more convenient form. For example, the geometry column can be discarded:

Then, we can rely on the fact that row and column order matches the vector layer feature and raster layer order, respectively, to assign names. For example, rainfall1 rows correspond to the selected station names:

Columns in rainfall1 are already named according to image dates. Where did these names come from?

The analogous operation in ArcGIS is the “Extract Multi Values to Points” tool (Figures 10.26–10.27).

"Extract Multi Values to Points" tool in ArcGIS

Figure 10.26: “Extract Multi Values to Points” tool in ArcGIS

"Extract Multi Values to Points" tool in ArcGIS

Figure 10.27: “Extract Multi Values to Points” tool in ArcGIS

10.7.4 Extracting to polygons: single-band

Extracting raster values to polygons (or to lines) is different from extracting to points. When extracting raster values to polygons (or to lines), each vector feature is potentially associated with more than one pixel value (Figure 10.21). Moreover, the number of pixels may vary between features. For example, a large polygon may cover many more pixels than a small polygon. Therefore, it is often convenient to summarize the raster values per polygon feauture by applying a function on them and obtaining a single number, such as the average. This is done by aggregating the raster based on a vector layer using function aggregate52.

When summarizing raster values per polygon feature, the aggregate function accepts53:

  • x—The stars raster to summarize
  • by—An sf layer determining areas
  • FUN—The function to be applied on each set of extracted pixel values from x

For example, we can calculate the average (mean) NDVI value (r_avg) per administrative area (nafot) as follows:

The result is a special type of a stars object, with just one dimension for the vector geometries:

It can be converted to the more familiar sf structure using st_as_sf:

Now we have a polygon layer of the administrative areas, with an NDVI attribute containing the average NDVI value according to the r_avg raster:

The result can be visualized as follows (Figure 10.28):

Average NDVI per "Nafa"

Figure 10.28: Average NDVI per “Nafa”

10.7.5 Extracting to polygons: multi-band

Extracting raster values from a multi-band raster works the same way as extracting from a multi-band raster. The difference is that the resulting stars object has two dimensions: one for the geometries and one for the input raster bands.

For example, the following expression calculates the average r value for each nafot feature over five months starting from 2000-02-01, based on the first five layers of r:

When transformed to an sf layer, the values extracted from each raster layer are placed in separate attributes:

The result reflects the average NDVI, for the first five months in the MODIS NDVI series, for each “Nafa” (Figure 10.29):

NDVI per "Nafa" for five months

Figure 10.29: NDVI per “Nafa” for five months

  1. http://rpubs.com/etiennebr/visualraster↩

  2. http://rpubs.com/etiennebr/visualraster↩

  3. http://rpubs.com/etiennebr/visualraster↩

  4. In rare cases, the point may fall exactly on the border between two raster pixels, in which case it will be joined to two different raster values. Therefore, to be safe, the join result may be filtered to remove duplicated points.↩

  5. Recall that we already used aggregate for a different use case: dissolving vector layers by attribute (Section 8.4).↩

  6. Note that using aggregate to extract raster values by polygons is only appropriate for non-overlapping polygons. For polygons that are potentially overlapping, use raster::extract.↩