- Simple Features for R

- What is a feature?
- How simple features in R are organized
- sf: objects with simple features
- sfc: simple feature geometry list-column
- Mixed geometry types
- sfg: simple feature geometry
- Well-known text, well-known binary, precision
- Reading and writing
- Coordinate reference systems and transformations
- Conversion, including to and from sp
- Geometrical operations
- Non-valid geometries

- Units
- How attributes relate to geometries

Simple features or *simple feature access* refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers, with emphasis on the *spatial* geometry of these objects. It also describes how such objects can be stored in and retrieved from databases, and which geometrical operations should be defined for them.

The standard is widely implemented in spatial databases (such as PostGIS), commercial GIS (e.g., ESRI ArcGIS) and forms the vector data basis for libraries such as GDAL. A subset of simple features forms the GeoJSON standard.

R has well-supported classes for storing spatial data (sp) and interfacing to the above mentioned environments (rgdal, rgeos), but has so far lacked a complete implementation of simple features, making conversions at times convoluted, inefficient or incomplete. The package sf tries to fill this gap, and aims at succeeding sp in the long term.

This vignette:

- explains what is meant by features, and by simple features
- shows how they are implemented in R
- provides examples of how you can work with them
- shows how they can be read from and written to external files or resources
- discusses how they can be converted to and from sp objects
- shows how they can be used for meaningful spatial analysis

A feature is thought of as a thing, or an object in the real world, such as a building or a tree. As is the case with objects, they often consist of other objects. This is the case with features too: a set of features can form a single feature. A forest stand can be a feature, a forest can be a feature, a city can be a feature. A satellite image pixel can be a feature, a complete image can be a feature too.

Features have a *geometry* describing *where* on Earth the feature is located, and they have attributes, which describe other properties. The geometry of a tree can be the delineation of its crown, of its stem, or the point indicating its centre. Other properties may include its height, color, diameter at breast height at a particular date, and so on.

The standard says: “*A simple feature is defined by the OpenGIS Abstract specification to have both spatial and non-spatial attributes. Spatial attributes are geometry valued, and simple features are based on 2D geometry with linear interpolation between vertices.*” We will see soon that the same standard will extend its coverage beyond 2D and beyond linear interpolation. Here, we take simple features as the data structures and operations described in the standard.

All geometries are composed of points. Points are coordinates in a 2-, 3- or 4-dimensional space. All points in a geometry have the same dimensionality. In addition to X and Y coordinates, there are two optional additional dimensions:

- a Z coordinate, denoting altitude
- an M coordinate (rarely used), denoting some
*measure*that is associated with the point, rather than with the feature as a whole (in which case it would be a feature attribute); examples could be time of measurement, or measurement error of the coordinates

The four possible cases then are:

- two-dimensional points refer to x and y, easting and northing, or longitude and latitude, we refer to them as XY
- three-dimensional points as XYZ
- three-dimensional points as XYM
- four-dimensional points as XYZM (the third axis is Z, fourth M)

The following seven simple feature types are the most common, and are for instance the only ones used for GeoJSON:

type | description |
---|---|

`POINT` |
zero-dimensional geometry containing a single point |

`LINESTRING` |
sequence of points connected by straight, non-self intersecting line pieces; one-dimensional geometry |

`POLYGON` |
geometry with a positive area (two-dimensional); sequence of points form a closed, non-self intersecting ring; the first ring denotes the exterior ring, zero or more subsequent rings denote holes in this exterior ring |

`MULTIPOINT` |
set of points; a MULTIPOINT is simple if no two Points in the MULTIPOINT are equal |

`MULTILINESTRING` |
set of linestrings |

`MULTIPOLYGON` |
set of polygons |

`GEOMETRYCOLLECTION` |
set of geometries of any type except GEOMETRYCOLLECTION |

Each of the geometry types can also be a (typed) empty set, containing zero coordinates (for `POINT`

the standard is not clear how to represent the empty geometry). Empty geometries can be thought of being the analogue to missing (`NA`

) attributes, NULL values or empty lists.

