Vector data manipulation¶
This chapter illustrates some ways in which we can manipulate vector data. We start with an example SpatVector that we read from a shapefile.
library(terra)
## terra version 1.0.10
##
## Attaching package: 'terra'
## The following object is masked from 'package:knitr':
##
## spin
f <- system.file("ex/lux.shp", package="terra")
p <- vect(f)
p
## class : SpatVector
## geometry : polygons
## dimensions : 12, 5 (geometries, attributes)
## 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
## type : <num> <chr> <num> <chr> <num>
## values : 1 Diekirch 1 Clervaux 312
## 1 Diekirch 2 Diekirch 218
## 1 Diekirch 3 Redange 259
We can plot it in many ways
plot(p, "NAME_2")
Basics¶
Geometry and attributes¶
To extract the attributes (data.frame) from a SpatVector, 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
You can also extracting the geometry (rarely needed) as a a matrix.
g <- geom(p)
head(g)
## geom 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
Or as “well-known-text”
g <- geom(p, wkt=TRUE)
substr(g, 1, 50)
## [1] "POLYGON ((6.026519 50.17767, 6.031361 50.165627, 6"
## [2] "POLYGON ((6.178368 49.876823, 6.185479 49.870525, "
## [3] "POLYGON ((5.881378 49.870148, 5.881672 49.868866, "
## [4] "POLYGON ((6.131309 49.972565, 6.134291 49.972382, "
## [5] "POLYGON ((5.977929 50.026016, 5.982312 50.022949, "
## [6] "POLYGON ((6.385532 49.837029, 6.3886 49.833683, 6."
## [7] "POLYGON ((6.316665 49.623375, 6.31835 49.623157, 6"
## [8] "POLYGON ((6.425158 49.731644, 6.42657 49.73082, 6."
## [9] "POLYGON ((5.998312 49.699924, 5.998632 49.698559, "
## [10] "POLYGON ((6.039474 49.448261, 6.036906 49.448696, "
## [11] "POLYGON ((6.155963 49.685047, 6.159284 49.685036, "
## [12] "POLYGON ((6.067982 49.828465, 6.071922 49.825478, "
Variables¶
You can extracting a variable as you would do with a data.frame.
p$NAME_2
## [1] "Clervaux" "Diekirch" "Redange" "Vianden"
## [5] "Wiltz" "Echternach" "Remich" "Grevenmacher"
## [9] "Capellen" "Esch-sur-Alzette" "Luxembourg" "Mersch"
To sub-setting a SpatVector to one or more variables you can use the notation below. 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
## dimensions : 12, 1 (geometries, attributes)
## extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : NAME_2
## type : <chr>
## values : Clervaux
## Diekirch
## Redange
You can adding a new variable as with a data.frame.
set.seed(0)
p$lets <- sample(letters, nrow(p))
p
## class : SpatVector
## geometry : polygons
## dimensions : 12, 6 (geometries, attributes)
## 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 lets
## type : <num> <chr> <num> <chr> <num> <chr>
## values : 1 Diekirch 1 Clervaux 312 n
## 1 Diekirch 2 Diekirch 218 y
## 1 Diekirch 3 Redange 259 d
Note that to get the number of geometries of SpatVector p, you 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] 117100.12 93477.28 84502.45 44919.14 85032.61 74708.05 57991.42
## [8] 81203.93 74443.82 95564.74 80618.76 70810.67
Assigning a new value to an existing variable.
p$lets <- sample(LETTERS, nrow(p))
head(p)
## ID_1 NAME_1 ID_2 NAME_2 AREA lets
## 1 1 Diekirch 1 Clervaux 312 J
## 2 1 Diekirch 2 Diekirch 218 V
## 3 1 Diekirch 3 Redange 259 N
## 4 1 Diekirch 4 Vianden 76 Z
## 5 1 Diekirch 5 Wiltz 263 G
## 6 2 Grevenmacher 6 Echternach 188 I
To get rid of a variable, set it to NULL.
p$lets <- NULL
Merge¶
You can assign an attributes table (data.frame) to a SpatVector with
values<-. To add attributes to a SpatVector that already has
attributes use 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
## dimensions : 50, 6 (geometries, attributes)
## 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
## type : <chr> <chr> <num> <num> <num> <num>
## values : Diekirch Diekirch 1 1 312 406
## Diekirch Diekirch 1 1 312 534
## Diekirch Diekirch 1 1 312 640
head(pm)
## NAME_1 NAME_2 ID_1 ID_2 AREA Value
## 1 Diekirch Diekirch 1 1 312 406
## 2 Diekirch Diekirch 1 1 312 534
## 3 Diekirch Diekirch 1 1 312 640
## 4 Diekirch Diekirch 1 1 312 544
## 5 Diekirch Diekirch 1 1 312 268
## 6 Diekirch Diekirch 1 2 218 406
Note the new variable Value added to pm
Records¶
Selecting rows (records).
i <- which(p$NAME_1 == 'Grevenmacher')
g <- p[i,]
g
## class : SpatVector
## geometry : polygons
## dimensions : 3, 5 (geometries, attributes)
## extent : 6.169137, 6.528252, 49.46498, 49.85403 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : ID_1 NAME_1 ID_2 NAME_2 AREA
## type : <num> <chr> <num> <chr> <num>
## values : 2 Grevenmacher 6 Grevenmacher 188
## 2 Grevenmacher 7 Grevenmacher 129
## 2 Grevenmacher 12 Grevenmacher 210
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 geometries 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)
dim(z) <- c(2,2)
values(z) <- 1:4
names(z) <- 'Zone'
# coerce SpatRaster to SpatVector polygons
z <- as.polygons(z)
z
## class : SpatVector
## geometry : polygons
## dimensions : 4, 1 (geometries, attributes)
## extent : 5.74414, 6.529321, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : Zone
## type : <num>
## values : 1
## 2
## 3
z2 <- z[2,]
plot(p)
plot(z, add=TRUE, border='blue', lwd=5)
plot(z2, add=TRUE, border='red', lwd=2, col='red')
To append SpatVector objects of the same (vector) type you can use c
b <- c(p, z)
head(b)
## ID_1 NAME_1 ID_2 NAME_2 AREA Zone
## 1 1 Diekirch 1 Diekirch 312 NA
## 2 1 Diekirch 2 Diekirch 218 NA
## 3 1 Diekirch 3 Diekirch 259 NA
## 4 1 Diekirch 4 Diekirch 76 NA
## 5 1 Diekirch 5 Diekirch 263 NA
## 6 2 Grevenmacher 6 Grevenmacher 188 NA
tail(b)
## ID_1 NAME_1 ID_2 NAME_2 AREA Zone
## 11 3 Luxembourg 10 Luxembourg 237 NA
## 12 3 Luxembourg 11 Luxembourg 233 NA
## 13 NaN <NA> NaN <NA> NaN 1
## 14 NaN <NA> NaN <NA> NaN 2
## 15 NaN <NA> NaN <NA> NaN 3
## 16 NaN <NA> NaN <NA> NaN 4
Note how c allows you to append SpatVect objects with different
attribute names.
Aggregate¶
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
## class : SpatVector
## geometry : polygons
## dimensions : 1, 0 (geometries, attributes)
## extent : 5.74414, 6.529321, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
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
## class : SpatVector
## geometry : polygons
## dimensions : 4, 0 (geometries, attributes)
## extent : 5.74414, 6.529321, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
Overlay¶
Intersect¶
Intersect SpatVectors
i <- intersect(p, z2)
plot(i)
You can also intersect with a SpatExtent (rectangle).
e <- ext(6, 6.4, 49.7, 50)
pe <- crop(p, e)
plot(p)
plot(e, add=TRUE, lwd=3, col="red")
plot(pe, col='light blue', add=TRUE)
plot(e, add=TRUE, lwd=3, border="blue")
Union¶
Get the union of two SpatVectors.
u <- union(p, z)
u
## class : SpatVector
## geometry : polygons
## dimensions : 28, 7 (geometries, attributes)
## extent : 5.74414, 6.529321, 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 Zone aggregate_by_variable
## type : <num> <chr> <num> <chr> <num> <num> <num>
## values : 1 Diekirch 1 Diekirch 312 1 1
## 1 Diekirch 1 Diekirch 312 2 1
## 1 Diekirch 2 Diekirch 218 1 1
Note that there are many more polygons now. One for each unique combination of polygons (and attributes in this case).
set.seed(5)
plot(u, col=sample(rainbow(length(u))))
Cover¶
Cover is a combination of intersect and union
cov <- cover(p, z[c(1,4),])
cov
## class : SpatVector
## geometry : polygons
## dimensions : 11, 7 (geometries, attributes)
## extent : 5.74414, 6.529321, 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 Zone aggregate_by_variable
## type : <num> <chr> <num> <chr> <num> <num> <num>
## values : 1 Diekirch 1 Diekirch 312 NA NaN
## 1 Diekirch 2 Diekirch 218 NA NaN
## 1 Diekirch 3 Diekirch 259 NA NaN
plot(cov)
Difference¶
The symmetrical difference of two SpatVectors
dif <- symdif(z,p)
plot(dif, col=rainbow(length(dif)))
dif
## class : SpatVector
## geometry : polygons
## dimensions : 4, 2 (geometries, attributes)
## extent : 5.74414, 6.529321, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : Zone aggregate_by_variable
## type : <num> <num>
## values : 1 1
## 2 1
## 3 1
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)
extract is used for queries between SpatVector and SpatRaster
objects, and also for queries between SpatVectors.
extract(spts, p)
## id.x ID_1 NAME_1 ID_2 NAME_2 AREA
## 1 1 1 Diekirch 5 Diekirch 263
## 2 2 1 Diekirch 2 Diekirch 218
## 3 3 1 Diekirch 3 Diekirch 259
## 4 4 NA <NA> NA <NA> NA
## 5 5 NA <NA> NA <NA> NA
extract(spts, z)
## id.x Zone aggregate_by_variable
## 1 1 1 1
## 2 2 1 1
## 3 3 3 1
## 4 4 NA NA
## 5 5 4 1