# Mapping of census crime and safety data to Baltimore neighborhoods

# libraries
library(ggplot2)
library(ggthemes)
library(dplyr)
library(reshape2)
library(magrittr)
library(tidyr)
library(viridis)
library(scales)

# Latitudinal and longitudinal coordinates for neighborhods scraped from 
# Crime and Safety (2010-2013) KML file 
# (original file: https://data.baltimorecity.gov/Neighborhoods/Crime-Safety-2010-2013-Shape/y3aa-iwrk).
# Square region without data is designated as "Unassigned -- Jail" in KML file.

# Note: fill argument for read.table determines row length from only first five 
# rows of data, thereby mistakenly slicing longer rows later in file; moved the 
# three longer rows being sliced to the beginning of the file as a work-around

census.map <- read.table("census_basemaps/Bmore_neighborhoods_from_Crime_and_Safety_kml.txt", 
                         header=FALSE, 
                         sep=" ", 
                         quote="\"", 
                         fill=TRUE, 
                         row.names=1) %>% 
                t() %>% 
                data.frame(check.names=FALSE) %>% 
                gather(key=Neighborhood, value=coords)

# remove empty coords generated by fill argument
census.map <- census.map[census.map$coords != "", ]

# separte coords into latitute and longitude and convert to numeric
census.map %<>% separate(col=coords, into=c("longitude", "latitude"), sep=",")
census.map[[2]] %<>% as.numeric()
census.map[[3]] %<>% as.numeric()

# Create column with unique IDs used for mapping.  Mapping without unique IDs 
# results in erroneous connecting lines, or "tears," for neighborhoods composed 
# of  multiple polygons.  
census.map$unique.mapping <- census.map$Neighborhood

# Format neighborhood IDs to match those in the census in order to merge data.
census.map$Neighborhood <- gsub(pattern="\\*", 
                                replacement=" ", 
                                x=census.map$Neighborhood)

for (i in 1:3){
  rep1 <- c("Brooklyn/Curtis Bay/Hawkins Point1", 
            "Brooklyn/Curtis Bay/Hawkins Point2", 
            "Brooklyn/Curtis Bay/Hawkins Point3")
  
  census.map$Neighborhood %<>% gsub(pattern=rep1[i], 
                                    replacement="Brooklyn/Curtis Bay/Hawkins Point", .)
}

for (i in 1:2){
  rep2 <- c("South Baltimore1", 
            "South Baltimore2")
  
  census.map$Neighborhood %<>% gsub(pattern=rep2[i], 
                                    replacement="South Baltimore", .)
}

# read in census data 
crime <- read.csv("Crime___Safety__2010-2013_.csv")

### Definitions (all per 1000 residents)
# violent: homicide, rape, aggravated assault, and robbery
# prop: burglary and auto theft

### dot plot of average violent and property crime from 
### 2010-2013 and 2011-2013, respectively.
crime.dot <- crime %>% select(Neighborhood, 
                              prop11, prop12, prop13, 
                              viol10, viol11, viol12, viol13) %>% 
  mutate(prop.av=(prop11 + prop12 + prop13)/3) %>%
  mutate(viol.av=(viol10 + viol11 + viol12 + viol13)/4) %>%
  select(Neighborhood, prop.av, viol.av) %>% 
  melt(value.name="freq") 

# for appropriate order and labeling of facets
levels(crime.dot$variable) <- c("Property Crime", "Violent Crime")
crime.dot$variable <- relevel(crime.dot$variable, "Violent Crime")

# reorder neighborhoods so plotted in descending alphabetical order
crime.dot$Neighborhood <- factor(crime.dot$Neighborhood, levels(crime.dot$Neighborhood)[rev(crime.dot$Neighborhood)])

# remove city-wise statistic ("Baltimore City"), leaving only neighborhood statistics
crime.dot <- crime.dot[crime.dot$Neighborhood != "Baltimore City", ]

# facetted plot
crime.dotplot <- ggplot(crime.dot, 
                        aes(x = freq, y = Neighborhood, color = variable)) +
                  geom_point(size = 12) + 
                  scale_color_manual(values = c("grey30", "#2D2D83")) +
                  facet_grid( . ~ variable, scales = "free") +
                  labs(y = NULL, x = "Incidences per 1000 residents") +
                  theme_bw() + 
                  theme(legend.position="none", 
                        axis.text.y=element_text(size=25), 
                        axis.text.x=element_text(size=25), 
                        axis.title.x=element_text(size=25),
                        strip.text.x = element_text(size = 25))
# plot
png("baltimore_census_crime_by_neighborhood_dotplot.png", height = 2000, width = 1152)
crime.dotplot
dev.off()

### Violent crime for all neighborhoods
viol <- crime %>% 
  select(Neighborhood, viol10, viol11, viol12, viol13) %>% 
  mutate(viol.av=(viol10 + viol11 + viol12 + viol13)/4) %>%
  select(Neighborhood, viol.av) %>% 
  melt(value.name="viol") %>% 
  right_join(census.map, by = "Neighborhood") %>%
  ggplot(mapping = aes(x = longitude, y = latitude, group = unique.mapping)) + 
  geom_polygon(aes(fill = viol), color = "white") + 
  theme_map() + coord_map("mercator") + 
  scale_fill_gradientn(colours = viridis(256), 
                       limits = c(0, 82), 
                       breaks = seq(from = 0, to = 80, by = 20), 
                       guide = guide_colorbar(barwidth = 1, 
                                              barheight = 3.7, 
                                              title = NULL, 
                                              label.position = "left", 
                                              label.theme = element_text(size = 12, 
                                                                         angle = 0))) +
  labs(title = "Violent crimes per 1000 residents \n (2010-2013 average)") + 
  theme(plot.title = element_text(size=12))


