The R Script associated with this page is available here. Download this file and open it (or copy-paste into a new script) with RStudio so you can follow along.

Setup

You may need to install a few of these packages with install.packages() or using the GUI. Note, much of the material below was taken from the sf vignettes available here.

Load packages

library(rgdal)
library(rgeos)
## rgeos version: 0.3-23, (SVN revision 546)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.2-4 
##  Polygon checking: TRUE
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(maptools)

Background

There are currently two main approaches in R to handle geographic vector data.

The sp package

The first package to provide classes and methods for spatial data types in R is called sp1. Development of the sp package began in the early 2000s in an attempt to standardize how spatial data would be treated in R and to allow for better interoperability between different analysis packages that use spatial data. The package (first release on CRAN in 2005) provides classes and methods to create points, lines, polygons, and grids and to operate on them. About 350 of the spatial analysis packages use the spatial data types that are implemented in sp i.e. they “depend” on the sp package and many more are indirectly dependent.

The foundational structure for any spatial object in sp is the Spatial class. It has two “slots” (new-style S4 class objects in R have pre-defined components called slots):

  • a bounding box

  • a CRS class object to define the Coordinate Reference System

This basic structure is then extended, depending on the characteristics of the spatial object (point, line, polygon).

To build up a spatial object in sp we could follow these steps:

I. Create geometric objects (topology)

Points (which may have 2 or 3 dimensions) are the most basic spatial data objects. They are generated out of either a single coordinate or a set of coordinates, like a two-column matrix or a dataframe with a column for latitude and one for longitude.
Lines are generated out of Line objects. A Line object is a spaghetti collection of 2D coordinates2 and is generated out of a two-column matrix or a dataframe with a column for latitude and one for longitude. A Lines object is a list of one or more Line objects, for example all the contours at a single elevation.
Polygons are generated out of Polygon objects. A Polygon object is a spaghetti collection of 2D coordinates with equal first and last coordinates and is generated out of a two-column matrix or a dataframe with a column for latitude and one for longitude. A Polygons object is a list of one or more Polygon objects, for example islands belonging to the same country.

See here for a very simple example for how to create a Line object:

ln <- Line(matrix(runif(6), ncol=2))
str(ln)
## Formal class 'Line' [package "sp"] with 1 slot
##   ..@ coords: num [1:3, 1:2] 0.815 0.636 0.878 0.783 0.179 ...

See here for a very simple example for how to create a Lines object:

lns <- Lines(list(ln), ID = "a") # this contains just one Line!
str(lns)
## Formal class 'Lines' [package "sp"] with 2 slots
##   ..@ Lines:List of 1
##   .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
##   .. .. .. ..@ coords: num [1:3, 1:2] 0.815 0.636 0.878 0.783 0.179 ...
##   ..@ ID   : chr "a"

II. Create spatial objects Spatial* object (* stands for Points, Lines, or Polygons).

This step adds the bounding box (automatically) and the slot for the Coordinate Reference System or CRS (which needs to be filled with a value manually). SpatialPoints can be directly generated out of the coordinates. SpatialLines and SpatialPolygons objects are generated using lists of Lines or Polygons objects respectively (more below).

See here for how to create a SpatialLines object:

sp_lns <- SpatialLines(list(lns))
str(sp_lns)
## Formal class 'SpatialLines' [package "sp"] with 3 slots
##   ..@ lines      :List of 1
##   .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
##   .. .. .. ..@ Lines:List of 1
##   .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
##   .. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 0.815 0.636 0.878 0.783 0.179 ...
##   .. .. .. ..@ ID   : chr "a"
##   ..@ bbox       : num [1:2, 1:2] 0.636 0.179 0.878 0.783
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

III. Add attributes (Optional:)

Add a data frame with attribute data, which will turn your Spatial* object into a Spatial*DataFrame object. The points in a SpatialPoints object may be associated with a row of attributes to create a SpatialPointsDataFrame object. The coordinates and attributes may, but do not have to be keyed to each other using ID values.
SpatialLinesDataFrame and SpatialPolygonsDataFrame objects are defined using SpatialLines and SpatialPolygons objects and data frames. The ID fields are here required to match the data frame row names.

See here for how to create a SpatialLinesDataframe:

dfr <- data.frame(id = "a", use = "road", cars_per_hour = 10) # note how we use the ID from above!
sp_lns_dfr <- SpatialLinesDataFrame(sp_lns, dfr, match.ID = "id")
str(sp_lns_dfr)
## Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  1 obs. of  3 variables:
##   .. ..$ id           : Factor w/ 1 level "a": 1
##   .. ..$ use          : Factor w/ 1 level "road": 1
##   .. ..$ cars_per_hour: num 10
##   ..@ lines      :List of 1
##   .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
##   .. .. .. ..@ Lines:List of 1
##   .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
##   .. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 0.815 0.636 0.878 0.783 0.179 ...
##   .. .. .. ..@ ID   : chr "a"
##   ..@ bbox       : num [1:2, 1:2] 0.636 0.179 0.878 0.783
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

A number of spatial methods are available for the classes in sp. Common ones include:

function and what it does
bbox() returns the bounding box coordinates
proj4string() sets or retrieves projection attributes using the CRS object.
CRS() creates an object of class of coordinate reference system arguments
spplot() plots a separate map of all the attributes unless specified otherwise
coordinates() set or retrieve the spatial coordinates. For spatial polygons it returns the centroids.
over(a, b) used for example to retrieve the polygon or grid indices on a set of points
spsample() sampling of spatial points within the spatial extent of objects

The sf package

The second package, first released on CRAN in late October 2016, is called sf3. It implements a formal standard called “Simple Features” that specifies a storage and access model of spatial geometries (point, line, polygon). A feature geometry is called simple when it consists of points connected by straight line pieces, and does not intersect itself. This standard has been adopted widely, not only by spatial databases such as PostGIS, but also more recent standards such as GeoJSON.

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.

If you work with PostGis or GeoJSON you may have come across the WKT (well-known text) format, for example like these:

