Working with spatial and geometric operations in `{sf}`

Session 5

2023-09-27

1 Getting started

library(tidyverse)
library(sf)
library(tigris)
options(tigris_use_cache = TRUE)

options("tigris_use_cache" = TRUE)
us_states <- states()
us_highways <- primary_roads()
maryland <- filter(us_states, NAME == "Maryland")
lower48 <- filter(us_states, !(NAME %in% c("United States Virgin Islands", "Commonwealth of the Northern Mariana Islands", "Guam", "Alaska", "American Samoa", "Puerto Rico", "Hawaii")))

2 Converting objects between types

sfg <- st_point(c(-76.712838, 39.253383))

sfg

st_as_sfc(list(sfg))

Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -76.71284 ymin: 39.25338 xmax: -76.71284 ymax: 39.25338
CRS:           NA

sfc <- st_as_sfc(list(sfg), crs = 4326)

sfc

Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -76.71284 ymin: 39.25338 xmax: -76.71284 ymax: 39.25338
Geodetic CRS:  WGS 84

st_as_sf(sfc)

Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -76.71284 ymin: 39.25338 xmax: -76.71284 ymax: 39.25338
Geodetic CRS:  WGS 84
                           x
1 POINT (-76.71284 39.25338)

sf_df <- data.frame(
  name = "Example data frame",
  lon = -76.712838,
  lat = 39.253383
)

sf <- st_as_sf(sf_df, coords = c("lon", "lat"), crs = 4326)

sf

Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -76.71284 ymin: 39.25338 xmax: -76.71284 ymax: 39.25338
Geodetic CRS:  WGS 84
                name                   geometry
1 Example data frame POINT (-76.71284 39.25338)

3 Types of spatial operations

Spatial filtering and topological relations
Spatial joins and non-overlapping joins
Spatial aggregation

Spatial filtering and topological relations

st_filter() efaults to intersecting geometries:

st_filter(us_states, maryland)

Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -83.67539 ymin: 36.54085 xmax: -74.68956 ymax: 42.51607
Geodetic CRS:  NAD83
  REGION DIVISION STATEFP  STATENS GEOID STUSPS                 NAME LSAD MTFCC
1      3        5      54 01779805    54     WV        West Virginia   00 G4000
2      3        5      24 01714934    24     MD             Maryland   00 G4000
3      3        5      10 01779781    10     DE             Delaware   00 G4000
4      1        2      42 01779798    42     PA         Pennsylvania   00 G4000
5      3        5      51 01779803    51     VA             Virginia   00 G4000
6      3        5      11 01702382    11     DC District of Columbia   00 G4000
  FUNCSTAT        ALAND     AWATER    INTPTLAT     INTPTLON
1        A  62266456923  489045863 +38.6472854 -080.6183274
2        A  25151771744 6979295311 +38.9466584 -076.6744939
3        A   5046703789 1399207462 +38.9985661 -075.4416440
4        A 115882119641 3397575101 +40.9046042 -077.8275233
5        A 102258178227 8528072639 +37.5222512 -078.6681938
6        A    158316184   18709787 +38.9042474 -077.0165167
                        geometry
