Chapter 5 Manipulating Geometries
Simple feature geometries can be queried for properties, combined into new geometries, and combinations of geometries can be queried for properties. This chapter will give an overview of the operations offered by sf
, entirely focusing on geometrical properties. The next chapter, 6, focuses on the analysis of nongeometrical feature properties, in relationship to their geometries. Some of the material in this chapter also appeared as (Pebesma 2018).
Several of the concepts of geometric manipulations were introduced in chapter 3. This chapter gives a complete listing of all geometries permitted on geometries, illustrating some of them.
We can categorise operations in terms of what they take as input, and what they give as output. In terms of output we have operations that give one or more
 predicates: a logical asserting a certain property is
TRUE
,  measures: a value (e.g. a numeric value with measurement unit), or
 geometries
and in terms of what they operate on, we distinguish operations that work on
 a single geometry (unary operations)
 pairs of geometries (binary operations)
 sets of geometries (nary operations)
Before we will go through all combinations, we make two observations:
 most functions are implemented as methods, and operate equally on single geometry objects (
sfg
), geometry set objects (sfc
) or simple feature (sf
) objects.  also for binary and nary operations,
sfg
orsf
objects are accepted as input, and taken as a set of geometries.
5.1 Predicates
Predicates return a logical, TRUE
or FALSE
value, or a set of those.
5.1.1 Unary predicates
st_is_simple
returns whether a geometry is simple:
st_is_simple(st_sfc(
st_point(c(0,1)),
st_linestring(rbind(c(0,0), c(1,1), c(0,1), c(1,0))))) # selfintersects
#> [1] TRUE FALSE
st_is_valid
returns whether a geometry is valid
st_is_valid(st_sfc(
st_linestring(rbind(c(1,1), c(1,2))),
st_linestring(rbind(c(1,1), c(1,1))))) # zerolength
#> [1] TRUE FALSE
st_is_empty
returns whether a geometry is empty
st_is_empty(st_point())
#> [1] TRUE
st_is_longlat
returns whether the coordinate reference system is geographic (chapter 2, 7):
demo(nc, ask = FALSE, echo = FALSE)
#> Reading layer `nc.gpkg' from data source `/home/edzer/R/x86_64pclinuxgnulibrary/3.6/sf/gpkg/nc.gpkg' using driver `GPKG'
#> Simple feature collection with 100 features and 14 fields
#> Attributegeometry relationship: 0 constant, 8 aggregate, 6 identity
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: 84.3 ymin: 33.9 xmax: 75.5 ymax: 36.6
#> geographic CRS: NAD27
st_is_longlat(nc)
#> [1] TRUE
nc2 < st_transform(nc, 3857) # to web Mercator
st_is_longlat(nc2)
#> [1] FALSE
st_is_longlat(st_point(0:1))
#> [1] NA
st_is
is an easy way to check for the simple feature geometry type:
st_is(st_point(0:1), "POINT")
#> [1] TRUE
all(st_is(nc, "POLYGON"))
#> [1] FALSE
all(st_is(nc, "MULTIPOLYGON"))
#> [1] TRUE
Equality and inequality of geometries can be checked by ==
or !=
; it uses geometric equality, and is insensitive to the order of traversal of nodes:
st_sfc(st_point(0:1), st_point(1:2)) == st_sfc(st_point(0:1))
#> [1] TRUE FALSE
st_linestring(rbind(c(0,0), c(1,1))) == st_linestring(rbind(c(1,1), c(0,0)))
#> [1] TRUE
Under the hood, it uses st_equals
, discussed by the binary predicates.
5.1.2 Binary predicates
Binary predicates result in a TRUE
or FALSE
value for every pair of inputs. For two sets of inputs with \(n\) and \(m\) geometries respectively, this results in an \(n \times m\) logical matrix. Because \(n\) and/or \(m\) may be very large and the predicate matrix typically contains mostly FALSE
values, a sparse representation of it, a sparse geometry binary predicate (sgbp
) object, is returned by all binary predicate functions (which are shown below). sgbp
objects are lists of indices of the TRUE
values in each row:
(r < st_touches(nc2[1:2,], nc2))
#> Sparse geometry binary predicate list of length 2, where the predicate was `touches'
#> 1: 2, 18, 19
#> 2: 1, 3, 18
str(r)
#> List of 2
#> $ : int [1:3] 2 18 19
#> $ : int [1:3] 1 3 18
#>  attr(*, "predicate")= chr "touches"
#>  attr(*, "region.id")= chr [1:2] "1" "2"
#>  attr(*, "ncol")= int 100
#>  attr(*, "class")= chr "sgbp"
sgbp
objects have the following methods:
methods(class = 'sgbp')
#> [1] as.data.frame as.matrix dim Ops print
#> [6] t
#> see '?methods' for accessing help and source code
For understanding predicates, the dimensionally extended 9intersection model (DE9IM) is adopted, which is explained in more detail on Wikipedia (Clementini, Di Felice, and Oosterom 1993, Egenhofer and Franzosa (1991)). Briefly, DE9IM considers that every geometry has an interior, a boundary and an exterior. For polygons this is trivial, for points the boundary is an empty set, for linestrings the boundary is formed by the end points and the interior by all non end points. Also, any geometry has a dimension of 0 (points), 1 (lines) or 2 (polygons) or nonexistent in the case of an empty geometry.
A relationship between two geometries A and B is expressed by the dimension of the overlap (intersections) of 9 intersections, formed by the 9 pairs from the interior, boundary and exterior of A, and the interior, boundary and exterior of B. We can query this relation by using st_relate
B = st_linestring(rbind(c(0,0), c(1,0)))
A = st_point(c(0.5, 0)) # halfway the line
st_relate(A, B)
#> [,1]
#> [1,] "0FFFFF102"
In the relationship found, 0FFFFF102
, F
indicates empty geometries, and we see from
0FF
that the (interior of the) point has 0dimensional overlap with the interior of line (i.e., the overlap is a point), and no overlap with the boundary (end points) or the exterior of the line,FFF
that the (empty) border of the point has nothing in common with the line, and102
that the exterior of the point (all points except this one) have a 1dimensional overlap with the interior of the line, a 0dimensional overlap with the boundary of the line (its end points), and a 2dimensional overlap with the exterior of the line.
We can query whether a particular relationship holds by giving st_relate
a pattern. To check for instance whether point A overlaps with an end point of linestring B, we can use
st_relate(A, B, pattern = "F0FFFFFFF") %>% as.matrix()
#> [,1]
#> [1,] FALSE
In these patterns, *
can be used for anything, and T
for nonempty (0, 1 or 2). The standard relationships below are all expressed as particular query patterns, the Wikipedia page gives details on the patterns used.
The binary predicates provided by package sf
are
predicate  value  inverse of 

