Introduction

Original Proposal: This past year, there have been a rise of thefts and burglaries in Washington, especially auto thefts. According to a Fox13 Seattle News article, auto thefts spiked 88% in 2022 compared to 2023 (Fox 13 News Staff, 2022). For the past months, there have been reports made of auto thefts on campus as well. As such, I plan to a point pattern analysis of burglaries and auto thefts within Thurston County.

Actual Project: Due to difficulties of finding crime data in Thurston County, I was able to find one in Seattle. Therefore, this project is a point pattern analysis of robberies within the city of Seattle .

Import the necessary libraries

library(sp) #handles spatial data
library(janitor) #data cleaning and format tasks

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(sf) # provides classes and functions for spatial operations

## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(ggplot2) # creating visualizations and graphics 
library(dplyr) #data manipulation,transformation, and data wrangling

## 
## 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(spatstat) #analyzes spatial point patterns

## Loading required package: spatstat.data

## Loading required package: spatstat.geom

## spatstat.geom 3.2-1

## Loading required package: spatstat.random

## spatstat.random 3.1-4

## Loading required package: spatstat.explore

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## spatstat.explore 3.1-0

## Loading required package: spatstat.model

## Loading required package: rpart

## spatstat.model 3.2-3

## Loading required package: spatstat.linnet

## spatstat.linnet 3.1-0

## 
## spatstat 3.0-5 
## For an introduction to spatstat, type 'beginner'

library(raster)# provides functions for reading, writing, and manipulating raster files

## 
## Attaching package: 'raster'

## The following object is masked from 'package:nlme':
## 
##     getData

## The following object is masked from 'package:dplyr':
## 
##     select

library(leaflet)# creates interactive maps and visualizations using leaflet.js
library(rgdal)#provides bindings to the GDAL (Geospatial Data Abstraction Library) for reading and writing geospatial data formats

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

library(tidyverse) 

## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──

## ✔ tibble  3.1.8     ✔ purrr   1.0.1
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.3     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ nlme::collapse() masks dplyr::collapse()
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ raster::select() masks dplyr::select()

Read a shape file of Seattle

city <- st_read("~/Downloads/WSDOT_-_City_Limits.geojson")

## Reading layer `WSDOT_-_City_Limits' from data source 
##   `/Users/tedyheaven-litabarber/Downloads/WSDOT_-_City_Limits.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 281 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.4184 ymin: 45.55897 xmax: -117.0234 ymax: 49.00237
## Geodetic CRS:  WGS 84

Create a new object that only has the city of Seattle and Cchange the coordinate system.

df1 <- city %>%
  filter(CityName == "Seattle")
seattle <- st_transform(df1, 4326)
#Get rid of objects we no longer need
rm(city)
rm(df1)

Read Seattle’s police data from Seattle’s Open Data and filter by robbery. The data downloaded was already filter to have data from January 1, 2023 - Present

crimes <- read_csv("~/Downloads/SPD_Crime_Data__2008-Present (2).csv") %>% clean_names()

## Rows: 12853 Columns: 17
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (14): Report Number, Offense Start DateTime, Offense End DateTime, Repor...
## dbl  (3): Offense ID, Longitude, Latitude
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

robbery <- filter(crimes, offense_parent_group == "ROBBERY")

Transform the data frame with crime information into a sf object, storing our other spacial objects.

robbery_spatial <-  st_as_sf(robbery, coords = c("longitude", "latitude"), 
                 crs = 4326, agr = "constant")

Select only robberies that take place within the space defined by the boundaries by intersection

rob_seattle <- st_intersects(seattle, robbery_spatial)

# Subsetting
rob_seattle <- robbery_spatial[unlist(rob_seattle),]

# Again, we remove things we don't need
rm(crimes)
rm(robbery)

Now that the data is cleaned and our files are prepared. Let’s see the results.

ggplot() + 
  geom_sf(data = seattle) + 
  geom_sf(data = rob_seattle) + 
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_blank()) 

In our spatial object one of these marks is the type of crime (although in this case it’s of little interest since we have filtered on it).

Using Spatstat

Transform the CRS of our Seattle polygon into projected coordinates

seattle_proj <- st_transform(seattle, 2285)

Define the window

window <- as.owin(seattle_proj)
class(window)

## [1] "owin"

Extract the coordinates from our sf point data into a matrix

rob_seattle <- st_transform(rob_seattle, 2285) #we must transform these too in order to match our window
sf_rob_sea_coords <- matrix(unlist(rob_seattle$geometry), ncol = 2, byrow = T)

Create the object using the information from our matrix and the window that we created.

rob_ppp <- ppp(x = sf_rob_sea_coords[,1], y = sf_rob_sea_coords[,2],
                   window = window, check = T)

## Warning: data contain duplicated points

plot(rob_ppp)

“When the data has coincidence points, some statistical procedures will be severely affected. So it is always strongly advisable to check for duplicate points and to decide on a strategy for dealing with them if they are present” - Baddeley et al., 2016: 60

any(duplicated(rob_ppp))

## [1] TRUE

Count the number of coincidence points

multiplicity(rob_ppp)

