3. Manipulating Simple Feature Geometries

This vignetted describes how simple feature geometries can be manipulated, where manipulations include

• type transformations (e.g., POLYGON to MULTIPOLYGON)
• affine transformation (shift, scale, rotate)
• transformation into a different coordinate reference system
• geometrical operations, e.g. finding the centroid of a polygon, detecting whether pairs of feature geometries intersect, or find the union (overlap) of two polygons.

Type transformations

This sections discusses how simple feature geometries of one type can be converted to another. For converting lines to polygons, see also st_polygonize below.

For single geometries

For single geometries, st_cast will

1. convert from XX to MULTIXX, e.g. LINESTRING to MULTILINESTRING
2. convert from MULTIXX to XX if MULTIXX has length one (else, it will still convert but warn about loss of information)
3. convert from MULTIXX to XX if MULTIXX does not have length one, but it will warn about the loss of information
4. convert GEOMETRYCOLLECTION of length one to its component if

Examples of the first three types are

library(sf)
suppressPackageStartupMessages(library(dplyr))
st_point(c(1,1)) %>% st_cast("MULTIPOINT")
## MULTIPOINT (1 1)
st_multipoint(rbind(c(1,1))) %>% st_cast("POINT")
## Warning in st_cast.MULTIPOINT(., "POINT"): point from first coordinate only
## POINT (1 1)
st_multipoint(rbind(c(1,1),c(2,2))) %>% st_cast("POINT")
## Warning in st_cast.MULTIPOINT(., "POINT"): point from first coordinate only
## POINT (1 1)

Examples of the fourth type are

st_geometrycollection(list(st_point(c(1,1)))) %>% st_cast("POINT")
## POINT (1 1)

For collections of geometry (sfc) and simple feature collections (sf)

It should be noted here that when reading geometries using st_read, the type argument can be used to control the class of the returned geometry:

shp = system.file("shape/nc.shp", package="sf")
class(st_geometry(st_read(shp, quiet = TRUE)))
## [1] "sfc_MULTIPOLYGON" "sfc"
class(st_geometry(st_read(shp, quiet = TRUE, type = 3)))
## [1] "sfc_POLYGON" "sfc"
class(st_geometry(st_read(shp, quiet = TRUE, type = 1)))
## [1] "sfc_GEOMETRY" "sfc"

This option is handled by the GDAL library; in case of failure to convert to the target type, the original types are returned, which in this case is a mix of POLYGON and MULTIPOLYGON geometries, leading to a GEOMETRY as superclass. When we try to read multipolygons as polygons, all secondary rings of multipolygons get lost.

When functions return objects with mixed geometry type (GEOMETRY), downstream functions such as st_write may have difficulty handling them. For some of these cases, st_cast may help modifying their type. For sets of geometry objects (sfc) and simple feature sets (sf),st_cast` can be used by specifying the target type, or without specifying it.

ls <- st_linestring(rbind(c(0,0),c(1,1),c(2,1)))
mls <- st_multilinestring(list(rbind(c(2,2),c(1,3)), rbind(c(0,0),c(1,1),c(2,1))))
(sfc <- st_sfc(ls,mls))
## Geometry set for 2 features 
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 2 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
## LINESTRING (0 0, 1 1, 2 1)
## MULTILINESTRING ((2 2, 1 3), (0 0, 1 1, 2 1))
st_cast(sfc, "MULTILINESTRING")
## Geometry set for 2 features 
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 2 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
## MULTILINESTRING ((0 0, 1 1, 2 1))
## MULTILINESTRING ((2 2, 1 3), (0 0, 1 1, 2 1))
sf <- st_sf(a = 5:4, geom = sfc)
st_cast(sf, "MULTILINESTRING")
## Simple feature collection with 2 features and 1 field
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 2 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##   a                       geometry
## 1 5 MULTILINESTRING ((0 0, 1 1,...
## 2 4 MULTILINESTRING ((2 2, 1 3)...

When no target type is given, st_cast tries to be smart for two cases:

1. if the class of the object is GEOMETRY, and all elements are of identical type, and
2. if all elements are length-one GEOMETRYCOLLECTION objects, in which case GEOMETRYCOLLECTION objects are replaced by their content (which may be a GEOMETRY mix again)

Examples are:

ls <- st_linestring(rbind(c(0,0),c(1,1),c(2,1)))
mls1 <- st_multilinestring(list(rbind(c(2,2),c(1,3)), rbind(c(0,0),c(1,1),c(2,1))))
mls2 <- st_multilinestring(list(rbind(c(4,4),c(4,3)), rbind(c(2,2),c(2,1),c(3,1))))
(sfc <- st_sfc(ls,mls1,mls2))
## Geometry set for 3 features 
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 4
## epsg (SRID):    NA
## proj4string:    NA
## LINESTRING (0 0, 1 1, 2 1)
## MULTILINESTRING ((2 2, 1 3), (0 0, 1 1, 2 1))
## MULTILINESTRING ((4 4, 4 3), (2 2, 2 1, 3 1))
class(sfc[2:3])
## [1] "sfc_MULTILINESTRING" "sfc"
class(st_cast(sfc[2:3]))
## [1] "sfc_MULTILINESTRING" "sfc"

gc1 <- st_geometrycollection(list(st_linestring(rbind(c(0,0),c(1,1),c(2,1)))))
gc2 <- st_geometrycollection(list(st_multilinestring(list(rbind(c(2,2),c(1,3)), rbind(c(0,0),c(1,1),c(2,1))))))
gc3 <- st_geometrycollection(list(st_multilinestring(list(rbind(c(4,4),c(4,3)), rbind(c(2,2),c(2,1),c(3,1))))))
(sfc <- st_sfc(gc1,gc2,gc3))
## Geometry set for 3 features 
## geometry type:  GEOMETRYCOLLECTION
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 4
## epsg (SRID):    NA
## proj4string:    NA
## GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1, 2 1))
## GEOMETRYCOLLECTION (MULTILINESTRING ((2 2, 1 3)...
## GEOMETRYCOLLECTION (MULTILINESTRING ((4 4, 4 3)...
class(st_cast(sfc))
## [1] "sfc_GEOMETRY" "sfc"
class(st_cast(st_cast(sfc), "MULTILINESTRING"))
## [1] "sfc_MULTILINESTRING" "sfc"

Affine transformations

Affine transformations are transformations of the type \(f(x) = xA + b\), where matrix \(A\) is used to flatten, scale and/or rotate, and \(b\) to translate \(x\). Low-level examples are:

(p = st_point(c(0,2)))
## POINT (0 2)
p + 1
## POINT (1 3)
p + c(1,2)
## POINT (1 4)
p + p
## POINT (0 4)
p * p
## POINT (0 4)
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
p * rot(pi/4)
## POINT (1.414214 1.414214)
p * rot(pi/2)
## POINT (2 1.224647e-16)
p * rot(pi)
## POINT (2.449294e-16 -2)

Just to make the point, we can for instance rotate the counties of North Carolina 90 degrees clockwise around their centroid, and shrink them to 75% of their original size:

nc = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
ncg = st_geometry(nc)
plot(ncg, border = 'grey')
cntrd = st_centroid(ncg)
## Warning in st_centroid.sfc(ncg): st_centroid does not give correct
## centroids for longitude/latitude data
ncg2 = (ncg - cntrd) * rot(pi/2) * .75 + cntrd
plot(ncg2, add = TRUE)
plot(cntrd, col = 'red', add = TRUE, cex = .5)