Chapter 10 Combining rasters and vector layers

Last updated: 2021-01-16 22:27:42

Aims

Our aims in this chapter are to learn how to:

  • Crop and mask a raster according to a vector layer
  • Switch from vector to raster representation, and vice versa
  • Calculate a raster of distances to nearest point
  • 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 raster (http://rpubs.com/etiennebr/visualraster)

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 demonstrates the masking algorithm: pixels whose centroid does not intersect with the polygon are converted to NA (Figure 10.4).

A small portion of the raster `r_avg`, brefore (left) after (right) cropping and maskingA small portion of the raster `r_avg`, brefore (left) after (right) cropping and masking

Figure 10.4: A small portion of the raster r_avg, brefore (left) after (right) cropping and masking

10.1.3 Masking-only

In case we need the output to keep the same extent as the input, we can mask a raster without cropping. 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 The st_rasterize function

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). For pixels that do overlap with the vector layer, the value of the (first) vector layer attribute is “burned” into the pixels. The default meaning of the term “overlap” in the context of point, line, or polygon geometries, is:

The additional parameter options="ALL_TOUCHED=TRUE" can be passed to st_rasterize so that all pixels intersecting with the geometry are considered “overlapping” (Figure 10.6).

Rasterizing points, lines and polygons to raster with `st_rasterize`, using the default algorithms (top) and `options="ALL_TOUCHED=TRUE"` (bottom).

Figure 10.6: Rasterizing points, lines and polygons to raster with st_rasterize, using the default algorithms (top) and options="ALL_TOUCHED=TRUE" (bottom).

By default, 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 attribute value of the last feature “wins”. Alternatively, we can add up all values of coinciding geometries, using the argument options="MERGE_ALG=ADD" (Section 10.2.3).

10.2.2 Rasterizing polygon attributes

For an example of rasterizing a polygon layer, let’s convert the nafot vector layer to a raster, using r_avg as the 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 in 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 Id attribute value.

Rasterizing the `Id` of the `nafot` administrative areas layer into the average NDVI rasterRasterizing the `Id` of the `nafot` administrative areas layer into the average NDVI rasterRasterizing 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

Instead of retaining the original values, typically we want to start with an empty template, so that the only values in the resulting raster are those coming 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 templateRasterizing into an empty templateRasterizing 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. These Shapefiles contain a point layer of rare plants observations and a polygon with the borders of Negev Mountains Nature Reserve, the largest nature reserve in Israel, respectively:

We will also reproject both layers to UTM, so that all pixels correspond to equal area size, which is more appropriate for density calculations:

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, where all pixels within the area of interest—the nature reserve—are given an initial value of zero. Note that we are hereby using r_avg as template, but, if necessary, we can prepare a new template with resolution of our choice, using st_as_stars (Section 9.2.1).

Figure 10.10 shows the template:

Template for calculating point density

Figure 10.10: Template for calculating point density

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" (Section 10.2.1). 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 values reflect the sum of the count attribute from all overlapping points per pixel, i.e., the number of rare plant observations per pixel (Figure 10.11):

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

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

The raster coloring is almost uniform, because the distribution is highly skewed. Namely, there are a lot of pixels with zero count, and few pixels with high count (try running hist(s[[1]]) to see that). We can use a logarithmic scale to visualize the density pattern more clearly (Figure 10.12). Now we can see which parts of the reserve were characterized by higher occurence of rare plants.

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

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

10.2.4 Standardizing density units

In the observation density raster s (Section 10.2.3), the units of measurement are plants per pixel. Since a pixel happens to have an area of 0.86 \(km^2\) (see below), raster values are, in fact, counts per 0.86 \(km^2\).

How do we know that the area of a pixel is 0.86 \(km^2\)? We can apply the st_area function to calculate pixel area sizes. Recall that, when applied on a vector layer, st_area returns a vectot of units values, with the corresponding areas of all geometries (Section 8.3.2.2). When st_area is applied on a raster, it returns a raster of area values per pixel:

As usual, the area values are given in CRS units (in this case, \(m^2\)). We can transform the matrix with the entire raster values (a[[1]]) from \(m^2\) to \(km^2\) units, using set_units (Section 8.3.2.2):

Finally, we can examine one of the raster values. It doesn’t matter which one, since all pixels have exactly the same area. For example, the following expression gives the top-left corner value from the stars values matrix:

Let’s return to the density raster s. It is more convenient to work with counts per standard area unit, such as plants per 1 \(km^2\). To do that, we can divide the plants count raster s, by pixel area raster a:

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

The result is shown in Figure 10.13:

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

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

We now move on to the opposite operation, converting a raster to a vector layer. As we will see, there are three variants of this conversion:

  • Raster to polygons (Section 10.3)
  • Raster to points (Section 10.4)
  • Raster to contour (Section 10.6)

10.3 Raster → Polygons

10.3.1 Raster to polygons conversion

