R in Space - Geometry manipulation

April 8, 2018

  R in Space R Spatial

  raster rgeos sf sp

David Beauchesne   Marie-Hélène Brice   Nicolas Casajus   Kevin Cazelles   Elliot Dreujou   Steve Vissault  


To manipulate geometries, there are few options. For spatial objects defined with sp, there is a specific package to manipulate vector objects: rgeos. The equivalent of most common functions in sp are found in sf. Finally, there are a set of interesting functions in raster to handle the manipulation of rasters with vector objects.

Let’s first load the packages required and download a simple shapefile of Belgium.

library(sp)
library(sf)
library(rgeos)
library(raster)
dir.create("data")
bel2 <- getData(name="GADM", country="BEL", path="data", level = 2)
plot(bel2)

sp objects

Subset

A quick subset:

## Create a subset (A sub-shapefile)
slc <- c(8,11,12)
bel_slc <- bel2[slc,]
class(bel_slc)
#R> [1] "SpatialPolygonsDataFrame"
#R> attr(,"package")
#R> [1] "sp"

Let’s plot our selection:

plot(bel_slc)

Unions

Let’s do the union of a selection of spatial polygons:

bel_south <- gUnionCascaded(bel_slc)
bel_north <- gUnionCascaded(bel2[-slc,])
bel_one <- gUnionCascaded(bel2)
par(mfrow=c(1,3), mar=c(1,1,1,1))
plot(bel_south)
plot(bel_north)
plot(bel_one)

Let’s combine them on one plot:

plot(bel2, lwd = .1)
plot(bel_south, add=T, col=2)
plot(bel_north, add=T, col=4)
plot(bel2, lwd = .5, add = T)
plot(bel_one, lwd = 2, add = T, border = "grey55")

Buffers

Buffer must be done on planar coordinates.

plot(bel2)
bufs <- gBuffer(bel_south, width=0.4)
bufn <- gBuffer(bel_north, width=0.1)
plot(bufs, add=T, lwd=2, lty=2)
plot(bufn, add=T, lwd=3)

Intersections

Intersections between bel_south and bel_north:

par(mfrow = c(1,2))
plot(bel_one)
plot(gIntersection(bel_south, bel_north), col=5, add=T, lwd=2)
plot(bel_one)
plot(gIntersection(bufs, bufn), col=5, add=T, lwd=2)

Differences

bel_diff <- gDifference(bufs, bufn)
plot(bel_one)
plot(bel_diff, add=T, lwd=2, lty=3, col=2)
plot(bel_south, add=T, col=2)

Overlays

One very useful tool is the over() function which provides a consistent spatial overlay. For instance, here, I create a set of random points and look which one are in bel2, bel_north and bel_south:

pts <- SpatialPoints(
  coords = cbind(
    runif(50, bel2@bbox[1, 1], bel2@bbox[1, 2]),
    runif(50, bel2@bbox[2, 1], bel2@bbox[2, 2])
  ),
  proj4string = bel2@proj4string
)
# Make a SpatialPolygons
bel2_geom <- SpatialPolygons(bel2@polygons, proj4string = bel2@proj4string)
# plots
plot(bel2_geom)
points(pts, col = is.na(over(pts, bel2_geom))+2, pch =19)

sf objects

Have a look at the dedicated vignette: https://cran.r-project.org/web/packages/sf/vignettes/sf3.html

raster objects

Let’s now download a simple raster of elevation.

bel_elv <- getData(name="alt", country="BEL", path="data", level = 2)
plot(bel_elv)

Extract values

val_south <- extract(bel_elv, bel_south)
val_north <- extract(bel_elv, bel_north)
#
par(mfrow = c(1,2))
hist(val_south[[1]])
hist(val_north[[1]])

Mask

elv_south <- rasterize(bel_south, bel_elv, mask = T)
elv_north <- rasterize(bel_north, bel_elv, mask = T)
par(mfrow = c(1,2))
plot(elv_south)
plot(elv_north)