Geog 380A Week 5

Welcome to Week 5! This week, we’ll be moving beyond mapping to start learning some basic spatial data operations.

#load in the libraries
library(sf)     

Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

library(dplyr)   

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(spData) 

To access larger datasets in this package, install the
spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

library(urbnmapr)
library(ggplot2)

We’ll start by loading in shapefile data.

setwd("~/Binghamton/geog380/NYS_Place_Points_SHP")

Warning: The working directory was changed to C:/Users/mhaller/Documents/Binghamton/geog380/NYS_Place_Points_SHP inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

place_shape <- st_read("NYS_Place_Points.shp")

Reading layer `NYS_Place_Points' from data source 
  `C:\Users\mhaller\Documents\Binghamton\geog380\NYS_Place_Points_SHP\NYS_Place_Points.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 4571 features and 8 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 108752.5 ymin: 4484510 xmax: 756007.4 ymax: 4984062
Projected CRS: NAD83 / UTM zone 18N

place_shape <- place_shape %>% st_transform("EPSG:32116")

#load in NYS county shape file and set crs
counties <- get_urbn_map("counties", sf = TRUE)

#filter the data to get just NYS
counties <- counties %>% 
  filter(state_abbv == "NY") %>% 
  st_transform("EPSG:32116")

old-style crs object detected; please recreate object with a recent sf::st_crs()

#Let's map them now!
ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = place_shape,
          mapping = aes(), color = "green")+
  theme_minimal()

Now, what if we just wanted to get points of interest for Broome County?

broome <- counties %>% 
  filter(county_name == "Broome County")

#Map Broome county to inspect that it's correct:

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue", color = "black")+
  theme_minimal()

Looks good! Now, we will filter our place data based on the Broome County layer, like this:

#By default, this will behave like the function st_intersects, which returns all features that touch, cross or are within the source ‘subsetting’ object

broome_points <- place_shape[broome,]

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue", color = "black")+
  geom_sf(data = broome_points,
          mapping = aes(), color = "forestgreen")+
  theme_minimal()

That’s a good start! But I’d like to zoom in on Broome county a bit more. How can I do that?

#Let's check the dimensions of the Broome county shape
broome

Simple feature collection with 1 feature and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 282328.4 ymin: 223493.3 xmax: 342746.8 ymax: 269994.9
Projected CRS: NAD83 / New York Central
  county_fips state_abbv state_fips   county_name fips_class
1       36007         NY         36 Broome County         H1
  state_name                       geometry
