require(scales)
require(dplyr)
require(gtools)
require(ggplot2)
require(rgdal)
require(ggmap)
require(Cairo)
require(gpclib)
require(maptools)
require(reshape)
library(devtools)
library(stringr)
library(raster)
library(sp)
library(lubridate)
## install_github("trendct/ctnamecleaner")
library(ctnamecleaner)
library(leaflet)
library(tidyr)
## Quickly mapping the locations
mega <- read.csv("mega.csv", stringsAsFactors=FALSE)
mega_no_na <- mega[!is.na(mega$InterventionLocationLatitude),]
m <- leaflet(mega_no_na) %>% addTiles('http://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png',
attribution='Map tiles by <a href="http://stamen.com">Stamen Design</a>, <a href="http://creativecommons.org/licenses/by/3.0">CC BY 3.0</a> — Map data © <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>')
m %>% setView(-72.690940, 41.651426, zoom = 8)
m %>% addCircles(~InterventionLocationLongitude, ~InterventionLocationLatitude, popup=mega$DepartmentName, weight = 3, radius=40,
color="#ffa500", stroke = TRUE, fillOpacity = 0.8)
## Stops within town
gpclibPermit()
## [1] TRUE
gpclibPermitStatus()
## [1] TRUE
towntracts40 <- readOGR(dsn="maps", layer="buffer60")
## OGR data source with driver: ESRI Shapefile
## Source: "maps", layer: "buffer60"
## with 169 features
## It has 20 fields
towntracts_only40 <- towntracts40
towntracts40 <- fortify(towntracts40, region="NAME10")
coords <- mega_no_na[c("InterventionLocationLongitude", "InterventionLocationLatitude")]
coords <- coords[complete.cases(coords),]
sp <- SpatialPoints(coords)
proj4string(sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(towntracts_only40)
plot(sp, col="red" , add=TRUE)
res40 <- over(sp, towntracts_only40)
sixty <- res40 %>%
group_by(NAME10) %>%
summarise(sixty=n())
sixty <- sixty[!is.na(sixty$NAME10),]
colnames(sixty) <- c("id", "sixty")
## Stops by town border
gpclibPermit()
## [1] TRUE
gpclibPermitStatus()
## [1] TRUE
towntracts60 <- readOGR(dsn="maps", layer="cut60")
## OGR data source with driver: ESRI Shapefile
## Source: "maps", layer: "cut60"
## with 169 features
## It has 20 fields
towntracts_only60 <- towntracts60
towntracts60 <- fortify(towntracts60, region="NAME10")
coords <- mega_no_na[c("InterventionLocationLongitude", "InterventionLocationLatitude")]
coords <- coords[complete.cases(coords),]
sp <- SpatialPoints(coords)
proj4string(sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(towntracts_only60)
plot(sp, col="red" , add=TRUE)
res60 <- over(sp, towntracts_only60)
forty <- res60 %>%
group_by(NAME10) %>%
summarise(forty=n())
forty <- forty[!is.na(forty$NAME10),]
colnames(forty) <- c("id", "forty")
forty_map <- left_join(towntracts60, forty)
sixty_map <- left_join(towntracts40, sixty)
## Choropleth totals
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur", fill="")
dtm12
## Choropleths by percent
town_total <- mega_no_na %>%
group_by(DepartmentName) %>%
summarise(total=n())
town_total$DepartmentName <- gsub(" Town", "", town_total$DepartmentName)
sixty2 <- sixty
colnames(sixty2) <- c("DepartmentName", "sixty")
sixty2$DepartmentName <- as.character(sixty2$DepartmentName)
town_total <- left_join(town_total, sixty2)
town_total$forty <- town_total$total - town_total$sixty
town_total$sixty_per <- round((town_total$sixty/town_total$total)*100,2)
town_total$forty_per <- round((town_total$forty/town_total$total)*100,2)
names(town_total)[names(town_total) == 'DepartmentName'] <- 'id'
forty_map <- left_join(towntracts60, town_total)
sixty_map <- left_join(towntracts40, town_total)
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur", fill="")
dtm12
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty_per), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty_per), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (percent by town)", fill="")
dtm12
### White analysis
gpclibPermit()
## [1] TRUE
gpclibPermitStatus()
## [1] TRUE
towntracts40 <- readOGR(dsn="maps", layer="buffer60")
## OGR data source with driver: ESRI Shapefile
## Source: "maps", layer: "buffer60"
## with 169 features
## It has 20 fields
towntracts_only40 <- towntracts40
towntracts40 <- fortify(towntracts40, region="NAME10")
coords_white <- subset(mega_no_na, SubjectRaceCode=="W")
coords_white <- coords_white[c("InterventionLocationLongitude", "InterventionLocationLatitude")]
coords_white <- coords_white[complete.cases(coords_white),]
sp <- SpatialPoints(coords_white)
proj4string(sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(towntracts_only40)
plot(sp, col="red" , add=TRUE)
res40 <- over(sp, towntracts_only40)
sixty <- res40 %>%
group_by(NAME10) %>%
summarise(sixty=n())
sixty <- sixty[!is.na(sixty$NAME10),]
colnames(sixty) <- c("id", "sixty")
## White stops by town border
gpclibPermit()
## [1] TRUE
gpclibPermitStatus()
## [1] TRUE
towntracts60 <- readOGR(dsn="maps", layer="cut60")
## OGR data source with driver: ESRI Shapefile
## Source: "maps", layer: "cut60"
## with 169 features
## It has 20 fields
towntracts_only60 <- towntracts60
towntracts60 <- fortify(towntracts60, region="NAME10")
coords_white <- subset(mega_no_na, SubjectRaceCode=="W")
coords_white <- coords_white[c("InterventionLocationLongitude", "InterventionLocationLatitude")]
coords_white <- coords_white[complete.cases(coords_white),]
sp <- SpatialPoints(coords_white)
proj4string(sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(towntracts_only60)
plot(sp, col="red" , add=TRUE)
res60 <- over(sp, towntracts_only60)
forty <- res60 %>%
group_by(NAME10) %>%
summarise(forty=n())
forty <- forty[!is.na(forty$NAME10),]
colnames(forty) <- c("id", "forty")
forty_map <- left_join(towntracts60, forty)
sixty_map <- left_join(towntracts40, sixty)
## Choropleth totals (White)
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (White drivers)", fill="")
dtm12
## Choropleths by percent (White)
town_total <- mega_no_na %>%
group_by(DepartmentName, SubjectRaceCode) %>%
summarise(total=n()) %>%
spread(SubjectRaceCode, total) %>%
mutate(total=sum(A,B,I,W, na.rm=TRUE))
town_total$DepartmentName <- gsub(" Town", "", town_total$DepartmentName)
sixty2 <- sixty
colnames(sixty2) <- c("DepartmentName", "sixty")
sixty2$DepartmentName <- as.character(sixty2$DepartmentName)
town_total <- left_join(town_total, sixty2)
town_total$forty <- town_total$W - town_total$sixty
town_total$sixty_per <- round((town_total$sixty/town_total$W)*100,2)
town_total$forty_per <- round((town_total$forty/town_total$W)*100,2)
names(town_total)[names(town_total) == 'DepartmentName'] <- 'id'
forty_map <- left_join(towntracts60, town_total)
sixty_map <- left_join(towntracts40, town_total)
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (White drivers)", fill="")
dtm12
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty_per), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty_per), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (White drivers - percent by town)", fill="")
dtm12
### Black analysis
gpclibPermit()
## [1] TRUE
gpclibPermitStatus()
## [1] TRUE
towntracts40 <- readOGR(dsn="maps", layer="buffer60")
## OGR data source with driver: ESRI Shapefile
## Source: "maps", layer: "buffer60"
## with 169 features
## It has 20 fields
towntracts_only40 <- towntracts40
towntracts40 <- fortify(towntracts40, region="NAME10")
coords_black <- subset(mega_no_na, SubjectRaceCode=="B")
coords_black <- coords_black[c("InterventionLocationLongitude", "InterventionLocationLatitude")]
coords_black <- coords_black[complete.cases(coords_black),]
sp <- SpatialPoints(coords_black)
proj4string(sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(towntracts_only40)
plot(sp, col="red" , add=TRUE)
res40 <- over(sp, towntracts_only40)
sixty <- res40 %>%
group_by(NAME10) %>%
summarise(sixty=n())
sixty <- sixty[!is.na(sixty$NAME10),]
colnames(sixty) <- c("id", "sixty")
## Black stops by town border
gpclibPermit()
## [1] TRUE
gpclibPermitStatus()
## [1] TRUE
towntracts60 <- readOGR(dsn="maps", layer="cut60")
## OGR data source with driver: ESRI Shapefile
## Source: "maps", layer: "cut60"
## with 169 features
## It has 20 fields
towntracts_only60 <- towntracts60
towntracts60 <- fortify(towntracts60, region="NAME10")
coords_black <- subset(mega_no_na, SubjectRaceCode=="B")
coords_black <- coords_black[c("InterventionLocationLongitude", "InterventionLocationLatitude")]
coords_black <- coords_black[complete.cases(coords_black),]
sp <- SpatialPoints(coords_black)
proj4string(sp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(sp)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(towntracts_only60)
plot(sp, col="red" , add=TRUE)
res60 <- over(sp, towntracts_only60)
forty <- res60 %>%
group_by(NAME10) %>%
summarise(forty=n())
forty <- forty[!is.na(forty$NAME10),]
colnames(forty) <- c("id", "forty")
forty_map <- left_join(towntracts60, forty)
sixty_map <- left_join(towntracts40, sixty)
## Choropleth totals (Black)
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (Black drivers)", fill="")
dtm12
## Choropleths by percent (Black)
town_total <- mega_no_na %>%
group_by(DepartmentName, SubjectRaceCode) %>%
summarise(total=n()) %>%
spread(SubjectRaceCode, total) %>%
mutate(total=sum(A,B,I,W, na.rm=TRUE))
town_total$DepartmentName <- gsub(" Town", "", town_total$DepartmentName)
sixty2 <- sixty
colnames(sixty2) <- c("DepartmentName", "sixty")
sixty2$DepartmentName <- as.character(sixty2$DepartmentName)
town_total <- left_join(town_total, sixty2)
town_total$forty <- town_total$B - town_total$sixty
town_total$sixty_per <- round((town_total$sixty/town_total$B)*100,2)
town_total$forty_per <- round((town_total$forty/town_total$B)*100,2)
names(town_total)[names(town_total) == 'DepartmentName'] <- 'id'
forty_map <- left_join(towntracts60, town_total)
sixty_map <- left_join(towntracts40, town_total)
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (Black drivers)", fill="")
dtm12
dtm12 <- ggplot() +
geom_polygon(data = forty_map, aes(x=long, y=lat, group=group, fill=forty_per), color = "black", size=0.2) +
geom_polygon(data = sixty_map, aes(x=long, y=lat, group=group, fill=sixty_per), color = "black", size=0.2) +
coord_map() +
scale_fill_distiller(type="seq", trans="reverse", palette = "Reds", breaks=pretty_breaks(n=10)) +
theme_nothing(legend=TRUE) +
labs(title="Where traffic stops occur (Black drivers - percent by town)", fill="")
dtm12