Vector data manipulation

Example SpatialPolygons

library(terra)
f <- system.file("exdata/lux.shp", package="terra")
p <- vect(f)
p
## class       : SpatVector
## geometry    : polygons
## elements    : 12
## extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : ID_1, NAME_1, ID_2, NAME_2, AREA
par(mai=c(0,0,0,0))
plot(p)

image0

Basics

Basic operations are pretty much like working with a data.frame.

Geometry and attributes

To extract the attributes (data.frame) from a Spatial object, use:

d <- as.data.frame(p)
head(d)
##   ID_1       NAME_1 ID_2     NAME_2 AREA
## 1    1     Diekirch    1   Clervaux  312
## 2    1     Diekirch    2   Diekirch  218
## 3    1     Diekirch    3    Redange  259
## 4    1     Diekirch    4    Vianden   76
## 5    1     Diekirch    5      Wiltz  263
## 6    2 Grevenmacher    6 Echternach  188

Extracting geometry (rarely needed).

g <- geom(p)
head(g)
##   id part        x        y hole
## 1  1    1 6.026519 50.17767    0
## 2  1    1 6.031361 50.16563    0
## 3  1    1 6.035646 50.16410    0
## 4  1    1 6.042747 50.16157    0
## 5  1    1 6.043894 50.16116    0
## 6  1    1 6.048243 50.16008    0

Variables

Extracting a variable.

p$NAME_2
##  [1] "Clervaux"         "Diekirch"         "Redange"          "Vianden"
##  [5] "Wiltz"            "Echternach"       "Remich"           "Grevenmacher"
##  [9] "Capellen"         "Esch-sur-Alzette" "Luxembourg"       "Mersch"

Sub-setting by variable. Note how this is different from the above example. Above a vector of values is returned. With the approach below you get a new SpatVector with only one variable.

p[, 'NAME_2']
## class       : SpatVector
## geometry    : polygons
## elements    : 12
## extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : NAME_2

Adding a new variable.

set.seed(0)
p$new <- sample(letters, nrow(p))
p
## class       : SpatVector
## geometry    : polygons
## elements    : 12
## extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : ID_1, NAME_1, ID_2, NAME_2, AREA, new

Note that to get the number of spatial object (“features”), of SpatVector p, you need can use nrow(p), or size(p). You can also do perimeter(p) to get the “length” of the spatial objects (zero for points, the length of the lines, or the perimeter of the polygons).

perimeter(p)
##  [1] 1.3413234 1.0656032 0.9548349 0.4965615 0.9641204 0.8743230 0.6520187
##  [8] 0.9584838 0.8589857 1.1414055 0.9473256 0.8280279

Assigning a new value to an existing variable.

p$new <- sample(LETTERS, nrow(p))
p
## class       : SpatVector
## geometry    : polygons
## elements    : 12
## extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : ID_1, NAME_1, ID_2, NAME_2, AREA, new

To get rid of a variable.

p$new <- NULL

Merge

You can join a table (data.frame) with a SpatialVector object with merge.

dfr <- data.frame(District=p$NAME_1, Canton=p$NAME_2, Value=round(runif(length(p), 100, 1000)))
dfr <- dfr[order(dfr$Canton), ]
pm <- merge(p, dfr, by.x=c('NAME_1', 'NAME_2'), by.y=c('District', 'Canton'))
pm
## class       : SpatVector
## geometry    : polygons
## elements    : 12
## extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names       : NAME_1, NAME_2, ID_1, ID_2, AREA, Value
head(pm)
##         NAME_1     NAME_2 ID_1 ID_2 AREA Value
## 1     Diekirch   Clervaux    1    1  312   406
## 2     Diekirch   Diekirch    1    2  218   534
## 3     Diekirch    Redange    1    3  259   640
## 4     Diekirch    Vianden    1    4   76   544
## 5     Diekirch      Wiltz    1    5  263   268
## 6 Grevenmacher Echternach    2    6  188   845

Note that new variable Value added to pm

Records

Selecting rows (records).

i <- which(p$NAME_1 == 'Grevenmacher')
g <- p[i,]
g
## class       : SpatVector
## geometry    : polygons
## elements    : 3
## extent      : 6.169137, 6.528252, 49.46498, 49.85403  (xmin, xmax, ymin, ymax)
## coord. ref. :
## names       : ID_1, NAME_1, ID_2, NAME_2, AREA

It is also possible to interactively select and query records by clicking on a plotted dataset. That is difficult to show here. See ?select for interactively selecting spatial features and ?click to identify attributes by clicking on a plot (map).

Append and aggregate

Append

More example data. Object z, consisting of four polygons, and z2 which is one of these four polygons.

z <- rast(p, nrow=2, ncol=2, vals=1:4)
names(z) <- 'Zone'
# coerce SpatRaster to SpatVector polygons
z <- as.polygons(z)
z
## class       : SpatVector
## geometry    : polygons
## elements    : 4
## extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
z2 <- z[2,]
plot(p)
plot(z, add=TRUE, border='blue', lwd=5)
plot(z2, add=TRUE, border='red', lwd=2, col='red')

image1

To append SpatVector objects of the same (vector) type you can use bind

# not implemented yet
#b <- bind(p, z)
#head(b)
#tail(b)

Note how bind allows you to append SpatVect objects with different attribute names.

Aggregate

#not implemented yet
#pa <- aggregate(p, by='NAME_1')
#za <- aggregate(z)
#plot(za, col='light gray', border='light gray', lwd=5)
#plot(pa, add=TRUE, col=rainbow(3), lwd=3, border='white')

Aggregate without dissolve

#zag <- aggregate(z, dissolve=FALSE)
#zag
#plot(zag, col="light gray")

This is a structure that is similar to what you may get for an archipelago: multiple polygons represented as one entity (one row). Use disaggregate to split these up into their parts.

#zd <- disaggregate(zag)
#zd

Overlay

Erase

Erase a part of a SpatialPolygons object

#e <- erase(p, z2)

This is equivalent to

#e <- p - z2
#plot(e)

Intersect

Intersect SpatialPolygons

#i <- intersect(p, z2)
#plot(i)

This is equivalent to

#i <- p * z2

You can also intersect with a SpatExtent (rectangle).

#e <- ext(6, 6.4, 49.7, 50)
#pe <- crop(p, e)
#plot(p)
#plot(pe, col='light blue', add=TRUE)
#plot(e, add=TRUE, lwd=3, col='red')

Union

Get the union of two SpatialPolygon* objects.

#u <- union(p, z)

This is equivalent to

#u <- p + z

Note that there are many more polygons now. One for each unique combination of polygons (and attributes in this case).

#u
#set.seed(5)
#plot(u, col=sample(rainbow(length(u))))

Cover

Cover is a combination of intersect and union

#cov <- cover(p, z)
#cov
#plot(cov)

Difference

The symmetrical difference of two SpatialPolygons* objects

#dif <- symdif(z,p)
#plot(dif, col=rainbow(length(dif)))
#dif

Spatial queries

Query polygons with points.

#pts <- matrix(c(6, 6.1, 5.9, 5.7, 6.4, 50, 49.9, 49.8, 49.7, 49.5), ncol=2)
#spts <- vect(pts, proj4string=crs(p))
#plot(z, col='light blue', lwd=2)
#points(spts, col='light gray', pch=20, cex=6)
#text(spts, 1:nrow(pts), col='red', font=2, cex=1.5)
#lines(p, col='blue', lwd=2)

Use over for queries between SpatVector objects

#over(spts, p)
#over(spts, z)

extract is generally used for queries between SpatVector and SpatRaster objects, but it can also be used here.

#extract(z, pts)