25  Advanced Spatial Analysis

Today we will be focusing on the practice of fancy geospatial data analysis and visualization.

25.1 Load libraries

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

25.2 Data

The school database from California School Campus Database. I have posted a cleaned up version for the Inland Empire to this class github.

schoolURL <- 'https://raw.githubusercontent.com/RadicalResearchLLC/WarehouseAssemblyBill/main/IEschools.geojson'
schools <- sf::read_sf(dsn = schoolURL) %>% 
  st_transform(crs = 4326)

head(schools)
Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -117.6022 ymin: 33.89027 xmax: -117.4624 ymax: 34.11768
Geodetic CRS:  WGS 84
# A tibble: 6 × 6
  School             District City  GradesServed Level                  geometry
  <chr>              <chr>    <chr> <chr>        <chr>        <MULTIPOLYGON [°]>
1 Cucamonga Middle   Central… Ranc… 6-8          Inte… (((-117.5996 34.11482, -…
2 Allan Orrenmaa El… Alvord … Rive… K-5          Elem… (((-117.4766 33.89399, -…
3 Alvord Continuati… Alvord … Rive… 10-12        High… (((-117.4885 33.89154, -…
4 Arizona Middle     Alvord … Rive… 6-8          Inte… (((-117.4637 33.89473, -…
5 Arlanza Elementary Alvord … Rive… K-5          Elem… (((-117.4684 33.94243, -…
6 Collett Elementary Alvord … Rive… K-5          Elem… (((-117.479 33.90745, -1…

fig-Schools would show the location of school campuses in the Inland Empire if I wasn’t running out of space.

palGrades <- colorFactor(palette = 'Dark2', domain = schools$Level)

leaflet() %>% 
  addTiles() %>% 
  setView(lat = 33.9, lng = -117.4, zoom = 9) %>% 
  addProviderTiles(providers$CartoDB.PositronNoLabels) %>% 
  addPolygons(data = schools,
               color = ~palGrades(Level),
               fillOpacity = 0.8,
               weight = 2,
               group = 'Schools',
               label = ~htmlEscape(School)) %>% 
  addLegend(data = schools,
            pal = palGrades,
            values = ~Level,
            position = 'bottomleft')

Let’s import the planned Warehouse layer for comparison.

I’ve shown the quiet = TRUE argument in st_read to remove the message import from st_read()

plannedWH.url <- 'https://raw.githubusercontent.com/RadicalResearchLLC/PlannedWarehouses/main/plannedWarehouses.geojson'
plannedWarehouses <- st_read(plannedWH.url, quiet = TRUE) %>% 
  st_transform("+proj=longlat +ellps=WGS84 +datum=WGS84")

Add the planned warehouses next to the schools and zoom in a bit to the City of Fontana.

fig-plannedWHs shows the data for schools in the Inland Empire with planned warehouse projects, or it would if I wasn’t running out of space.

palGrades <- colorFactor(palette = 'Dark2', domain = schools$Level)

leaflet() %>% 
  addTiles() %>% 
  setView(lat = 34, lng = -117.5, zoom = 11) %>% 
  addProviderTiles(providers$CartoDB.PositronNoLabels) %>% 
  addPolygons(data = schools,
               color = ~palGrades(Level),
               fillOpacity = 0.8,
               weight = 2,
               group = 'Schools',
               label = ~htmlEscape(School)) %>% 
  addLegend(data = schools,
            pal = palGrades,
            values = ~Level,
            position = 'bottomleft') %>% 
  addPolygons(data = plannedWarehouses,
              color = 'black',
              fillOpacity = 0.3,
              weight = 2,
              group = 'Planned and Approved Warehouses')

Figure 25.1: Inland Empire public schools

25.2.1 Analysis

Create a 1,000 foot buffer around planned warehouses. The goal is to identify schools within a set distance of warehouses.

The st_buffer() function can be used to create shape specific boundaries around polygons.

buffWH1000 <- plannedWarehouses %>% 
  #default units are meters, 304 m = 1000 ft
  st_buffer(dist = 304)

Figure 25.2 shows the schools with their new buffer as a gray overlay.

leaflet() %>% 
  addTiles() %>% 
  addProviderTiles(provider = providers$CartoDB.PositronNoLabels) %>% 
  setView(lat = 34, lng = -117.5, zoom = 12) %>% 
  addPolygons(data = buffWH1000,
              fillOpacity = 0.4,
              color = 'gray',
              weight = 1,
              group = '1000 foot buffer') %>% 
    addPolygons(data = plannedWarehouses,
              color = 'black',
              fillOpacity = 0.3,
              weight = 2,
              group = 'Planned and Approved Warehouses') %>% 
  addPolygons(data = schools,
               color = ~palGrades(Level),
               fillOpacity = 0.8,
               weight = 2,
               group = 'Schools',
               label = ~htmlEscape(School)) %>% 
  addLegend(data = schools,
            pal = palGrades,
            values = ~Level,
            position = 'bottomleft') #%>%

Figure 25.2: Testing buffer around schools to make sure the calculation and distance looks right

Alright, we have the two pieces - a buffered set of schools and a set of planned warehouses layer.

The next function we’ll use is st_join(). st_join() is a spatial filter that joins records of one geospatial dataset (x) with another geospatial dataset (y) through spatial intersections.

If the two spatial layers intersect (i.e., touch or overlap), the records will be returned.

nearWH1kSchools <- buffWH1000 %>% 
  # join by things that intersect
  st_join(schools, left = FALSE) %>% 
  # remove doubled geometries
  st_set_geometry(value = NULL) %>% 
  # re
  left_join(schools) %>% 
  st_as_sf() %>% 
  distinct()
Joining with `by = join_by(School, District, City, GradesServed, Level)`

Figure 25.3 shows the schools within 1,000 feet of planned warehouses.

leaflet() %>% 
  leaflet() %>% 
  addTiles() %>% 
  addProviderTiles(provider = providers$CartoDB.PositronNoLabels) %>% 
  setView(lat = 34, lng = -117.5, zoom = 10) %>% 
  addPolygons(data = nearWH1kSchools,
               color = ~palGrades(Level),
               fillOpacity = 0.4,
               weight = 2,
               group = 'Schools',
               label = ~htmlEscape(School)) %>% 
  addLegend(data = schools,
            pal = palGrades,
            values = ~Level,
            position = 'bottomleft') %>% 
  addPolygons(data = plannedWarehouses,
            color = 'black',
            fillOpacity = 0.3,
            weight = 1,
            group = 'Planned and Approved Warehouses') #%>%

Figure 25.3: Schools within 1,000 feet of planned Warehouses in the Inland Empire

The table below shows a sortable table of schools within 1,000 feet of planned warehouses.

I am using the DT package to make an interactive and sortable table.

narrowNearSchools <- nearWH1kSchools %>% 
  st_set_geometry(value = NULL) %>% 
  select(School, Level, District, City, name) %>% 
  distinct() %>% 
  rename(WarehouseName = name)
DT::datatable(narrowNearSchools,
               rownames = FALSE,
               options = list(dom = 'fp',
                              pageLength = 15))