### Seasonal Variation in Spatial Distribution of Burglary Incidents in Philadelphia, 2014-2016, at the street segment level ###
### Yinuo Yin ###
### February 20, 2018 ###

# set working directory
setwd("G:/Yinuo Yin/Study/Upenn/2018 Spring/MUSA 620/Week5 Spatial Queries")

# load libraries
library(RPostgreSQL)
library(sf)
library(postGIStools)
library(tidyverse)
library(viridis)
library(classInt)
library(lubridate)

# connect to database
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname = "spatialdb",
                 host = "127.0.0.1", port = 7437,
                 user = "postgres", password = '321qaz')

# Step 1: Import data + data cleaning
# add Philly shapefile to the database
phillySF <- st_read('census-tracts-philly.shp', stringsAsFactors = FALSE)

# we want lowercase fields
phillySF <- rename(phillySF, gisjoin = GISJOIN)

# transform to Mercator CRS
phillySF <- st_transform(phillySF, 3785)

# import to the database, overwriting the previous one
st_write_db(con, phillySF, "phillysf", drop = TRUE)
dbGetQuery(con, "SELECT UpdateGeometrySRID('phillysf','wkb_geometry',3785)")

# add Philly's road network - same procedure: read, rename, tranform, write
phillyStreets <- st_read('Street_Centerline.shp', stringsAsFactors = FALSE)
phillyStreets <- rename(phillyStreets,seg_id=SEG_ID, stname = STNAME)
phillyStreets <- st_transform(phillyStreets, 3785)
st_write_db(con, phillyStreets, "phillystreets", drop = TRUE)

# add Philly's crime data
rawcrime <- read.csv("crime.csv")

# remove rows of data with unknown longitude and latitude
cleancrime <- subset(rawcrime, !is.na(Lon))
cleancrime <- subset(cleancrime, !is.na(Lat))

# transform data
crimeSF <- st_as_sf(cleancrime, coords = c("Lon", "Lat"), crs = 4326)
crimeSF <- st_transform(crimeSF, 3785)

# filter raw crime data
# filter time (2014 - 2016)
time <- ymd_hms(crimeSF$Dispatch_Date_Time)
crime1 <-
  crimeSF %>%
  mutate(Year = year(time), Month = month(time))
crime1 <- filter(crime1, Year == 2014 | Year == 2015 | Year == 2016)

# filter burglary crime incidents
crime2 <- filter(crime1, Text_General_Code == "Burglary Non-Residential" | Text_General_Code == "Burglary Residential")

# add a column: season
crime.final <- crime2 %>%
  mutate(Month = as.numeric(Month))
crime.final$season <- cut(crime.final$Month, breaks = c(0, 2, 5, 8, 11, 12),
              labels = c("Winter", "Spring", "Summer", "Fall", "Win"))
crime.final$season[crime.final$season == "Win"] <- "Winter"
crime.final <- arrange(crime.final, Month)

# add filtered crime data to the database
st_write_db(con, crime.final, "burglary14to16", drop = TRUE)

# add spatial index
dbGetQuery(con, "CREATE INDEX phillystreets_gix ON phillystreets USING GIST (wkb_geometry)")
dbGetQuery(con, "CREATE INDEX phillysf_gix ON phillysf USING GIST (wkb_geometry)")
dbGetQuery(con, "CREATE INDEX burglary14to16_gix ON burglary14to16 USING GIST (wkb_geometry)")

# associate burglary incidents with street segments based on this criteria:
# the incident is within 100 m of the street segment
startTime <- Sys.time()
spatialQuery <- paste0("SELECT DISTINCT ON (b.wkb_geometry) b.*, p.seg_id, ",
                       "ST_Distance(b.wkb_geometry, p.wkb_geometry) AS distance ",
                       "FROM burglary14to16 AS b, phillystreets AS p ",
                       "WHERE ST_Distance(b.wkb_geometry, p.wkb_geometry) < 100 ",
                       "ORDER BY b.wkb_geometry, ST_Distance(b.wkb_geometry, p.wkb_geometry) ASC")
StreetJoinBurglary <- st_read_db(con, query=spatialQuery, geom_column='wkb_geometry')
endTime <- Sys.time()
endTime - startTime # 9.8 min # 20951

# just summer
startTime <- Sys.time()
spatialQuery <- paste0("SELECT DISTINCT ON (b.wkb_geometry) b.*, p.seg_id, ",
                       "ST_Distance(b.wkb_geometry, p.wkb_geometry) AS distance ",
                       "FROM burglary14to16 AS b, phillystreets AS p ",
                       "WHERE b.season = 'Summer' ",
                       "AND ST_Distance(b.wkb_geometry, p.wkb_geometry) < 100 ",
                       "ORDER BY b.wkb_geometry, ST_Distance(b.wkb_geometry, p.wkb_geometry) ASC")
