5. Plotting Simple Features

Edzer Pebesma

This vignette describes the functions insf that canhelp to plot simple features. It tries to be complete about the plotmethodssf 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 classsfc, obtained bythest_geometry method) only show the geometry:

library(sf)demo(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 ifaxes = TRUE.

Geometry with attributes:`sf`

The default plot of ansf object is a multi-plot of allattributes, 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. 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 formax.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; fornc 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 andnbreaks.nbreaksspecifies the number of breaks;breaks is either a vectorwith break values:

plot(nc["AREA"], breaks = c(0,.05,.1,.15,.2,.25))

orbreaks is used to indicate a breaks-finding methodthat is passed as thestyle argument toclassInt::classIntervals(). Its default value,pretty, results in rounded class breaks, and has as a sideeffect thatnbreaks 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?

Packagesf 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 anequirectangularprojection (also calledequidistant 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 assf 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 ofusa), 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 acrs object tograticuleto obtain a graticule in a datum different from the default (WGS84).st_graticule() takes parameters, and we can pass an objectreturned by it to thegraticule 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`

Packagesf 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 code

which convert simple simple feature objects intogrob(“graphics objects”) objects;grobs are the graphicprimitives of thegrid plotting package. These methods canbe used by plotting packages that build ongrid, such asggplot2 (which uses them ingeom_sf()) andtmap. In addition,st_viewport() can be usedto set up a grid viewport from ansf object, with an aspectratio similar to that ofplot.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.

library(ggplot2)ggplot() + geom_sf(data = usa)

Polygons can be colored usingaes:

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

library(dplyr)library(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

Packagemapview creates interactive maps in html pages,using packageleaflet as a workhorse. Extensive examplesare foundhere.

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

Packagetmap 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 bookElegant and informative maps withtmap by Martijn Tennekes and Jakub Nowosad is found athttps://r-tmap.github.io/

Movatterモバイル変換