The sf package encodes spatial vector data according to the Simple Features standard described by the Open Geospatial Consortium (OGC) and International Organization for Standardization (ISO). The sf package provides bindings to GDAL for reading and writing data, to GEOS for geometrical operations, and to proj.4 for projection conversions and datum transformations.

library(sf)

The sfpackage represents simple features as native R objects, using simple data structures (S3 classes, lists, matrices, vectors). All functions and methods in sf that operate on spatial data are prefixed by st_, which refers to spatial and temporal.

In the sf package three classes are used to represent simple features:

Note that the sf will gradually replace the well known sp package for spatial analysis. For more details on the sf package refer to the documentation on CRAN, and visit the vignettes (here, here, and here).


sf: objects with simple features

A sf object is a data set consisting of sets of features with attributes.

To showcase the functionality of the sf package we download a polygon data set, an ESRI shapefile, and load it into memory with the st_read() function.

The shapefile stores the geometry of districts in Germany and additional information on the population size. The shapefile was downloaded from the Federal Agency for Cartography and Geodesy on June 25, 2017. The data is found here, and the file description is found here (GeoBasis-DE/BKG 2017). Additional information with respect to the population size was obtained from the The Federal Statistical Office on June 26, 2017 (here).

## build a temporary folder on disk
temp <- tempfile()
download.url <- "https://userpage.fu-berlin.de/soga/300/30100_data_sets/spatial/"
zipfile <- "germany_districts.zip"
## download the file
download.file(paste0(download.url,zipfile),temp, mode="wb")
## unzip the file(s)
unzip(temp)
## close file connection
unlink(temp)
## read the data
dists <- st_read("./germany_districts.shp")
## Reading layer `germany_districts' from data source `/Users/jokr/Documents/soga/germany_districts.shp' using driver `ESRI Shapefile'
## Simple feature collection with 401 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 280371.1 ymin: 5235856 xmax: 921292.4 ymax: 6101444
## epsg (SRID):    NA
## proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs

The short report printed gives the file name, the driver (ESRI Shapefile), shows that there are 401 features (records, represented as rows) and 6 fields (attributes, represented as columns) in the data set.

The class() function shows that the sf object class extends the data.frame object class.

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

The column with name geometry holds the geometries, and the other columns refer to the additional information available for each records.

colnames(dists)
## [1] "district" "area"     "pop"      "pop_m"    "pop_w"    "pop_dens"
## [7] "geometry"

We can easily access the first 20 rows of the column that contains the names of the districts of Germany.

head(dists$district, 20)
##  [1] Ahrweiler                 Aichach-Friedberg        
##  [3] Alb-Donau-Kreis           Altenburger Land         
##  [5] Altenkirchen (Westerwald) Altmarkkreis Salzwedel   
##  [7] Altötting                 Alzey-Worms              
##  [9] Amberg                    Amberg-Sulzbach          
## [11] Ammerland                 Anhalt-Bitterfeld        
## [13] Ansbach                   Ansbach                  
## [15] Aschaffenburg             Aschaffenburg            
## [17] Augsburg                  Augsburg                 
## [19] Aurich                    Bad Dürkheim             
## 379 Levels: Ahrweiler Aichach-Friedberg ... Zwickau

The sf object class offers a vast range of functions and spatial data operations to be performed.

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

Let us try the plot() function on the sf object dists:

plot(dists)