The remaining geometries 10 are more rare, but increasingly find implementations:

type | description |
---|---|

`CIRCULARSTRING` |
The CIRCULARSTRING is the basic curve type, similar to a LINESTRING in the linear world. A single segment requires three points, the start and end points (first and third) and any other point on the arc. The exception to this is for a closed circle, where the start and end points are the same. In this case the second point MUST be the center of the arc, ie the opposite side of the circle. To chain arcs together, the last point of the previous arc becomes the first point of the next arc, just like in LINESTRING. This means that a valid circular string must have an odd number of points greated than 1. |

`COMPOUNDCURVE` |
A compound curve is a single, continuous curve that has both curved (circular) segments and linear segments. That means that in addition to having well-formed components, the end point of every component (except the last) must be coincident with the start point of the following component. |

`CURVEPOLYGON` |
Example compound curve in a curve polygon: CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,2 0, 2 1, 2 3, 4 3),(4 3, 4 5, 1 4, 0 0)), CIRCULARSTRING(1.7 1, 1.4 0.4, 1.6 0.4, 1.6 0.5, 1.7 1) ) |

`MULTICURVE` |
A MultiCurve is a 1-dimensional GeometryCollection whose elements are Curves, it can include linear strings, circular strings or compound strings. |

`MULTISURFACE` |
A MultiSurface is a 2-dimensional GeometryCollection whose elements are Surfaces, all using coordinates from the same coordinate reference system. |

`CURVE` |
A Curve is a 1-dimensional geometric object usually stored as a sequence of Points, with the subtype of Curve specifying the form of the interpolation between Points |

`SURFACE` |
A Surface is a 2-dimensional geometric object |

`POLYHEDRALSURFACE` |
A PolyhedralSurface is a contiguous collection of polygons, which share common boundary segments |

`TIN` |
A TIN (triangulated irregular network) is a PolyhedralSurface consisting only of Triangle patches. |

`TRIANGLE` |
A Triangle is a polygon with 3 distinct, non-collinear vertices and no interior boundary |

Note that `CIRCULASTRING`

, `COMPOUNDCURVE`

and `CURVEPOLYGON`

are not described in the SFA standard, but in the SQL-MM part 3 standard. The descriptions above were copied from the PostGIS manual.

Coordinates can only be placed on the Earth’s surface when their coordinate reference system (CRS) is known; this may be an spheroid CRS such as WGS84, a projected, two-dimensional (Cartesian) CRS such as a UTM zone or Web Mercator, or a CRS in three-dimensions, or including time. Similarly, M-coordinates need an attribute reference system, e.g. a measurement unit.

Package `sf`

represents simple features as native R objects. Similar to PostGIS, all functions and methods in `sf`

that operate on spatial data are prefixed by `st_`

, which refers to *spatial and temporal*; this makes them easily findable by command-line completion. Simple features are implemented as R native data, using simple data structures (S3 classes, lists, matrix, vector). Typical use involves reading, manipulating and writing of sets of features, with attributes and geometries.

As attributes are typically stored in `data.frame`

objects (or the very similar `tbl_df`

), we will also store feature geometries in a `data.frame`

column. Since geometries are not single-valued, they are put in a list-column, a list of length equal to the number of records in the `data.frame`

, with each list element holding the simple feature geometry of that feature. The three classes used to represent simple features are:

`sf`

, the table (`data.frame`

) with feature attributes and feature geometries, which contains`sfc`

, the list-column with the geometries for each feature (record), which is composed of`sfg`

, the feature geometry of an individual simple feature.

We will now discuss each of these three classes.

As we usually do not work with geometries of single simple features, but with datasets consisting of sets of features with attributes, the two are put together in `sf`

(simple feature) objects. The following command reads the `nc`

dataset from a file that is contained in the `sf`

package:

