This vignette describes the functions in
sf that can help to plot simple features. It tries to be complete about the plot methods
sf provides, and give examples and pointers to options to plot simple feature objects with other packages (mapview, tmap, ggplot2).
Geometry list-columns (objects of class
sfc, obtained by the
st_geometry method) only show the geometry:
library(sf) ## Loading required package: methods ## Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3 demo(nc, ask = FALSE, echo = FALSE) ## Reading layer `nc.gpkg' from data source `/home/edzer/R/x86_64-pc-linux-gnu-library/3.4/sf/gpkg/nc.gpkg' using driver `GPKG' ## Simple feature collection with 100 features and 14 fields ## Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965 ## epsg (SRID): 4267 ## proj4string: +proj=longlat +datum=NAD27 +no_defs plot(st_geometry(nc))
which can be further annotated with colors, symbols, etc., as the usual 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 in st_centroid.sfc(st_geometry(x), of_largest_polygon = ## of_largest_polygon): st_centroid does not give correct centroids for ## longitude/latitude data
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 in case of longitude/latitude coordinates, degree symbols and orientation are added if
axes = TRUE.
The default plot of an
sf object is a multi-plot of all attributes, up to a reasonable maximum:
plot(nc) ## Warning: plotting the first 10 out of 14 attributes; use max.plot = 14 to ## plot all
with a warning when not all attributes can be reasonably plotted. One can 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 is maximally filled. The default value for
max.plot can be controlled, e.g. by setting the global option
In case a single attribute is selected, by default a color key is given on the right side,
which can be moved to other sides
plot(nc["AREA"], key.pos = 1)
and the size of which can be controlled:
plot(nc["AREA"], key.pos = 1, axes = TRUE, key.size = lcm(1.3))
Keys for factor variables are a bit different, as we typically don’t want to rotate text for them:
nc$f = cut(nc$AREA, 10) plot(nc["f"], axes = TRUE, pal = sf.colors(10), key.size = lcm(4.5))
color breaks (class intervals) can be controlled by plot arguments
nbreaks specifies the number of breaks;
breaks is either a vector with break values:
plot(nc["AREA"], breaks = c(0,.05,.1,.15,.2,.25))
breaks is used to indicate a breaks finding method that is passed as the
style argument to
classInt::classIntervals. Its default value,
pretty, results in rounded class breaks, and has as a side effect that
nbreaks may be honoured only approximately. Other methods include
"equal" to define the data range into
"nbreaks" equal classes,
"quantile" to use quantiles as class breaks, and
"jenks", used in other software.
plot(nc["AREA"], breaks = "jenks")
sf plots projected maps in their native projection, meaning that easting and northing are mapped linearly to the x and y axis, keeping an aspect ratio of 1 (one unit east equals one unit north). For geographic data, where coordinates constitute degrees longitude and latitude, it chooses an equirectangular projection (also called equidistant circular), where at the center of the plot (or of the bounding box) one unit north equals one unit east.
Proj.4 also lets you project data to this projection, and the plot of
plot(st_geometry(nc), axes = TRUE)
should, apart from the values along axes, be otherwise identical to
lat_ts = mean(st_bbox(nc)[c(2,4)]) # latitude of true scale eqc = st_transform(nc, paste0("+proj=eqc +lat_ts=", lat_ts)) plot(st_geometry(eqc), axes = TRUE)
Graticules are grid lines along equal longitude (meridians) or latitude (parallels) that, depending on the projection used, often plot as curved lines on a map, giving it reference in terms of longitude and latitude. The
st_graticule tries to create a graticule grid for arbitrary maps. As there are infinitely many projections, there are 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 area usa <- 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 the axes show projected coordinates.
When we compute the graticule within the plotting function, we know the plotting region and can compute it up to the plot margins, and add axes in graticule units:
plot(usa, graticule = TRUE, key.pos = NULL, axes = TRUE)
We can also pass a
crs object to
graticule to obtain a graticule in a datum different from the default (WGS84).
st_graticule takes parameters, and we can pass an object returned by it to the
graticule parameter of
plot, 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 plotting region of base plots is not easy.
sf provides a number of methods for
methods(st_as_grob) ##  st_as_grob.CIRCULARSTRING* st_as_grob.COMPOUNDCURVE* ##  st_as_grob.CURVEPOLYGON* st_as_grob.GEOMETRYCOLLECTION* ##  st_as_grob.LINESTRING* st_as_grob.MULTILINESTRING* ##  st_as_grob.MULTIPOINT* st_as_grob.MULTIPOLYGON* ##  st_as_grob.MULTISURFACE* st_as_grob.POINT* ##  st_as_grob.POLYGON* ## see '?methods' for accessing help and source code
which convert simple simple feature objects into
grob (“graphics objects”) objects;
grobs are the graphic primitives of the
grid plotting package. These methods can be used by plotting packages that build on
grid, such as
ggplot2 (which uses them in
tmap. In addition,
st_viewport can be used to set up a grid viewport from an
sf object, with an aspect ratio similar to that of
ggplot2 version to be released after 2.2.1 (on CRAN), which currently has to be installed by
contains a geom specially for simple feature objects, with support for graticule white lines in the background using
sf::st_graticule. Support is currently good for polygons; for lines or points your milage may vary.
polygons can be colored using
if (utils::packageVersion("ggplot2") > "2.2.1") 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 the
sf object, e.g. by
library(dplyr) ## ## Attaching package: 'dplyr' ## The following objects are masked from 'package:stats': ## ## filter, lag ## The following objects are masked from 'package:base': ## ## intersect, setdiff, setequal, union library(tidyr) nc2 <- nc %>% select(SID74, SID79, geom) %>% gather(VAR, SID, -geom) if (utils::packageVersion("ggplot2") > "2.2.1") ggplot() + geom_sf(data = nc2, aes(fill = SID)) + facet_wrap(~VAR, ncol = 1) + scale_y_continuous(breaks = 34:36)
mapview creates interactive maps in html pages, using package
leaflet as a work horse. Extensive examples are found here.
gives a map which is interactive: you can zoom and pan, and query features by clicking on them.
tmap is another package for plotting maps, with emphasis on production-ready maps.
tmap also has interactive leaflet maps:
replotting the last map in non-interactive mode is as simple as:
tmap currently converts
sf objects first to
Spatial* objects using
as(obj, "Spatial"). For sf objects that do not have an equivalent representation in
sp, this fails. Examples of this are objects with
GEOMETRYCOLLECTION geometries, or mixed geometries such as obtained by
st_sf(a = 1:2, geom = st_sfc(st_point(0:1), st_linestring(rbind(c(1,1), c(2,2))))) ## Simple feature collection with 2 features and 1 field ## geometry type: GEOMETRY ## dimension: XY ## bbox: xmin: 0 ymin: 1 xmax: 2 ymax: 2 ## epsg (SRID): NA ## proj4string: NA ## a geom ## 1 1 POINT (0 1) ## 2 2 LINESTRING (1 1, 2 2)