POINT (30 10)
LINESTRING (30 10, 10 30, 40 40)
POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))

sf implements this standard natively in R. Data are structured and conceptualized very differently from the sp approach.

In sf spatial objects are stored as a simple data frame with a special column that contains the information for the geographic coordinates. That special column is a list with the same length as the number of rows in the data frame. Each of the individual list elements then can be of any length needed to hold the coordinates that correspond to an individual feature.

To create a spatial object manually the basic steps would be:

I. Create geometric objects (topology)

Geometric objects (simple features) can be created from a numeric vector, matrix or a list with the coordinates. They are called sfg objects for Simple Feature Geometry.

See here for an example of how a LINESTRING sfg object is created:

lnstr_sfg <- st_linestring(matrix(runif(6), ncol=2)) 
class(lnstr_sfg)
## [1] "XY"         "LINESTRING" "sfg"

II. Combine all individual single feature objects for the special column.

In order to work our way towards a data frame for all features we create what is called an sfc object with all individual features, which stands for Simple Feature Collection. The sfc object also holds the bounding box and the projection information.

See here for an example of how a sfc object is created:

(lnstr_sfc <- st_sfc(lnstr_sfg)) # just one feature here
## Geometry set for 1 feature 
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 0.04144419 ymin: 0.03059071 xmax: 0.9553763 ymax: 0.2406219
## epsg (SRID):    NA
## proj4string:    NA
## LINESTRING (0.955376301659271 0.03059071325697,...
class(lnstr_sfc) 
## [1] "sfc_LINESTRING" "sfc"

III. Add attributes.

We now combine the dataframe with the attributes and the simple feature collection. See here how its done.

(lnstr_sf <- st_sf(dfr , lnstr_sfc))
## Simple feature collection with 1 feature and 3 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 0.04144419 ymin: 0.03059071 xmax: 0.9553763 ymax: 0.2406219
## epsg (SRID):    NA
## proj4string:    NA
##   id  use cars_per_hour                      lnstr_sfc
## 1  a road            10 LINESTRING (0.9553763016592...
class(lnstr_sf)
## [1] "sf"         "data.frame"

There are many methods available in the sf package, to find out use methods(class="sp")

Here are some of the other highlights of sf you might be interested in:

  • provides fast I/O, particularly relevant for large files
  • directly reads from and writes to spatial databases such as PostGIS
  • stay tuned for a new ggplot release that will be able to read and plot the sf format without the need of conversion to a data frame, like the sp format

Note that sp and sf are not the only way spatial objects are conceptualized in R. Other spatial packages may use their own class definitions for spatial data (for example spatstat). Usually you can find functions that convert sp and increasingly sf objects to and from these formats.

sf: objects with simple features

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:

file=system.file("shape/nc.shp", package="sf")
file
## [1] "/Library/Frameworks/R.framework/Versions/3.4/Resources/library/sf/shape/nc.shp"
nc <- st_read(file)
## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.4/Resources/library/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 you will rarely use system.file but instead 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           anti_join           arrange            
##  [7] cbind               coerce              distinct           
## [10] filter              full_join           gather             
## [13] group_by            initialize          inner_join         
## [16] left_join           merge               mutate             
## [19] nest                plot                print              
## [22] rbind               rename              right_join         
## [25] sample_frac         sample_n            select             
## [28] semi_join           separate            show               
## [31] slice               slotsFromS3         spread             
## [34] st_agr              st_agr<-            st_as_sf           
## [37] st_bbox             st_boundary         st_buffer          
## [40] st_cast             st_centroid         st_convex_hull     
## [43] st_coordinates      st_crs              st_crs<-           
## [46] st_difference       st_geometry         st_geometry<-      
## [49] st_intersection     st_is               st_line_merge      
## [52] st_make_valid       st_point_on_surface st_polygonize      
## [55] st_precision        st_segmentize       st_set_precision   
## [58] st_simplify         st_split            st_sym_difference  
## [61] st_transform        st_triangulate      st_union           
## [64] st_voronoi          st_zm               summarise          
## [67] transmute           ungroup             unite              
## 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

sfc: simple feature geometry list-column

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)

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.4727554321289 36.2343559265137, -81.5408401489258 36.2725067138672, -81.5619812011719 36.2735939025879, -81.6330642700195 36.3406867980957, -81.7410736083984 36.3917846679688, -81.6982803344727 36.4717788696289, -81.7027969360352 36.5193405151367, -81.6699981689453 36.5896492004395, -81.3452987670898 36.5728645324707, -81.347541809082 36.537914276123, -81.3247756958008 36.5136795043945, -81.3133239746094 36.4806976318359, -81.2662353515625 36.4372062683105, -81.2628402709961 36.4050407409668, -81.2406921386719 36.3794174194336, -81.2398910522461 36.365364074707, -81.2642440795898 36.3524131774902, -81.3289947509766 36.3635025024414, -81.3613739013672 36.3531608581543, -81.3656921386719 36.3390502929688, -81.354133605957 36.2997169494629, -81.3674545288086 36.2786979675293, -81.4063873291016 36.2850532531738, -81.4123306274414 36.2672920227051, -81.431037902832 36.2607192993164, -81.4528884887695 36.2395858764648, -81.4727554321289 36.2343559265137)))

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:

plot(nc[1])
plot(nc[1,1], col = 'grey', add = TRUE)

but some of the polygons in this dataset have multiple exterior rings; they can be identified by

w <- which(sapply(nc_geom, length) > 1)
plot(nc[w,1], col = 2:7)

Following the MULTIPOLYGON datastructure, in R we have a list of lists of lists of matrices. For instance, we get the coordinate pairs of the second exterior ring (first ring is always exterior) for the geometry of feature 4 by

nc_geom[[4]][[2]][[1]]
##           [,1]     [,2]
## [1,] -76.02717 36.55672
## [2,] -75.99866 36.55665
## [3,] -75.91192 36.54253
## [4,] -75.92480 36.47398
## [5,] -75.97728 36.47802
## [6,] -75.97629 36.51793
## [7,] -76.02717 36.55672