st_contains 
None of the points of A are outside B  st_within 
st_contains_properly 
A contains B and B has no points in common with the boundary of A  
st_covers 
No points of B lie in the exterior of A  st_covered_by 
st_covered_by 
inverse of st_covers 

st_crosses 
A and B have some but not all interior points in common  
st_disjoint 
A and B have no points in common  st_intersects 
st_equals 
A and B are geometrically equal; node order number of nodes may differ; identical to A contains B AND A within B  
st_equals_exact 
A and B are geometrically equal, and have identical node order  
st_intersects 
A and B are not disjoint  st_disjoint 
st_is_within_distance 
A is closer to B than a given distance  
st_within 
None of the points of B are outside A  st_contains 
st_touches 
A and B have at least one boundary point in common, but no interior points  
st_overlaps 
A and B have some points in common; the dimension of these is identical to that of A and B  
st_relate 
given a pattern, returns whether A and B adhere to this pattern 
5.1.3 Nary
Higherorder predicates are not supported by special functions.
5.2 Measures
5.2.1 Unary
Unary measures return a single value that describes a property of the geometry:
function 
returns 

st_dimension 
0 for points, 1 for linear, 2 for polygons, NA for empty geometries 
st_area 
the area for geometries 
st_length 
the lengths of linear geometries 
lwgeom::st_geohash 
the geohash for geometries 
st_geometry_type 
the types of a set of geometries 
5.2.2 Binary
st_distance
returns the distances between pairs of geometries, either as a vector with distances between the two first, the two second, … pairs, or as a matrix with all pairwise distances. The result is numeric, or is of class units
(E. Pebesma, Mailund, and Hiebert 2016a) when distance units can be derived from the coordinate reference system (chapter 7):
st_distance(nc[1:3,], nc[2:4,], by_element = TRUE) %>% setNames(NULL)
#> Units: [m]
#> [1] 0 0 367505
st_distance(nc[1:3,], nc[2:4,])
#> Units: [m]
#> [,1] [,2] [,3]
#> [1,] 0 25650 440513
#> [2,] 0 0 409370
#> [3,] 0 0 367505
st_relate
returns the relation pattern, as explained in section 5.1.2, or an sgbp
object when given a pattern template to match to.
5.2.3 Nary
No higherorder functions returning a measure are available.
5.3 Geometry generating functions
5.3.1 Unary
Unary operations work on a pergeometry basis, and for each geometry return a new geometry. None of these functions operate on more than one feature geometry. Most functions are implemented as (S3) generic, with methods for sfg
, sfc
and sf
; their output is of the same class as their input:
 for
