Author

Carolyn Koehn

Load libraries

Code 
library(sf)
library(tidyverse)

Centroids vs Point on Surface

Code 
id.counties <- tigris::counties(state = "ID", progress_bar=FALSE)
id.centroid <- st_centroid(id.counties)
id.pointonsurf <- st_point_on_surface(id.counties)
Code 
plot(st_geometry(id.counties))
plot(st_geometry(id.centroid), col="blue", add=TRUE)
plot(st_geometry(id.pointonsurf), col="red", add=TRUE)

Practice Example 1: Distance on Points and Polygons

Code 
system.time(poly_dist <- st_distance(id.counties))
   user  system elapsed 
   4.32    0.03    4.36
Code 
system.time(cent_dist <- st_distance(id.centroid))
   user  system elapsed 
   0.00    0.00    0.01
Code 
system.time(pos_dist <- st_distance(id.pointonsurf))
   user  system elapsed 
   0.02    0.00    0.00

Practice Example 2: Intersections and Buffers

Code 
# get roads data
roads <- tigris::primary_secondary_roads("ID", progress_bar=FALSE)

# get a polygon of Ada county
ada.cty <- filter(id.counties, NAME == "Ada")

# find all road sections within Ada county
ada.roads <- st_intersection(roads, ada.cty)
Code 
# plot result
plot(st_geometry(ada.cty))
plot(st_geometry(ada.roads), col="purple", add=TRUE)

Code 
# check the units of the CRS
st_crs(id.centroid)
Coordinate Reference System:
  User input: NAD83 
  wkt:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]
Code 
# create a 50km buffer around the centroid of Ada county
ada.cent <- filter(id.centroid, NAME == "Ada")
ada.buff <- st_buffer(ada.cent, dist = 50000)
Code 
# get roads within buffer zone
roads.buff <- st_intersection(roads, ada.buff)
Code 
# plot result
plot(st_geometry(roads.buff))

If we want to plot county boundaries, st_intersection won’t work. We’ll only get parts of each county polygon:

Code 
cty.buff <- st_intersection(id.counties, ada.buff)
plot(st_geometry(cty.buff))

Back to predicates!

Code 
# find counties that intersect with buffer
cty.50 <- st_intersects(id.counties, ada.buff, sparse = FALSE)
# convert to vector so filter is happy
cty.50 <- as.vector(cty.50)

# get counties that intersect buffer
map.counties <- filter(id.counties, cty.50)
Code 
# plot result
plot(map.counties$geometry)
plot(st_geometry(roads.buff), col="purple", add=TRUE)