# png
png("baltimore_census_violent_crime_2010_to_2013.png", height = 1152/4, width = 1152/4)
viol
dev.off()

### Property crimes for all neighborhoods
prop <- crime %>% 
  select(Neighborhood, prop11, prop12, prop13) %>% 
  mutate(prop.av=(prop11 + prop12 + prop13)/3) %>%
  select(Neighborhood, prop.av) %>% 
  melt(value.name="prop") %>% 
  right_join(census.map, by = "Neighborhood") %>%
  ggplot(mapping = aes(x = longitude, y = latitude, group = unique.mapping)) + 
  geom_polygon(aes(fill = prop), color = "white") + 
  theme_map() + 
  coord_map("mercator") + 
  scale_fill_gradientn(colours = viridis(256), 
                       limits = c(0, 271), 
                       breaks = seq(from = 0, to = 250, by = 50), 
                       guide = guide_colorbar(barwidth = 1, 
                                              barheight = 3.8, 
                                              title = NULL, 
                                              label.position = "left", 
                                              label.theme = element_text(size = 11, 
                                                                         angle = 0))) +
  labs(title = "Property crimes per 1000 residents \n (2011-2013 average)") + 
  theme(plot.title = element_text(size=12))

# png
png("baltimore_census_property_crime_2011_to_2013.png", height = 1152/4, width = 1152/4)
prop
dev.off()

### Violent crime excluding Downtown/Seton Hill
viol.wo.dt <- crime %>% 
  select(Neighborhood, viol10, viol11, viol12, viol13) %>% 
  filter(Neighborhood != "Downtown/Seton Hill") %>%
  mutate(viol.av=(viol10 + viol11 + viol12 + viol13)/4) %>%
  select(Neighborhood, viol.av) %>% 
  melt(value.name="viol") %>% 
  right_join(census.map, by = "Neighborhood") %>%
  ggplot(mapping = aes(x = longitude, y = latitude, group = unique.mapping)) + 
  geom_polygon(aes(fill = viol), color = "white") + 
  theme_map() + 
  coord_map("mercator") + 
  scale_fill_gradientn(colours = viridis(256), 
                       limits = c(0, 30), 
                       breaks = seq(from = 0, to = 30, by = 5), 
                       guide = guide_colorbar(barwidth = 5, 
                                              barheight = 23, 
                                              title = NULL, 
                                              label.position = "left", 
                                              label.theme = element_text(size = 30, 
                                                                         angle = 0))) +
  labs(title = "Violent crimes per 1000 residents \n (2010-2013 average)") + 
  theme(plot.title = element_text(size = 40))

# png
png("baltimore_census_violent_crime_wo_dt_2010_to_2013.png", height = 1152, width = 1152)
viol.wo.dt
dev.off()

### Property crime excluding Downtown/Seton Hill
prop.wo.dt <- crime %>% 
  select(Neighborhood, prop11, prop12, prop13) %>% 
  filter(Neighborhood != "Downtown/Seton Hill") %>%
  mutate(prop.av=(prop11 + prop12 + prop13)/3) %>%
  select(Neighborhood, prop.av) %>% 
  melt(value.name="prop") %>% 
  right_join(census.map, by = "Neighborhood") %>%
  ggplot(mapping = aes(x = longitude, y = latitude, group = unique.mapping)) + 
  geom_polygon(aes(fill = prop), color = "white") + 
  theme_map() + 
  coord_map("mercator") + 
  scale_fill_gradientn(colours = viridis(256), 
                       limits = c(0, 125), 
                       breaks = seq(from = 0, to = 125, by = 25),
                       guide = guide_colorbar(barwidth = 5, 
                                              barheight = 23, 
                                              title = NULL, 
                                              label.position = "left", 
                                              label.theme = element_text(size = 30, 
                                                                         angle = 0))) +
  labs(title = "Property crimes per 1000 residents \n (2011-2013 average)") + 
  theme(plot.title = element_text(size = 40))

# png
png("baltimore_census_property_crime_wo_dt_2011_to_2013.png", height = 1152, width = 1152)
prop.wo.dt
dev.off()

### Homicide
gunhom <- crime %>% select(Neighborhood, gunhom11, gunhom12, gunhom13) %>% 
  mutate(gunhom.av=(gunhom11+gunhom12+gunhom13)/3) %>%
  select(Neighborhood, gunhom.av) %>% 
  melt(value.name="gunhom") %>% 
  right_join(census.map, by = "Neighborhood") %>%
  ggplot(mapping = aes(x = longitude, y = latitude, group = unique.mapping)) + 
  geom_polygon(aes(fill = gunhom), color = "white") + 
  theme_map() + 
  coord_map("mercator") + 
  scale_fill_gradientn(colours = viridis(256), 
                       limits = c(0, 1), 
                       breaks = seq(from = 0, to = 1.0, by = 0.2),
                       guide = guide_colorbar(barwidth = 5, 
                                              barheight = 23, 
                                              title = NULL, 
                                              label.position = "left", 
                                              label.theme = element_text(size = 30, 
                                                                         angle = 0))) +
  labs(title = "Homicides by firearm per 1000 residents \n (2011-2013 average)") + 
  theme(plot.title = element_text(size=40))


# png
png("baltimore_census_average_gun_homicide_2011_to_2013.png", height = 1152, width = 1152)
gunhom
dev.off()