Chapter 9: Mapping Point PatternS

9.1 Lab Goals

This chapter explores mapping of point patterns:

• Crating grids covering study area
• Point pattern visualization
• Kernel Density Estimation Mapping

9.2 Good Practice

9.2.1 Organizing Folders & Sub-folders

Under the course folder, please create a folder called “lab9”. Next, in the lab9 folder, please create two sub-folders that one is called “data” and another one is “plot”.

9.2.2 Data

This chapter uses:

• Crash data reported to CSPD in 2024.
• 2023 Tiger Line Road Data from US Census Bureau
• 2023 Cartographic Boundary County Subdivision from US Census Bureau

Download the lab’s data:

• Go to https://github.com/fuzhen-yin/uccs_geoviz/blob/main/data/lab9_data.zip
• Download the file “lab9_data.zip”
• Unzip folder “lab9_data.zip”
• Move all files from the “lab9_data” folder to the “data” folder under “lab9” see Step 9.2.1

If there you have any questions about the above-mentioned steps, please refer to Chapter 3.2.3 for detailed instructions.

9.2.3 Launching R Studio

Again, we would like to start a new project from scratch with a clean R Script. Please do the following steps. If you have any questions about these steps, please refer to the relevant chapters for help.

• Step 1: Make sure all existing R projects are properly closed.
  • If not, please close it by going to File –> Close Project –> Save changes (see Chapter 2.5).
• Step 2: Create a New Project using Existing Directory, navigate to lab9, click open, then Create Project. (see Chapter 1.3).
• Step 3: Create a New Script by go to File –> New File –> R Script. Save the script by giving it a proper name.

9.3 Libraries & Data

9.3.1 Load Libraries

Please install necessary packages as needed.

library(dplyr)    # mutate, left_join
library(sf)
library(geojsonio)  ## import geojson
library(tmap)       ## mapping
library(spatstat)   ## kernel density map

9.3.2 Data

Simple Features (sf) - spatial data

# colorado state: county subdivision (cartographic boundary)
sf_co_cousub <- st_read("data/cb_2023_us_cousub/cb_2023_us_cousub_500k.shp") %>% 
  filter(STATE_NAME == "Colorado")

## Reading layer `cb_2023_us_cousub_500k' from data source 
##   `D:\OneDrive\UCCS\Teaching\GES4070\Labs\01_git_wk\uccs_geoviz\data\cb_2023_us_cousub\cb_2023_us_cousub_500k.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 36367 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1467 ymin: -14.5487 xmax: 179.7785 ymax: 71.38782
## Geodetic CRS:  NAD83

# crash data
sf_crash <- st_read("data/crash_data/cos_crash_2024.shp")

## Reading layer `cos_crash_2024' from data source 
##   `D:\OneDrive\UCCS\Teaching\GES4070\Labs\01_git_wk\uccs_geoviz\data\crash_data\cos_crash_2024.shp' 
##   using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 6017 features and 29 fields (with 160 geometries empty)
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -105.3126 ymin: 38.74474 xmax: -104.6317 ymax: 39.03464
## Geodetic CRS:  WGS 84

# road data 
sf_road <- st_read("data/tl_2023_08041_roads/tl_2023_08041_roads.shp")

## Reading layer `tl_2023_08041_roads' from data source 
##   `D:\OneDrive\UCCS\Teaching\GES4070\Labs\01_git_wk\uccs_geoviz\data\tl_2023_08041_roads\tl_2023_08041_roads.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 29689 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -105.0685 ymin: 38.51945 xmax: -104.0516 ymax: 39.13004
## Geodetic CRS:  NAD83

9.4 Data Processing

9.4.1 Spatial Boundaries

Project to CRS NAD83 / Colorado North (ftUS). CRS code: 2231

# project to CRS NAD83 / Colorado North (ftUS) 
sf_co_cousub <- st_transform(sf_co_cousub, 2231)

# el paso county: county subdivision 
sf_ep_cousub <- sf_co_cousub %>% filter(NAMELSADCO == "El Paso County")
# olorado springs 
sf_cos_cousub <- sf_ep_cousub %>% filter(NAME == "Colorado Springs")

9.4.2 Crash Data

Examine the crash data.

# top 5 rows 
head(sf_crash, 5)
# data type 
str(sf_crash)

Filter data

# extract records with non-empty geometry 
sf_crash_24 <- sf_crash[!st_is_empty(sf_crash),]