StreetJoinBurglary_summer <- st_read_db(con, query=spatialQuery, geom_column='wkb_geometry')
endTime <- Sys.time()
endTime - startTime # 2.66 min # 6262

# just winter
startTime <- Sys.time()
spatialQuery <- paste0("SELECT DISTINCT ON (b.wkb_geometry) b.*, p.seg_id, ",
                       "ST_Distance(b.wkb_geometry, p.wkb_geometry) AS distance ",
                       "FROM burglary14to16 AS b, phillystreets AS p ",
                       "WHERE b.season = 'Winter' ",
                       "AND ST_Distance(b.wkb_geometry, p.wkb_geometry) < 100 ",
                       "ORDER BY b.wkb_geometry, ST_Distance(b.wkb_geometry, p.wkb_geometry) ASC")
StreetJoinBurglary_winter <- st_read_db(con, query=spatialQuery, geom_column='wkb_geometry')
endTime <- Sys.time()
endTime - startTime #  2.07 min # 5153

# write these files to database
StreetJoinBurglary <- rename(StreetJoinBurglary_summer ,oldgeom=wkb_geometry)
st_write_db(con, StreetJoinBurglary, "streetjoinburglary",drop=TRUE)
dbGetQuery(con, "CREATE INDEX streetjoinburglary_gix ON streetjoinburglary USING GIST (wkb_geometry)")

StreetJoinBurglary_summer <- rename(StreetJoinBurglary_summer ,oldgeom=wkb_geometry)
st_write_db(con, StreetJoinBurglary_summer, "streetjoinburglary_summer",drop=TRUE)
dbGetQuery(con, "CREATE INDEX streetjoinburglary_summer_gix ON streetjoinburglary_summer USING GIST (wkb_geometry)")

StreetJoinBurglary_winter <- rename(StreetJoinBurglary_winter ,oldgeom=wkb_geometry)
st_write_db(con, StreetJoinBurglary_winter, "streetjoinburglary_winter",drop=TRUE)
dbGetQuery(con, "CREATE INDEX streetjoinburglary_winter_gix ON streetjoinburglary_winter USING GIST (wkb_geometry)")

# count the number of burglary at each street segment
# summer
spatialQuery <- paste0("SELECT p.seg_id, p.wkb_geometry AS geom, COUNT(b.*) AS crime_count ",
                     "FROM phillystreets AS p ",
                     "JOIN streetjoinburglary_summer AS b ",
                     "ON p.seg_id = b.seg_id ",
                     "GROUP BY (p.seg_id, p.wkb_geometry) ",
                     "ORDER BY crime_count DESC")
summercount <- st_read_db(con, query=spatialQuery, geom_column='geom')

# winter
spatialQuery <- paste0("SELECT p.seg_id, p.wkb_geometry AS geom, COUNT(b.*) AS crime_count ",
                       "FROM phillystreets AS p ",
                       "JOIN streetjoinburglary_winter AS b ",
                       "ON p.seg_id = b.seg_id ",
                       "GROUP BY (p.seg_id, p.wkb_geometry) ",
                       "ORDER BY crime_count DESC")
wintercount <- st_read_db(con, query=spatialQuery, geom_column='geom')

# get natural breaks of my data
natural <- classIntervals(summercount$crime_count, n=4, style="jenks")
print(natural, unique=TRUE)

natural2 <- classIntervals(wintercount$crime_count, n=4, style="jenks")
print(natural2, unique=TRUE)

# join back to street layer
summerjoin <- left_join(phillyStreets, as.data.frame(summercount), by = "seg_id")
winterjoin <- left_join(phillyStreets, as.data.frame(wintercount), by = "seg_id")
summerjoin$crime_count[is.na(summerjoin$crime_count)] <- 0
winterjoin$crime_count[is.na(winterjoin$crime_count)] <- 0

summerjoin$crime_level <- factor(
  cut(summerjoin$crime_count, c(-1, 0, 2, 4, 10)),
  labels = c("0", "1 or 2", "3 or 4", "> 5"))

winterjoin$crime_level <- factor(
  cut(winterjoin$crime_count, c(-1, 0, 2, 4, 10)),
  labels = c("0", "1 or 2", "3 or 4", "> 5"))

# > 3 crime in summer & winter
summerfew <-filter(summerjoin, crime_level == "1 or 2")
summerok <- filter(summerjoin, crime_level == "3 or 4")
summermore <- filter(summerjoin, crime_level == "> 5")
winterfew <-filter(winterjoin, crime_level == "1 or 2")
winterok <- filter(summerjoin, crime_level == "3 or 4")
wintermore <- filter(winterjoin, crime_level == "> 5")