Geometry columns have their own class,

class(nc_geom)
## [1] "sfc_MULTIPOLYGON" "sfc"

Methods for geometry list-columns include

methods(class = 'sfc')
##  [1] [                   [<-                 as.data.frame      
##  [4] c                   coerce              format             
##  [7] initialize          Ops                 print              
## [10] show                slotsFromS3         st_as_binary       
## [13] st_as_text          st_bbox             st_boundary        
## [16] st_buffer           st_cast             st_centroid        
## [19] st_convex_hull      st_coordinates      st_crs             
## [22] st_crs<-            st_difference       st_geometry        
## [25] st_intersection     st_is               st_line_merge      
## [28] st_make_valid       st_point_on_surface st_polygonize      
## [31] st_precision        st_segmentize       st_set_precision   
## [34] st_simplify         st_split            st_sym_difference  
## [37] st_transform        st_triangulate      st_union           
## [40] st_voronoi          st_zm               str                
## [43] summary             type_sum           
## see '?methods' for accessing help and source code

Coordinate reference systems (st_crs and st_transform) are discussed in the section on coordinate reference systems. st_as_wkb and st_as_text convert geometry list-columns into well-known-binary or well-known-text, explained below. st_bbox retrieves the coordinate bounding box.

Attributes include

attributes(nc_geom)
## $n_empty
## [1] 0
## 
## $crs
## $epsg
## [1] 4267
## 
## $proj4string
## [1] "+proj=longlat +datum=NAD27 +no_defs"
## 
## attr(,"class")
## [1] "crs"
## 
## $class
## [1] "sfc_MULTIPOLYGON" "sfc"             
## 
## $precision
## [1] 0
## 
## $bbox
##      xmin      ymin      xmax      ymax 
## -84.32385  33.88199 -75.45698  36.58965

Mixed geometry types

The class of nc_geom is c("sfc_MULTIPOLYGON", "sfc"): sfc is shared with all geometry types, and sfc_TYPE with TYPE indicating the type of the particular geometry at hand.

There are two “special” types: GEOMETRYCOLLECTION, and GEOMETRY. GEOMETRYCOLLECTION indicates that each of the geometries may contain a mix of geometry types, as in

mix <- st_sfc(st_geometrycollection(list(st_point(1:2))),
    st_geometrycollection(list(st_linestring(matrix(1:4,2)))))
class(mix)
## [1] "sfc_GEOMETRYCOLLECTION" "sfc"

Still, the geometries are here of a single type.

The second GEOMETRY, indicates that the geometries in the geometry list-column are of varying type:

mix <- st_sfc(st_point(1:2), st_linestring(matrix(1:4,2)))
class(mix)
## [1] "sfc_GEOMETRY" "sfc"

These two are fundamentally different: GEOMETRY is a superclass without instances, GEOMETRYCOLLECTION is a geometry instance. GEOMETRY list-columns occur when we read in a data source with a mix of geometry types. GEOMETRYCOLLECTION is a single feature’s geometry: the intersection of two feature polygons may consist of points, lines and polygons, see the example below.

sfg: simple feature geometry

Simple feature geometry (sfg) objects carry the geometry for a single feature, e.g. a point, linestring or polygon.

Simple feature geometries are implemented as R native data, using the following rules

  1. a single POINT is a numeric vector
  2. a set of points, e.g. in a LINESTRING or ring of a POLYGON is a matrix, each row containing a point
  3. any other set is a list

Creator functions are rarely used in practice, since we typically bulk read and write spatial data. They are useful for illustration:

(x <- st_point(c(1,2)))
## POINT (1 2)
str(x)
## Classes 'XY', 'POINT', 'sfg'  num [1:2] 1 2
(x <- st_point(c(1,2,3)))
## POINT Z (1 2 3)
str(x)
## Classes 'XYZ', 'POINT', 'sfg'  num [1:3] 1 2 3
(x <- st_point(c(1,2,3), "XYM"))
## POINT M (1 2 3)
str(x)
## Classes 'XYM', 'POINT', 'sfg'  num [1:3] 1 2 3
(x <- st_point(c(1,2,3,4)))
## POINT ZM (1 2 3 4)
str(x)
## Classes 'XYZM', 'POINT', 'sfg'  num [1:4] 1 2 3 4
st_zm(x, drop = TRUE, what = "ZM")
## POINT (1 2)

This means that we can represent 2-, 3- or 4-dimensional coordinates. All geometry objects inherit from sfg (simple feature geometry), but also have a type (e.g. POINT), and a dimension (e.g. XYM) class name. A figure illustrates six of the seven most common types.

With the exception of the POINT which has a single point as geometry, the remaining six common single simple feature geometry types that correspond to single features (single records, or rows in a data.frame) are created like this

p <- rbind(c(3.2,4), c(3,4.6), c(3.8,4.4), c(3.5,3.8), c(3.4,3.6), c(3.9,4.5))
(mp <- st_multipoint(p))
## MULTIPOINT (3.2 4, 3 4.6, 3.8 4.4, 3.5 3.8, 3.4 3.6, 3.9 4.5)
s1 <- rbind(c(0,3),c(0,4),c(1,5),c(2,5))
(ls <- st_linestring(s1))
## LINESTRING (0 3, 0 4, 1 5, 2 5)
s2 <- rbind(c(0.2,3), c(0.2,4), c(1,4.8), c(2,4.8))
s3 <- rbind(c(0,4.4), c(0.6,5))
(mls <- st_multilinestring(list(s1,s2,s3)))
## MULTILINESTRING ((0 3, 0 4, 1 5, 2 5), (0.2 3, 0.2 4, 1 4.8, 2 4.8), (0 4.4, 0.6 5))
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
(mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5))))
## MULTIPOLYGON (((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3)))
(gc <- st_geometrycollection(list(mp, mpol, ls)))
## GEOMETRYCOLLECTION (MULTIPOINT (3.2 4, 3 4.6, 3.8 4.4, 3.5 3.8, 3.4 3.6, 3.9 4.5), MULTIPOLYGON (((0 0, 1 0, 3 2, 2 4, 1 4, 0 0), (1 1, 1 2, 2 2, 1 1)), ((3 0, 4 0, 4 1, 3 1, 3 0), (3.3 0.3, 3.3 0.8, 3.8 0.8, 3.8 0.3, 3.3 0.3)), ((3 3, 4 2, 4 3, 3 3))), LINESTRING (0 3, 0 4, 1 5, 2 5))

