2.2 Vector data
Take care when using the word ‘vector’ as it can have two meanings in this book:geographic vector data and the vector
class (note the monospace
font) in R.The former is a data model, the latter is an R class just like data.frame
and matrix
.Still, there is a link between the two: the spatial coordinates which are at the heart of the geographic vector data model can be represented in R using vector
objects.
The geographic vector model is based on points located within a coordinate reference system (CRS).Points can represent self-standing features (e.g., the location of a bus stop) or they can be linked together to form more complex geometries such as lines and polygons.Most point geometries contain only two dimensions (3-dimensional CRSs contain an additional (z) value, typically representing height above sea level).
In this system London, for example, can be represented by the coordinates c(-0.1, 51.5)
.This means that its location is -0.1 degrees east and 51.5 degrees north of the origin.The origin in this case is at 0 degrees longitude (the Prime Meridian) and 0 degree latitude (the Equator) in a geographic (‘lon/lat’) CRS (Figure 2.1, left panel).The same point could also be approximated in a projected CRS with ‘Easting/Northing’ values of c(530000, 180000)
in the British National Grid, meaning that London is located 530 km East and 180 km North of the (origin) of the CRS.This can be verified visually: slightly more than 5 ‘boxes’ — square areas bounded by the gray grid lines 100 km in width — separate the point representing London from the origin (Figure 2.1, right panel).
The location of National Grid’s origin, in the sea beyond South West Peninsular, ensures that most locations in the UK have positive Easting and Northing values.6There is more to CRSs, as described in Sections 2.4 and 6 but, for the purposes of this section, it is sufficient to know that coordinates consist of two numbers representing distance from an origin, usually in (x) then (y) dimensions.
Figure 2.1: Illustration of vector (point) data in which location of London (the red X) is represented with reference to an origin (the blue circle). The left plot represents a geographic CRS with an origin at 0° longitude and latitude. The right plot represents a projected CRS with an origin located in the sea west of the South West Peninsula.
sf is a package providing a class system for geographic vector data.Not only does sf supersede sp, it also provides a consistent command-line interface to GEOS and GDAL, superseding rgeos and rgdal (described in Section 1.5).This section introduces sf classes in preparation for subsequent chapters (Chapters 5 and 7 cover the GEOS and GDAL interface, respectively).
2.2.1 An introduction to simple features
Simple features is an open standard developed and endorsed by the Open Geospatial Consortium (OGC), a not-for-profit organization whose activities we will revisit in a later chapter (in Section 7.5).Simple Features is a hierarchical data model that represents a wide range of geometry types.Of 17 geometry types supported by the specification, only 7 are used in the vast majority of geographic research (see Figure 2.2);these core geometry types are fully supported by the R package sf(E. Pebesma 2018).7
Figure 2.2: Simple feature types fully supported by sf.
sf can represent all common vector geometry types (raster data classes are not supported by sf): points, lines, polygons and their respective ‘multi’ versions (which group together features of the same type into a single feature).sf also supports geometry collections, which can contain multiple geometry types in a single object.sf largely supersedes the sp ecosystem, which comprises sp(E. Pebesma and Bivand 2018), rgdal for data read/write (Bivand, Keitt, and Rowlingson 2018) and rgeos for spatial operations (Bivand and Rundel 2018).The package is well documented, as can be seen on its website and in 6 vignettes, which can be loaded as follows:
vignette(package = "sf") # see which vignettes are available
vignette("sf1") # an introduction to the package
As the first vignette explains, simple feature objects in R are stored in a data frame, with geographic data occupying a special column, usually named ‘geom’ or ‘geometry’.We will use the world
dataset provided by the spData, loaded at the beginning of this chapter (see nowosad.github.io/spData for a list of datasets loaded by the package).world
is a spatial object containing spatial and attribute columns, the names of which are returned by the function names()
(the last column contains the geographic information):
names(world)
#> [1] "iso_a2" "name_long" "continent" "region_un" "subregion"
#> [6] "type" "area_km2" "pop" "lifeExp" "gdpPercap"
#> [11] "geom"
The contents of this geom
column give sf
objects their spatial powers: world$geom
is a ‘list column’ that contains all the coordinates of the country polygons.The sf package provides a plot()
method for visualizing geographic data:the following command creates Figure 2.3.
plot(world)
Figure 2.3: A spatial plot of the world using the sf package, with a facet for each attribute.
Note that instead of creating a single map, as most GIS programs would, the plot()
command has created multiple maps, one for each variable in the world
datasets.This behavior can be useful for exploring the spatial distribution of different variables and is discussed further in Section 2.2.3 below.
Being able to treat spatial objects as regular data frames with spatial powers has many advantages, especially if you are already used to working with data frames.The commonly used summary()
function, for example, provides a useful overview of the variables within the world
object.
summary(world["lifeExp"])
#> lifeExp geom
#> Min. :50.6 MULTIPOLYGON :177
#> 1st Qu.:65.0 epsg:4326 : 0
#> Median :72.9 +proj=long...: 0
#> Mean :70.9
#> 3rd Qu.:76.8
#> Max. :83.6
#> NA's :10
Although we have only selected one variable for the summary
command, it also outputs a report on the geometry.This demonstrates the ‘sticky’ behavior of the geometry columns of sf objects, meaning the geometry is kept unless the user deliberately removes them, as we’ll see in Section 3.2.The result provides a quick summary of both the non-spatial and spatial data contained in world
: the mean average life expectancy is 71 years (ranging from less than 51 to more than 83 years with a median of 73 years) across all countries.
The word MULTIPOLYGON
in the summary output above refers to the geometry type of features (countries) in the world
object.This representation is necessary for countries with islands such as Indonesia and Greece.Other geometry types are described in Section 2.2.5.
It is worth taking a deeper look at the basic behavior and contents of this simple feature object, which can usefully be thought of as a ‘spatial data frame’.
sf
objects are easy to subset.The code below shows its first two rows and three columns.The output shows two major differences compared with a regular data.frame
: the inclusion of additional geographic data (geometry type
, dimension
, bbox
and CRS information - epsg (SRID)
, proj4string
), and the presence of a geometry
column, here named geom
:
world_mini = world[1:2, 1:3]
world_mini
#> Simple feature collection with 2 features and 3 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: -180 ymin: -18.3 xmax: 180 ymax: -0.95
#> epsg (SRID): 4326
#> proj4string: +proj=longlat +datum=WGS84 +no_defs
#> iso_a2 name_long continent geom
#> 1 FJ Fiji Oceania MULTIPOLYGON (((180 -16.1, ...
#> 2 TZ Tanzania Africa MULTIPOLYGON (((33.9 -0.95,...
All this may seem rather complex, especially for a class system that is supposed to be simple.However, there are good reasons for organizing things this way and using sf.
Before describing each geometry type that the sf package supports, it is worth taking a step back to understand the building blocks of sf
objects.Section 2.2.8 shows how simple features objects are data frames, with special geometry columns.These spatial columns are often called geom
or geometry
: world$geom
refers to the spatial element of the world
object described above.These geometry columns are ‘list columns’ of class sfc
(see Section 2.2.7).In turn, sfc
objects are composed of one or more objects of class sfg
: simple feature geometries that we describe in Section 2.2.6.
To understand how the spatial components of simple features work, it is vital to understand simple feature geometries.For this reason we cover each currently supported simple features geometry type in Section 2.2.5 before moving on to describe how these can be represented in R using sfg
objects, which form the basis of sfc
and eventually full sf
objects.
The preceding code chunk uses =
to create a new object called world_mini
in the command world_mini = world[1:2, 1:3]
.This is called assignment.An equivalent command to achieve the same result is world_mini <- world[1:2, 1:3]
.Although ‘arrow assigment’ is more commonly used, we use ‘equals assignment’ because it’s slightly faster to type and easier to teach due to compatibility with commonly used languages such as Python and JavaScript.Which to use is largely a matter of preference as long as you’re consistent (packages such as styler can be used to change style).
2.2.2 Why simple features?
Simple features is a widely supported data model that underlies data structures in many GIS applications including QGIS and PostGIS.A major advantage of this is that using the data model ensures your work is cross-transferable to other set-ups, for example importing from and exporting to spatial databases.
A more specific question from an R perspective is “why use the sf package when sp is already tried and tested”?There are many reasons (linked to the advantages of the simple features model) including:
- Fast reading and writing of data.
- Enhanced plotting performance.
- sf objects can be treated as data frames in most operations.
- sf functions can be combined using
%>%
operator and works well with the tidyverse collection of R packages. - sf function names are relatively consistent and intuitive (all begin with
st_
). Due to such advantages, some spatial packages (including tmap, mapview and tidycensus) have added support for sf.However, it will take many years for most packages to transition and some will never switch.Fortunately, these can still be used in a workflow based onsf
objects, by converting them to theSpatial
class used in sp:
library(sp)
world_sp = as(world, Class = "Spatial")
# sp functions ...
Spatial
objects can be converted back to sf
in the same way or with st_as_sf()
:
world_sf = st_as_sf(world_sp, "sf")
2.2.3 Basic map making
Basic maps are created in sf with plot()
.By default this creates a multi-panel plot (like sp’s spplot()
), one sub-plot for each variable of the object, as illustrated in the left-hand panel in Figure 2.4.A legend or ‘key’ with a continuous color is produced if the object to be plotted has a single variable (see the right-hand panel).Colors can also be set with col =
, although this will not create a continuous palette or a legend.
plot(world[3:6])
plot(world["pop"])
Figure 2.4: Plotting with sf, with multiple variables (left) and a single variable (right).
Plots are added as layers to existing images by setting add = TRUE
.8To demonstrate this, and to provide a taster of content covered in Chapters 3 and 4 on attribute and spatial data operations, the subsequent code chunk combines countries in Asia:
world_asia = world[world$continent == "Asia", ]
asia = st_union(world_asia)
We can now plot the Asian continent over a map of the world.Note that the first plot must only have one facet for add = TRUE
to work.If the first plot has a key, reset = FALSE
must be used (result not shown):
plot(world["pop"], reset = FALSE)
plot(asia, add = TRUE, col = "red")
Adding layers in this way can be used to verify the geographic correspondence between layers:the plot()
function is fast to execute and requires few lines of code, but does not create interactive maps with a wide range of options.For more advanced map making we recommend using dedicated visualization packages such as tmap (see Chapter 8).
2.2.4 Base plot arguments
There are various ways to modify maps with sf’s plot()
method.Because sf extends base R plotting methods plot()
’s arguments such as main =
(which specifies the title of the map) work with sf
objects (see ?graphics::plot
and ?par
).9
Figure 2.5 illustrates this flexibility by overlaying circles, whose diameters (set with cex =
) represent country populations, on a map of the world.A basic version of the map can be created with the following commands (see exercises at the end of this chapter and the script 02-contplot.R
to create Figure 2.5):
plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)
plot(st_geometry(world_cents), add = TRUE, cex = cex)
Figure 2.5: Country continents (represented by fill color) and 2015 populations (represented by circles, with area proportional to population).
The code above uses the function st_centroid()
to convert one geometry type (polygons) to another (points) (see Chapter 5), the aesthetics of which are varied with the cex
argument.
sf’s plot method also has arguments specific to geographic data. expandBB
, for example, can be used plot an sf
object in context:it takes a numeric vector of length four that expands the bounding box of the plot relative to zero in the following order: bottom, left, top, right.This is used to plot India in the context of its giant Asian neighbors, with an emphasis on China to the east, in the following code chunk, which generates Figure 2.6 (see exercises below on adding text to plots):
india = world[world$name_long == "India", ]
plot(st_geometry(india), expandBB = c(0, 0.2, 0.1, 1), col = "gray", lwd = 3)
plot(world_asia[0], add = TRUE)
Figure 2.6: India in context, demonstrating the expandBB argument.
Note the use of [0]
to keep only the geometry column and lwd
to emphasize India.See Section 8.6 for other visualization techniques for representing a range of geometry types, the subject of the next section.
2.2.5 Geometry types
Geometries are the basic building blocks of simple features.Simple features in R can take on one of the 17 geometry types supported by the sf package.In this chapter we will focus on the seven most commonly used types: POINT
, LINESTRING
, POLYGON
, MULTIPOINT
, MULTILINESTRING
, MULTIPOLYGON
and GEOMETRYCOLLECTION
.Find the whole list of possible feature types in the PostGIS manual.
Generally, well-known binary (WKB) or well-known text (WKT) are the standard encoding for simple feature geometries.WKB representations are usually hexadecimal strings easily readable for computers.This is why GIS and spatial databases use WKB to transfer and store geometry objects.WKT, on the other hand, is a human-readable text markup description of simple features.Both formats are exchangeable, and if we present one, we will naturally choose the WKT representation.
The basis for each geometry type is the point.A point is simply a coordinate in 2D, 3D or 4D space (see vignette("sf1")
for more information) such as (see left panel in Figure 2.7):
POINT (5 2)
A linestring is a sequence of points with a straight line connecting the points, for example (see middle panel in Figure 2.7):LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)
A polygon is a sequence of points that form a closed, non-intersecting ring.Closed means that the first and the last point of a polygon have the same coordinates (see right panel in Figure 2.7).10Polygon without a hole:
POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))
Figure 2.7: Illustration of point, linestring and polygon geometries.
So far we have created geometries with only one geometric entity per feature.However, sf also allows multiple geometries to exist within a single feature (hence the term ‘geometry collection’) using “multi” version of each geometry type:
- Multipoint:
MULTIPOINT (5 2, 1 3, 3 4, 3 2)
- Multilinestring:
MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
- Multipolygon:
MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5), (0 2, 1 2, 1 3, 0 3, 0 2)))
Figure 2.8: Illustration of multi* geometries.
Finally, a geometry collection can contain any combination of geometries including (multi)points and linestrings (see Figure 2.9):
- Geometry collection:
GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2), LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))
Figure 2.9: Illustration of a geometry collection.
2.2.6 Simple feature geometries (sfg)
The sfg
class represents the different simple feature geometry types in R: point, linestring, polygon (and their ‘multi’ equivalents, such as multipoints) or geometry collection.
Usually you are spared the tedious task of creating geometries on your own since you can simply import an already existing spatial file.However, there are a set of functions to create simple feature geometry objects (sfg
) from scratch if needed.The names of these functions are simple and consistent, as they all start with the st_
prefix and end with the name of the geometry type in lowercase letters:
- A point:
st_point()
- A linestring:
st_linestring()
- A polygon:
st_polygon()
- A multipoint:
st_multipoint()
- A multilinestring:
st_multilinestring()
- A multipolygon:
st_multipolygon()
A geometry collection:
st_geometrycollection()
sfg
objects can be created from three base R data types:A numeric vector: a single point
- A matrix: a set of points, where each row represents a point, a multipoint or linestring
- A list: a collection of objects such as matrices, multilinestrings or geometry collections
The function
st_point()
creates single points from numeric vectors:
st_point(c(5, 2)) # XY point
#> POINT (5 2)
st_point(c(5, 2, 3)) # XYZ point
#> POINT Z (5 2 3)
st_point(c(5, 2, 1), dim = "XYM") # XYM point
#> POINT M (5 2 1)
st_point(c(5, 2, 3, 1)) # XYZM point
#> POINT ZM (5 2 3 1)
The results show that XY (2D coordinates), XYZ (3D coordinates) and XYZM (3D with an additional variable, typically measurement accuracy) point types are created from vectors of length 2, 3, and 4, respectively.The XYM type must be specified using the dim
argument (which is short for dimension).
By contrast, use matrices in the case of multipoint (st_multipoint()
) and linestring (st_linestring()
) objects:
# the rbind function simplifies the creation of matrices
## MULTIPOINT
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
st_multipoint(multipoint_matrix)
#> MULTIPOINT (5 2, 1 3, 3 4, 3 2)
## LINESTRING
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix)
#> LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)
Finally, use lists for the creation of multilinestrings, (multi-)polygons and geometry collections:
## POLYGON
polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
st_polygon(polygon_list)
#> POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))
## POLYGON with a hole
polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list)
#> POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5), (2 4, 3 4, 3 3, 2 3, 2 4))
## MULTILINESTRING
multilinestring_list = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)),
rbind(c(1, 2), c(2, 4)))
st_multilinestring((multilinestring_list))
#> MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
## MULTIPOLYGON
multipolygon_list = list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
st_multipolygon(multipolygon_list)
#> MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5)), ((0 2, 1 2, 1 3, 0 3, 0 2)))
## GEOMETRYCOLLECTION
gemetrycollection_list = list(st_multipoint(multipoint_matrix),
st_linestring(linestring_matrix))
st_geometrycollection(gemetrycollection_list)
#> GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2),
#> LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))
2.2.7 Simple feature columns (sfc)
One sfg
object contains only a single simple feature geometry.A simple feature geometry column (sfc
) is a list of sfg
objects, which is additionally able to contain information about the coordinate reference system in use.For instance, to combine two simple features into one object with two features, we can use the st_sfc()
function.This is important since sfc
represents the geometry column in sf data frames:
# sfc POINT
point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc
#> Geometry set for 2 features
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 1 ymin: 2 xmax: 5 ymax: 3
#> epsg (SRID): NA
#> proj4string: NA
#> POINT (5 2)
#> POINT (1 3)
In most cases, an sfc
object contains objects of the same geometry type.Therefore, when we convert sfg
objects of type polygon into a simple feature geometry column, we would also end up with an sfc
object of type polygon, which can be verified with st_geometry_type()
.Equally, a geometry column of multilinestrings would result in an sfc
object of type multilinestring:
# sfc POLYGON
polygon_list1 = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
polygon1 = st_polygon(polygon_list1)
polygon_list2 = list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2)))
polygon2 = st_polygon(polygon_list2)
polygon_sfc = st_sfc(polygon1, polygon2)
st_geometry_type(polygon_sfc)
#> [1] POLYGON POLYGON
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON ... TRIANGLE
# sfc MULTILINESTRING
multilinestring_list1 = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)),
rbind(c(1, 2), c(2, 4)))
multilinestring1 = st_multilinestring((multilinestring_list1))
multilinestring_list2 = list(rbind(c(2, 9), c(7, 9), c(5, 6), c(4, 7), c(2, 7)),
rbind(c(1, 7), c(3, 8)))
multilinestring2 = st_multilinestring((multilinestring_list2))
multilinestring_sfc = st_sfc(multilinestring1, multilinestring2)
st_geometry_type(multilinestring_sfc)
#> [1] MULTILINESTRING MULTILINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON ... TRIANGLE
It is also possible to create an sfc
object from sfg
objects with different geometry types:
# sfc GEOMETRY
point_multilinestring_sfc = st_sfc(point1, multilinestring1)
st_geometry_type(point_multilinestring_sfc)
#> [1] POINT MULTILINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON ... TRIANGLE
As mentioned before, sfc
objects can additionally store information on the coordinate reference systems (CRS).To specify a certain CRS, we can use the epsg (SRID)
or proj4string
attributes of an sfc
object.The default value of epsg (SRID)
and proj4string
is NA
(Not Available), as can be verified with st_crs()
:
st_crs(points_sfc)
#> Coordinate Reference System: NA
All geometries in an sfc
object must have the same CRS.We can add coordinate reference system as a crs
argument of st_sfc()
.This argument accepts an integer with the epsg
code such as 4326
, which automatically adds the ‘proj4string’ (see Section 2.4):
# EPSG definition
points_sfc_wgs = st_sfc(point1, point2, crs = 4326)
st_crs(points_sfc_wgs)
#> Coordinate Reference System:
#> EPSG: 4326
#> proj4string: "+proj=longlat +datum=WGS84 +no_defs"
It also accepts a raw proj4string (result not shown):
# PROJ4STRING definition
st_sfc(point1, point2, crs = "+proj=longlat +datum=WGS84 +no_defs")
Sometimes st_crs()
will return a proj4string
but not an epsg
code.This is because there is no general method to convert from proj4string
to epsg
(see Chapter 6).
2.2.8 The sf class
Sections 2.2.5 to 2.2.7 deal with purely geometric objects, ‘sf geometry’ and ‘sf column’ objects, respectively.These are geographic building blocks of geographic vector data represented as simple features.The final building block is non-geographic attributes, representing the name of the feature or other attributes such as measured values, groups, and other things.
To illustrate attributes, we will represent a temperature of 25°C in London on June 21st, 2017.This example contains a geometry (the coordinates), and three attributes with three different classes (place name, temperature and date).11Objects of class sf
represent such data by combining the attributes (data.frame
) with the simple feature geometry column (sfc
).They are created with st_sf()
as illustrated below, which creates the London example described above:
lnd_point = st_point(c(0.1, 51.5)) # sfg object
lnd_geom = st_sfc(lnd_point, crs = 4326) # sfc object
lnd_attrib = data.frame( # data.frame object
name = "London",
temperature = 25,
date = as.Date("2017-06-21")
)
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom) # sf object
What just happened? First, the coordinates were used to create the simple feature geometry (sfg
).Second, the geometry was converted into a simple feature geometry column (sfc
), with a CRS.Third, attributes were stored in a data.frame
, which was combined with the sfc
object with st_sf()
.This results in an sf
object, as demonstrated below (some output is ommited):
lnd_sf
#> Simple feature collection with 1 features and 3 fields
#> ...
#> name temperature date geometry
#> 1 London 25 2017-06-21 POINT (0.1 51.5)
class(lnd_sf)
#> [1] "sf" "data.frame"
The result shows that sf
objects actually have two classes, sf
and data.frame
.Simple features are simply data frames (square tables), but with spatial attributes stored in a list column, usually called geometry
, as described in Section 2.2.1.This duality is central to the concept of simple features:most of the time a sf
can be treated as and behaves like a data.frame
.Simple features are, in essence, data frames with a spatial extension.