This vignette describes the functions in sf that canhelp to plot simple features. It tries to be complete about the plotmethods sf provides, and give examples and pointers tooptions to plot simple feature objects with other packages (mapview,tmap, ggplot2).
Plot methods for sf and sfc objects
Geometry only: sfc
Geometry list-columns (objects of class sfc, obtained bythe st_geometry method) only show the geometry:
library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUEdemo(nc, ask = FALSE, echo = FALSE)plot(st_geometry(nc))
which can be further annotated with colors, symbols, etc., as theusual base plots, e.g.points are added to a polygon plot by:
plot(st_geometry(nc), col = sf.colors(12, categorical = TRUE), border = 'grey', axes = TRUE)plot(st_geometry(st_centroid(nc)), pch = 3, col = 'red', add = TRUE)## Warning: st_centroid assumes attributes are constant over geometries
and legends, titles and so on can be added afterwards.border = NA removes the polygon borders.
As can be seen, the axes plotted are sensitive to the CRS, and incase of longitude/latitude coordinates, degree symbols and orientationare added if axes = TRUE.
Geometry with attributes: sf
The default plot of an sf object is a multi-plot of allattributes, up to a reasonable maximum:
plot(nc)## Warning: plotting the first 9 out of 14 attributes; use max.plot = 14 to plot## all
with a warning when not all attributes can be reasonably plotted. Onecan increase the maximum number of maps to be plotted by
plot(nc, max.plot = 14)
The row/column layout is chosen such that the plotting area ismaximally filled. The default value for max.plot can becontrolled, e.g.by setting the global optionsf_max.plot:
options(sf_max.plot=1)plot(nc)
Color key place and size
In case a single attribute is selected, by default a color key isgiven the side of the plot where it leaves as much as possible room forthe plotted map; for nc this is below:
plot(nc["AREA"])
but this can be controlled, and set to a particular side (1=below,2=left, 3=above and 4=right):
plot(nc["AREA"], key.pos = 4)
The size of a color key can be controlled, using either relativeunits (a number between 0 and 1) or absolute units (likelcm(2) for 2 cm):
plot(nc["AREA"], key.pos = 1, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)
Keys for factor variables are a bit different, as we typically don’twant to rotate text for them:
nc$f = cut(nc$AREA, 10)plot(nc["f"], axes = TRUE, key.pos = 4, pal = sf.colors(10), key.width = lcm(5))
Class intervals
Color breaks (class intervals) can be controlled by plot argumentsbreaks and nbreaks. nbreaksspecifies the number of breaks; breaks is either a vectorwith break values:
plot(nc["AREA"], breaks = c(0,.05,.1,.15,.2,.25))
or breaks is used to indicate a breaks-finding methodthat is passed as the style argument toclassInt::classIntervals(). Its default value,pretty, results in rounded class breaks, and has as a sideeffect that nbreaks may be honoured only approximately.Other methods include "equal" to break the data range into"nbreaks" equal classes, "quantile" to usequantiles as class breaks, and "jenks", used in othersoftware.
plot(nc["AREA"], breaks = "jenks")
How does sf project geographic coordinates?
Package sf plots projected maps in their nativeprojection, meaning that easting and northing are mapped linearly to thex and y axis, keeping an aspect ratio of 1 (one unit east equals oneunit north). For geographic data, where coordinates constitute degreeslongitude and latitude, it chooses an equirectangularprojection (also called equidistant circular), where at thecenter of the plot (or of the bounding box) one unit north equals oneunit east.
Proj.4 also lets you project data to this projection, and the plotof
plot(st_geometry(nc), axes = TRUE)
should, apart from the values along axes, be otherwise identicalto
lat_ts = mean(st_bbox(nc)[c(2,4)]) # latitude of true scaleeqc = st_transform(nc, paste0("+proj=eqc +lat_ts=", lat_ts))plot(st_geometry(eqc), axes = TRUE)
Graticules
Graticules are grid lines along equal longitude (meridians) orlatitude (parallels) that, depending on the projection used, often plotas curved lines on a map, giving it reference in terms of longitude andlatitude. sf::st_graticule() tries to create a graticulegrid for arbitrary maps. As there are infinitely many projections, thereare most likely many cases where it does not succeed in doing this well,and examples of these are welcomed as sf issues.
The following plot shows a graticule geometry on itself,
library(maps)usa = st_as_sf(map('usa', plot = FALSE, fill = TRUE))laea = st_crs("+proj=laea +lat_0=30 +lon_0=-95") # Lambert equal areausa <- st_transform(usa, laea)g = st_graticule(usa)plot(st_geometry(g), axes = TRUE)
where we see that the graticule does not reach the plot boundaries(but is cut off at the bounding box of usa), and that theaxes show projected coordinates.
When we compute the graticule within the plotting function, we knowthe plotting region and can compute it up to the plot margins, and addaxes in graticule units:
plot(usa, graticule = TRUE, key.pos = NULL, axes = TRUE)
We can also pass a crs object to graticuleto obtain a graticule in a datum different from the default (WGS84).st_graticule() takes parameters, and we can pass an objectreturned by it to the graticule parameter ofplot, to get finer control:
g = st_graticule(usa, lon = seq(-130,-65,5))plot(usa, graticule = g, key.pos = NULL, axes = TRUE, xlim = st_bbox(usa)[c(1,3)], ylim = st_bbox(usa)[c(2,4)], xaxs = "i", yaxs = "i")
which still doesn’t look great – completely controlling the plottingregion of base plots is not easy.
Plotting sf objects with other packages
grid: st_as_grob
Package sf provides a number of methods forst_as_grob():
methods(st_as_grob)## [1] st_as_grob.CIRCULARSTRING* st_as_grob.COMPOUNDCURVE* ## [3] st_as_grob.CURVEPOLYGON* st_as_grob.GEOMETRYCOLLECTION* ## [5] st_as_grob.LINESTRING* st_as_grob.MULTILINESTRING* ## [7] st_as_grob.MULTIPOINT* st_as_grob.MULTIPOLYGON* ## [9] st_as_grob.MULTISURFACE* st_as_grob.POINT* ## [11] st_as_grob.POLYGON* st_as_grob.sfc* ## [13] st_as_grob.sfc_CIRCULARSTRING* st_as_grob.sfc_LINESTRING* ## [15] st_as_grob.sfc_MULTILINESTRING* st_as_grob.sfc_MULTIPOINT* ## [17] st_as_grob.sfc_MULTIPOLYGON* st_as_grob.sfc_POINT* ## [19] st_as_grob.sfc_POLYGON* ## see '?methods' for accessing help and source codewhich convert simple simple feature objects into grob(“graphics objects”) objects; grobs are the graphicprimitives of the grid plotting package. These methods canbe used by plotting packages that build on grid, such asggplot2 (which uses them in geom_sf()) andtmap. In addition, st_viewport() can be usedto set up a grid viewport from an sf object, with an aspectratio similar to that of plot.sf().
ggplot2
contains a geom specially for simple feature objects, with supportfor graticule white lines in the background usingsf::st_graticule(). Support is currently good for polygons;for lines or points, your mileage may vary.
Polygons can be colored using aes:
ggplot() + geom_sf(data = nc, aes(fill = BIR74)) + scale_y_continuous(breaks = 34:36)
and sets of maps can be plotted as facet plots after rearranging thesf object, e.g.by
## ## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':## ## filter, lag## The following objects are masked from 'package:base':## ## intersect, setdiff, setequal, unionlibrary(tidyr)nc2 <- nc %>% select(SID74, SID79, geom) %>% gather(VAR, SID, -geom)ggplot() + geom_sf(data = nc2, aes(fill = SID)) + facet_wrap(~VAR, ncol = 1) + scale_y_continuous(breaks = 34:36)
mapview
Package mapview creates interactive maps in html pages,using package leaflet as a workhorse. Extensive examplesare found here.
An example is obtained by
library(mapview)mapviewOptions(fgb = FALSE) # needed when creating web pagesmapview(nc["BIR74"], col.regions = sf.colors(10), fgb = FALSE)gives a map which is interactive: you can zoom and pan, and queryfeatures by clicking on them.
tmap
Package tmap is another package for plotting maps, withemphasis on production-ready maps.
library(tmap)qtm(nc)
tmap also has interactive leaflet maps:
tmap_mode("view")tm_shape(nc) + tm_fill("BIR74", palette = sf.colors(5))Replotting the last map in non-interactive mode is as simple as:
ttm()tmap_last()A draft version of the book Elegant and informative maps withtmap by Martijn Tennekes and Jakub Nowosad is found at https://r-tmap.github.io/