The objects created are shown here:

Geometries can also be empty, as in

x <- st_geometrycollection()
length(x)
## [1] 0

Well-known text, well-known binary, precision

WKT and WKB

Well-known text (WKT) and well-known binary (WKB) are two encodings for simple feature geometries. Well-known text, e.g. seen in

x <- st_linestring(matrix(10:1,5))
st_as_text(x)
## [1] "LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)"

(but without the leading ## [1] and quotes), is human-readable. Coordinates are usually floating point numbers, and moving large amounts of information as text is slow and imprecise. For that reason, we use well-known binary (WKB) encoding

st_as_binary(x)
##  [1] 01 02 00 00 00 05 00 00 00 00 00 00 00 00 00 24 40 00 00 00 00 00 00
## [24] 14 40 00 00 00 00 00 00 22 40 00 00 00 00 00 00 10 40 00 00 00 00 00
## [47] 00 20 40 00 00 00 00 00 00 08 40 00 00 00 00 00 00 1c 40 00 00 00 00
## [70] 00 00 00 40 00 00 00 00 00 00 18 40 00 00 00 00 00 00 f0 3f

WKT and WKB can both be transformed back into R native objects by

st_as_sfc("LINESTRING(10 5, 9 4, 8 3, 7 2, 6 1)")[[1]]
## LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)
st_as_sfc(structure(list(st_as_binary(x)), class = "WKB"))[[1]]
## LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)

GDAL, GEOS, spatial databases and GIS read and write WKB which is fast and precise. Conversion between R native objects and WKB is done by package sf in compiled (C++/Rcpp) code, making this a reusable and fast route for I/O of simple feature geometries in R.

Reading and writing

As we’ve seen above, reading spatial data from an external file can be done by

filename <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filename)
## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.4/Resources/library/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

we can suppress the output by adding argument quiet=TRUE or by using the otherwise nearly identical but more quiet

nc <- read_sf(filename)

Writing takes place in the same fashion, using st_write:

st_write(nc, "nc.shp")
## Writing layer `nc' to data source `nc.shp' using driver `ESRI Shapefile'
## features:       100
## fields:         14
## geometry type:  Multi Polygon

If we repeat this, we get an error message that the file already exists, and we can overwrite by

st_write(nc, "nc.shp", delete_layer = TRUE)
## Deleting layer `nc' using driver `ESRI Shapefile'
## Writing layer `nc' to data source `/Users/adamw/Documents/repos/RDataScience/nc.shp' using driver `ESRI Shapefile'
## features:       100
## fields:         14
## geometry type:  Multi Polygon

or its quiet alternative that does this by default,

write_sf(nc, "nc.shp") # silently overwrites

Driver-specific options

The dsn and layer arguments to st_read and st_write denote a data source name and optionally a layer name. Their exact interpretation as well as the options they support vary per driver, the GDAL driver documentation is best consulted for this. For instance, a PostGIS table in database postgis might be read by

meuse <- st_read("PG:dbname=postgis", "meuse")

where the PG: string indicates this concerns the PostGIS driver, followed by database name, and possibly port and user credentials. When the layer and driver arguments are not specified, st_read tries to guess them from the datasource, or else simply reads the first layer, giving a warning in case there are more.

st_read typically reads the coordinate reference system as proj4string, but not the EPSG (SRID). GDAL cannot retrieve SRID (EPSG code) from proj4string strings, and, when needed, it has to be set by the user. See also the section on crs.

st_drivers() returns a data.frame listing available drivers, and their metadata: names, whether a driver can write, and whether it is a raster and/or vector driver. All drivers can read. Reading of some common data formats is illustrated below:

st_layers(dsn) lists the layers present in data source dsn, and gives the number of fields, features and geometry type for each layer:

st_layers(system.file("osm/overpass.osm", package="sf"))
## Driver: OSM 
## Available layers:
##         layer_name       geometry_type features fields
## 1           points               Point       NA     10
## 2            lines         Line String       NA      9
## 3 multilinestrings   Multi Line String       NA      4
## 4    multipolygons       Multi Polygon       NA     25
## 5  other_relations Geometry Collection       NA      4

we see that in this case, the number of features is NA because for this xml file the whole file needs to be read, which may be costly for large files. We can force counting by

Sys.setenv(OSM_USE_CUSTOM_INDEXING="NO")
st_layers(system.file("osm/overpass.osm", package="sf"), do_count = TRUE)
## Driver: OSM 
## Available layers:
##         layer_name       geometry_type features fields
## 1           points               Point        1     10
## 2            lines         Line String        0      9
## 3 multilinestrings   Multi Line String        0      4
## 4    multipolygons       Multi Polygon        0     25
## 5  other_relations Geometry Collection        0      4

Another example of reading kml and kmz files is:

# Download .shp data
u_shp <- "http://coagisweb.cabq.gov/datadownload/biketrails.zip"
download.file(u_shp, "biketrails.zip")
unzip("biketrails.zip")
u_kmz <- "http://coagisweb.cabq.gov/datadownload/BikePaths.kmz"
download.file(u_kmz, "BikePaths.kmz")
# Read file formats
biketrails_shp <- st_read("biketrails.shp")
if(Sys.info()[1] == "Linux") # may not work if not Linux
  biketrails_kmz <- st_read("BikePaths.kmz")
