Spatial Data II

Computing on geometries with sf

Last Time

1. A coordinate system locates points from an origin.
2. A datum models the shape of the earth.
3. A projection flattens the globe and distorts it.
4. A simple feature pairs a geometry with attributes.

Key Ideas for Today

1. An sf object is a data frame plus geometries.
2. Unary functions measure or reshape a geometry.
3. Binary predicates test how two geometries relate.
4. A spatial join matches rows by location.

Concept Acquisition

1. Spatial data come from R packages and web APIs, shapefiles, GeoJSON, and KML.
2. An sf object is a data frame of attributes plus a geometry column (an sfc) built from individual st_ geometries and a CRS.
3. dplyr verbs work on the attributes of an sf object and carry the geometry along; summarize() also combines the geometries.
4. Unary measures (dimension, area, length) return a number for each geometry; unary transformers (centroid, jitter, buffer) return a new geometry.
5. A binary predicate such as within tests how two geometries relate; a spatial join uses one to match rows, keeping (left) or dropping (inner) unmatched rows.

Tool Acquisition

1. st_read(), st_geometry(), and base plot() / points().
2. st_point(), st_polygon(), st_sfc(), and st_sf().
3. st_dimension(), st_area(), st_length(), st_centroid(), st_jitter(), and st_buffer().
4. st_within() and st_join() with join = and left =.

Concept Application

1. Predict the rows and columns of a left join with st_within (timed, 1:00).
2. PS 5: map the centroid of California (and, as a challenge, find its county with a spatial join); build lines and a polygon from scratch, compute their measures, split the polygon with the lines, and join the pieces to a triangle.

(Last Time) Simple Feature Object:

Features: “things” or objects that have a spatial location or extent; they may be physical objects like a building, or social conventions like a political state

Feature geometry: the spatial properties (location or extent) of a feature, and can be described by a point, a point set, a polygon, etc.

Feature attributes: other non-spatial properties measured on the feature (e.g. color, name, landtype)

Agenda

1. Reading in Spatial Data
2. Building simple features from Scratch
3. Computing on Attributes
4. Computing on Geometries
5. Spatial Joins

Reading in Spatial Data

Sources of sf data

1. R packages
2. web APIs
3. Shapefiles

We saw R packages last time. web APIs is coming soon.

Shapefile

A bundle of files that describe spatial data. Developed by Esri for its ArcView GIS; the format’s specification is public, so most GIS software can read it.

Google: CA shapefile

CA Shapefile

library(sf)
ca <- st_read("data/CA_State.shp")

Reading layer `CA_State' from data source 
  `/workspace/24-spatial-2/data/CA_State.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -13857270 ymin: 3832931 xmax: -12705030 ymax: 5162406
Projected CRS: Popular Visualisation CRS / Mercator

CA Shapefile

ca

Simple feature collection with 1 feature and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -13857270 ymin: 3832931 xmax: -12705030 ymax: 5162406
Projected CRS: Popular Visualisation CRS / Mercator
  OBJECTID REGION DIVISION STATEFP  STATENS GEOID STUSPS       NAME LSAD MTFCC
1       56      4        9      06 01779778    06     CA California   00 G4000
  FUNCSTAT        ALAND      AWATER    INTPTLAT     INTPTLON Shape_Leng
1        A 403467168329 20499799756 +37.1551773 -119.5434183    5258041
  Shape_Le_1   Shape_Area                       geometry