```
library(sf)
## Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3
nc <- st_read(system.file("shape/nc.shp", package="sf"))
## Reading layer `nc' from data source `/tmp/Rtmpc6p3fp/Rinst12262f308bb5/sf/shape/nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID): 4267
## proj4string: +proj=longlat +datum=NAD27 +no_defs
```

(Note that users will not use `system.file`

but give a `filename`

directly, and that shapefiles consist of more than one file, all with identical basename, which reside in the same directory.) The short report printed gives the file name, the driver (ESRI Shapefile), mentions that there are 100 features (records, represented as rows) and 14 fields (attributes, represented as columns). This object is of class

```
class(nc)
## [1] "sf" "data.frame"
```

meaning it extends (and “is” a) `data.frame`

, but with a single list-column with geometries, which is held in the column with name

```
attr(nc, "sf_column")
## [1] "geometry"
```

If we print the first three features, we see their attribute values and an abridged version of the geometry

`print(nc[9:15], n = 3)`

which would give the following output:

In the output we see:

- in green a simple feature: a single record, or
`data.frame`

row, consisting of attributes and geometry - in blue a single simple feature geometry (an object of class
`sfg`

) - in red a simple feature list-column (an object of class
`sfc`

, which is a column in the`data.frame`

) - that although geometries are native R objects, they are printed as well-known text

Methods for `sf`

objects are

```
methods(class = "sf")
## [1] $<- [ [[<-
## [4] aggregate cbind coerce
## [7] identify initialize merge
## [10] plot print rbind
## [13] show slotsFromS3 st_agr
## [16] st_agr<- st_as_sf st_bbox
## [19] st_boundary st_buffer st_cast
## [22] st_centroid st_collection_extract st_convex_hull
## [25] st_coordinates st_crs st_crs<-
## [28] st_difference st_geometry st_geometry<-
## [31] st_intersection st_is st_line_merge
## [34] st_node st_point_on_surface st_polygonize
## [37] st_precision st_segmentize st_set_precision
## [40] st_simplify st_snap st_sym_difference
## [43] st_transform st_triangulate st_union
## [46] st_voronoi st_wrap_dateline st_write
## [49] st_zm
## see '?methods' for accessing help and source code
```

It is also possible to create `data.frame`

objects with geometry list-columns that are not of class `sf`

, e.g. by

```
nc.no_sf <- as.data.frame(nc)
class(nc.no_sf)
## [1] "data.frame"
```

However, such objects:

- no longer register which column is the geometry list-column
- no longer have a plot method, and
- lack all of the other dedicated methods listed above for class
`sf`

The column in the `sf`

data.frame that contains the geometries is a list, of class `sfc`

. We can retrieve the geometry list-column in this case by `nc$geom`

or `nc[[15]]`

, but the more general way uses `st_geometry`

:

```
(nc_geom <- st_geometry(nc))
## Geometry set for 100 features
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID): 4267
## proj4string: +proj=longlat +datum=NAD27 +no_defs
## First 5 geometries:
## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...
## MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...
## MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...
## MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...
## MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...
```

Geometries are printed in abbreviated form, but we can can view a complete geometry by selecting it, e.g. the first one by

```
nc_geom[[1]]
## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))
```

The way this is printed is called *well-known text*, and is part of the standards. The word `MULTIPOLYGON`

is followed by three parenthesis, because it can consist of multiple polygons, in the form of `MULTIPOLYGON(POL1,POL2)`

, where `POL1`

might consist of an exterior ring and zero or more interior rings, as of `(EXT1,HOLE1,HOLE2)`

. Sets of coordinates are held together with parenthesis, so we get `((crds_ext)(crds_hole1)(crds_hole2))`

where `crds_`

is a comma-separated set of coordinates of a ring. This leads to the case above, where `MULTIPOLYGON(((crds_ext)))`

refers to the exterior ring (1), without holes (2), of the first polygon (3) - hence three parentheses.

We can see there is a single polygon with no rings:

```
par(mar = c(0,0,1,0))
plot(nc[1])
plot(nc[1,1], col = 'grey', add = TRUE)
```