# project to CRS NAD83 / Colorado North (ftUS)
sf_crash_24 <- st_transform(sf_crash_24, 2231)

9.4.3 Exploring Analysis

Visualization: quick interactive map.

# set tmap mode for interactive mapping 
tmap_mode("view")

# plot crash data and cos county subdivision
tm_shape(sf_cos_cousub) + tm_polygons(alpha = 0, lwd = 5, border.col = "blue") + 
  tm_shape(sf_crash_24) + tm_dots(alpha = 0.5)

Here are some issues with the crash data, (1) many dots overlap with each other; (2) there are crash dots outside of cos city boundary. Let’s extract only crash within OS city boundary

# extract only crash within OS city boundary 
sf_crash_24_cos <- sf_crash_24[lengths(st_intersects(sf_crash_24, sf_cos_cousub))>0, ]

9.5 Viz: Grid Point Pattern

9.5.1 Grids

# create a grid cover the entire city by defining the cell size (each grid will be 2000 ft * 2000 ft * )
sf_grid_cos <- st_make_grid(sf_cos_cousub, cellsize =  2000, square = TRUE ) %>% st_as_sf() 

# visualize the grid 
tm_shape(sf_grid_cos) + tm_polygons() + 
  tm_shape(sf_cos_cousub) + tm_polygons(alpha = 0, border.col = "red")

# remove grid outside of COS boundary 
sf_grid_cos <- sf_grid_cos[lengths(st_intersects(sf_grid_cos, sf_cos_cousub)) > 0, ]

# view 
tm_shape(sf_grid_cos) + tm_polygons() + 
  tm_shape(sf_cos_cousub) + tm_polygons(alpha = 0, border.col = "red")

# count crashes in each grid 
sf_grid_cos$n_crash <- lengths(st_intersects(sf_grid_cos, sf_crash_24_cos))

# view 
tm_shape(sf_grid_cos) + tm_polygons(col = "n_crash", border.alpha = 0.2) + 
  tm_shape(sf_cos_cousub) + tm_polygons(alpha = 0, border.col = "red") 

9.5.2 Roads

# projection 
sf_road <- st_transform(sf_road, 2231)

# extract main roads
# MTFCC code: https://www2.census.gov/geo/pdfs/reference/mtfccs2022.pdf
sf_r_main <- sf_road %>% filter(MTFCC %in% c("S1100", "S1200")) 
sf_r_main <- sf_r_main[lengths(st_intersects(sf_r_main, sf_grid_cos)) > 0, ]

# prepare labels 
sf_r_main_lbl <- sf_r_main[!duplicated(sf_r_main$FULLNAME),]

## Add main roads to the grid point map
tm_shape(sf_grid_cos) + tm_polygons(col = "n_crash", border.alpha = 0.2) + 
  tm_shape(sf_cos_cousub) + tm_polygons(alpha = 0, border.col = "red") + 
  tm_shape(sf_r_main) + tm_lines() + 
  tm_shape(sf_r_main_lbl) + tm_text(text = "FULLNAME", size = 0.7, remove.overlap = T)

9.5.3 Crash Points by Vehicle Types

# vehicle types 
lt_v_type <- unique(sf_crash_24_cos$vehicletyp) 

# remove NA
lt_v_type <- lt_v_type[!is.na(lt_v_type)]

# lbl 
lt_v_type_lbl <- c("auto", "motor", "truck", "bus", "bike", "trailor")

print(lt_v_type)

## [1] "Automobile"           "Motor, Trike etc...." "Trucks"              
## [4] "Bus"                  "Bicycle"              "Semi-Truck/Trailor"

print(lt_v_type_lbl)

## [1] "auto"    "motor"   "truck"   "bus"     "bike"    "trailor"

# county the number of crash by vehicle types in each grid  
for (i in 1:length(lt_v_type)){
  print(i)
  print(lt_v_type[[i]])
  
  data <- sf_crash_24_cos %>% filter(vehicletyp == lt_v_type[[i]] )
  sf_grid_cos[as.character(lt_v_type_lbl[[i]])] <- lengths(st_intersects(sf_grid_cos, 
                                                                         data))
}

## [1] 1
## [1] "Automobile"
## [1] 2
## [1] "Motor, Trike etc...."
## [1] 3
## [1] "Trucks"
## [1] 4
## [1] "Bus"
## [1] 5
## [1] "Bicycle"
## [1] 6
## [1] "Semi-Truck/Trailor"