# now ready to plot maps!
# get color palette: viridis
library(colormap)
scales::show_col(colormap(), labels = T)

myTheme <- function() {
  theme_void() +
    theme(
      text = element_text(size = 7),
      plot.title = element_text(size = 20, color = "white", hjust = 0.5, vjust = 0, face = "bold"),
      plot.subtitle = element_text(size = 16, color = "#cccccc", hjust = 0.5, vjust = 0),
      axis.ticks = element_blank(),
      panel.grid.major = element_line(colour = "black"),
      panel.background = element_rect(fill = "black"),
      plot.background = element_rect(fill = "black"),
      legend.direction = "horizontal",
      legend.position = "bottom",
      plot.margin = margin(0.5, 0.5, 0.5, 0.5, 'cm'),
      legend.key.height = unit(0.4, "cm"), legend.key.width = unit(0.8, "cm"),
      legend.title = element_text(size = 14, color = "white", hjust = 0.5, vjust = 0, face = "bold"),
      legend.text = element_text(size = 12, color = "#cccccc", hjust = 0.5, vjust = 0)
    )
}

# summer
ggplot() +
  geom_sf(data = summerjoin, aes(color = crime_level)) +
  geom_sf(data = summerfew, fill = NA, color = "#218f8d", size = 1) +
  geom_sf(data = summerok, fill = NA, color = "#5cc863", size = 2) +
  geom_sf(data = summermore, fill = NA, color = "#fde725", size = 4) +
  # scale_color_viridis(discrete = TRUE, direction = 1, option="viridis", name = "# of Burglary Incident") +
  scale_colour_manual(values = c("#472775", "#218f8d", "#5cc863", "#fde725"),
                      name="Number of Burglary Incidents") +
  labs(title = 'BURGLARY INCIDENTS, SUMMER 2014 - 2016',
       subtitle = "By Street Segment, Philadelphia, PA") +
  guides(colour = guide_legend(override.aes = list(colour = "black", fill=c("#472775", "#218f8d", "#5cc863", "#fde725")))) +
  myTheme()

# winter
ggplot() +
  geom_sf(data = winterjoin, aes(color = crime_level)) +
  geom_sf(data = winterfew, fill = NA, color = "#218f8d", size = 1) +
  geom_sf(data = winterok, fill = NA, color = "#5cc863", size = 2) +
  geom_sf(data = wintermore, fill = NA, color = "#fde725", size = 4) +
  # scale_color_viridis(discrete = TRUE, direction = 1, option="viridis", name = "# of Burglary Incident") +
  scale_colour_manual(values = c("#472775", "#218f8d", "#5cc863", "#fde725"),
                      name="Number of Burglary Incidents") +
  labs(title = 'BURGLARY INCIDENTS, WINTER 2014 - 2016',
       subtitle = "By Street Segment, Philadelphia, PA") +
  guides(colour = guide_legend(override.aes = list(colour = "black", fill=c("#472775", "#218f8d", "#5cc863", "#fde725")))) +
  myTheme()

# add bar plots for better comparison
# spring and fall
spatialQuery <- paste0("SELECT p.seg_id, p.wkb_geometry AS geom, COUNT(b.*) AS crime_count ",
                       "FROM phillystreets AS p ",
                       "JOIN streetjoinburglary AS b ",
                       "ON p.seg_id = b.seg_id ",
                       "WHERE season = 'Spring' ",
                       "GROUP BY (p.seg_id, p.wkb_geometry) ",
                       "ORDER BY crime_count DESC")
springcount <- st_read_db(con, query=spatialQuery, geom_column='geom')

spatialQuery <- paste0("SELECT p.seg_id, p.wkb_geometry AS geom, COUNT(b.*) AS crime_count ",
                       "FROM phillystreets AS p ",
                       "JOIN streetjoinburglary AS b ",
                       "ON p.seg_id = b.seg_id ",
                       "WHERE season = 'Fall' ",
                       "GROUP BY (p.seg_id, p.wkb_geometry) ",
                       "ORDER BY crime_count DESC")
fallcount <- st_read_db(con, query=spatialQuery, geom_column='geom')

# join back to street layer
springjoin <- left_join(phillyStreets, as.data.frame(springcount), by = "seg_id")
falljoin <- left_join(phillyStreets, as.data.frame(fallcount), by = "seg_id")
springjoin$crime_count[is.na(springjoin$crime_count)] <- 0
falljoin$crime_count[is.na(falljoin$crime_count)] <- 0

springjoin$crime_level <- factor(
  cut(springjoin$crime_count, c(-1, 0, 2, 4, 10)),
  labels = c("0", "1 or 2", "3 or 4", "> 5"))

