You can also switch to street view
Let's get started. Make sure you have installed the required libraries. The new library needed is googleVis which can be installed by launching R in sudo mode and entering the following command
install.packages("googleVis")
q()
First step is to download the data. This previous post shows how to obtain the spatial data used in this example from the NYC website. Now launch R in regular mode and enter the following code to create the map.
First we need to load the required libraries
# Load the needed libraries library(sp) # classes for spatial data library(maptools) library(rgeos) # and their dependencies library(rgdal) library(ggplot2) library(ggmap) library(scales) library(gmodels) library(googleVis)
Next we need to set the working directory where our scripts are located and declare a variable to folder where our data is located.
# Set default paths for data
setwd("~/Work/Transportation/NYC/scripts")
shpfilepath="../data/shpfiles"
Next we will read the shape files with data. We have two shape files from NYC. First contains only fatal accidents, and second contains injury accidents. We will later combine these files.
# Read injury and fatality statistics ogrInfo(shpfilepath,"fatality_yearly") ogrInfo(shpfilepath,"injury_yearly") ftlty<-readOGR(shpfilepath, "fatality_yearly") injry<-readOGR(shpfilepath, "injury_yearly") names(ftlty); dim(ftlty) names(injry); dim(injry)
Next we need to rename the column names across the two datasets so that they are identical. Following step shows us how to do that. We also will add a new variable that classifies the incident as an injury or a fatality incident.
#Classify severity and rename columns to be identical in fatalities and injury datasets ftlty$Severity<-"Fatality" names(ftlty)[names(ftlty)=="Fatalities"]<-"Total" names(ftlty)[names(ftlty)=="PedFatalit"]<-"Ped" names(ftlty)[names(ftlty)=="BikeFatali"]<-"Bike" names(ftlty)[names(ftlty)=="MVOFatalit"]<-"MVO" injry$Severity<-"Injury" names(injry)[names(injry)=="Injuries"]<-"Total" names(injry)[names(injry)=="PedInjurie"]<-"Ped" names(injry)[names(injry)=="BikeInjuri"]<-"Bike" names(injry)[names(injry)=="MVOInjurie"]<-"MVO" names(injry);names(ftlty)
Now we can combine the two datasets to create one single dataset called acc.
#Combine datasets acc<- rbind(ftlty, injry) dim(acc);names(acc);head(acc)
We need to use the fortify command to make sure the coordinate information can be accessed by spatial functions.
#Now fortify the accident data acc.f<-fortify(as.data.frame(acc)) dim(acc.f);names(acc.f);head(acc.f)
We only want to plot 2015 data. The following command shows us how to filter only 2015 records from the data.
#Subsetting only 2015 accident records acc_2015<-acc.f[which(acc.f$YR==2015),]
Now we are ready to plot this on the map. The gvisMap function expects a single field with latitude and longitude information in lat:long format, that is separated by a colon. We don't need more than 4 digit precision so we will round the lat long information and use the paste command to concatenate everything into a single string.
#Add a single column "LatLong" that is of the format lat:long for the accident location acc_2015$LatLong<-paste(round(acc_2015$coords.x2,4),round(acc_2015$coords.x1,4),sep=":")
The gVisMap function also expects a data field which will be displayed as a map tip. We need to concatenate all available fields with html break symbol to create one single field that can pop up when we hover or click on the icon.
#Add a single column "Tip" for Accident descriptions separated by html breaks ("<br/>")
acc_2015$Tip<-paste("Type: ",acc_2015$Severity,"Total: ",acc_2015$Total,"Ped: ",acc_2015$Ped,"Bike: ",acc_2015$Bike,"MVO: ",acc_2015$MVO,sep="<br/>")
head(acc_2015)
Finally, we call the gVisMap function passing the dataset, name of concatenated lat long column, name of the Map tip column and set several defaults for the Google map. Calling the plot function launches an embedded http server and launches the default browser on your machine to display the generated html page.
#Plot 2015 accidents on Google
map1<-gvisMap(acc_2015, "LatLong", "Tip",
options=list(showTip=TRUE, showLine=F, enableScrollWheel=TRUE,
mapType='roads', useMapTypeControl=TRUE, width=800,height=1024))
plot(map1)
That's it... We are done. Here is the output
Here is the complete source code
# Load the needed libraries
library(sp) # classes for spatial data
library(maptools)
library(rgeos) # and their dependencies
library(rgdal)
library(ggplot2)
library(ggmap)
library(scales)
library(gmodels)
library(googleVis)
# Set default paths for data
setwd("~/Work/Transportation/NYC/scripts")
shpfilepath="../data/shpfiles"
# Read injury and fatality statistics
ogrInfo(shpfilepath,"fatality_yearly")
ogrInfo(shpfilepath,"injury_yearly")
ftlty<-readOGR(shpfilepath, "fatality_yearly")
injry<-readOGR(shpfilepath, "injury_yearly")
names(ftlty); dim(ftlty)
names(injry); dim(injry)
#Classify severity and rename columns to be identical in fatalities and injury datasets
ftlty$Severity<-"Fatality"
names(ftlty)[names(ftlty)=="Fatalities"]<-"Total"
names(ftlty)[names(ftlty)=="PedFatalit"]<-"Ped"
names(ftlty)[names(ftlty)=="BikeFatali"]<-"Bike"
names(ftlty)[names(ftlty)=="MVOFatalit"]<-"MVO"
injry$Severity<-"Injury"
names(injry)[names(injry)=="Injuries"]<-"Total"
names(injry)[names(injry)=="PedInjurie"]<-"Ped"
names(injry)[names(injry)=="BikeInjuri"]<-"Bike"
names(injry)[names(injry)=="MVOInjurie"]<-"MVO"
names(injry);names(ftlty)
#Combine datasets
acc<- rbind(ftlty, injry)
dim(acc);names(acc);head(acc)
#Now fortify the accident data
acc.f<-fortify(as.data.frame(acc))
dim(acc.f);names(acc.f);head(acc.f)
#Subsetting only 2015 accident records
acc_2015<-acc.f[which(acc.f$YR==2015),]
#Add a single column "LatLong" that is of the format lat:long for the accident location
acc_2015$LatLong<-paste(round(acc_2015$coords.x2,4),round(acc_2015$coords.x1,4),sep=":")
#Add a single column "Tip" for Accident descriptions separated by html breaks ("<br/>")
acc_2015$Tip<-paste("Type: ",acc_2015$Severity,"Total: ",acc_2015$Total,"Ped: ",acc_2015$Ped,"Bike: ",acc_2015$Bike,"MVO: ",acc_2015$MVO,sep="<br/>")
head(acc_2015)
#Plot 2015 accidents on Google
map1<-gvisMap(acc_2015, "LatLong", "Tip",
options=list(showTip=TRUE, showLine=F, enableScrollWheel=TRUE,
mapType='roads', useMapTypeControl=TRUE, width=800,height=1024))
plot(map1)
No comments:
Post a Comment