The st_as_sf function makes the raster to polygons conversion, when using the (default) as_points=FALSE argument. The st_as_sf function ignores pixels that have a NA values in all layers35. 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 small subset of the raster r:

We will replace some of the pixel values with NA:

and round the values to one digit:

The resulting small raster u is visualized in Figure 10.14. 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.14: Sample raster

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

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

The resulting layer p, which has eight features and two attributes, is shown in Figure 10.15:

Polygon layer created from a raster

Figure 10.15: 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 to transform each of the layers separately:

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

Figure 10.16 displays both p1 and p2:

Polygon layer created from a rasterPolygon layer created from a raster

Figure 10.16: Polygon layer created from a raster

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

10.3.2 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.17) shows the segments on top of the NDVI raster:

Segments with NDVI>0.2

Figure 10.17: 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. Just like in a conversion to polygons (Section 10.3), the attribute table of the resulting point layer contains the raster values.

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

Once again (Section 10.3), the resulting point layer p has eight points, even though the raster has nine pixels, because one of the pixels had NA in all layers and therefore was not converted to a point36:

The point layer p is shown in Figure 10.18:

Point layer created from a raster

Figure 10.18: 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 the nearest meteorological stations in Israel, to evaluate where coverage is too sparse and reliable meteorological data are missing.

We are already familiar with the st_distance function for calculating distances (Section 8.3.2.3). However, st_distance expects two vector layers as input, not a vector layer and a raster. Therefore, to calculate a raster of distances from nearest point, we first need to convert the “template” raster to a point later (Section 10.4). In this example, our template for the distances raster, hereby named grid, is the \(926\times926\) \(m^2\) raster r_avg (though we could use any other raster template). The NDVI attribute is discarded, since it is irrelevant for the distance calculation:

Let’s also import the point layer of meteorological stations (Section 7.4):

Now we can use 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. This is a very large matrix, due to the large number of grid points:

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

It is also convenient to convert the distances from from \(m\) to \(km\). Since apply returns numeric, thus “losing” the units of measurement, we first need to convert the distances to units and only then make the conversion:

The grid point layer now has a distance attribute, with the distance to nearest meteorological station, in \(km\):

The result is shown in Figure 10.19:

Distance to nearest meteorological station (point grid)

Figure 10.19: 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:

The final raster with distances to the nearest stations is shown in Figure 10.20. 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.20: Distance to nearest meteorological station (raster)

10.6 Raster → Lines (contours)

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

To illustrate contour lines calculation, we will reconstruct the Haifa DEM at 90 \(m\) resolution (Section 9.3):

Raster contours can be calculated 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.21:

Contour lines

Figure 10.21: 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.22). 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.22: Extracting raster values (http://rpubs.com/etiennebr/visualraster)

When extracting values to liner or polygons, it is common to summarize the values that the geometry intersects with, using a function such as mean. That way, the values can “fit” into a new attribute or attributes. Extracting to points is simpler, since each geometry corresponds to one pixel, so there is nothing to summarize.

Extracting to points can be accompished with the st_extract function, while extracting to lines or to polygons can be accomplished with the aggregate function. In the next few sections, we will see examples of the most common extract scenarios:

  • 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 using the st_extract function. The function accepts:

  • x—a stars raster
  • pts—a sf point layer

and returns a new sf layer, with an additional attribute containing the values from x from the matching pixel for each point in pts.

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

The extracted values can be attached to the rainfall layer, using assignment:

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

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

Raster values extracted to points

Figure 10.23: 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.24):

Average NDVI as function of average rainfall in December

Figure 10.24: 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.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:

Now, let’s extract the NDVI values from a subset with the first five bands of r, into rainfall1:

This time, since extraction from a multi-band raster took place, the result is a stars object and not an sf object. Moreover, this is a special type of a stars object, with one dimension for the vector geometries (geometry) and another dimension for the raster layers (time):

To work with the data, it is often more convenient to transform the stars object to an sf layer. This can be done with st_as_sf:

In case we are interested only in the extracted values, not the geometries, the latter can be discarded with st_drop_geometry:

Finally, 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, therefore we can assign the station names into rownames:

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

"Extract Multi Values to Points" tool in ArcGIS

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

"Extract Multi Values to Points" tool in ArcGIS

Figure 10.26: “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.22). 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 aggregate37.

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

  • 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. Any additional arguments, such as na.rm=TRUE, are passed to the function, in this case so that mean excludes NA values:

The result is, again (Section 10.7.3), 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.27):

Average NDVI per "Nafa"

Figure 10.27: 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.28):

NDVI per "Nafa" for five months

Figure 10.28: NDVI per “Nafa” for five months


  1. If necessary, we can convert those pixels that have NA in all layers to polygons, too, by specifying na.rm=FALSE.

  2. If necessary, we can convert those pixels that have NA in all layers to points, too, by specifying na.rm=FALSE.

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

  4. 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.