# check the output data 
head(sf_grid_cos, 10)

## Simple feature collection with 10 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 3211006 ymin: 779273.4 xmax: 3221006 ymax: 787273.4
## Projected CRS: NAD83 / Colorado North (ftUS)
##                                 x n_crash auto motor truck bus bike trailor
## 1  POLYGON ((3217006 779273.4,...       0    0     0     0   0    0       0
## 2  POLYGON ((3219006 779273.4,...       0    0     0     0   0    0       0
## 3  POLYGON ((3215006 781273.4,...       0    0     0     0   0    0       0
## 4  POLYGON ((3217006 781273.4,...       0    0     0     0   0    0       0
## 5  POLYGON ((3219006 781273.4,...       0    0     0     0   0    0       0
## 6  POLYGON ((3213006 783273.4,...       0    0     0     0   0    0       0
## 7  POLYGON ((3215006 783273.4,...       0    0     0     0   0    0       0
## 8  POLYGON ((3217006 783273.4,...       0    0     0     0   0    0       0
## 9  POLYGON ((3219006 783273.4,...       0    0     0     0   0    0       0
## 10 POLYGON ((3211006 785273.4,...       0    0     0     0   0    0       0

Map crash involving motors.

tmap_mode("plot")

tm_shape(sf_grid_cos) + tm_polygons(col = "motor", border.alpha = 0.2, palette = "Blues") + 
  tm_shape(sf_cos_cousub) + tm_polygons(alpha = 0, border.col = "red") + 
  tm_shape(sf_r_main) + tm_lines() + 
  tm_shape(sf_r_main_lbl) + tm_text(text = "FULLNAME", size = 0.5, remove.overlap = T) + 
  tm_compass(north = 0, size = 1, position = c("right", "top")) + 
  tm_scale_bar(position = c("right", "bottom"), width = 0.15)  + 
  tm_layout(legend.outside = T)

Q1

Please adapt the code chunk above to produce the following maps, and paste the screenshots to the report.

• crashes involving automobiles using the color palette “PuBuGn”
• crashes involving trailers using the color palette “YlOrRd”
• crashes involving buses using the color palette “Greens”

9.6 Viz: Kernel Density Map

Hint You would need to re-write this section (9.6) to answer Q2.

Q2: Please produce a kernel density map for crash involving trucks, and paste the screenshot to the report.

9.6.1 Data Preparation

# data
# spatial polygon of city boundary: sf_cos_cousub 
# spatial point of crash data: sf_crash_24_cos

## check different vehicle types 
print(lt_v_type)

## [1] "Automobile"           "Motor, Trike etc...." "Trucks"              
## [4] "Bus"                  "Bicycle"              "Semi-Truck/Trailor"

## filter crash of automobile 
data <- sf_crash_24_cos %>% filter(vehicletyp == "Automobile")

# check the class of "data"
class(data)

## [1] "sf"         "data.frame"

# convert from simple features to spatial point 
sp_data <- as(data, "Spatial")

# check class "sp_data"
class(sp_data)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

# coordinates top 10 rows 
head(sp_data@coords, 10)

##       coords.x1 coords.x2
##  [1,]   3193120  811891.1
##  [2,]   3221835  840870.2
##  [3,]   3198922  859368.4
##  [4,]   3193682  830805.8
##  [5,]   3192054  833678.2
##  [6,]   3196537  833468.2
##  [7,]   3194730  833535.2
##  [8,]   3191926  825034.4
##  [9,]   3214185  832822.7
## [10,]   3222087  834037.1

9.6.2 KDE Mapping

# create observation window (owin)
wndw_cos <- as.owin(sf_cos_cousub)
plot(wndw_cos)

#create a ppp (point pattern) object
ppp_crash_auto_24_cos <- ppp(x=sp_data@coords[,1],
                        y=sp_data@coords[,2],
                        window=wndw_cos)

# kernel density map
ppp_crash_auto_24_cos %>%
  density(., sigma=500) %>%
  plot()

# play with sigma 
ppp_crash_auto_24_cos %>%
  density(., sigma=2000) %>%
  plot()

9.7 Close & Exit

Congratulations!! You have completed the entire tutorial and learnt the intro to network visualization!!

Please submit your report on time, otherwise, a late penalty (5 pts / day) will be applied.

Please go “File”–> “Close Project” – a pop window asking “Do you want to save these changes” –> “Yes”.

Don’t forget to submit the lab9 report and your script to Canvas.