u_kml = "http://www.northeastraces.com/oxonraces.com/nearme/safe/6.kml"
download.file(u_kml, "bikeraces.kml")
bikraces <- st_read("bikeraces.kml")

Create, read, update and delete

GDAL provides the crud (create, read, update, delete) functions to persistent storage. st_read (or read_sf) are used for reading. st_write (or write_sf) creates, and has the following arguments to control update and delete:

  • update=TRUE causes an existing data source to be updated, if it exists; this options is by default TRUE for all database drivers, where the database is updated by adding a table.
  • delete_layer=TRUE causes st_write try to open the the data source and delete the layer; no errors are given if the data source is not present, or the layer does not exist in the data source.
  • delete_dsn=TRUE causes st_write to delete the data source when present, before writing the layer in a newly created data source. No error is given when the data source does not exist. This option should be handled with care, as it may wipe complete directories or databases.

Directly reading and writing to spatial databases

Two further functions, st_read_db and st_write_db attempt to read and write from spatial databases, directly reading WKB or WKT without using GDAL. The advantage over st_read may be that instead of a complete table, the result of a (spatial) query may be fetched, limiting the amount of data that is read into R, and potentially benefiting from the spatial index of the database. Although intended to use the DBI interface, current use and testing of these functions are limited to PostGIS.

Coordinate reference systems and transformations

Coordinate reference systems (CRS) are like measurement units for coordinates: they specify which location on Earth a particular coordinate pair refers to. We saw above that sfc objects (geometry list-columns) have two attributes to store a CRS: epsg and proj4string. This implies that all geometries in a geometry list-column must have the same CRS. Both may be NA, e.g. in case the CRS is unknown, or when we work with local coordinate systems (e.g. inside a building, a body, or an abstract space).

proj4string is a generic, string-based description of a CRS, understood by the PROJ.4 library. It defines projection types and (often) defines parameter values for particular projections, and hence can cover an infinite amount of different projections. This library (also used by GDAL) provides functions to convert or transform between different CRS. epsg is the integer ID for a particular, known CRS that can be resolved into a proj4string. There is no (known, simple and general) way to resolve proj4string values into epsg IDs.

The importance of having epsg values stored with data besides proj4string values is that the epsg refers to particular, well-known CRS, whose parameters may change (improve) over time; fixing only the proj4string may remove the possibility to benefit from such improvements, and limit the provenance of datasets.

Coordinate reference system transformations can be carried out using st_transform, e.g. converting longitudes/latitudes in NAD27 to web mercator (EPSG:3857) can be done by

nc.web_mercator <- st_transform(nc, 3857)
st_geometry(nc.web_mercator)[[4]][[2]][[1]][1:3,]
##          [,1]    [,2]
## [1,] -8463306 4377498
## [2,] -8460133 4377489
## [3,] -8450476 4375532

Conversion, including to and from sp

sf objects and objects deriving from Spatial (package sp) can be coerced both ways:

# anticipate that sp::CRS will expand proj4strings:
p4s <- "+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat"
st_crs(nc) <- p4s
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
# anticipate geometry column name changes:
names(nc)[15] = "geometry"
attr(nc, "sf_column") = "geometry"
nc.sp <- as(nc, "Spatial")
class(nc.sp)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
nc2 <- st_as_sf(nc.sp)
#all.equal(nc, nc2)

As the Spatial* objects only support MULTILINESTRING and MULTIPOLYGON, LINESTRING and POLYGON geometries are automatically coerced into their MULTI form. When converting Spatial* into sf, if all geometries consist of a single POLYGON (possibly with holes), a POLYGON and otherwise all geometries are returned as MULTIPOLYGON: a mix of POLYGON and MULTIPOLYGON (such as common in shapefiles) is not created. Argument forceMulti=TRUE will override this, and create MULTIPOLYGONs in all cases. For LINES the situation is identical.

Geometrical operations

The standard for simple feature access defines a number of geometrical operations.

st_is_valid and st_is_simple return a boolean indicating whether a geometry is valid or simple.

st_is_valid(nc[1:2,])
## [1] TRUE TRUE

st_distance returns a dense numeric matrix with distances between geometries. st_relate returns a character matrix with the DE9-IM values for each pair of geometries:

x = st_transform(nc, 32119)
st_distance(x[c(1,4,22),], x[c(1, 33,55,56),])
## Units: m
##           [,1]     [,2]      [,3]     [,4]
## [1,]      0.00 312176.2 128338.51 475608.8
## [2,] 440548.35 114938.1 590417.79      0.0
## [3,]  18943.74 352708.6  78754.75 517511.6
st_relate(nc[1:5,], nc[1:4,])
## although coordinates are longitude/latitude, it is assumed that they are planar
##      [,1]        [,2]        [,3]        [,4]       
## [1,] "2FFF1FFF2" "FF2F11212" "FF2FF1212" "FF2FF1212"
## [2,] "FF2F11212" "2FFF1FFF2" "FF2F11212" "FF2FF1212"
## [3,] "FF2FF1212" "FF2F11212" "2FFF1FFF2" "FF2FF1212"
## [4,] "FF2FF1212" "FF2FF1212" "FF2FF1212" "2FFF1FFF2"
## [5,] "FF2FF1212" "FF2FF1212" "FF2FF1212" "FF2FF1212"

The commands st_intersects, st_disjoint, st_touches, st_crosses, st_within, st_contains, st_overlaps, st_equals, st_covers, st_covered_by, st_equals_exact and st_is_within_distance return a sparse matrix with matching (TRUE) indexes, or a full logical matrix:

st_intersects(nc[1:5,], nc[1:4,])
## although coordinates are longitude/latitude, it is assumed that they are planar
## [[1]]
## [1] 1 2
## 
## [[2]]
## [1] 1 2 3
## 
## [[3]]
## [1] 2 3
## 
## [[4]]
## [1] 4
## 
## [[5]]
## integer(0)
st_intersects(nc[1:5,], nc[1:4,], sparse = FALSE)
## although coordinates are longitude/latitude, it is assumed that they are planar
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE  TRUE FALSE FALSE
## [2,]  TRUE  TRUE  TRUE FALSE
## [3,] FALSE  TRUE  TRUE FALSE
## [4,] FALSE FALSE FALSE  TRUE
## [5,] FALSE FALSE FALSE FALSE