1    5258041 671892661618 MULTIPOLYGON (((-13060109 3...

Annotate the header: geometry type, bounding box, CRS. Then box the attribute columns vs. the geometry column.

CA Shapefile

plot(st_geometry(ca))

This is a good teaching moment of a “getter” like we saw with lubridate (st_geometry()) as well as a plot method for the geometry class.

Sources of sf data

1. R packages
2. web APIs
3. Shapefiles

4. GeoJSON: Web-friendly, single file
5. KML: Web-friendly, designed for Google Earth/Maps

Building sf from scratch

Building a point set

st_: “spatial type”

sfc: “simple feature (geometry) column”

p1 <- st_point(c(7.35, 52.31))
p2 <- st_point(c(7.22, 52.18))
p3 <- st_point(c(7.44, 52.19))
point_sfc <- st_sfc(list(p1, p2, p3), crs = 'OGC:CRS84')
point_sfc

Geometry set for 3 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)

Check typeof and class as you run each line.

On the pinned board: start the sf table with a “build” column (st_point, st_sfc, st_sf). Add to it at each section.

Building a point set

my_point_sf <- st_sf(elev = c(33.2, 52.1, 81.2), 
                     marker = c("Id01", "Id02", "Id03"),
                     geom = point_sfc)
my_point_sf

Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev marker               geom
1 33.2   Id01 POINT (7.35 52.31)
2 52.1   Id02 POINT (7.22 52.18)
3 81.2   Id03 POINT (7.44 52.19)

Draw:

Quick plotting

plot(st_geometry(my_point_sf), pch = 16)

Building a polygon

# record coordinates
square <- matrix(c(7.17, 52.11,   # lower left  
                   7.40, 52.11,   # lower right  
                   7.40, 52.35,   # upper right  
                   7.17, 52.35,   # upper left  
                   7.17, 52.11),  # back to start
                   ncol = 2, byrow = TRUE)

# make polygon and simple feature
square <- st_polygon(list(square))
square_sfc <- st_sfc(square, crs = "OGC:CRS84")
my_square_sf <- st_sf(elev = 45, 
                      marker = "Id04",
                      geom = square_sfc)

Number each row of the matrix and sketch the square with matching corner numbers; circle the repeated first point.

Building a polygon

my_square_sf

Simple feature collection with 1 feature and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 7.17 ymin: 52.11 xmax: 7.4 ymax: 52.35
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                           geom
1   45   Id04 POLYGON ((7.17 52.11, 7.4 5...

Building a polygon

plot(st_geometry(my_square_sf), lwd = 5)

Building a polygon

plot(st_geometry(my_square_sf), lwd = 5)
points(st_geometry(my_point_sf), pch = 16)

Successive base R plot commands draw on top of one another.

Computing on Attributes

Computing on Attributes

library(dplyr)
my_point_sf |>
  slice(1:3)

Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev marker               geom
1 33.2   Id01 POINT (7.35 52.31)
2 52.1   Id02 POINT (7.22 52.18)
3 81.2   Id03 POINT (7.44 52.19)

Computing on Attributes

my_point_sf |>
  arrange(desc(elev))

Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev marker               geom
1 81.2   Id03 POINT (7.44 52.19)
2 52.1   Id02 POINT (7.22 52.18)
3 33.2   Id01 POINT (7.35 52.31)

Computing on Attributes

my_point_sf |>
  summarize(sum(elev))

Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  sum(elev)                           geom
1     166.5 MULTIPOINT ((7.35 52.31), (...

The image is a lead-in to “Computing on Geometries” and shows the different types that exist.

Computing on Geometries

I’ll be showing just a subset of the functions that can be used for computations on geometries. The important part is not to memorize them exhaustively, but think in a structured manner about the types of functions that should be runnable on different types of geometries.

Unary Measures

1. Dimension

my_point_sf |>
  mutate(dim = st_dimension(geom))

Simple feature collection with 3 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev marker               geom dim
1 33.2   Id01 POINT (7.35 52.31)   0
2 52.1   Id02 POINT (7.22 52.18)   0
3 81.2   Id03 POINT (7.44 52.19)   0

Unary Measures

1. Dimension

my_square_sf |>
  mutate(dim = st_dimension(geom))

Simple feature collection with 1 feature and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 7.17 ymin: 52.11 xmax: 7.4 ymax: 52.35
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                           geom dim
1   45   Id04 POLYGON ((7.17 52.11, 7.4 5...   2

Unary Measures

1. Dimension
2. Area

my_square_sf |>
  mutate(area = st_area(geom))

Simple feature collection with 1 feature and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 7.17 ymin: 52.11 xmax: 7.4 ymax: 52.35
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                           geom            area
1   45   Id04 POLYGON ((7.17 52.11, 7.4 5... 418033285 [m^2]

Unary Measures

1. Dimension
2. Area
3. Length

my_square_sf |>
  mutate(len = st_length(geom))

Simple feature collection with 1 feature and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 7.17 ymin: 52.11 xmax: 7.4 ymax: 52.35
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                           geom   len
1   45   Id04 POLYGON ((7.17 52.11, 7.4 5... 0 [m]

The objective here is to show that not every measure is computable for every geometry. Can show that st_perimeter() is the appropriate measure for a polygon.

Unary Transformers

1. Centroid

my_centroid <- my_square_sf |>
  st_centroid()
my_centroid

Simple feature collection with 1 feature and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.285 ymin: 52.22995 xmax: 7.285 ymax: 52.22995
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                   geom
1   45   Id04 POINT (7.285 52.22995)

Unary Transformers

1. Centroid

my_centroid <- my_square_sf |>
  st_centroid()
plot(st_geometry(my_square_sf), lwd = 5)

Unary Transformers

1. Centroid

my_centroid <- my_square_sf |>
  st_centroid()
plot(st_geometry(my_square_sf), lwd = 5)
points(st_geometry(my_centroid), pch = 16)

Unary Transformers

1. Centroid
2. Jitter

my_jittered <- my_point_sf |>
  st_jitter(factor = .1)
my_jittered

Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.218624 ymin: 52.17505 xmax: 7.432613 ymax: 52.33405
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                      geom
1 33.2   Id01 POINT (7.359401 52.33405)
2 52.1   Id02 POINT (7.218624 52.17505)
3 81.2   Id03 POINT (7.432613 52.21221)

Unary Transformers

1. Centroid
2. Jitter

my_jittered <- my_point_sf |>
  st_jitter(factor = .1)
plot(st_geometry(my_jittered), pch = 16)

Unary Transformers

1. Centroid
2. Jitter

my_jittered <- my_point_sf |>
  st_jitter(factor = .1)
plot(st_geometry(my_jittered), pch = 16)
points(st_geometry(my_point_sf), pch = 13)

Unary Transformers

1. Centroid
2. Jitter
3. Buffer

my_buff <- my_square_sf |>
  st_buffer(dist = 2000)
my_buff

Simple feature collection with 1 feature and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 7.137086 ymin: 52.0894 xmax: 7.433105 ymax: 52.37073
Geodetic CRS:  WGS 84 (CRS84)
  elev marker                           geom
1   45   Id04 POLYGON ((7.429425 52.1882,...

Unary Transformers

1. Centroid
2. Jitter
3. Buffer

my_buff <- my_square_sf |>
  st_buffer(dist = 2000)
plot(st_geometry(my_buff), lwd = 5)

Unary Transformers

1. Centroid
2. Jitter
3. Buffer

my_buff <- my_square_sf |>
  st_buffer(dist = 2000)
plot(st_geometry(my_buff), lwd = 5)
plot(st_geometry(my_square_sf), lwd = 5, add = TRUE)

Unary Transformers

1. Centroid
2. Jitter
3. Buffer
4. and Many More

Binary Predicates

1. Within

wi <- my_point_sf |>
  st_within(my_square_sf)
wi

Sparse geometry binary predicate list of length 3, where the predicate
was `within'
 1: 1
 2: 1
 3: (empty)

Binary Predicates

1. Within
2. and Many More

wi <- my_point_sf |>
  st_within(my_square_sf)
wi

Sparse geometry binary predicate list of length 3, where the predicate
was `within'
 1: 1
 2: 1
 3: (empty)

Spatial Joins

Joins

Table joins:

• Join records of two tables when keys match

Joins

Table joins:

• Join records of two tables when keys match

Spatial joins:

• Join records of two tables based on a spatial predicate.

my_point_sf

Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev marker               geom
1 33.2   Id01 POINT (7.35 52.31)
2 52.1   Id02 POINT (7.22 52.18)
3 81.2   Id03 POINT (7.44 52.19)

Your Turn

Question: How many rows and columns will st_join(my_point_sf, my_square_sf, join = st_within) return? What goes in the square’s columns for Id03?

01:00

Debrief by asking for dimensions, then reveal on the next slide.

Left join within

my_point_sf |>
  st_join(my_square_sf, join = st_within)

Simple feature collection with 3 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev.x marker.x elev.y marker.y               geom
1   33.2     Id01     45     Id04 POINT (7.35 52.31)
2   52.1     Id02     45     Id04 POINT (7.22 52.18)
3   81.2     Id03     NA     <NA> POINT (7.44 52.19)

Inner join within

my_point_sf |>
  st_join(my_square_sf, join = st_within, left = FALSE)

Simple feature collection with 2 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.35 ymax: 52.31
Geodetic CRS:  WGS 84 (CRS84)
  elev.x marker.x elev.y marker.y               geom
1   33.2     Id01     45     Id04 POINT (7.35 52.31)
2   52.1     Id02     45     Id04 POINT (7.22 52.18)

Guess the Metaphor

A spatial join is like the pirate islands: the islands are polygons and the ships are points. To learn whose waters each ship sails in, you match them by location, not by a shared key, and a ship out in open water finds no island at all.

Key Ideas for Today

1. An sf object is a data frame plus geometries.
2. Unary functions measure or reshape a geometry.
3. Binary predicates test how two geometries relate.
4. A spatial join matches rows by location.

For Next Time

1. Sign up for Quiz 5.
2. Work on PS 5, due Saturday 10/31 at 11:59 pm.