R for GIS
BayGeo, Spring
2024
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_bgs24\exercises\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" "UNITNAME" "GROUPCODE" "REGIONCODE" "ISEXTANT" "MAPMETHOD"
## [13] "MAPSOURCE" "SRCESCALE" "SOURCEDATE" "XYERROR" "LOCATIONID" "ASSETID" "NOTES" "DISTRIBUTE" "RESTRICTIO" "CREATEDATE" "CREATEUSER" "EDITDATE"
## [25] "EDITUSER" "FEATUREID" "GEOMETRYID" "PLACESID" "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.
Now that the trails are in geographic coordinates, we can overlay them on the boundary and historical points.
yose_hp_ll <- st_read("./data/yose_historic_pts.kml", layer="yose_historic_places", quiet=TRUE)
## 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_ll |> 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_ll |> st_geometry(), col="gray30", pch=16, add=TRUE)
plot(yose_bnd_ll |> st_geometry(), col=NA, border="chartreuse4", lwd=3, add=TRUE)## [1] TRUE
To turn a data frame into a sf object, use
st_as_sf(df, coords, crs)
where:
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).
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:
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
## Checking rgeos availability: TRUEdata(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"
See also: