Example: Workflow in R

This example provides a schematic workflow for processing vector and raster data in R.

Get raster data

Firstly, we import some raster data into our working environment Therefore, we need to load a package to handle raster data in R, preferable raster. If the package is not available, we need to install it first with install.packages("raster").

library("raster")

We now can use the function getData() to download some raster data: In this example a global map of precipitation values at 10 minutes spatial resolution.

prec <- getData("worldclim", var="prec", res=10)

Fortunately, the downloaded data already have a correct CRS:

## class      : RasterStack 
## dimensions : 900, 2160, 1944000, 12  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names      : prec1, prec2, prec3, prec4, prec5, prec6, prec7, prec8, prec9, prec10, prec11, prec12 
## min values :     0,     0,     0,     0,     0,     0,     0,     0,     0,      0,      0,      0 
## max values :   885,   736,   719,   820,   955,  1850,  2088,  1386,   904,    980,    893,    914 

… and can be quickly and simply visualized with plot(). Note that the object type is a RasterStack with 12 layers, one for each month of the year.

plot(prec$prec1)

Get vector data

Secondly, we add some vector data to our working environment. For example the administrative boundaries of France at the country level:

fra <- getData('GADM', country='FRA', level=0)

Fortunately, also these downloaded data already have a CRS, defined by a proj4 string:

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -5.143751, 9.560416, 41.33375, 51.0894  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : GID_0, NAME_0 
## value       :   FRA, France 

Note that the object type here is a SpatialPolygonsDataFrame (defined in package sp) with one feature (i.e. with a single polygon), a certain extent (which can also be extracted with extent(fra)), a CRS (which can be extracted with CRS(fra)), and two variables with some values (in this case country abbreviation and name).

Also vector data can quickly and simply be visualized with plot()

plot(fra)

Set extent

We can use the extent of one spatial object to crop (i.e. to cut out) another spatial object. In this example, we will crop the raster map(s) with the extent defined in the vector object. This is going to work because both objects have the same CRS. Note that crop() processes all layers of the input raster stack.

cropped_prec <- crop(prec, extent(fra))

For function arguments see ?crop(). Now we have precipitation maps of France:

plot(cropped_prec$prec1)

Vector operations

A simple vector operation would be to clip a spatial object not by the extent of another spatial object but by features or polygons of any shape.

We will now use the country boundary of France for clipping of the already cropped precipitation maps with the mask() function of the raster package:

clipped_prec <- mask(cropped_prec, fra)

The result is a RasterBrick object with 12 layers, one for each month of the year (the difference to a RasterStack as obove is only how memory is allocated, what does not matter here).

Again, the result can quickly and simply be visualized with plot()

plot(clipped_prec$prec1)

Raster operations

The created RasterBrick now contains precipitation values of France on a monthly basis. What if we want to have a single precipitation layer with annual precipitation values? We would need to sum up the values of all 12 precipitation layers for each pixel location. This can be done by:

clipped_prec_sum <- sum(clipped_prec)

An alternative would be to use stackApply() (see ?stackApply() for details. Note the na.rm argument).

clipped_prec_sum_2 <- raster::stackApply(clipped_prec, rep(1,12), sum, na.rm=FALSE)

The resulting raster map looks like this:

plot(clipped_prec_sum)

Note the different value range, which now stands for the annual amount of precipitation (in mm).

Mapping

We now combine all the above created spatial objects to create a single simple map:

plot(clipped_prec_sum)
plot(fra, add=T)
points(2.349014, 48.864716, pch=8, cex=2) # roughly the location of Paris

Write out

The above created overlay of maps can be written to file as an ordinary image with e.g.

jpeg("FirstSimpleMap.jpg")
plot(clipped_prec_sum)
plot(fra, add=T)
points(2.349014, 48.864716, pch=8, cex=2)
dev.off()

Note that this is a raster image without geographic information. If you want to write out the spatial objects with geographic information, use e.g. raster::writeRaster() or raster::shapefile().

Other important functions

  • Reading in raster data from file: raster::raster()
  • Reading in vector data from file: rgdal::readOGR()

More information

For more details see www.rspatial.org and Geocomputation with R

Note that the packages rgdal, rgeos and maptools, upon which many other spatial packages depend, will retire by the end of 2023. So it will be wise to switch to the packages sf, stars, and terra in due time. For details, see this blog.