Some packages for working with vector and raster data in R

In this tutorial we will use the sf package. This package focuses solely on vector data. It provides a standardized encoding of vector data and uses GDAL to read and write data, GEOS for geometrical operations and PROJ for projection conversions and datum transformations.

The GDAL library is well-documented (http://gdal.org/), but with a catch for R and Python programmers. The GDAL (and associated OGR) library and command line tools are all written in C and C++. Bindings are available that allow access from a variety of other languages including R and Python but the documentation is all written for the C++ version of the libraries. This can make reading the documentation rather challenging. Fortunately, the sf package, providing GDAL bindings in R, is also well documented with lots of examples. The same is valid for the Python libaries.

Thus, functionality that you commonly find in expensive GIS software is also available within R, using free but very powerful software libraries. Here is handy ‘cheatsheet’ for spatial operations with sf. The functions of the sf package are prefixed by st_, short for ‘spatial type’.

The possibilities are huge. In this course we can only scratch the surface with some essentials, which hopefully invite you to experiment further and use them in your research. Details can be found in the book Applied Spatial Data Analysis with R and several vignettes authored by Roger Bivand, Edzer Pebesma and Virgilio Gomez-Rubio.

Raster and vector integration and conversion

Historically, the first package handling vectors in R was the sp package, and the first package handling rasters was the raster package. They had cross-integration so that you could perform operations such as cropping a raster by a vector or extracting raster information over a vector location. However, the sp package got deprecated and sf was made as its successor, with much easier handling of the data (as a regular data.frame). Many of the packages that previously handled sp objects, including raster, got updated to also handle sf objects. Next, the stars package was developed by the creators of the sf package as a means of having multidimensional rasters. However, the older raster package was still the go-to solution for raster handling in R. Finally, as you noticed in a previous tutorial, in 2020 the raster package was deprecated and replaced with the terra package, which is a much faster C++ version of raster, but it also includes its own definition of vector data.

The result is that some packages provide or work with sp, sf and terra vectors, and some provide or work with raster, stars and terra rasters. This can be quite confusing. Therefore, in the course we currently use only the sf package for handling vectors and the terra package for handling rasters. The workflow is then to convert any object that is not sf or terra into sf or terra, do processing in sf and/or terra, and optionally convert the objects back if integration with some other package is necessary (e.g. sf to SpatVector for use in terra in some cases). Below is a matrix showing how to convert the various objects.

The matrix above shows lossless conversions, also known as casting or coercing data types. You can also convert data from raster to vector format and vice-versa. However, whenever you start converting between rasters and vectors, you should wonder whether you are taking the right approach to solve your problem. An approach that does not involve converting your data from vector to raster or the opposite should almost always be preferred.

As a result, because these functions are only useful for some very particular situations, we only give a brief description of them below.

• Vector to raster:
  • There is one function that allows to convert a vector object to a raster object. It is the rasterize() function.
• Raster to vector:
  • Three functions allow to convert raster data to vector: as.points(), as.contour() and as.polygons() produce SpatVector (terra) objects. The conversion to polygons can be useful to convert the result of a classification. In that case, set dissolve = to TRUE so that the polygons with the same attribute value will be dissolved into multi-polygon regions.

Geometric operations

Raw raster data do not usually conform to any notion of administrative or geographical boundaries. Vector data (and extents) can be used to mask or crop data to a desired region of interest.

• Crop:
  • Cropping consists in reducing the extent of a spatial object to a smaller extent. As a result, the output of crop() will automatically be rectangular and will not consider any features such as polygons to perform the subsetting. It is often useful to crop data as tightly as possible to the area under investigation to reduce the amount of data and have a more focused view when visualizing the data.
  • crop() uses objects of class SpatExtent to define the new extent, or any object that can be coerced to an extent (see ?ext for more info on this). This means that practically all spatial objects (raster or vector) can be used directly in crop. Considering two rasters r1 and r2 with r2 smaller than r1, you can simply use crop(r1, r2) in order to crop r1 to the extent of r2.
  • You can easily define an extent interactively (by clicking) thanks to the draw() function.
• Mask:
  • mask() can be used with almost all spatial objects to mask (= set to NA) values of a raster object. When used with a polygon object, mask() will keep values of the raster overlayed by polygons and mask the values outside of polygons.
  • Note the very useful inverse = argument of mask(), which allows to mask the inverse of the area covered by the features. We will use this feature of mask later in the tutorial to exclude water areas of a raster, defined in an independent polygon object.
• Buffer:
  • buffer() Calculate a buffer around all cells that are not NA in a SpatRaster, or around the geometries of a SpatVector (currently only implemented for points).
  • Note that the distance unit of the buffer width parameter is meters if the CRS is (+proj=longlat), and in map units (typically also meters) if not.
  • See ?buffer for more info.
• Extract:
  • The most common operation when combining vector and raster data is the extraction, using the function extract(). It simply consists in extracting the values of a raster object for locations specified by a vector object (e.g. points).
  • When using extract() with polygons or lines, individual features of the vector layer may overlay or intersect several pixels. In that case a function (fun =) can be used to summarize the values into one. Note that although most often the functions min, max, mean or median are used for the spatial aggregation, also any custom-made function can be used. The result is a matrix of raster values at the given locations. See ?extract() example section.