7.1 Lab Goals

This chapter explores mapping of point patterns:

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

7.2 Good Practice

7.2.1 Organizing Folders & Sub-folders

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

7.2.2 Data

This chapter uses:

Download the lab’s data:

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

7.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 lab7, 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.

7.3 Libraries & Data

7.3.1 Load Libraries

Please install necessary packages as needed.

library(dplyr)      ## mutate, left_join
library(sf)         ## simple feature manipulation and processing
library(tmap)       ## mapping
library(spatstat)   ## kernel density map

7.3.2 Data

Simple Features (sf) - spatial data

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

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

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

7.4 Data Processing

7.4.1 Spatial Boundaries

Check the coordinate reference system of “sf_co_cousub” - county subdivisions. sf_co_cousub is in a Geodetic Coordinate Reference System “NAD83”.

# check coordinate reference systems
sf_co_cousub$geometry
## Geometry set for 209 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -109.0603 ymin: 36.99243 xmax: -102.0415 ymax: 41.00344
## Geodetic CRS:  NAD83
## First 5 geometries:

the following line project the Geodetic CRS to a planar CRS “NAD83 / Colorado North (ftUS)”. CRS code: 2231

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

Filter county subdivision based on county & city names

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

7.4.2 Crash Data

Examine the crash data.

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

The crash data has some records with empty geometry, e.g., row 10. Let’s have a look at an example ‘empty point’ in the 10th row.

sf_crash[10, ]$geometry

Data filtering: We can use the function st_is_empty() identifies empty geometry.

sf_crash[st_is_empty(sf_crash),]

The logical NOT operator ! can be used to reverse the condition to extract NON-EMPTY geometry

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

7.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(fill_alpha = 0, lwd = 3, col = "blue") + 
  tm_shape(sf_crash_24) + tm_dots(fill_alpha  = 0.5)

The above visualization has some issues with the crash data, e.g.,

    1. many dots overlap with each other;
    1. there are crash dots outside of cos city boundary.

Let’s extract only crash points that fall within the COS city boundary by using the function st_intersects() and lengths().

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

7.5 Viz: Grid Point Pattern

As we have seen in the Exploratory Analysis, many points data overlaping with each other make it hard to analyze. This section explores ways to aggregate points into a higher spatial level (e.g., using grid).

7.5.1 Grids

First, let’s create a grid cover the entire city by defining the cell size (each grid will be 2000 ft * 2000 ft * ).

sf_grid_cos(): create a square or hexagonal grid covering the bounding box of the geometry of an sf or sfc object. See documentation.

# 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(fill_alpha = 0, col = "red")

There are some grids are outside of the COS boundary. Let’s remove them.

# 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(fill_alpha = 0, col = "red")

Count the number of crashes in each grid and save the number as a new column n_crash into the table sf_grid_cos.

# 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(fill = "n_crash", lwd = 0.1, col="grey") + 
  tm_shape(sf_cos_cousub) + tm_polygons(fill_alpha = 0, col = "red")

7.5.2 Roads

Remember, we have also loaded the road data at the beginning of the tutorial. Let’s checking on it.

# check coordinate reference systems
sf_road$geometry
## Geometry set for 29689 features 
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -105.0685 ymin: 38.51945 xmax: -104.0516 ymax: 39.13004
## Geodetic CRS:  NAD83
## First 5 geometries:
# projection 
sf_road <- st_transform(sf_road, 2231)

# have a look at the data 
head(sf_road)
## Simple feature collection with 6 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 3183686 ymin: 838568.6 xmax: 3204586 ymax: 872959.3
## Projected CRS: NAD83 / Colorado North (ftUS)
##        LINEARID                      FULLNAME RTTYP MTFCC
## 1 1105633035042           N Union Blvd on Rmp     M S1400
## 2 1105633035229           N Union Blvd on Rmp     M S1400
## 3 1105633034991          N Union Blvd Off Rmp     M S1400
## 4 1105633035073 Austin Bluffs Parkway Off Rmp     M S1400
## 5 1105633035301  Austin Bluffs Parkway on Rmp     M S1400
## 6 1105602055298                     Sumac Exd     M S1400
##                         geometry
## 1 LINESTRING (3203796 839097....
## 2 LINESTRING (3203358 838669,...
## 3 LINESTRING (3203102 839060....
## 4 LINESTRING (3203438 839023....
## 5 LINESTRING (3203384 839022....
## 6 LINESTRING (3183844 871602....

The MTFCC code denotes different road categories, we can use it to extract main roads.

# extract main roads
sf_r_main <- sf_road %>% filter(MTFCC %in% c("S1100", "S1200")) 

# main roads that are overlapping with COS city boundary 
sf_r_main <- sf_r_main[lengths(st_intersects(sf_r_main, sf_cos_cousub)) > 0, ]

Add a layer of road to the grid map.

## Add main roads to the grid point map
tm_shape(sf_grid_cos) + tm_polygons(fill = "n_crash", lwd = 0.1, col="grey") + 
  tm_shape(sf_cos_cousub) + tm_polygons(fill_alpha = 0, col = "red") + 
  tm_shape(sf_r_main) + tm_lines()

7.5.3 Crash Points by Vehicle Types

It would be interesting to plot crashes by vehicle types. Let’s check unique vehicle types involved in crashes and save it as a new list lt_v_type.

# vehicle types 
lt_v_type <- unique(sf_crash_24_cos$vehicletyp) 
print(lt_v_type)
## [1] "Automobile"           NA                     "Motor, Trike etc...."
## [4] "Trucks"               "Bus"                  "Bicycle"             
## [7] "Semi-Truck/Trailor"
# remove NA
lt_v_type <- lt_v_type[!is.na(lt_v_type)]
print(lt_v_type)
## [1] "Automobile"           "Motor, Trike etc...." "Trucks"              
## [4] "Bus"                  "Bicycle"              "Semi-Truck/Trailor"

Give those vehicle types a new label.

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

create a “for” loop to count different crashes by vehicle types in “lt_v_type_lbl”

# 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.

# Static mapping 
tmap_mode("plot")

# motor
tm_shape(sf_grid_cos) + tm_polygons(fill = "motor", col_alpha = 0.2, 
                                    fill.scale = tm_scale_intervals(values = "brewer.blues")) + 
  tm_shape(sf_cos_cousub) + tm_polygons(fill_alpha = 0, col = "red") + 
  tm_shape(sf_r_main) + tm_lines() + 
  tm_compass(north = 0, size = 1, position = c("right", "top")) + 
  tm_scalebar(position = c("right", "top"), width = 5)  + 
  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 “brewer.greens”
  • crashes involving trailers using the color palette “brewer.purples”
  • crashes involving buses using the color palette “brewer.yl_gn”

7.6 Viz: Kernel Density Map

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

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

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

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

7.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 lab7 report and your script to Canvas.