##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   2   1   1   1   1   2   1   2   4   5   1   1   1   1   1   1   1   1   1   4 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   1   3   1   1   5   1   1   1   2   1   1   7   3   7   1   2   1   1   1   1 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   1   1   2   1   1   1   3   4   1   1   1   2   1   1   1   2   4   2   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   2   2   1   1   1   1   1   5   5   1   2   1   1   1   1   3   1   1   4   1 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   1   2   1   1   1   1   1   1   1   1   1   4   1   1   2   1   4   1   1   1 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   1   1   1   1   1   1   7   1   3   1   2   1   1   1   3   1   4   1   1   1 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   1   7   1   1   2   2   1   3   1   1   1   1   2   1   1   1   1   1   1   2 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   1   5   1   1   1   1   2   1   1   1   1   1   2   1   4   1   5   1   1   1 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   2   5   1   1   1   3   1   2   1   1   1   1   7   2   1   4   1   2   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   1   1   1   1   3   1   2   4   1   2   7   2   1   1   1   1   1   1   1   1 
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 
##   2   1   2   7   4   1   1   1   1   2   1   1   2   1   1   3   3   1   1   4 
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 
##   4   5   2   1   1   1   1   2   1   2   1   2   1   2   1   1   1   1   1   1 
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 
##   1   1   1   1   1   4   1   1   1   1   1   1   1   1   5   1   2   1   1   2 
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 
##   1   7   1   2   1   1   2   1   1   1   1   1   1   1   1   1   1   3   2   1 
## 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 
##   1   1   1   2   7   1   1   1   1   2   1   1   2   1   1   2   4   4   1   1 
## 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 
##   3   1   1   1   1   1   1   1   2   1   1   1   3   7   4   2   1   1   1   1 
## 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 
##   1   4   1   1   2   7   3   1   1   1   1   1   7   1   2   2   5   1   2   1 
## 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 
##   1   1   1   1   2   1   7   1   1   1   1   1   1   1   1   1   3   1   1   1 
## 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 
##   1   3   1   1   7   1   3   1   1   1   1   1   2   4   1   1

How many locations have more than one event?

sum(multiplicity(rob_ppp) > 1)

## [1] 116

ggplot() + 
  geom_sf(data = seattle) + 
  geom_sf(data = rob_seattle, alpha = 0.4) + 
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_blank()) 

Jitter

It will give the duplicates a slight displacement so they do not occupy the same space.

jitter_rob <- rjitter(rob_ppp, retry=TRUE, nsim=1, drop=TRUE)
plot(jitter_rob)

One could divide the window of observation into quadrants and count the number of points into each of these quadrants.

Q <- quadratcount(jitter_rob, nx = 7, ny = 10)
plot(jitter_rob)
plot(Q, add = TRUE, cex =.2, col = "red")

Kernel Density Estimates

Generate a smooth continuous surface aiming to represent the density or volume of crimes across the target area.

ds <- density(jitter_rob)
class(ds)

## [1] "im"

plot(ds, main='Robbery density in Seattle')

Define the bandwidth of the density estimation to select an appropriate bandwidth.

bw.diggle(jitter_rob)

##    sigma 
## 198.2672

bw.ppl(jitter_rob)

##    sigma 
## 1384.009

bw.scott(jitter_rob)

##  sigma.x  sigma.y 
## 2753.166 7504.833

par(mfrow=c(2,2))
plot(density.ppp(jitter_rob, sigma = bw.diggle(jitter_rob),edge=T),
     main = paste("diggle"))

plot(density.ppp(jitter_rob, sigma = bw.ppl(jitter_rob),edge=T),
     main=paste("likelihood cross-validation"))

plot(density.ppp(jitter_rob, sigma = bw.scott(jitter_rob)[2],edge=T),
     main=paste("scott 1"))

plot(density.ppp(jitter_rob, sigma = bw.scott(jitter_rob)[1],edge=T),
     main=paste("scott 2"))

par(mfrow=c(2,2))
plot(density.ppp(jitter_rob, sigma = bw.ppl(jitter_rob),edge=T),
     main=paste("Gaussian"))
plot(density.ppp(jitter_rob, kernel = "epanechnikov", sigma = bw.ppl(jitter_rob),edge=T),
     main=paste("Epanechnikov"))
plot(density.ppp(jitter_rob, kernel = "quartic", sigma = bw.ppl(jitter_rob),edge=T),
     main=paste("Quartic"))
plot(density.ppp(jitter_rob, kernel = "disc", sigma = bw.ppl(jitter_rob),edge=T),
     main=paste("Disc"))

K-Cross Function

Are there any clustering of two sets of patterns?

#A ppp class is created with marks by using as.factor() around the column of interest.
mpp <- ppp(sf_rob_sea_coords[,1], sf_rob_sea_coords[,2], window = window, marks=as.factor(rob_seattle$precinct))

## Warning: data contain duplicated points

#Using the envelope() function, we will examine clustering between south precinct and east precinct point patterns.
ekc <- envelope(mpp, Kcross, nsim = 5, i = "S", j = "E")

## Generating 5 simulations of CSR  ...
## 1, 2, 3, 4,  5.
## 
## Done.

plot(ekc, main = "S and E")

Interpretation

The dashed red line in our KCross plot represents the reference line for complete spatial randomness between the two point patterns.

The grey line around the red line is the randomization envelope.

The black line represent the Kcross function for these two point patterns.

So our black line is quite a bit higher than the red dashed line. This gives us evidence that there is clustering of the two point patterns in this example.