author: Janet Prumachuk

setwd(".")

library("dplyr")
## 
## Attaching package: 'dplyr'
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("ggmap")
## Warning: package 'ggmap' was built under R version 3.1.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.1.3
library("ggplot2")
library("ggthemes")
## Warning: package 'ggthemes' was built under R version 3.1.3
library("gpclib")
## General Polygon Clipper Library for R (version 1.5-5)
##  Type 'class ? gpc.poly' for help
library("graphics")
library("grDevices")
library("grid")
library("gridExtra")
library("RColorBrewer")
library("lubridate")
library("maptools")
## Loading required package: sp
## Checking rgeos availability: TRUE
library("reshape2")
library("rgdal")
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj
library("rgeos")
## rgeos version: 0.3-8, (SVN revision 460)
##  GEOS runtime version: 3.3.3-CAPI-1.7.4 
##  Polygon checking: TRUE 
## 
## 
## Attaching package: 'rgeos'
## 
## The following objects are masked from 'package:gpclib':
## 
##     append.poly, area.poly, get.bbox, get.pts, read.polyfile,
##     scale.poly, triangulate, tristrip, write.polyfile
library("sp")

rawdata <- read.csv("./DATA/NYPD_Motor_Vehicle_Collisions-2014.csv", header=TRUE)
names(rawdata)[4] <- "ZIP"
rawdata$DATE <- as.character(rawdata$DATE)

data <- group_by( rawdata, ZIP)
data <- summarize( data, 
                   Collisions = n(), 
                   Fatalities = sum(`NUMBER.OF.PERSONS.KILLED`),
                   BicycleFatalities = sum(`NUMBER.OF.CYCLIST.KILLED`),
                   BicycleInjuries = sum(`NUMBER.OF.CYCLIST.INJURED`))

rawdata$PDATE <- mdy(rawdata$DATE)
rawdata$DAY   <- wday(rawdata$PDATE, label=TRUE, abbr=TRUE)

ftotal <- summarize(rawdata, sum(`NUMBER.OF.PERSONS.KILLED`))

data_by_day <- group_by(rawdata, DAY)
data_by_day <- summarize (data_by_day, Fatalities = sum(`NUMBER.OF.PERSONS.KILLED`))
data_by_day$Pct <- (data_by_day$Fatalities / as.numeric(ftotal)) * 100

Get map

setwd("./zip_codes_shp")
zipmap <- readOGR(dsn=".","PostalBoundary")
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "PostalBoundary"
## with 2332 features and 6 fields
## Feature type: wkbPolygon with 2 dimensions
zipmap_nyc <- zipmap[zipmap$POSTAL %in% factor(data$ZIP),]
zipmap_nyc <- spTransform(zipmap_nyc, CRS("+proj=longlat + datum=WGS84"))
setwd("../")

Merge map with data by ZIP

zipmap_nyc_df <- fortify(zipmap_nyc)
## Regions defined for each Polygons
map_data <- data.frame(id=rownames(zipmap_nyc@data), ZIP=zipmap_nyc@data$POSTAL, ZIPNAME=zipmap_nyc@data$NAME)
map_data <- merge(map_data, data, by="ZIP")
map_df   <- merge(zipmap_nyc_df, map_data, by="id")

map_df$CollisionRange <- cut(map_df$Collisions, breaks=quantile(map_df$Collisions), dig.lab=4, include.lowest=TRUE)

Make plots

p1 <- qmap('new york, ny', zoom=11, maptype='roadmap', color='bw', legend='topleft') +
        geom_polygon(aes(long, lat, group=group, fill=CollisionRange), 
                     data=map_df, 
                     color='white',
                     alpha=.6, 
                     size=.3) + 
        ggtitle("Collisions Reported in New York City - 2014") +
        scale_fill_brewer("Collision Range")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=new+york,+ny&zoom=11&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=new+york,+ny&sensor=false
p2 <- qmap('new york, ny', zoom=11, maptype='roadmap', color='bw', legend='topleft') +
        geom_polygon(aes(long, lat, group=group, fill=factor(Fatalities)),
                     data=map_df,
                     color='white',
                     alpha=.6, 
                     size=.3) + 
        ggtitle("Collisions Reported in New York City - 2014") +
        scale_fill_brewer("Fatalities")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=new+york,+ny&zoom=11&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
p3 <- qmap('new york, ny', zoom=11, maptype='roadmap', color='bw', legend='topleft') +
        geom_polygon(aes(long, lat, group=group, fill=factor(BicycleFatalities)),
                     data=map_df,
                     color='white',
                     alpha=.6, 
                     size=.3) + 
        ggtitle("Collisions Reported in New York City - 2014") +
        scale_fill_brewer("Bicycle Fatalities")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=new+york,+ny&zoom=11&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
p4 <- ggplot(data_by_day, aes(x=DAY, y=Pct, fill=Pct)) + 
        geom_bar(stat="identity") + 
        scale_fill_gradient(low="#aaaaaa", high="#333333") + 
        ggtitle("NYC Fatal Collisions by Day of the Week - 2014") + 
        theme_fivethirtyeight() + 
        xlab(NULL) + ylab("% of Total") + 
        guides(fill=FALSE) + 
        scale_color_fivethirtyeight()

Show plots and save as images