1 MULTIPOLYGON (((-77.75438 3...
2 MULTIPOLYGON (((-75.756 39....
3 MULTIPOLYGON (((-75.50949 3...
4 MULTIPOLYGON (((-75.5931 39...
5 MULTIPOLYGON (((-76.4915 36...
6 MULTIPOLYGON (((-77.11975 3...

st_filter() can also support alternate topological relations:

st_filter(us_states, maryland, .predicate = st_disjoint)

Simple feature collection with 50 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
Geodetic CRS:  NAD83
First 10 features:
   REGION DIVISION STATEFP  STATENS GEOID STUSPS           NAME LSAD MTFCC
1       3        5      12 00294478    12     FL        Florida   00 G4000
2       2        3      17 01779784    17     IL       Illinois   00 G4000
3       2        4      27 00662849    27     MN      Minnesota   00 G4000
4       1        1      44 01219835    44     RI   Rhode Island   00 G4000
5       4        8      16 01779783    16     ID          Idaho   00 G4000
6       1        1      33 01779794    33     NH  New Hampshire   00 G4000
7       3        5      37 01027616    37     NC North Carolina   00 G4000
8       1        1      50 01779802    50     VT        Vermont   00 G4000
9       1        1      09 01779780    09     CT    Connecticut   00 G4000
10      4        8      35 00897535    35     NM     New Mexico   00 G4000
   FUNCSTAT        ALAND      AWATER    INTPTLAT     INTPTLON
1         A 138962819934 45971472526 +28.3989775 -082.5143005
2         A 143778515726  6216539665 +40.1028754 -089.1526108
3         A 206244837557 18937184315 +46.3159573 -094.1996043
4         A   2677763373  1323686975 +41.5964850 -071.5264901
5         A 214049908397  2391592787 +44.3484222 -114.5588538
6         A  23190126218  1025960758 +43.6726907 -071.5843145
7         A 125935585728 13453835222 +35.5397100 -079.1308636
8         A  23872569964  1030754610 +44.0589536 -072.6710173
9         A  12541666998  1816447681 +41.5798637 -072.7466572
10        A 314198573403   726463825 +34.4346843 -106.1316181
                         geometry
1  MULTIPOLYGON (((-83.10874 2...
2  MULTIPOLYGON (((-87.89243 3...
3  MULTIPOLYGON (((-95.31989 4...
4  MULTIPOLYGON (((-71.67881 4...
5  MULTIPOLYGON (((-116.3584 4...
6  MULTIPOLYGON (((-70.83887 4...
7  MULTIPOLYGON (((-77.89977 3...
8  MULTIPOLYGON (((-72.04187 4...
9  MULTIPOLYGON (((-72.5279 41...
10 MULTIPOLYGON (((-103.0647 3...

Spatial joins and non-overlapping joins

st_join(maryland, us_highways)

Simple feature collection with 470 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.48765 ymin: 37.88661 xmax: -74.98628 ymax: 39.72304
Geodetic CRS:  NAD83
First 10 features:
    REGION DIVISION STATEFP  STATENS GEOID STUSPS     NAME LSAD MTFCC.x
1        3        5      24 01714934    24     MD Maryland   00   G4000
1.1      3        5      24 01714934    24     MD Maryland   00   G4000
1.2      3        5      24 01714934    24     MD Maryland   00   G4000
1.3      3        5      24 01714934    24     MD Maryland   00   G4000
1.4      3        5      24 01714934    24     MD Maryland   00   G4000
1.5      3        5      24 01714934    24     MD Maryland   00   G4000
1.6      3        5      24 01714934    24     MD Maryland   00   G4000
1.7      3        5      24 01714934    24     MD Maryland   00   G4000
1.8      3        5      24 01714934    24     MD Maryland   00   G4000
1.9      3        5      24 01714934    24     MD Maryland   00   G4000
    FUNCSTAT       ALAND     AWATER    INTPTLAT     INTPTLON      LINEARID
1          A 25151771744 6979295311 +38.9466584 -076.6744939 1104258795090
1.1        A 25151771744 6979295311 +38.9466584 -076.6744939 1103075320064
1.2        A 25151771744 6979295311 +38.9466584 -076.6744939 1104258795072
1.3        A 25151771744 6979295311 +38.9466584 -076.6744939 1104258845568
1.4        A 25151771744 6979295311 +38.9466584 -076.6744939  110191599447
1.5        A 25151771744 6979295311 +38.9466584 -076.6744939  110191599448
1.6        A 25151771744 6979295311 +38.9466584 -076.6744939 1102149806067
1.7        A 25151771744 6979295311 +38.9466584 -076.6744939 1102149813008
1.8        A 25151771744 6979295311 +38.9466584 -076.6744939 1102149806121
1.9        A 25151771744 6979295311 +38.9466584 -076.6744939 1102166635559
         FULLNAME RTTYP MTFCC.y                       geometry
1   Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.1 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.2 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.3 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.4 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.5 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.6 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.7 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.8 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....
1.9 Salisbury Byp     M   S1100 MULTIPOLYGON (((-75.756 39....

st_is_within_distance() is an example of a non-overlapping join:

st_is_within_distance(us_highways, maryland, dist = units::set_units(100, "mi"))

Sparse geometry binary predicate list of length 17389, where the
predicate was `is_within_distance'
first 10 elements:
 1: (empty)
 2: (empty)
 3: (empty)
 4: (empty)
 5: (empty)
 6: (empty)
 7: (empty)
 8: (empty)
 9: (empty)
 10: (empty)

You can use these joins for spatial filters using the base R syntax:

is_hwy_within_dist <- st_is_within_distance(us_highways, maryland, dist = units::set_units(100, "mi"), sparse = FALSE)
us_highways[is_hwy_within_dist, ]

Simple feature collection with 2241 features and 4 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -81.72573 ymin: 36.57507 xmax: -74.04969 ymax: 41.34626
Geodetic CRS:  NAD83
First 10 features:
          LINEARID      FULLNAME RTTYP MTFCC                       geometry
17  11013853809147 State Rte 611     S S1100 LINESTRING (-75.1366 40.350...
18  11013853809270 State Rte 611     S S1100 LINESTRING (-75.12634 40.28...
19  11013853809430 State Rte 611     S S1100 LINESTRING (-75.12603 40.28...
77  11013853864377  State Rte 29     S S1100 LINESTRING (-74.76432 40.20...
224  1105598256809  State Rte 26     S S1100 LINESTRING (-77.73622 40.89...
226  1108296536943  State Rte 26     S S1100 LINESTRING (-74.76814 40.20...
227  1108296536969  State Rte 26     S S1100 LINESTRING (-74.71174 40.27...
229  1104486135566  State Rte 26     S S1100 LINESTRING (-77.73356 40.89...
244  1108296538090 State Rte 129     S S1100 LINESTRING (-74.72239 40.18...
245  1108296538091 State Rte 129     S S1100 LINESTRING (-74.7224 40.183...

Or using filter() if you use the name of the geomtery column (and convert the matrix output into a logical vector):

filter(us_highways, as.logical(st_is_within_distance(geometry, maryland, dist = units::set_units(100, "mi"), sparse = FALSE)))

Simple feature collection with 2241 features and 4 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -81.72573 ymin: 36.57507 xmax: -74.04969 ymax: 41.34626
Geodetic CRS:  NAD83
First 10 features:
         LINEARID      FULLNAME RTTYP MTFCC                       geometry
1  11013853809147 State Rte 611     S S1100 LINESTRING (-75.1366 40.350...
2  11013853809270 State Rte 611     S S1100 LINESTRING (-75.12634 40.28...
3  11013853809430 State Rte 611     S S1100 LINESTRING (-75.12603 40.28...
4  11013853864377  State Rte 29     S S1100 LINESTRING (-74.76432 40.20...
5   1105598256809  State Rte 26     S S1100 LINESTRING (-77.73622 40.89...
6   1108296536943  State Rte 26     S S1100 LINESTRING (-74.76814 40.20...
7   1108296536969  State Rte 26     S S1100 LINESTRING (-74.71174 40.27...
8   1104486135566  State Rte 26     S S1100 LINESTRING (-77.73356 40.89...
9   1108296538090 State Rte 129     S S1100 LINESTRING (-74.72239 40.18...
10  1108296538091 State Rte 129     S S1100 LINESTRING (-74.7224 40.183...

Spatial aggregation

us_states_combine <- us_states |>
  group_by(DIVISION) |>
  summarise(
    geometry = st_combine(geometry)
  )

plot(us_states_combine[, 1])

4 Types of geometric operations

Simplification
Centroids
Buffers
Clipping (with or without subsetting)
Unions
Type transformations

Simplification

plot(st_simplify(lower48, dTolerance = 100000)[, 1])

Centroids

plot(st_centroid(lower48)[, 1])

Buffers

plot(st_buffer(lower48, dist = 100)[, 1])

plot(st_buffer(lower48, dist = units::set_units(100, "mi"))[, 1])

Clipping (with or without subsetting)

plot(st_intersection(us_highways, maryland)[, 1])

Unions

us_highways_union <- us_highways |>
  group_by(RTTYP) |>
  summarise(
    geometry = st_union(geometry)
  )

plot(us_highways_union[, 1])

Type transformations

plot(st_cast(us_highways[1, ], to = "POINT"))

Working with spatial and geometric operations in {sf}