The commands st_buffer, st_boundary, st_convexhull, st_union_cascaded, st_simplify, st_triangulate, st_polygonize, st_centroid, st_segmentize, and st_union return new geometries, e.g.:

sel <- c(1,5,14)
geom = st_geometry(nc.web_mercator[sel,])
buf <- st_buffer(geom, dist = 30000)
plot(buf, border = 'red')
plot(geom, add = TRUE)
plot(st_buffer(geom, -5000), add = TRUE, border = 'blue')

Commands st_intersection, st_union, st_difference, st_sym_difference return new geometries that are a function of pairs of geometries:

par(mar = rep(0,4))
u <- st_union(nc)
plot(u)

The following code shows how computing an intersection between two polygons may yield a GEOMETRYCOLLECTION with a point, line and polygon:

opar <- par(mfrow = c(1, 2))
a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0))))
b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5,-0.5,-0.5,1,1))))
plot(a, ylim = c(-1,1))
title("intersecting two polygons:")
plot(b, add = TRUE, border = 'red')
(i <- st_intersection(a,b))
## GEOMETRYCOLLECTION (POINT (1 0), LINESTRING (4 0, 3 0), POLYGON ((5.5 0, 7 0, 7 -0.5, 6 -0.5, 5.5 0)))
plot(a, ylim = c(-1,1))
title("GEOMETRYCOLLECTION")
plot(b, add = TRUE, border = 'red')
plot(i, add = TRUE, col = 'green', lwd = 2)

par(opar)

Non-valid geometries

Invalid geometries are for instance self-intersecting lines (left) or polygons with slivers (middle) or self-intersections (right).

x1 <- st_linestring(cbind(c(0,1,0,1),c(0,1,1,0)))
x2 <- st_polygon(list(cbind(c(0,1,1,1,0,0),c(0,0,1,0.6,1,0))))
x3 <- st_polygon(list(cbind(c(0,1,0,1,0),c(0,1,1,0,0))))
st_is_simple(st_sfc(x1))
## [1] FALSE
st_is_valid(st_sfc(x2,x3))
## Warning in evalq((function (..., call. = TRUE, immediate. =
## FALSE, noBreaks. = FALSE, : Self-intersection at or near point 1
## 0.59999999999999998
## Warning in evalq((function (..., call. = TRUE, immediate. = FALSE,
## noBreaks. = FALSE, : Self-intersection at or near point 0.5 0.5
## [1] FALSE FALSE

Units

Where possible geometric operations such as st_distance(), st_length() and st_area() report results with a units attribute appropriate for the CRS:

a <- st_area(nc[1,])
attributes(a)
## $units
## $numerator
## [1] "m" "m"
## 
## $denominator
## character(0)
## 
## attr(,"class")
## [1] "symbolic_units"
## 
## $class
## [1] "units"

The units package can be used to convert between units:

units::set_units(a, km^2) # result in square kilometers
## 1137.389 km^2
units::set_units(a, ha) # result in hectares
## 113738.9 ha

The result can be stripped of their attributes if needs be:

as.numeric(a)
## [1] 1137388604

How attributes relate to geometries

(This will eventually be the topic of a new vignette; now here to explain the last attribute of sf objects)

The standard documents about simple features are very detailed about the geometric aspects of features, but say nearly nothing about attributes, except that their values should be understood in another reference system (their units of measurement, e.g. as implemented in the package units). But there is more to it. For variables like air temperature, interpolation usually makes sense, for others like human body temperature it doesn’t. The difference is that air temperature is a field, which continues between sensors, where body temperature is an object property that doesn’t extend beyond the body – in spatial statistics bodies would be called a point pattern, their temperature the point marks. For geometries that have a non-zero size (positive length or area), attribute values may refer to the every sub-geometry (every point), or may summarize the geometry. For example, a state’s population density summarizes the whole state, and is not a meaningful estimate of population density for a give point inside the state without the context of the state. On the other hand, land use or geological maps give polygons with constant land use or geology, every point inside the polygon is of that class. Some properties are spatially extensive, meaning that attributes would summed up when two geometries are merged: population is an example. Other properties are spatially intensive, and should be averaged, with population density the example.

Simple feature objects of class sf have an agr attribute that points to the attribute-geometry-relationship, how attributes relate to their geometry. It can be defined at creation time:

nc <- st_read(system.file("shape/nc.shp", package="sf"),
    agr = c(AREA = "aggregate", PERIMETER = "aggregate", CNTY_ = "identity",
        CNTY_ID = "identity", NAME = "identity", FIPS = "identity", FIPSNO = "identity",
        CRESS_ID = "identity", BIR74 = "aggregate", SID74 = "aggregate", NWBIR74 = "aggregate",
        BIR79 = "aggregate", SID79 = "aggregate", NWBIR79 = "aggregate"))
## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.4/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
## 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.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID):    4267
## proj4string:    +proj=longlat +datum=NAD27 +no_defs
st_agr(nc)
##      AREA PERIMETER     CNTY_   CNTY_ID      NAME      FIPS    FIPSNO 
## aggregate aggregate  identity  identity  identity  identity  identity 
##  CRESS_ID     BIR74     SID74   NWBIR74     BIR79     SID79   NWBIR79 
##  identity aggregate aggregate aggregate aggregate aggregate aggregate 
## Levels: constant aggregate identity
data(meuse, package = "sp")
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, agr = "constant")
st_agr(meuse_sf)
##  cadmium   copper     lead     zinc     elev     dist       om    ffreq 
## constant constant constant constant constant constant constant constant 
##     soil     lime  landuse   dist.m 
## constant constant constant constant 
## Levels: constant aggregate identity

When not specified, this field is filled with NA values, but if non-NA, it has one of three possibilities

