2.6 Geometry types

Do you remember then?

They can be in well-known binary (WKB) or well-known text (WKT) standard. (the first one is a bit hard to read!). st_as_binary provide a way to convert to binary and it provides intersesting examples regarding precision.

The standard look like that (more here):

POINT (5 2)
LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)
POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5)) # w/o hole
POLYGON ((0 0 0,4 0 0,4 4 0,0 4 0,0 0 0),(1 1 0,2 1 0,2 2 0,1 2 0,1 1 0)) # with hole / interior boundaries
MULTIPOINT (5 2, 1 3, 3 4, 3 2)
MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5), (0 2, 1 2, 1 3, 0 3, 0 2)))
GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2), LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))

2.6.1 In sf: the workflow:

knitr::include_graphics("https://raw.githubusercontent.com/geocompx/geocompr/main/figures/02-sfdiagram.png")

lnd_point = st_point(c(0.1, 51.5))                 # sfg object
lnd_geom = st_sfc(lnd_point, crs = "EPSG:4326")    # sfc object
lnd_attrib = data.frame(                           # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
  )
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object
lnd_sf

## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
## Geodetic CRS:  WGS 84
##     name temperature       date         geometry
## 1 London          25 2017-06-21 POINT (0.1 51.5)

2.6.2 First is building geometry (sfg)

st_point(c(5, 2, 3))              # XYZ point

## POINT Z (5 2 3)

st_point(c(5, 2, 3, 1))           # XYZM point

## POINT ZM (5 2 3 1)

# stop at the M diemsion
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix)

## LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)

polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
st_polygon(polygon_list)

## POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))

## POLYGON with a hole
polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list)

## POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5), (2 4, 3 4, 3 3, 2 3, 2 4))

It can be a bit tedious but this is usefull when you want to build some small/shareable example.

2.6.3 Simple feature columns (sfc)

We move from one feature to more:

# sfc POINT
point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc

## Geometry set for 2 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1 ymin: 2 xmax: 5 ymax: 3
## CRS:           NA

## POINT (5 2)

## POINT (1 3)

st_crs()

## Coordinate Reference System: NA

All geometries in sfc objects must have the same CRS

it can be provided when building sfc

# Set the CRS with an identifier referring to an 'EPSG' CRS code:
points_sfc_wgs = st_sfc(point1, point2, crs = "EPSG:4326")
st_crs(points_sfc_wgs)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

# TODO erase EPSG and see what we have

2.6.4 sfheaders package

fully compatible with sf
every function come with a c++ implementation (header -> hence the name) and can be linked in your package
it is also a bit less “verbose”

# matrices
m = matrix(1:8, ncol = 2)
sfheaders::sfg_linestring(obj = m)

## LINESTRING (1 5, 2 6, 3 7, 4 8)

#> LINESTRING (1 5, 2 6, 3 7, 4 8)
# data.frames
df = data.frame(x = 1:4, y = 4:1)
sfheaders::sfg_polygon(obj = df)

## POLYGON ((1 4, 2 3, 3 2, 4 1, 1 4))

#> POLYGON ((1 4, 2 3, 3 2, 4 1, 1 4))

sfheaders is also very good at “casting”: ie converting between types (cf later!)

2.6.5 Spherical geometry with s2

It is an implementation of a Discrete Global Grid System (other example is H3)

sf_use_s2()

## [1] TRUE

india_buffer_with_s2 = st_buffer(india, 1)
sf_use_s2(FALSE)

## Spherical geometry (s2) switched off

#> Spherical geometry (s2) switched off
india_buffer_without_s2 = st_buffer(india, 1)

## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
## endCapStyle, : st_buffer does not correctly buffer longitude/latitude data

## dist is assumed to be in decimal degrees (arc_degrees).

#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).