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 non-geometrical 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 (n-ary 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 n-ary operations, sfg or sf 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_linestring(rbind(c(0,0), c(1,1), c(0,1), c(1,0))))) # self-intersects
#> [1]  TRUE FALSE

st_is_valid returns whether a geometry is valid

    st_linestring(rbind(c(1,1), c(1,2))),
    st_linestring(rbind(c(1,1), c(1,1))))) # zero-length
#> [1]  TRUE FALSE

st_is_empty returns whether a geometry is empty

#> [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_64-pc-linux-gnu-library/3.6/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.3 ymin: 33.9 xmax: -75.5 ymax: 36.6
#> geographic CRS: NAD27
#> [1] TRUE
nc2 <- st_transform(nc, 3857) # to web Mercator
#> [1] FALSE
#> [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
#> List of 2
#>  $ : int [1:3] 2 18 19
#>  $ : int [1:3] 1 3 18
#>  - attr(*, "predicate")= chr "touches"
#>  - attr(*, "")= chr [1:2] "1" "2"
#>  - attr(*, "ncol")= int 100
#>  - attr(*, "class")= chr "sgbp"

sgbp objects have the following methods:

methods(class = 'sgbp')
#> [1] as.matrix     dim           Ops           print        
#> [6] t            
#> see '?methods' for accessing help and source code

For understanding predicates, the dimensionally extended 9-intersection model (DE-9IM) is adopted, which is explained in more detail on Wikipedia (Clementini, Di Felice, and Oosterom 1993, Egenhofer and Franzosa (1991)). Briefly, DE-9IM 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 non-existent 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 0-dimensional 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, and
  • 102 that the exterior of the point (all points except this one) have a 1-dimensional overlap with the interior of the line, a 0-dimensional overlap with the boundary of the line (its end points), and a 2-dimensional 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 non-empty (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 N-ary

Higher-order 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 N-ary

No higher-order functions returning a measure are available.

5.3 Geometry generating functions

5.3.1 Unary

Unary operations work on a per-geometry 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, an sfg value is returned
  • for sfc input, a new set of geometries is returned as sfc
  • for sf objects, the same sf 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 LINESTRINGs.
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')
mp = st_multipoint(matrix(runif(20), 10))
plot(st_voronoi(mp), add = TRUE, col = NA, border = 'red')
left: convex hull (red) around a polygon (black); right: voronoi diagram (red) from a `MULTIPOINT` (black)

Figure 5.1: left: convex hull (red) around a polygon (black); right: voronoi diagram (red) from a MULTIPOINT (black)

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)
#> 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 non-empty 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 n-ary operation, explained in the next section.

5.3.3 N-ary operations returning a geometry 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))
#> [1] TRUE
st_union(x) %>% st_is_valid()
#> [1] TRUE N-ary intersection and difference

N-ary 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: 1-1, 1-1, 1-3, 2-1, 2-2, 2-3, 3-1, 3-2, 3-3:

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 list-column with the indexes of contributing feature geometries.

left: three overlapping boxes -- how do we identify the small box where all three overlap? right: unique, non-overlapping n-ary intersections

Figure 5.2: left: three overlapping boxes – how do we identify the small box where all three overlap? right: unique, non-overlapping n-ary intersections

Similarly, one can compute n-ary differences from a set \(\{s_1, s_2, s_3, ...\}\) by creating differences \(\{s_1, s_2-s_1, s_3-s_2-s_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.

difference between subsequent boxes, left: in original order; right: in reverse order

Figure 5.3: difference between subsequent boxes, left: in original order; right: in reverse order

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 small-scale (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 8-byte 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 well-known 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, 8-byte doubles get converted to 4-byte 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 / Pseudo-Mercator
#> 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 non-noded 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))
#> [1] TRUE
n-ary intersection that may lead to invalid geometries

Figure 5.4: n-ary intersection that may lead to invalid geometries

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.


Clementini, Eliseo, Paolino Di Felice, and Peter van Oosterom. 1993. “A Small Set of Formal Topological Relationships Suitable for End-User 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. “Point-Set 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/RJ-2018-009.

Pebesma, Edzer, Thomas Mailund, and James Hiebert. 2016a. “Measurement Units in R.” The R Journal 8 (2): 486–94. doi:10.32614/RJ-2016-061.