knitr::opts_chunk$set(collapse = TRUE)
This vignetted describes how simple features, i.e. records that come with a geometry, can be manipulated, for the case where these manipulations involve geometries. Manipulations include:
Features are represented by records in an
sf object, and have feature attributes (all non-geometry fields) and feature geometry. Since
sf objects are a subclass of
tbl_df, operations on feature attributes work identical to how they work on
library(sf) nc <- st_read(system.file("shape/nc.shp", package="sf")) nc <- st_transform(nc, 2264) nc[1,]
prints the first record.
Many of the tidyverse/dplyr verbs have methods for
sf objects. This means that given that both
dplyr are loaded, manipulations such as selecting a single attribute will return an
library(dplyr) nc %>% select(NWBIR74) %>% head(2)
which implies that the geometry is sticky, and gets added automatically. If we want to drop geometry, we can coerce to
data.frame first, this drops geometry list-columns:
nc %>% as.data.frame %>% select(NWBIR74) %>% head(2)
We can subset feature sets by using the square bracket notation
and use the
drop argument to drop geometries:
nc[1, "NWBIR74", drop = TRUE]
but we can also use a spatial object as the row selector, to select features that intersect with another spatial feature:
Ashe = nc[nc$NAME == "Ashe",] class(Ashe) nc[Ashe,]
we see that in the result set
Ashe is included, as the default value for argument
Ashe intersects with itself. We could exclude self-intersection by using predicate
st_touches (overlapping features don't touch):
Ashe = nc[nc$NAME == "Ashe",] nc[Ashe, op = st_touches]
dplyr, we can do the same by calling the predicate directly:
nc %>% filter(lengths(st_touches(., Ashe)) > 0)
Suppose we want to compare the 1974 fraction of SID (sudden infant death) of the counties that intersect with
Ashe to the remaining ones. We can do this by
a <- aggregate(nc[, c("SID74", "BIR74")], list(Ashe_nb = lengths(st_intersects(nc, Ashe)) > 0), sum) (a <- a %>% mutate(frac74 = SID74 / BIR74) %>% select(frac74)) plot(a, col = c(grey(.8), grey(.5))) plot(st_geometry(Ashe), border = '#ff8888', add = TRUE, lwd = 2)
The usual join verbs base R (
merge) and of dplyr (
left_join, etc) work for
sf objects as well; the joining takes place on attributes (ignoring geometries). In case of no matching geometry, an empty geometry is substituted. The second argument should be a
data.frame (or similar), not an
x = st_sf(a = 1:2, geom = st_sfc(st_point(c(0,0)), st_point(c(1,1)))) y = data.frame(a = 2:3) merge(x, y) merge(x, y, all = TRUE) right_join(x, y)
For joining based on spatial intersections (of any kind),
st_join is used:
x = st_sf(a = 1:3, geom = st_sfc(st_point(c(1,1)), st_point(c(2,2)), st_point(c(3,3)))) y = st_buffer(x, 0.1) x = x[1:2,] y = y[2:3,] plot(st_geometry(x), xlim = c(.5, 3.5)) plot(st_geometry(y), add = TRUE)
The join method is a left join, retaining all records of the first attribute:
st_join(x, y) st_join(y, x)
and the geometry retained is that of the first argument.
The spatial join predicate can be controlled with any function compatible to
st_intersects (the default), e.g.
st_join(x, y, join = st_covers) # no matching y records: points don't cover circles st_join(y, x, join = st_covers) # matches for those circles covering a point
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.