Building a multi-layer map in R using ggplot2 and spatial data

In this post, we will render spatial point and polygon data on the same map (as separate layers). The finished sample is the following:

This exercise involves using several libraries from R, as well as transforming data from one coordinate system to another. Here are the prerequisites

1. Install required spatial libraries into R: Please refer to this post, on how to do so.
2. Download spatial data for NYC pertaining to safety: Please refer to this post on how to do so.

So let's get started.

First step is to load all the libraries

# Load required libraries
library(rgdal)
library(ggplot2)
library(ggmap)
library(rgeos)
library(maptools)
library(scales)

Each of the above have a specific function in the process. Next step is to set the working directory for the data files. We use the setwd command for that and I have also set a variable for setting the shape file path.

# Set default paths for data
setwd("~/Work/Transportation/NYC/scripts")
shpfilepath="../data/shpfiles"

Next we read the shape files that contain annual fatality and injury data as well as data on borough boundaries for New York Let's start with the injury and fatality data. ogrInfo command provides high level summary of the data in each shape file.

# Read injury and fatality statistics
ogrInfo(shpfilepath,"fatality_yearly")
ogrInfo(shpfilepath,"injury_yearly")

Next, we read the shape file data into R and assign it to an R object variable.

ftlty_yr<-readOGR(shpfilepath, "fatality_yearly")
injry_yr<-readOGR(shpfilepath, "injury_yearly")
names(ftlty_yr); dim(ftlty_yr)
names(injry_yr); dim(injry_yr)

At this point the geometric data is not bound to the data frame, as we can see.

Also, lets create our first graph with ggplot. We read the data frame from the fatalities data. As we can see there is no geometry for us to plot at this point, so we will simply plot Fatalities against the year. GGPlot2 library distinguishes between data, aesthetics, from the graph. In the first example, we simply add the data to the base layer and ask ggplot to render it as a point.

# Extract the data frame and use ggplot
ftldata<-ftlty_yr@data
names(ftldata)
p<-ggplot(ftlty_yr@data, aes(Fatalities, YR))
p + geom_point()
p

This is the graph that was plotted...

This is our first attempt at understanding ggplot, but not what we wanted. The problem is that when we read the shapefile, the spatial data was not brought along.

We need to use the fortify command to add the spatial data. In the following steps, we can see how to do that. Once we fortify the data, we can extract the x and y coordinates as well.

ftlty.f <-fortify(as.data.frame(ftlty_yr))
names(ftlty.f);dim(ftlty.f)
fx<-ftlty.f$coords.x1
fy<-ftlty.f$coords.x2

Now, since we have the x and y, we can plot this data using ggplot

# Plotting the data using ggplot
Map <- ggplot(data=ftlty.f, aes(fx,fy,colour=factor(Fatalities)))
Map <- Map + layer(geom="point")
Map

Here is the plot output

Next step is to read the boroughs data. We use the familiar ogrInfo and readOGR functions to read the shapefile data into R

#Load boroughs data
ogrInfo(shpfilepath,"nybbwi")
bbnd <- readOGR(shpfilepath,"nybbwi")

The Boroughs data is in a different coordinate system than the accident data. While the accident data is in latitude longitude, the borough data is in Lambert conformal conic coordinate system. If we want to have the two map layers sit on top of one another, we need to have them in the same coordinate system. We use the spTransform function to transform the borough boundaries to latlong, as shown below...

#Transform the coordinates from Lambert Conformal to Lat Long
CRS.latlong<-CRS("+proj=longlat +datum=WGS84 +no_defs")
bbnd<-spTransform(bbnd, CRS.latlong)

In this case, again the spatial data and non-spatial data is not together. We need one dataset with both, so we need to fortify and merge the spatial and non-spatial datasets.

# Add the spatial data back to the boundary data.frame
bbnd.f <-fortify(bbnd, region="BoroCode")
names(bbnd.f);fix(bbnd.f);dim(bbnd.f)

# Combine the geometry and data back together using common keys
bbnd.f<-merge(bbnd.f,bbnd@data, by.x="id", by.y="BoroCode")
names(bbnd.f);dim(bbnd.f)

As we can see in the following window, using fortify alone brought in the 7 geometric attributes and merge brought the total to 10 attributes with attributes from the data frame. So, we need both fortify and merge commands to bring in the geometry and associated non-spatial attributes such as area and length together.

Next step is to build the map one layer at a time.