1   New York MULTIPOLYGON (((282328.4 26...

#We can use xlim and ylim to set the dimensions of the plot (just larger than the dimensions of Broome County)

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue", color = "black")+
  geom_sf(data = broome_points,
          mapping = aes(), color = "forestgreen")+
  xlim(280000, 350000)+
  ylim(220000, 270000)+
  theme_minimal()

What if we wanted to subset our point data by counties that are not Broome county (in other words, all other counties)? We can easily do that by specifying the op= argument (“op” stands for operator) within our subsetting:

#The st_disjoint argument tells R to return all objects that do not intersect with our subsetting object
not_broome <- place_shape[broome, , op = st_disjoint]

#Let's map it now
ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue", color = "black")+
  geom_sf(data = not_broome,
          mapping = aes(), color = "forestgreen")+
  theme_minimal()

Below are some spatial operations that you may use on your data:

st_intersects: Does object 1 touch, cross, or contain object 2? st_disjoint: Does object 1 have no spatial relationship with object 2? st_contains: Is object 2 fully enclosed or contained by object 1? st_within: Is object 1 fully within object 2? (*this is the reverse of st_contains) st_touches: Do object 1 and object 2 have at least one point in common (but their interiors do not intersect)? st_overlap: Do object 1 and object 2 share space without completely containing each other? st_equals: Are object 1 and object 2 the same?

I’ll illustrate some examples of these:

#load in roads data
setwd("~/Binghamton/geog380/SimplifiedStreets_SHP")

Warning: The working directory was changed to C:/Users/mhaller/Documents/Binghamton/geog380/SimplifiedStreets_SHP inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

street_shape <- st_read("SimplifiedStreetSegmentQrt.shp")

Reading layer `SimplifiedStreetSegmentQrt' from data source 
  `C:\Users\mhaller\Documents\Binghamton\geog380\SimplifiedStreets_shp\SimplifiedStreetSegmentQrt.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 236888 features and 9 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 105570 ymin: 4483142 xmax: 764038.9 ymax: 4985389
Projected CRS: NAD83 / UTM zone 18N

street_shape <- street_shape %>% st_transform("EPSG:32116")

broome_streets <- street_shape[broome, , op = st_intersects]

#let's visualize the intersection of the streets first as a comparison
ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue", color = "black")+
  geom_sf(data = broome_streets,
          mapping = aes(), color = "forestgreen")+
  xlim(280000, 350000)+
  ylim(220000, 270000)+
  theme_minimal()

#This looks cool, but I have streets that go outside of the county! Let's subset by streets that are fully contained within Broome county to clean it up:
#Note that I use st_within and now st_contains - can anyone tell me why? Which object is Object 1 and which is Object 2?

broome_streets <- street_shape[broome, , op = st_within]

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue", color = "black")+
  geom_sf(data = broome_streets,
          mapping = aes(), color = "forestgreen")+
  xlim(280000, 350000)+
  ylim(220000, 270000)+
  theme_minimal()

This looks much nicer! What does it look like if I subset by roads that aren’t in Broome County? That way, I can show all of the roads on the map, but I can make roads inside and outside of Broome county different colors

not_broome_streets <- street_shape[broome, , op = st_disjoint]
broome_streets <- street_shape[broome, , op = st_within]

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightgray", color = "black")+
  geom_sf(data = broome_streets,
          mapping = aes(), color = "gray")+
  geom_sf(data = not_broome_streets,
          mapping = aes(), color = "gray")+
  xlim(280000, 350000)+
  ylim(220000, 270000)+
  theme_minimal()

What if I wanted to map counties that border Broome county? I can do that pretty easily, too! Does anyone have any guesses as to what spatial operation I should use?

broome_borders <- counties[broome, , op = st_touches]


ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = broome,
          mapping = aes(), fill = "lightblue4", color = "black")+
  geom_sf(data = broome_borders,
          mapping = aes(), fill = "lightblue1")+
  theme_minimal()

What if I wanted to filter my spatial data by proximity? We can do this too! Let’s start by calculating the distance between different counties

#Start by finding the central-most point for each county
county_centroid = st_centroid(counties)

Warning: st_centroid assumes attributes are constant over geometries of x

#Calculate distances between each county
head(st_distance(county_centroid))

Units: [m]
          1        2         3         4        5         6         7
1      0.00 126400.6 308765.42 288206.25 266049.3  72229.08 115750.75
2 126400.58      0.0 305160.53 303561.87 383249.5 185354.56 220623.91
3 308765.42 305160.5      0.00  46650.44 358775.6 279252.41 259752.50
4 288206.25 303561.9  46650.44      0.00 313640.3 249781.30 225157.97
5 266049.25 383249.5 358775.65 313640.29      0.0 198225.95 164666.98
6  72229.08 185354.6 279252.41 249781.30 198226.0      0.00  44206.46
         8         9        10        11       12       13       14
1 225655.5 361779.43  98169.69 234614.17 391670.7 166851.4 273365.8
2 315128.2 354279.18 197645.36 249172.98 417732.7 103818.7 357753.7
3 230617.8  53014.16 251164.16  77519.45 129461.3 408241.5 230351.4
4 185496.2  86480.92 219674.42  55568.35 114545.0 404106.2 183711.5
5 128183.6 397390.62 188242.28 295482.67 348221.8 430726.0 143047.9
6 155546.3 331167.85  32461.60 201787.07 342566.0 238853.7 203821.4
         15        16        17       18       19        20       21
1 277678.53 379938.06 267122.29 430799.3 171936.2 382438.64 415453.5
2 400756.90 387775.37 273727.88 450673.1 281651.3 400761.71 405202.6
3 406564.80  84175.63  44092.31 152052.9 281389.9 105055.88 106688.5
4 362173.59  91782.06  32682.21 147657.4 240512.1  97714.03 136100.1
5  52320.21 379125.63 320690.00 388228.7 104753.8 359013.40 439183.4
6 216584.08 340347.68 235351.59 383459.2 100160.4 337431.65 383933.4
         22        23       24        25       26        27        28
1 346714.38 166846.57 242909.1 220484.32 182390.3 172744.17 156826.43
2 375153.50  42511.41 307086.4 220797.19 238292.1  52741.14  36175.78
3 101032.06 329063.36 163930.2  89354.70 157601.3 343931.94 332623.27
4  74625.26 332360.77 117998.9  83404.56 121254.8 347061.90 334169.34
5 315970.43 425722.59 196055.0 314279.68 216661.4 434268.31 417654.92
6 298135.80 227732.97 182147.3 197961.39 131763.4 236052.35 219457.20
         29        30        31        32        33        34
1 343868.69 109638.42  47962.85 307846.27 207978.41  40521.77
2 356266.54  36355.88 161090.74 330136.99 330421.99 125316.05
3  60446.45 272576.34 282314.57  65440.21 357787.49 268295.50
4  55763.72 268798.77 256332.83  29108.68 315927.64 247889.64
5 347882.92 354814.62 223132.19 308938.20  67591.85 257937.76
6 303578.80 158888.82  25269.56 264375.61 146268.13  60515.68
        35        36        37        38       39        40        41
1 180517.2 169165.91 122628.34 119361.63 164878.1  54755.47 258006.63
2 186606.0 263444.90 184675.13 245272.91 198829.3  72070.55 252507.08
3 129328.4 235220.27 197123.85 342923.19 147706.2 305902.63  53241.28
4 117932.7 194835.63 169557.01 308445.39 123339.1 293324.72  61652.31
5 296922.2 140905.42 225612.20 158501.77 256154.4 317703.23 337279.01
6 161952.1  99090.08  82249.17  75025.84 132101.7 120026.91 234059.14
         42        43       44        45       46        47        48
1 203830.82  97212.29 110128.3 168311.18 421960.0 250490.98  52265.88
2 233331.24 217875.79 140190.2  58984.98 428484.2 286610.62 178574.92
3 113763.58 302525.02 199095.7 359017.49 123550.0  99112.23 321502.66
4  84798.08 269445.32 182088.8 360266.88 133758.5  54933.59 294156.19
5 263803.10 168840.55 271630.9 432704.30 410225.4 259663.60 218721.93
6 166315.84  35834.00 101590.2 234753.33 381466.5 203738.03  46387.72
         49        50        51        52        53        54
1 121530.81  92131.53  65583.59 181038.60 295350.75 275409.38
2  22466.29  38189.77  84115.56  57525.46 333171.64 381464.36
3 325198.77 315186.60 261341.12 339342.24 102992.99 310553.96
4 321875.31 307932.36 248503.63 344058.23  58485.68 264322.79
5 384136.95 354935.37 300768.92 440664.13 270710.17  61590.59
6 185979.87 156902.32 104808.28 242616.73 244598.14 203459.57
         55        56       57        58        59        60       61
1 268149.95  76388.40 183981.1 190633.83 146863.26 225867.76 472157.6
2 293788.85 163390.83 259780.8  64391.06 268346.78 273036.22 459468.3
3  73201.56 244069.31 192422.6 327989.59 322061.12 127302.78 163392.2
4  28244.69 218042.51 151830.0 334950.06 284113.96  85078.54 191086.3
5 285822.50 223157.20 175515.8 446181.04 120330.48 232766.24 486302.2
6 225669.44  40505.17 121450.8 248906.94  84386.61 175228.51 440029.2
         62
1 111170.61
2  88167.08
3 224678.48
4 217949.81
5 322106.44
6 135860.86

It looks funky, but what’s actually happening here is we’re creating a distance matrix - this is a matrix with 62 rows and 62 columns - each row/column pair tells us the distance between two counties.

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = county_centroid,
          mapping = aes(), color = "black")+
  theme_minimal()

NA
NA

#Let's get the distance between Broome county and all other counties
broome_centroid <- st_centroid(broome)

Warning: st_centroid assumes attributes are constant over geometries of x

#get the distance between Broome and every other county
broome_dist <- as.data.frame(st_distance(county_centroid, broome_centroid))

#Turn those distances into a dataframe with a distance and county column
#County names come from the counties sf object
broome_dist <- cbind(broome_dist, counties$county_name)
colnames(broome_dist) <- c("distance", "county")

#Arrange in ascending order by distance and grab ten closest counties
broome_dist <- broome_dist %>% 
  arrange(distance) %>% 
  head(n = 10)

#Create new county object that contains only the ten closest counties to Broome
close_counties <- counties %>% 
  filter(county_name %in% broome_dist$county)

#Now we can map it
ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(data = close_counties,
          mapping = aes(), fill = "aquamarine1")+
  geom_sf(data = broome,
          mapping = aes(), fill = "aquamarine4", color = "black")+
  geom_sf(data = county_centroid,
          mapping = aes(), color = "green")+
  theme_minimal()

Why is it important to use centroids, and not just distances from county polygons? Let’s compare the output from each:

broome_dist <- as.data.frame(st_distance(counties, broome))

broome_dist <- cbind(broome_dist, counties$county_name)
colnames(broome_dist) <- c("distance", "county")

broome_dist <- broome_dist %>% 
  arrange(distance) %>% 
  head(n = 10)

close_polygons <- counties %>% 
  filter(county_name %in% broome_dist$county)

ggplot() +
  geom_sf(data = counties,
    mapping = aes(), color = "lightgray")+
  geom_sf(close_polygons, 
          mapping = aes(), fill = "red")+
  geom_sf(data = close_counties,
          mapping = aes(), fill = "aquamarine1", alpha = 0.8)+
  geom_sf(data = broome,
          mapping = aes(), fill = "aquamarine4", color = "black")+
  geom_sf(data = county_centroid,
          mapping = aes(), color = "green")+
  theme_minimal()

The red counties are counties that are close (top ten) based on polygon distance, but not centroid distance! Most of the counties are the same in both versions, but we can see that that there are some differences. It’s important to think carefully about how you specify distance between objects, as it can change the final product that you end up with.

Resources

Lovelace, R., Nowosad, J., Muenchow. J. (2019). Geocomputation with R. Retrieved from: https://geocompr.robinlovelace.net/.

NYS GIS Clearinghouse (2022). NYS Streets. [Data Set]. Retrieved from: http://gis.ny.gov/gisdata/.

NYS GIS Clearinghouse (2022). NYS Place Points. [Data Set]. Retrieved from: http://gis.ny.gov/gisdata/.

Soltoff, B. (2022). Importing Spatial Data Files Using SF. Computing for the Social Sciences. Retrieved from: https://cfss.uchicago.edu/notes/simple-features/.