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)
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')
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))))
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)