`sf` Package

sf is a relatively new package for working with vector data. It combines the functionality of the older sp package, as well as select functionality from rgdal and rgeos.

library(sf)

sf stands for ‘simple features’, which is a standard (developed by the Open Geospatial Consortium) for storing various geometry types in a hierarchical data model.

A ‘feature’ is just a representation of a thing in the real world (e.g. a building, a city…).

sf offers numerous specific benefits, including:

uses established open source libraries: GDAL, GEOS, PROJ
fast IO (Input and Output)
enhanced plotting
excellent integration with common R data structures (using data frames)
integration with the tidyverse packages and %>% piping syntax
consistent function names (all starting with st_)
objects are supported by an increasing number of other geospatial packages
alignment with other software that uses the simple features data model (e.g. QGIS, PostGIS), or an analogous one (e.g. Python’s Geopandas)

Importing Vector Data

To import spatial data, you have to specify a source and layer.

The Source can be a:

folder: "./gis_data"

geodatabase (which is really a folder): "./gis_data/yose_roads.gdb"

file: "trail_heads.kml", "cell_towers.geojson"

database connection string: "PG:dbname=postgis"

In sf functions, the argument where you provide the source is often named dsn (data source name)

The Layer can be a:

Shapefile name (minus the shp extension)

a named layer in the database / file

Import Functions

The two main functions for importing vector data are:

st_layers(source) - returns the names of available layers in a source

st_read(source, layer) - import into R

Most of the functions in the sf package start with st_, which stands for ‘space time’, and matches the names of similar functions in PostGIS.

To view the metadata of a layer before bringing it into R, use rgdal::ogrInfo()

Spatial Data File Formats

Shapefiles

Old format dating back to the early 90s
Originally from ESRI
Well-supported by nearly all GIS software
Actually requires 4-9 files with extensions *.shp, *.shx, *.prj, *.dbf, etc.
Column names in attribute table limited to 10 characters
Better options available today
st_read():

source - shp filename
layer - can be omitted

Example:

yose_bnd_ll <- st_read(dsn="./data", layer="yose_boundary")
yose_bnd_ll

## Reading layer `yose_boundary' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -119.8864 ymin: 37.4947 xmax: -119.1964 ymax: 38.18515
## Geodetic CRS:  North_American_Datum_1983
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -119.8864 ymin: 37.4947 xmax: -119.1964 ymax: 38.18515
## Geodetic CRS:  North_American_Datum_1983
##   UNIT_CODE                                                                           GIS_NOTES              UNIT_NAME  DATE_EDIT STATE REGION GNIS_ID
## 1      YOSE Lands - http://landsnet.nps.gov/tractsnet/documents/YOSE/METADATA/yose_metadata.xml Yosemite National Park 2016-01-27    CA     PW  255923
##       UNIT_TYPE CREATED_BY                                               METADATA PARKNAME                       geometry
## 1 National Park      Lands http://nrdata.nps.gov/programs/Lands/YOSE_METADATA.xml Yosemite POLYGON ((-119.8456 37.8327...

KML & KMZ

made popular with Google Earth
data and symbology are stored in one file
kmz files are simply zipped kml files
kml & kmz files can have more than one layer
limitations:

kmls support geographic coordinates only
you can not import a kmz file with st_read() (at least not on Windows). Unzip it first then import the doc.kml file
st_read() can not currently import the attribute table from a kml file on Windows (more info)
☞ Stay away from kml/kmz files if possible!

st_read():

source - kml filename
layer - name of a layer (required)

GeoJSON

GeoJSON is a good alternative to the Shapefile and KML formats
Supports all major geometry types, projections, attribute fields, etc.
GeoJSON typically stores one vector layer in one file that ends in .geojson (thus not very compact)
Format is plain text; can open and read in a text editor
st_read():

source - geojson filename
layer - named layer (optional if only 1 layer exists)

File Geodatabase

Proprietary spatial database format from ESRI
On disk is a folder with a .gdb extension
Supports all major feature types, projections, attribute tables, etc.
You can import (but not write to) a file geodatabase using the sf package
You can write to a geodatabase with the R-ArcGIS Bridge, if you have licensed copy of ArcGIS installed.
st_read():

source - folder containing the geodatabase
layer - named layer (required)

GeoPackage

GeoPackage is a relatively new open format for geospatial data. It is similar to a file geodatabase in that:

it can contain many layers
supports both raster and vector data
under the hood data are stored in a spatial database

