Spatial II

Computing on geometries with sf

(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 / web APIs
  2. Shapefiles

Shapefiles

A shapefile is a bundle of files that describe spatial data. Developed for ESRI ArcView GIS system (proprietary).

Google: CA shapefile

Shapefiles

A shapefile is a bundle of files that describe spatial data. Developed for ESRI ArcView GIS system (proprietary).

library(sf)
ca <- st_read("data/CA_State.shp")
Reading layer `CA_State' from data source 
  `/workspace/27-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
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...

Shapefiles

A shapefile is a bundle of files that describe spatial data. Developed for ESRI ArcView GIS system (proprietary).

plot(st_geometry(ca))

Sources of sf data

  1. R packages / web APIs
  2. Shapefiles
  1. GeoJSON: Web-friendly, single file
  2. KLM: 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)

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)

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)

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)

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), (...

Computing on 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]

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 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)

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)

gr
Simple feature collection with 70 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 406265 ymin: 347094.8 xmax: 3052877 ymax: 1044143
Projected CRS: NAD83 / North Carolina (ftUS)
First 10 features:
   label                           geom       col
31   G 1 POLYGON ((406265 347094.8, ... #8dd3c74d
32   G 2 POLYGON ((670926.1 347094.8... #ffffb34d
33   G 3 POLYGON ((935587.3 347094.8... #bebada4d
34   G 4 POLYGON ((1200248 347094.8,... #fb80724d
35   G 5 POLYGON ((1464910 347094.8,... #80b1d34d
36   G 6 POLYGON ((1729571 347094.8,... #fdb4624d
37   G 7 POLYGON ((1994232 347094.8,... #b3de694d
38   G 8 POLYGON ((2258893 347094.8,... #fccde54d
39   G 9 POLYGON ((2523554 347094.8,... #d9d9d94d
40   G10 POLYGON ((2788215 347094.8,... #bc80bd4d

nc
Simple feature collection with 100 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
Projected CRS: NAD83 / North Carolina (ftUS)
# A tibble: 100 × 15
    AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS  FIPSNO CRESS_ID BIR74 SID74 NWBIR74
 * <dbl>     <dbl> <dbl>   <dbl> <chr> <chr>  <dbl>    <int> <dbl> <dbl>   <dbl>
 1 0.114      1.44  1825    1825 Ashe  37009  37009        5  1091     1      10
 2 0.061      1.23  1827    1827 Alle… 37005  37005        3   487     0      10
 3 0.143      1.63  1828    1828 Surry 37171  37171       86  3188     5     208
 4 0.07       2.97  1831    1831 Curr… 37053  37053       27   508     1     123
 5 0.153      2.21  1832    1832 Nort… 37131  37131       66  1421     9    1066
 6 0.097      1.67  1833    1833 Hert… 37091  37091       46  1452     7     954
 7 0.062      1.55  1834    1834 Camd… 37029  37029       15   286     0     115
 8 0.091      1.28  1835    1835 Gates 37073  37073       37   420     0     254
 9 0.118      1.42  1836    1836 Warr… 37185  37185       93   968     4     748
10 0.124      1.43  1837    1837 Stok… 37169  37169       85  1612     1     160
# ℹ 90 more rows
# ℹ 4 more variables: BIR79 <dbl>, SID79 <dbl>, NWBIR79 <dbl>,
#   geom <MULTIPOLYGON [US_survey_foot]>

“Add to each county in NC the label and color of the grid cell that it most intersects.”

nc_j <- nc |>
  st_join(gr, join = st_intersects, largest = TRUE)
nc_j
Simple feature collection with 100 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
Projected CRS: NAD83 / North Carolina (ftUS)
# A tibble: 100 × 17
    AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS  FIPSNO CRESS_ID BIR74 SID74 NWBIR74
 * <dbl>     <dbl> <dbl>   <dbl> <chr> <chr>  <dbl>    <int> <dbl> <dbl>   <dbl>
 1 0.114      1.44  1825    1825 Ashe  37009  37009        5  1091     1      10
 2 0.061      1.23  1827    1827 Alle… 37005  37005        3   487     0      10
 3 0.143      1.63  1828    1828 Surry 37171  37171       86  3188     5     208
 4 0.07       2.97  1831    1831 Curr… 37053  37053       27   508     1     123
 5 0.153      2.21  1832    1832 Nort… 37131  37131       66  1421     9    1066
 6 0.097      1.67  1833    1833 Hert… 37091  37091       46  1452     7     954
 7 0.062      1.55  1834    1834 Camd… 37029  37029       15   286     0     115
 8 0.091      1.28  1835    1835 Gates 37073  37073       37   420     0     254
 9 0.118      1.42  1836    1836 Warr… 37185  37185       93   968     4     748
10 0.124      1.43  1837    1837 Stok… 37169  37169       85  1612     1     160
# ℹ 90 more rows
# ℹ 6 more variables: BIR79 <dbl>, SID79 <dbl>, NWBIR79 <dbl>,
#   geom <MULTIPOLYGON [US_survey_foot]>, label <chr>, col <chr>

What is the total value of the pink polygon on the right?

01:30

Spatial Interpolation

regions_sf <- st_sf(geometry = st_sfc(
    st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))),
    st_polygon(list(rbind(c(1,0), c(2,0), c(2,1), c(1,1), c(1,0))))),
  crs = 3857)
regions_sf
Simple feature collection with 2 features and 0 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
Projected CRS: WGS 84 / Pseudo-Mercator
                        geometry
1 POLYGON ((0 0, 1 0, 1 1, 0 ...
2 POLYGON ((1 0, 2 0, 2 1, 1 ...

Spatial Interpolation

# Make a grid that cuts across both
grid_sf <- st_make_grid(regions_sf, n = c(3, 1)) |> 
  st_sf()
grid_sf$value = c(10, 20, 30)
grid_sf
Simple feature collection with 3 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
Projected CRS: WGS 84 / Pseudo-Mercator
  st_make_grid.regions_sf..n...c.3..1.. value
1        POLYGON ((0 0, 0.6666667 0,...    10
2        POLYGON ((0.6666667 0, 1.33...    20
3        POLYGON ((1.333333 0, 2 0, ...    30

Spatial Interpolation

interp <- st_interpolate_aw(grid_sf["value"], regions_sf, extensive = TRUE)
interp
Simple feature collection with 2 features and 1 field
Attribute-geometry relationship: aggregate (1)
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
Projected CRS: WGS 84 / Pseudo-Mercator
  value                       geometry
1    20 POLYGON ((0 0, 1 0, 1 1, 0 ...
2    40 POLYGON ((1 0, 2 0, 2 1, 1 ...