sfg
input, ansfg
value is returned  for
sfc
input, a new set of geometries is returned assfc
 for
sf
objects, the samesf
object is returned which has geometries replaced with the new ones.
function  returns a geometry… 

st_centroid 
of type POINT with the geometry’s centroid 
st_buffer 
that is this larger (or smaller) than the input geometry, depending on the buffer size 
st_jitter 
that was moved in space a certain amount, using a bivariate uniform distribution 
st_wrap_dateline 
cut into pieces that do no longer cover the dateline 
st_boundary 
with the boundary of the input geometry 
st_convex_hull 
that forms the convex hull of the input geometry (figure 5.1) 
st_line_merge 
after merging connecting LINESTRING elements of a MULTILINESTRING into longer LINESTRING s. 
st_make_valid 
that is valid 
st_node 
with added nodes to linear geometries at intersections without a node; only works on individual linear geometries 
st_point_on_surface 
with a (arbitrary) point on a surface 
st_polygonize 
of type polygon, created from lines that form a closed ring 
st_segmentize 
a (linear) geometry with nodes at a given density or minimal distance 
st_simplify 
simplified by removing vertices/nodes (lines or polygons) 
lwgeom::st_split 
that has been split with a splitting linestring 
st_transform 
transformed to a new coordinate reference system (chapter 7) 
st_triangulate 
with triangulated polygon(s) 
st_voronoi 
with the voronoi tesselation of an input geometry (figure 5.1) 
st_zm 
with removed or added Z and/or M coordinates 
st_collection_extract 
with subgeometries from a GEOMETRYCOLLECTION of a particular type 
st_cast 
that is converted to another type 
par(mar = rep(0,4), mfrow = c(1, 2))
plot(st_geometry(nc)[1], col = NA, border = 'black')
plot(st_convex_hull(st_geometry(nc)[1]), add = TRUE, col = NA, border = 'red')
box()
set.seed(131)
mp = st_multipoint(matrix(runif(20), 10))
plot(mp)
plot(st_voronoi(mp), add = TRUE, col = NA, border = 'red')
box()
A number of operation can be applied directly to geometries
(A = st_point(c(1,2)))
#> POINT (1 2)
(B = st_linestring(rbind(c(2,2), c(3,4))))
#> LINESTRING (2 2, 3 4)
A
#> POINT (1 2)
B + A
#> LINESTRING (3 4, 4 6)
st_sfc(B + A) * matrix(c(1,0,0,2), 2, 2)
#> Geometry set for 1 feature
#> geometry type: LINESTRING
#> dimension: XY
#> bbox: xmin: 3 ymin: 8 xmax: 4 ymax: 12
#> CRS: NA
#> LINESTRING (3 8, 4 12)
st_sfc(A, B) * c(3, 5) # scale first by 3, second by 5:
#> Geometry set for 2 features
#> geometry type: GEOMETRY
#> dimension: XY
#> bbox: xmin: 3 ymin: 6 xmax: 15 ymax: 20
#> CRS: NA
#> POINT (3 6)
#> LINESTRING (10 10, 15 20)
5.3.2 Binary operations returning geometries
Binary functions that return a geometry include
function  returns  infix operator 