The GeoPackage is different from a file geodatabase in that GeoPackages are:

a single file disk with a .gpkg extension (as opposed to a directory)
completely open source following a OGC standard
based on SQLite instead of a proprietary database format
cross-platform

st_read():

source - .gpkg filename
layer - named layer (required)

Spatial Databases

PostGIS - see rpostgis

SpatiaLite - see RSQLite.spatialite

The workshop datasets descriptions includes sample import code for most of the layers.

Basic Plotting sf Objects

To plot just the geometry of a sf object (i.e., no symbology from the attribute table), use st_geometry():

yose_bnd_ll <- st_read(dsn="./data", layer="yose_boundary", quiet=TRUE)

## Plot the geometry (outline) of the Yosemite boundary
plot(yose_bnd_ll %>% st_geometry(), asp=1)

We add %>% st_geometry() to indicate we want to plot the polygon with a single color/fill (more on that later).

The asp=1 argument sets the aspect ratio. It is optional but generally a good idea for spatial data.

Overlaying Layers

To overlay more than one layer on a plot:

both layers must have the same CRS (plot() does not reproject on the fly)
include add=TRUE to the second (and all subsequent) plot() statements

## Import the 'yosem_historic_places' layer 
yose_hp_ll <- st_read("./data/yose_historic_pts.kml", layer="yose_historic_places", quiet=TRUE)

## Plot the boundary
plot(yose_bnd_ll %>% st_geometry(), asp=1)

## Add historic places
plot(yose_hp_ll %>% st_geometry(), add=TRUE)

Convert Plain Data Frames to sf Data Frames

To turn a data frame into a sf object, use

st_as_sf(df, coords, crs)

where:

df - the data frame
coords - a character vector containing the names of the columns containing x and y coordinates
crs - the coordinate reference system (an EPSG code or proj4 string)

Example: Turn the quakes Data Frame into a sf Object

quakes is a dataframe that comes with R (in the datasets package). Here we turn it into a sf data frame.

## Make a sf object from the built-in 'quakes' dataset
head(quakes)
quakes_sf <- st_as_sf(quakes, coords=c("long", "lat"), crs=4326)
plot(quakes_sf %>% st_geometry(), axes = TRUE, asp = 1)

##      lat   long depth mag stations
## 1 -20.42 181.62   562 4.8       41
## 2 -20.62 181.03   650 4.2       15
## 3 -26.00 184.10    42 5.4       43
## 4 -17.97 181.66   626 4.1       19
## 5 -20.42 181.96   649 4.0       11
## 6 -19.68 184.31   195 4.0       12

If you don’t know the datum for geographic coordinates, use WGS84 (epsg 4326).

R Notebook Exercise

nb_import-vector_scgis21.Rmd
Import all kinds of vector data into R

preview notebook | answer key

sp → sf

sf is relatively new, and many R packages still use the older data classes from the sp package.

Fortunately it’s relatively easy to convert back and forth:

sp → sf

st_as_sf()

sf → sp

as(x, “Spatial”)

Example: sp → sf

## Get a sample sp object from the maptools package
library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

data(wrld_simpl)
class(wrld_simpl)

## Convert sp to sf
wrld_simpl_sf <- sf::st_as_sf(wrld_simpl)
class(wrld_simpl_sf)

## Plot the sf object
plot(wrld_simpl_sf["REGION"])

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
## [1] "sf"         "data.frame"

Summary

In this section we saw how to:

import from Shapefiles, KMLs, GeoJson files, geodatabases, and GeoPackages
- st_layers() views the layer names
- st_read() imports the layer
use plot() for basic plotting
use plot(..., add=TRUE) to overlay additional layers with the same CRS
convert a plain data frame to a sf data frame with st_as_Sf()

Additional Resources

Importing spatial data with rgdal, by Zev Ross

Plotting Simple Features, sf package Vignette

Next: Visualization with tmap

Spatial Data Analysis with RSociety for Conservation GIS, July 2021

sf Package

Importing Vector Data

Import Functions

Spatial Data File Formats

Shapefiles

Example:

KML & KMZ

GeoJSON

File Geodatabase

GeoPackage

Spatial Databases

Basic Plotting sf Objects

Overlaying Layers

Convert Plain Data Frames to sf Data Frames

Example: Turn the quakes Data Frame into a sf Object

R Notebook Exercise

sp → sf

Example: sp → sf

Summary

Spatial Data Analysis with R
Society for Conservation GIS, July 2021

`sf` Package