ggplots in R

Table of Contents

Introduction

Even though the R base plot capabilities are more impressive than most other programming language, there are other packages to generate graphics available. Two of the more popular packages besides the base package is lattice and ggplot2. They are according to many superior to the base plot library when it comes to explanatory data analysis as they without much work generate trellis graphics eg. graphs that display a variable or the relationship between variables, conditioned on one or more other variables. Over the last years ggplot2 has become the standard plot library for many R users, especially as it keeps evolving and new features are added. In addition to being more convient for certain types of plots many feel that the defaults colors, axis types etc looks better on ggplot2 compared to the base R and lattice libraries.

Exercise

This exercise will not consist of any exercises for you to solve instead the code blocks below will read in data and generate map plots. Make sure that you can run the code and generate the plot. Once you have a working code go over the code line by line and try to obtain a fairly detailed understanding of what is going on. Read the help for some of the functions and try the plots with other arguments.

Before you start you should download the following files to your working directory.
Bolaget.csv
coop_Uppsala.kml
ica_Uppsala.kml

If you can not load the library install it either via the menu of Rstudio or by typing install.packages("ggmap") in R.

# Load the libraries necessary for working with maps
library(ggmap)
library(maptools)
library(sp)
library(rgdal)
library(deldir)
    
# Download map of Uppsala from stamen maps
google.map <- get_map(c(17.63,59.84), zoom=12, maptype = 'toner')
    
# Lad the list of shops
data <- read.table("Bolaget.csv", header=T, sep=";", quote="")
    
# Narrow down the list to Uppsala only
		data <- data[data$Address4 == "UPPSALA",]
		
# The description in the file says the coords are in RT90 datum.
# The map uses WGS84 thus we need a conversion:
latlonRT90 <- data[,c('RT90x', 'RT90y')]
colnames(latlonRT90) <- c('x','y')

# EPSG codes for RT90 and WGS84 are 3021 and 4326. 
# Here, we do the actual conversion
tmp <- data.frame(coords.x = latlonRT90$y, coords.y = latlonRT90$x)
coordinates(tmp)=~coords.x+coords.y
proj4string(tmp)=CRS("+init=epsg:3021") 
coords <- spTransform(tmp, CRS("+init=epsg:4326"))
		
# Create the data frame for ggplot2
coords <- data.frame(lat=coords@coords[,1], lon=coords@coords[,2])

# Do the Voronoi tesseleation
voronoi <- deldir(coords)

# Plot the map, shop density, shops as points and the tesselation lines.
map <- ggmap(google.map)
map +
stat_density2d(data = coords, 
	           aes(x = lat, y = lon,
			   fill = ..level..,
			   alpha = ..level..), 
			   size = 0.01, bins = 50, 
			   geom = "polygon") + 
			   scale_fill_gradient(low = "green", high = "olivedrab", guide = FALSE) + 
			   scale_alpha(range = c(0, 0.1), guide = FALSE) +
			   geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), 
			                size = .6, linetype=2, data = voronoi$dirsgs, color= "olivedrab") +
							geom_point(aes(x=lat, y=lon), data=coords, color='olivedrab')
# ICA & Coop
ica.kml <- getKMLcoordinates(kmlfile="ica_Uppsala.kml", ignoreAltitude=T)
tmp <- unlist(ica.kml)
ica.coords <- data.frame(lat=tmp[1:length(tmp) %% 2 == 0], lon=tmp[1:length(tmp) %% 2 == 1], type='ica')
coop.kml <- getKMLcoordinates(kmlfile="coop_Uppsala.kml", ignoreAltitude=T)
tmp <- unlist(coop.kml)
coop.coords <- data.frame(lat=tmp[1:length(tmp) %% 2 == 0], lon=tmp[1:length(tmp) %% 2 == 1], type='coop')
coords <- rbind(ica.coords, coop.coords)
google.map <- get_map(c(17.63,59.84), zoom=12)
voronoi <- deldir(coords)

map <- ggmap(google.map)
	         map +
			 scale_fill_gradient(low = "green", high = "olivedrab", guide = FALSE) + 
			 scale_alpha(range = c(0, 0.1), guide = FALSE) + 
			 geom_segment(aes(x = y1, y = x1, xend = y2, yend = x2), 
			 size = .4, linetype=1, data = voronoi$dirsgs, color= "olivedrab") + 
			 geom_point(aes(x=lon, y=lat, col=type), data=coords)