st_intersection 
the overlapping geometries for pair of geometries  & 
st_union 
the combination of the geometries; also removes duplicate points, nodes or line pieces   
st_difference 
the geometries of the first after removing the overlap with the second geometry  / 
st_sym_differenc 
the combinations of the geometries after removing where they overlap  %/% 
When operating on two sfg
, single geometries, it is clear what all these functions do: return a single geometry for this pair. When given two sets of geometries (sfc
or sf
objects), a new set of geometries is returned; for st_intersection
containing only the nonempty geometries, for all other operations the geometries from all pairwise evaluation. In case the arguments are of class sf
, the attributes of the objects are copied over to all intersections to which each of the features contributed:
a = st_sf(a = 1, geom = st_sfc(st_linestring(rbind(c(0,0), c(1,0)))))
b = st_sf(b = 1:3, geom = st_sfc(st_point(c(0,0)), st_point(c(1,0)), st_point(c(2,0))))
st_intersection(a, b)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
#> Simple feature collection with 2 features and 2 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 1 ymax: 0
#> CRS: NA
#> a b geom
#> 1 1 1 POINT (0 0)
#> 1.1 1 2 POINT (1 0)
When st_intersection
or st_difference
are called with a single set of geometries (an sfc
object), they perform an nary operation, explained in the next section.
5.3.3 Nary operations returning a geometry
5.3.3.1 Union, c, and combine
Calling st_union
with only a single argument leads either to computing the union of all geometries, or applying union to each of the individual geometries, depending on the setting of by_feature
:
st_union(b, by_feature = FALSE) # default
#> Geometry set for 1 feature
#> geometry type: MULTIPOINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 2 ymax: 0
#> CRS: NA
#> MULTIPOINT ((0 0), (1 0), (2 0))
st_union(b, by_feature = TRUE)
#> Simple feature collection with 3 features and 1 field
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 2 ymax: 0
#> CRS: NA
#> b geom
#> 1 1 POINT (0 0)
#> 2 2 POINT (1 0)
#> 3 3 POINT (2 0)
The c
method combines sets of geometries
bb = st_geometry(b)
c(bb, bb)
#> Geometry set for 6 features
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 2 ymax: 0
#> CRS: NA
#> First 5 geometries:
#> POINT (0 0)
#> POINT (1 0)
#> POINT (2 0)
#> POINT (0 0)
#> POINT (1 0)
or single geometries into single a new single geometry
c(st_point(0:1), st_point(1:2))
#> MULTIPOINT ((0 1), (1 2))
and st_combine
uses this to collapse features for different geometries into one:
st_combine(c(bb, bb))
#> Geometry set for 1 feature
#> geometry type: MULTIPOINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 2 ymax: 0
#> CRS: NA
#> MULTIPOINT ((0 0), (1 0), (2 0), (0 0), (1 0), ...
When using this on lines or polygons, it is easy to obtain invalid geometries, and one needs to use st_union
on the result.
(x = st_combine(st_sfc(st_linestring(rbind(c(0,0), c(1,1))), st_linestring(rbind(c(1,0),c(0,1))))))
#> Geometry set for 1 feature
#> geometry type: MULTILINESTRING
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 1 ymax: 1
#> CRS: NA
#> MULTILINESTRING ((0 0, 1 1), (1 0, 0 1))
st_is_valid(x)
#> [1] TRUE
st_union(x) %>% st_is_valid()
#> [1] TRUE
5.3.3.2 Nary intersection and difference
Nary st_intersection
and st_difference
take a single argument, but operate (sequentially) on all pairs, triples, quadruples etc. Consider the plot in figure 5.2: how do we identify the box where all three overlap? Using binary intersections, as of gives us intersections of all pairs, double since x
is passed twice: 11, 11, 13, 21, 22, 23, 31, 32, 33:
sq = function(pt, sz = 1) st_polygon(list(rbind(c(pt  sz),
c(pt[1] + sz, pt[2]  sz), c(pt + sz), c(pt[1]  sz, pt[2] + sz), c(pt  sz))))
x = st_sf(box = 1:3, st_sfc(sq(c(0,0)), sq(c(1.7, 0.5)), sq(c(0.5, 1))))
(ixx = st_intersection(x, x)) %>% nrow
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
#> [1] 9
lengths(st_overlaps(ixx, ixx))
#> [1] 4 5 5 5 4 5 5 5 4
When we use however
(i = st_intersection(x))
#> Simple feature collection with 7 features and 3 fields
#> geometry type: POLYGON
#> dimension: XY
#> bbox: xmin: 1 ymin: 1.5 xmax: 2.7 ymax: 2
#> CRS: NA
#> box n.overlaps origins st_sfc.sq.c.0..0....sq.c.1.7...0.5....sq.c.0.5..1...
#> 1 1 1 1 POLYGON ((0.7 1, 1 1, 1...
#> 1.1 1 2 1, 2 POLYGON ((1 0, 1 1, 0.7 1...
#> 2 2 1 2 POLYGON ((1.5 0.5, 2.7 0.5,...
#> 2.1 2 2 2, 3 POLYGON ((1 0.5, 1.5 0.5, 1...
#> 1.2 1 3 1, 2, 3 POLYGON ((1 0.5, 1 0, 0.7 0...
#> 1.3 1 2 1, 3 POLYGON ((0.5 1, 1 1, 1 0....
#> 3 3 1 3 POLYGON ((0.5 1, 0.5 2, 1...
we end up with a set of all seven distinct intersections, without overlaps.
lengths(st_overlaps(i, i))
#> [1] 0 0 0 0 0 0 0
When given an sf
object an sf
is returned, with two additional fields, one with the number of overlapping features, and a listcolumn with the indexes of contributing feature geometries.
Similarly, one can compute nary differences from a set \(\{s_1, s_2, s_3, ...\}\) by creating differences \(\{s_1, s_2s_1, s_3s_2s_1, ...\}\). This is done by
(xd = st_difference(x))
#> Simple feature collection with 3 features and 1 field
#> geometry type: POLYGON
#> dimension: XY
#> bbox: xmin: 1 ymin: 1.5 xmax: 2.7 ymax: 2
#> CRS: NA
#> box st_sfc.sq.c.0..0....sq.c.1.7...0.5....sq.c.0.5..1...
#> 1 1 POLYGON ((1 1, 1 1, 1 1,...
#> 2 2 POLYGON ((1 0.5, 2.7 0.5, 2...
#> 3 3 POLYGON ((0.5 1, 0.5 2, 1...
The result is shown in figure 5.3, for x
and for x[3:1]
, to make clear that the result here depends on order of the geometries.
Resulting geometries do not overlap:
lengths(st_overlaps(xd, xd))
#> [1] 0 0 0
5.3.4 Other geometry manipulators
st_make_grid
creates a grid of square or hexagonal polygons, based on an input bounding box and a grid cell size.
st_graticule
creates a set of graticules, lines of constant latitude or longitude, which can serve as a reference on smallscale (large area) maps.
5.4 Precision
Geometrical operations, such as finding out whether a certain point is on a line, may fail when coordinates are represented by highly precise floating point numbers, such as 8byte doubles in R. A remedy might be to limit the precision of the coordinates before the operation. For this, a precision model is adopted by sf
. It uses a precision value to round coordinates (X, Y, Z and M) right before they are encoded as wellknown binary, and passed on to the libraries where this may have an effect (GEOS, GDAL, liblwgeom). We demonstrate this by an R  WKB  R roundtrip.
Rounding can be done in two different ways. First, with a negative precision value, 8byte doubles get converted to 4byte floats and back again:
(p = st_sfc(st_point(c(1e6/3, 1e4/3)), crs = 3857))
#> Geometry set for 1 feature
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 333000 ymin: 3330 xmax: 333000 ymax: 3330
#> projected CRS: WGS 84 / PseudoMercator
#> POINT (333333 3333)
p %>% st_set_precision(1) %>% st_as_binary() %>% st_as_sfc() %>% `[[`(1) %>% print(digits = 16)
#> POINT (333333.34375 3333.333251953125)
Second, with a positive precision \(p\), each coordinate value \(c\) is replaced by \[c' = \mbox{round}(p \cdot c) / p\] This implies that for instance with a precision of 1000, the number of decimal places to round to is 1/1000, or to mm if the unit of coordinates is metre:
p %>% st_set_precision(1000) %>% st_as_binary() %>% st_as_sfc() %>% `[[`(1)
#> POINT (333333 3333)
With a precision of e.g. 0.001 or 0.05, rounding to the nearest 1/precision
, i.e. if the unit is m to the nearest 1000 m or 20 m, is obtained:
p %>% st_set_precision(0.001) %>% st_as_binary() %>% st_as_sfc() %>% `[[`(1) # to nearest 1000
#> POINT (333000 3000)
p %>% st_set_precision(0.05) %>% st_as_binary() %>% st_as_sfc() %>% `[[`(1) # to nearest 20
#> POINT (333340 3340)
As a convenience, precisions can also be specified as a units
object, with the unit to round to, e.g. to the nearest 5 cm:
p %>% st_set_precision(units::set_units(5, cm)) %>%
st_as_binary() %>%
st_as_sfc() %>% `[[`(1) %>% print(digits = 10)
#> POINT (333333.35 3333.35)
but this requires that the object, p
, has a coordinate reference system with known units.
In essence, these rounding methods bring the coordinates to points on a regular grid, which is beneficial for geometric computations. Of course, it also affects all computations like areas and distances. Which precision values are best for which application is often a matter of common sense combined with trial and error. A reproducible example illustrating the need for setting precision is found here.
5.5 Generating invalid geometries
It is rather easy to have st_intersection
generate invalid geometries, resulting in an error. Consider the graph constructed and shown in figure 5.4. In this case, not setting the precision (i.e., precision has value 0) would have led to the cryptic error message
Error in CPL_nary_intersection(x) :
Evaluation error: TopologyException: found nonnoded intersection between
LINESTRING (0.329035 0.0846201, 0.333671 0.0835073) and
LINESTRING (0.330465 0.0842769, 0.328225 0.0848146)
at 0.32965918719530368 0.084470389572422672.
Calls: st_intersection ... st_intersection > st_intersection.sfc > CPL_nary_intersection
However, with zero precision and a buf_size
of 0.7 we will not get this error.
n = 12 # n points, equally spread along unit circle:
pts = (1:n)/n * 2 * pi
xy = st_as_sf(data.frame(x = cos(pts), y = sin(pts)), coords = c("x", "y"))
buf_size = 0.8
precision = 1000
b = st_buffer(xy, buf_size)
i = st_intersection(st_set_precision(b, precision))
par(mar = rep(0, 4))
plot(i[1], col = sf.colors(nrow(i), categorical = TRUE))
all(st_is_valid(i))
#> [1] TRUE
5.6 Warnings for longitude/latitude geometries
When working on geodetic coordinates (degrees longitude/latitude), package sf
gives warnings when it makes the assumption that coordinates are Cartesian, e.g. in
i = st_intersects(nc[1,], nc[2,])
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
In many cases, making this assumption is not a problem. It might be a problem when we have polygons that cover very large areas, cover North or South pole, or have lines crossing or polygons covering the dateline.
References
Clementini, Eliseo, Paolino Di Felice, and Peter van Oosterom. 1993. “A Small Set of Formal Topological Relationships Suitable for EndUser Interaction.” In Advances in Spatial Databases, edited by David Abel and Beng Chin Ooi, 277–95. Berlin, Heidelberg: Springer Berlin Heidelberg.
Egenhofer, Max J., and Robert D. Franzosa. 1991. “PointSet Topological Spatial Relations.” International Journal of Geographical Information Systems 5 (2). Taylor & Francis: 161–74. doi:10.1080/02693799108927841.
Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. doi:10.32614/RJ2018009.
Pebesma, Edzer, Thomas Mailund, and James Hiebert. 2016a. “Measurement Units in R.” The R Journal 8 (2): 486–94. doi:10.32614/RJ2016061.