Spatial Data Analysis with R
Society for Conservation GIS, July 2021
Common Manipulations for sf Objects 
Common Manipulations for sf Objects
simple features is an ISO standard for storing and accessing spatial data. It is widely adopted in spatial databases, open source formats like GeoJSON, GIS software, tools, etc.
sf employs a standard called “Well Known Text” (WKT) to encode the geometry. WKT is a very simple written representation of the ‘connect-the-dots’ vector data model.
| Geometry | Sample WKT |
|---|---|
| point | POINT (2 4) |
| multipoint | MULTIPOINT (2 2, 3 3, 3 2) |
| linestring | LINESTRING (0 3, 1 4, 2 3) |
| polygon | POLYGON ((1 0, 3 4, 5 1, 1 0)) |
Import the Yosemite points-of-interest Shapefile.
library(sf)
## Import the Yosemite Points-of-Interest Shapefile
yose_poi_utm <- st_read(dsn="./data", layer="yose_poi")
## View column names in the attribute table
names(yose_poi_utm)
## Plot the points
plot(yose_poi_utm %>% st_geometry(), pch=16, asp=1, cex=0.5, axes=TRUE)## Reading layer `yose_poi' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data' using driver `ESRI Shapefile'
## Simple feature collection with 2720 features and 30 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 246416.2 ymin: 4153717 xmax: 301510.7 ymax: 4208419
## Projected CRS: NAD83 / UTM zone 11N
## [1] "OBJECTID" "POINAME" "POIALTNAME" "POILABEL" "POIFEATTYP" "POITYPE" "UNITCODE"
## [8] "UNITNAME" "GROUPCODE" "REGIONCODE" "ISEXTANT" "MAPMETHOD" "MAPSOURCE" "SRCESCALE"
## [15] "SOURCEDATE" "XYERROR" "LOCATIONID" "ASSETID" "NOTES" "DISTRIBUTE" "RESTRICTIO"
## [22] "CREATEDATE" "CREATEUSER" "EDITDATE" "EDITUSER" "FEATUREID" "GEOMETRYID" "PLACESID"
## [29] "GlobalID" "TAGS" "geometry"
Take a closer look at the geometry column:
## [1] "sfc_POINT" "sfc"
## Simple feature collection with 6 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 260859.4 ymin: 4153717 xmax: 270931.1 ymax: 4179771
## Projected CRS: NAD83 / UTM zone 11N
## geometry
## 1 POINT (260859.4 4178493)
## 2 POINT (264115.4 4158415)
## 3 POINT (268276.8 4153717)
## 4 POINT (268620.1 4178344)
## 5 POINT (269052.6 4178870)
## 6 POINT (270931.1 4179771)The geometry column of a sf object does not have to be called ‘geometry’. You can view the column name by running attr(x, "sf_column") where x is a sf object.
The features in a sf object can be a mix of points, lines, and polyons. The data class of the geometry column is therefore sfc (simple features collection).
You can project data from one CRS to another with:
st_transform(sf_object, new_crs)
Let’s ‘unproject’ the Yosemite trails layer from UTM to geographic coordinates.
## (Un)project the trails layer to geographic coordinates
yose_trails_ll <- st_transform(yose_trails, 4269)Now that the trails are in geographic coordinates, we can overlay them on the boundary and historical points.
## Plot the trails, then overlay the historic points and park boundary
plot(yose_trails_ll %>% st_geometry(), asp=1, col="bisque3", axes=T, main="Yosemite National Park")
plot(yose_hp %>% st_geometry(), col="gray30", pch=16, add=TRUE)
plot(yose_bnd_ll %>% st_geometry(), col=NA, border="chartreuse4", lwd=3, add=TRUE)
You can project one layer to match another without knowing the exact CRS. Just extract the CRS of the layer you’re trying to match with st_crs(), and use that as the second argument in st_transform().
Import the Yosemite roads and add those to your map. You may use the code below to get started.
[Solution]
## (Un)project the roads layer to geographic coordinates
yose_roads_ll <- st_transform(yose_roads_utm, 4269)
## Plot all four layers
plot(yose_trails_ll %>% st_geometry(), asp=1, col="bisque3", axes=T, main="Yosemite")
plot(yose_roads_ll %>% st_geometry(), col="navyblue", lwd=1, lty=5, add=TRUE)
plot(yose_hp %>% st_geometry(), col="gray30", pch=16, add=TRUE)
plot(yose_bnd_ll %>% st_geometry(), col=NA, border="chartreuse4", lwd=3, add=TRUE)## [1] TRUE