value meaning
constant a variable that has a constant value at every location over a spatial extent; examples: soil type, climate zone, land use
aggregate values are summary values (aggregates) over the geometry, e.g. population density, dominant land use
identity values identify the geometry: they refer to (the whole of) this and only this geometry

With this information (still to be done) we can for instance

  • either return missing values or generate warnings when a aggregate value at a point location inside a polygon is retrieved, or
  • list the implicit assumptions made when retrieving attribute values at points inside a polygon when relation_to_geometry is missing.
  • decide what to do with attributes when a geometry is split: do nothing in case the attribute is constant, give an error or warning in case it is an aggregate, change the relation_to_geometry to constant in case it was identity.

Simple feature geometries manipulation

Simple features can be manipulated including:

  • 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.

Coordinate reference systems conversion and transformation

Getting and setting coordinate reference systems of sf objects

The coordinate reference system of objects of class sf or sfc is obtained by st_crs, and replaced by st_crs<-:

geom = st_sfc(st_point(c(0,1)), st_point(c(11,12)))
s = st_sf(a = 15:16, geometry = geom)
st_crs(s)
## $epsg
## [1] NA
## 
## $proj4string
## [1] NA
## 
## attr(,"class")
## [1] "crs"
s1 = s
st_crs(s1) <- 4326
st_crs(s1)
## $epsg
## [1] 4326
## 
## $proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
## 
## attr(,"class")
## [1] "crs"
s2 = s
st_crs(s2) <- "+proj=longlat +datum=WGS84"
all.equal(s1, s2)
## [1] TRUE

an alternative, more pipe-friendly version of st_crs<- is

s1 %>% st_set_crs(4326)
## Simple feature collection with 2 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 0 ymin: 1 xmax: 11 ymax: 12
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##    a      geometry
## 1 15   POINT (0 1)
## 2 16 POINT (11 12)

Coordinate reference system transformations

If we change the coordinate reference system from one non-missing value into another non-missing value, the crs is is changed without modifying any coordinates, but a warning is issued that this did not reproject values:

s3 <- s1 %>% st_set_crs(4326) %>% st_set_crs(3857)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that

A cleaner way to do this that better expresses intention and does not generate this warning is to first wipe the CRS by assigning it a missing value, and then setting it to the intended value.

s3 <- s1  %>% st_set_crs(NA) %>% st_set_crs(3857)

To carry out a coordinate conversion or transformation, we use st_transform

s3 <- s1 %>% st_transform(3857)
s3
## Simple feature collection with 2 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 0 ymin: 111325.1 xmax: 1224514 ymax: 1345708
## epsg (SRID):    3857
## proj4string:    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
##    a                       geometry
## 1 15     POINT (0 111325.142866385)
## 2 16 POINT (1224514.39872601 134...

for which we see that coordinates are actually modified (projected).

Geometrical operations

All geometrical operations st_op(x) or or st_op2(x,y) work both for sf objects as well as sfc objects x and y; since the operations work on the geometries, the non-geometries parts of an sf object are simply discarded. Also, all binary operations st_op2(x,y) called with a single argument, as st_op2(x), are handled as st_op2(x,x).

We will illustrate the geometrical operations on a very simple dataset:

b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
b1 = b0 + 2
b2 = b0 + c(-0.2, 2)
x = st_sfc(b0, b1, b2)
a0 = b0 * 0.8
a1 = a0 * 0.5 + c(2, 0.7)
a2 = a0 + 1
a3 = b0 * 0.5 + c(2, -0.5)
y = st_sfc(a0,a1,a2,a3)
plot(x, border = 'red')
plot(y, border = 'green', add = TRUE)

Unary operations

st_is_valid returns whether polygon geometries are topologically valid:

b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
b1 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(0,-1), c(-1,-1))))
st_is_valid(st_sfc(b0,b1))
## Warning in evalq((function (..., call. = TRUE, immediate. = FALSE,
## noBreaks. = FALSE, : Self-intersection at or near point -1 -1
## [1]  TRUE FALSE

and st_is_simple whether line geometries are simple:

s = st_sfc(st_linestring(rbind(c(0,0), c(1,1))), 
    st_linestring(rbind(c(0,0), c(1,1),c(0,1),c(1,0))))
st_is_simple(s)
## [1]  TRUE FALSE

st_area returns the area of polygon geometries, st_length the length of line geometries:

st_area(x)
## [1] 4 4 4
st_area(st_sfc(st_point(c(0,0))))
## [1] 0
st_length(st_sfc(st_linestring(rbind(c(0,0),c(1,1),c(1,2))), st_linestring(rbind(c(0,0),c(1,0)))))
## [1] 2.414214 1.000000
st_length(st_sfc(st_multilinestring(list(rbind(c(0,0),c(1,1),c(1,2))),rbind(c(0,0),c(1,0))))) # ignores 2nd part!
## [1] 2.414214

Binary operations: distance and relate

st_distance computes the shortest distance matrix between geometries; this is a dense matrix:

st_distance(x,y)
##           [,1] [,2] [,3]     [,4]
## [1,] 0.0000000  0.6    0 0.500000
## [2,] 0.2828427  0.0    0 1.000000
## [3,] 0.2000000  0.8    0 1.220656

st_relate returns a dense character matrix with the DE9-IM relationships between each pair of geometries:

st_relate(x,y)
##      [,1]        [,2]        [,3]        [,4]       
## [1,] "212FF1FF2" "FF2FF1212" "212101212" "FF2FF1212"
## [2,] "FF2FF1212" "212101212" "212101212" "FF2FF1212"
## [3,] "FF2FF1212" "FF2FF1212" "212101212" "FF2FF1212"

element [i,j] of this matrix has nine characters, refering to relationship between x[i] and y[j], encoded as IxIy, IxBy, IxEy, BxIy, BxBy, BxEy, ExIy, ExBy, ExEy where I refers to interior, B to boundary, and E to exterior, and e.g. BxIy the dimensionality of the intersection of the the boundary Bx of x[i] and the interior Iy of y[j], which is one of {0,1,2,F}, indicating zero-, one-, two-dimension intersection, and (F) no intersection, respectively.