falljoin$crime_level <- factor(
  cut(falljoin$crime_count, c(-1, 0, 2, 4, 10)),
  labels = c("0", "1 or 2", "3 or 4", "> 5"))

# comparison
# create a dataset
# spring
tot1 <- nrow(springjoin)
few1 <- nrow(filter(springjoin, crime_level == "1 or 2"))
ok1 <- nrow(filter(springjoin, crime_level == "3 or 4"))
many1 <- nrow(filter(springjoin, crime_level == "> 5"))

# summer
tot2 <- nrow(summerjoin)
few2 <- nrow(filter(summerjoin, crime_level == "1 or 2"))
ok2 <- nrow(filter(summerjoin, crime_level == "3 or 4"))
many2 <- nrow(filter(summerjoin, crime_level == "> 5"))

#fall
tot3 <- nrow(falljoin)
few3 <- nrow(filter(falljoin, crime_level == "1 or 2"))
ok3 <- nrow(filter(falljoin, crime_level == "3 or 4"))
many3 <- nrow(filter(falljoin, crime_level == "> 5"))

# winter
tot4 <- nrow(winterjoin)
few4 <- nrow(filter(winterjoin, crime_level == "1 or 2"))
ok4 <- nrow(filter(winterjoin, crime_level == "3 or 4"))
many4 <- nrow(filter(winterjoin, crime_level == "> 5"))

# create a data frame
fewc <- data.frame(
  Season = factor(c("Spring","Summer", "Fall", "Winter"), levels=c("Spring","Summer", "Fall", "Winter")),
  Count = c(few1, few2, few3, few4))

okc <- data.frame(
  Season = factor(c("Spring","Summer", "Fall", "Winter"), levels=c("Spring","Summer", "Fall", "Winter")),
  Count = c(ok1,ok2,ok3,ok4))

manyc <- data.frame(
  Season = factor(c("Spring","Summer", "Fall", "Winter"), levels=c("Spring","Summer", "Fall", "Winter")),
  Count = c(many1,many2,many3,many4))

# 1 or 2
a <- ggplot(data=fewc, aes(x=Season, y=Count, fill=Season)) +
  geom_bar(stat="identity", position="dodge") +
  scale_fill_manual(values=c("#6FDEDC","#1A6F6E","#218F8D","#30CFCC")) +
  theme_minimal() +
  theme(legend.position="none",
        text = element_text(size = 12),
        axis.text=element_text(size=12, color = "#BEBEBE", face="bold"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "in"),
        plot.title = element_text(colour="white", size=18, face ="bold"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background = element_rect(fill = "black"),
        plot.background = element_rect(fill = 'black')) +
  ggtitle("# of Street Segments with 1 or 2 Burglary Incidents")

# 3 or 4
b<- ggplot(data=okc, aes(x=Season, y=Count, fill=Season)) +
  geom_bar(stat="identity", position="dodge") +
  scale_fill_manual(values=c("#36A53E","#297A2E","#66CB6D","#A6E1A9")) +
  theme_minimal() +
  theme(legend.position="none",
        text = element_text(size = 12),
        axis.text=element_text(size=12, color = "#BEBEBE", face="bold"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "in"),
        plot.title = element_text(colour="white", size=18, face ="bold"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background = element_rect(fill = "black"),
        plot.background = element_rect(fill = 'black')) +
  ggtitle("# of Street Segments with 3 or 4 Burglary Incidents")

# >5
c<- ggplot(data=manyc, aes(x=Season, y=Count, fill=Season)) +
  geom_bar(stat="identity", position="dodge") +
  scale_fill_manual(values=c("#FEEB4F","#F5DC03","#D9C200","#FEF188")) +
  theme_minimal() +
  theme(legend.position="none",
        text = element_text(size = 12),
        axis.text=element_text(size=12, color = "#BEBEBE", face="bold"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "in"),
        plot.title = element_text(colour="white", size=18, face ="bold"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background = element_rect(fill = "black"),
        plot.background = element_rect(fill = 'black')) +
  ggtitle("# of Street Segments with Over 5 Burglary Incidents")

# combine the three plot
library("cowplot")
plot_grid(a, b, c,
          labels = c("A", "B", "C"),
          ncol = 3, nrow = 1)

# did not use this
Season=c(rep("Spring" , 3) , rep("Summer" , 3) , rep("Fall" , 3) , rep("Winter" , 3) )
Condition=rep(c("1 or 2 Incidents" , "3 or 4 Incidents", "More Than 3 Incidents") , 4)
Value=c(few1, ok1, many1, few2, ok2, many2, few3, ok3, many3, few4, ok4, many4)
data=data.frame(Season,Condition,Value)

# disconnect from the database once you're done
dbDisconnect(con)
dbUnloadDriver(drv)