# Plot a map of the boroughs and accidents
baes <-aes(x=long, y=lat, group = group ,alpha=1/50,fill=id) # Aesthetics for Boroughs
bblyr <-geom_polygon(data=bbnd.f, baes) # Polygon layer for Bouroughs
faes <-aes(x=fx, y=fy, colour=factor(Fatalities)) # Aesthetics for Fatalities
flyr <-geom_point(data=ftlty.f,faes) # Point layer for Fatalities
Map <- ggplot() # Initialize the map
Map <-Map + bblyr # Add the boroughs layer
Map <-Map + flyr # Add the Fatalities layer
Map <- Map + labs(x = "Longitude", y = "Latitude", fill = "NYC Boroughs") # Add the labels
Map <- Map + ggtitle("NYC Boroughs with Accidents")
Map # Display the Map

Here is the console output.

Here is the map

That's it... For those who are interested, here is the complete source listing..

# Load required libraries
library(rgdal)
library(ggplot2)
library(ggmap)
library(rgeos)
library(maptools)
library(scales)

# 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_yr<-readOGR(shpfilepath, "fatality_yearly")
injry_yr<-readOGR(shpfilepath, "injury_yearly")
names(ftlty_yr); dim(ftlty_yr)
names(injry_yr); dim(injry_yr)

# Extract the data frame and use ggplot
ftldata<-ftlty_yr@data
names(ftldata)
p<-ggplot(ftlty_yr@data, aes(Fatalities, YR))
p + geom_point()
p

# Add the spatial data to the data frame
ftlty.f <-fortify(as.data.frame(ftlty_yr))
names(ftlty.f);dim(ftlty.f)
fx<-ftlty.f$coords.x1
fy<-ftlty.f$coords.x2

# Plotting the data using qplot
qplot(fx,fy,data=ftlty.f)
qplot(fx,fy,data=ftlty.f,colour=factor(Fatalities))
qplot(fx,fy, data=ftlty.f,colour=factor(Fatalities), facets= . ~ YR) + geom_point()

# Plotting the data using ggplot
Map <- ggplot(data=ftlty.f, aes(fx,fy,colour=factor(Fatalities)))
Map <- Map + layer(geom="point")
Map
 
#Load boroughs data
ogrInfo(shpfilepath,"nybbwi")
bbnd <- readOGR(shpfilepath,"nybbwi")

#Transform the coordinates from Lambert Conformal to Lat Long
CRS.latlong<-CRS("+proj=longlat +datum=WGS84 +no_defs")
bbnd<-spTransform(bbnd, CRS.latlong)

#Examine the data
bbndata <-bbnd@data
names(bbndata);fix(bbndata);dim(bbndata)

# Add the spatial data back to the boundary data.frame
bbnd.f <-fortify(bbnd, region="BoroCode")
names(bbnd.f);fix(bbnd.f);dim(bbnd.f)

# Combine the geometry and data back together using common keys
bbnd.f<-merge(bbnd.f,bbnd@data, by.x="id", by.y="BoroCode")
names(bbnd.f);dim(bbnd.f)

# Plot a map of the boroughs and accidents
baes<-aes(x=long, y=lat, group = group ,alpha=1/50,fill=id) # Aesthetics for Boroughs
bblyr<-geom_polygon(data=bbnd.f, baes) # Polygon layer for Bouroughs
faes<-aes(x=fx, y=fy, colour=factor(Fatalities)) # Aesthetics for Fatalities
flyr<-geom_point(data=ftlty.f,faes) # Point layer for Fatalities
Map <- ggplot() # Initialize the map
Map<-Map + bblyr # Add the boroughs layer
Map<-Map + flyr # Add the Fatalities layer
Map <- Map + labs(x = "Longitude", y = "Latitude", fill = "NYC Boroughs") # Add the labels
Map <- Map + ggtitle("NYC Boroughs with Accidents")
Map # Display the Map

Reading shape files in R

In the previous post, we saw how to read dbf files in R. I wanted to go a step further.

I was looking at ways to read shape file data in R. It turns out there are various mechanisms for doing just that. However, for the library to be installed correctly, your environment needs to be setup correctly.

Here are the steps I followed.

1. Install libgdal
2. Install libproj
3. Install rgdal package
4. Test loading a shape file

The following script should install these packages correctly for you.

sudo apt-get install apt-file
sudo apt-get update
sudo apt-get install libgdal1h
sudo apt-get install libgdal1-dev libproj-dev

Once these packages are setup correctly start R as sudo

sudo R

Enter the following commands on the R console

install.packages("rgdal")
install.packages("rgeos")
install.packages("ggmap")
install.packages("maptools")
q()

The following screen shows installation for rgeos

Next run the following commands in R (without sudo)

library(rgdal)
shpfilepath="../data/shpfiles"
ogrInfo(shpfilepath,"fatality_yearly")
ftlty_yr <-readOGR(shpfilepath, "fatality_yearly")

This shows the following:

Now, enter the following command in R

plot(ftlty_yr, axes=TRUE, border="gray")

Our R environment can now read shapefiles correctly