DE9-IM Reading from left-to-right and top-to-bottom, the DE-9IM(a,b) string code is ‘212101212’, the compact representation of IxIy = 2, IxBy = 1, IxEy = 2, BxIy = 1, BxBy = 0, BxEy = 1, ExIy = 2, ExBy = 1, ExEy = 2. Figure from here.

Binary logical operations:

Binary logical operations return either a sparse matrix

st_intersects(x,y)
## [[1]]
## [1] 1 3
## 
## [[2]]
## [1] 2 3
## 
## [[3]]
## [1] 3

or a dense matrix

st_intersects(x, x, sparse = FALSE)
##      [,1]  [,2]  [,3]
## [1,] TRUE  TRUE  TRUE
## [2,] TRUE  TRUE FALSE
## [3,] TRUE FALSE  TRUE
st_intersects(x, y, sparse = FALSE)
##       [,1]  [,2] [,3]  [,4]
## [1,]  TRUE FALSE TRUE FALSE
## [2,] FALSE  TRUE TRUE FALSE
## [3,] FALSE FALSE TRUE FALSE

where list element i of a sparse matrix contains the indices of the TRUE elements in row i of the the dense matrix. For large geometry sets, dense matrices take up a lot of memory and are mostly filled with FALSE values, hence the default is to return a sparse matrix.

st_intersects returns for every geometry pair whether they intersect (dense matrix), or which elements intersect (sparse).

Other binary predicates include (using sparse for readability):

st_disjoint(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3] [,4]
## [1,] FALSE  TRUE FALSE TRUE
## [2,]  TRUE FALSE FALSE TRUE
## [3,]  TRUE  TRUE FALSE TRUE
st_touches(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_crosses(s, s, sparse = FALSE)
##       [,1]  [,2]
## [1,] FALSE FALSE
## [2,] FALSE FALSE
st_within(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_contains(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_overlaps(x, y, sparse = FALSE)
##       [,1]  [,2] [,3]  [,4]
## [1,] FALSE FALSE TRUE FALSE
## [2,] FALSE  TRUE TRUE FALSE
## [3,] FALSE FALSE TRUE FALSE
st_equals(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_covers(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_covered_by(x, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE
st_covered_by(y, y, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,]  TRUE FALSE FALSE FALSE
## [2,] FALSE  TRUE FALSE FALSE
## [3,] FALSE FALSE  TRUE FALSE
## [4,] FALSE FALSE FALSE  TRUE
st_equals_exact(x, y,0.001, sparse = FALSE)
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE

Operations returning a geometry

u = st_union(x)
plot(u)

par(mfrow=c(1,2), mar = rep(0,4))
plot(st_buffer(u, 0.2))
plot(u, border = 'red', add = TRUE)
plot(st_buffer(u, 0.2), border = 'grey')
plot(u, border = 'red', add = TRUE)
plot(st_buffer(u, -0.2), add = TRUE)

plot(st_boundary(x))

par(mfrow = c(1:2))
plot(st_convex_hull(x))
plot(st_convex_hull(u))

par(mfrow = c(1,1))
par(mfrow=c(1,2))
plot(x)
plot(st_centroid(x), add = TRUE, col = 'red')
plot(x)
plot(st_centroid(u), add = TRUE, col = 'red')

The intersection of two geometries is the geometry covered by both; it is obtained by st_intersection:

plot(x)
plot(y, add = TRUE)
plot(st_intersection(st_union(x),st_union(y)), add = TRUE, col = 'red')

To get everything but the intersection, use st_difference or st_sym_difference`:

par(mfrow=c(2,2), mar = c(0,0,1,0))
plot(x, col = '#ff333388'); 
plot(y, add=TRUE, col='#33ff3388')
title("x: red, y: green")
plot(x, border = 'grey')
plot(st_difference(st_union(x),st_union(y)), col = 'lightblue', add = TRUE)
title("difference(x,y)")
plot(x, border = 'grey')
plot(st_difference(st_union(y),st_union(x)), col = 'lightblue', add = TRUE)
title("difference(y,x)")
plot(x, border = 'grey')
plot(st_sym_difference(st_union(y),st_union(x)), col = 'lightblue', add = TRUE)
title("sym_difference(x,y)")

Function st_segmentize adds points to straight line sections of a lines or polygon object:

par(mfrow=c(1,3),mar=c(1,1,0,0))
pts = rbind(c(0,0),c(1,0),c(2,1),c(3,1))
ls = st_linestring(pts)
plot(ls)
points(pts)
ls.seg = st_segmentize(ls, 0.3)
plot(ls.seg)
pts = ls.seg
points(pts)
pol = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0))))
pol.seg = st_segmentize(pol, 0.3)
plot(pol.seg, col = 'grey')
points(pol.seg[[1]])

Function st_polygonize polygonizes a multilinestring, as far as the points form a closed polygon:

par(mfrow=c(1,2),mar=c(0,0,1,0))
mls = st_multilinestring(list(matrix(c(0,0,0,1,1,1,0,0),,2,byrow=TRUE)))
x = st_polygonize(mls)
plot(mls, col = 'grey')
title("multilinestring")
plot(x, col = 'grey')
title("polygon")

Further reading:

  1. S. Scheider, B. Gräler, E. Pebesma, C. Stasch, 2016. Modelling spatio-temporal information generation. Int J of Geographic Information Science, 30 (10), 1980-2008. (pdf)
  2. Stasch, C., S. Scheider, E. Pebesma, W. Kuhn, 2014. Meaningful Spatial Prediction and Aggregation. Environmental Modelling & Software, 51, (149–165, open access).

  1. R Bivand (2011) Introduction to representing spatial objects in R

  2. Coordinates should be of type double and will be promoted if not.

  3. E. Pebesma & R. Bivand (2016)Spatial data in R: simple features and future perspectives



logo
Creative Commons License