The undiscovered country – a tutorial on plotting maps in R

The ability to handle maps and geospatial images is always a nice trick to have up your sleeve. Almost any sizeable report will contain a map – of a locality, of a country, or of the world. However, very few analysts have the ability to produce these plots for themselves and often resort to using Google Earth snapshots, images that they have drawn in paint, or – heaven forefend – a map that they found on Google Images. Making maps appears to be a skill held only by a few citizens of some undiscovered country, from whose bourn few analysts return.

In this post we chart the shores of the undiscovered country, and show that drawing maps, such as the one shown in Figure 1, is not so difficult. We will make use of R's powerful functions that allow plotting of maps, locations on those maps, and connecting lines that illustrate links between points of interest.

Figure 1: Map of Australia's state capitals. Connections have been selected arbitrarily for illustrative purposes.

We will also make use of Python, as I have found it far easier to use the Python bindings for the GoogleMaps API. Our approach consists of three steps:

1. Work out what you want to plot: the region, the nodes (ie, the locations), and the edges (ie, the connections between nodes).
2. Get the coordinates of the nodes.
3. Read and plot the data in R.

Step 1: What do you want to plot?

As a proud Australian, I am trying to increase the profile of my tiny country. All too often, the maps that I see on R bloggers are of the United States. Granted, the US has 300 million people and, despite the inexorable rise of China, remains the economic super-power of the world. Kudos to Uncle Sam. But Australia can boast a similar landmass, a bucket-load of natural resources, and more species of poisonous snakes and man-eating reptiles than any other country on the planet (a point that we have down-played since the unsuccessful tourism slogan "come to Australia and you might get eaten"). No, this article is going to be printed next to a map of the Land Downunder dammit!

To illustrate how one can plot points, label them, and connect them with lines, we will plot each of Australia's state capitals, and connect selected pairs with lines. These lines will be based on the geodesic between the two points that uses the geosphere package (see here). To complete this task, there are three pieces of input data that you will require:  a shape-file, a list of nodes, and an edge-list.

Shapefile
 I have used a shape-file located here: http://www.gadm.org/country. This website also contains shape files for other countries and is definitely worth a look.

List of nodes
The nodes could be anything from airports to electricity transmission hubs. In my case, the nodes are just the capital cities of each of the states and territories. The list, which I have stored in the CSV file "CapitalCities.csv", consists of two columns: the name and the state/territory of each city.

Edge-list
The edge-list contains two columns: the 'from' column and the 'to' column. Each entry represents a connection between two nodes. In this instance, the two columns are interchangeable as the network is not directed, however, it would be easy to incorporate this feature. The edge-list that I have used can be accessed here.

Step 2: Get the coordinates of the nodes

Getting coordinates of a specific location is a simple affair thanks to the wonders of Google Maps and Google Earth. The problem is not how to get the data, but instead how to get it programmatically. For ten nodes, you can just as easily do this manually. For 1000 nodes, you can either find some intern to sit in front of Google Earth for the next week, or accept that there must be a better way.

I am sure that there are solutions to this problem that use R, but I have not yet found anything that is as quick and easy as the Python wrappers for the Google Maps API. (If you are unfamiliar with Python, then I advise you to take the time to learn the basics – your time will not have been wasted.) The module can be downloaded at http://pypi.python.org/pypi/googlemaps/ or via easy_install. An API key is needed to interface with the GoogleMaps API, but this is easy enough to obtain – just search for 'Google Maps API key' and follow the instructions.

The following script retrieves the latitude and longitude of each city in the file 'CapitalCities.csv', and copies them into a new file 'CapCoordinates.csv'.

from googlemaps import GoogleMaps
import csv

ifile=csv.reader(open('CapitalCities.csv','rb'),delimiter=',')
ofile=open('CapCoordinates.csv','wb')

w=csv.writer(ofile,delimiter=',')
w.writerow(['Name','State','Latitude','Longitude'])

gmaps=GoogleMaps(API_KEY)

count=0

for row in ifile:
   if count!=0:
      address=str(row[0])+" "+str(row[1])
      lat, lng = gmaps.address_to_latlng(address)
      w.writerow([str(row[0]),str(row[1]),str(lat),str(lng)])
      print row[0],lat,lng

   count+=1

ofile.close() 

Step 3: Plotting the map in R

At last we arrive at the R portion of this article. We will need to load four packages: the maps, sp, maptools, and geosphere libraries. Below is the complete code.

# Load necessary packages
library(maps)
library(sp)
library(maptools)
library(geosphere)

# Function that returns coordinates of a given location
getCoords<-function(location,df){
  return(df[match(x=location,df[,"Name"]),][,c(4,3)])
}

# Plot a great circle from 'from' to 'to' names.
plotLine<-function(from,to,df){
  inter<-gcIntermediate(p1=getCoords(from,df),
                      p2=getCoords(to,df),
                      n=50,addStartEnd=TRUE)
  lines(inter, col="green",cex=2,lwd=2)
}

# Read in the Australia Shapefile
ozDATA<-readShapeSpatial("AUS_adm1")

# Read CSV file containing coordinates and edges
nodes<-read.csv(file="CapCoordinates.csv",sep=",",header=TRUE)
edges<-read.csv(file="EdgeListCapCities.csv",sep=",",header=TRUE)

# Plot the graph itself
plot(ozDATA,lwd=0.01,bg="#242424",col="#0000CD",ylim=c(-46,-10),xlim=c(125,145))

# Plot the nodes on the graph
points(nodes$Longitude,nodes$Latitude,pch=16,cex=1.5)

# Label the nodes
text(nodes$Longitude,nodes$Latitude,nodes$Name,cex=1,adj=0,pos=2,col="#C2C2C2")

# Plot each of the edges between nodes
for (j in 1:nrow(edges)){
   plotLine(edges[j,]$From,edges[j,]$To,nodes) 
}

Line-by-line walkthrough

We start by reading in the Australia shape-file using the readShapeSpatial function from the sp package.

ozDATA<-readShapeSpatial("AUS_adm1")

Then we read in our list of nodes and freshly obtained coordinates, as well as the edge list.

nodes<-read.csv(file="CapCoordinates.csv",sep=",",header=TRUE)
edges<-read.csv(file="EdgeListCapCities.csv",sep=",",header=TRUE)

Having loaded all the necessary data, we can start plotting.

plot(ozDATA,lwd=0.01,bg="#242424",col="#0000CD",ylim=c(-46,-10),xlim=c(125,145))

We plot the shape-file using the plot function with the following arguments:

• lwd – the line width;
• bg – the background colour;
• col – the colour of the regions (usually landmasses);
• ylim – a 2-tuple containing the latitude boundaries for the plot; and
• xlim – a 2-tuple containing the longitude boundaries for the plot.

Now we want to plot the cities (ie, the nodes) on this map. To do this we simply pass the coordinates to the points function. The parameters pch and cex set the type and size of the marker, respectively.

points(nodes$Longitude,nodes$Latitude,pch=16,cex=1.5)

Next we want to add labels to each of the nodes.

text(nodes$Longitude,nodes$Latitude,nodes$Name,cex=1,adj=0,pos=2,col="#C2C2C2")

The text function is used in almost identical fashion to the points function, but with an additional argument that contains the names of the nodes. We also add two additional arguments: adj and pos. The adj parameter adjusts the location of the labels. In this instance we have set it to zero, as it will be overwritten by our specification of pos. The pos parameter takes a value between 1 and 4 that specifies whether the label will be below, to the left, above, or to the right of the specified coordinates. (In this instance we have chosen to have labels to the left of the nodes.)

All that remains is to plot the lines between the respective cities. We loop through each of the edges and call the plotLine function.

for (j in 1:nrow(edges)){
   plotLine(edges[j,]$From,edges[j,]$To,nodes) 
}

getCoords<-function(location,df){
  return(df[match(x=location,df[,"Name"]),][,c(4,3)])
}

plotLine<-function(from,to,df){
  inter<-gcIntermediate(p1=getCoords(from,df),
                      p2=getCoords(to,df),
                      n=50,addStartEnd=TRUE)
  lines(inter, col="green",cex=2,lwd=2)
}

plotLine simply takes each pair of nodes in the edge-list, and uses the gcIntermediate function (as defined in the geosphere package) to join the two nodes with the shortest edge that lies upon a great circle (ie, the geodesic between the two nodes). To simplify this function, I have defined an additional function, getCoords, that takes as arguments a string (the node name) and a data frame (the nodes and their coordinates) and returns the coordinates of that particular node. After the for loop is completed, we arrive at the graph shown in Figure 1.

Conclusion
Nothing here was particularly complicated. We decided what we wanted to plot, collected the data, and used R to handle the data and produce the image. However, the applications for the tools that we have used are practically endless. To that end, I would be keen to hear from people who have been exploiting R's mapping functions across different applications. Ideally, I would like to showcase a range of applications in my next posting. Feel free to leave comments and links to your websites.

Until next time, enjoy being one of the few analysts who has visited the undiscovered country and can make maps in R with ease.

Web scraping with Python - the dark side of data

In searching for some information on web-scrapers, I found a great presentation given at Pycon in 2010 by Asheesh Laroia. I thought this might be a valuable resource for R users who are looking for ways to gather data from user-unfriendly websites. The presentation can be found here:

http://python.mirocommunity.org/video/1616/pycon-2010-scrape-the-web-stra

Highlights (at least from my perspective)

• Screen scraping is not about regular expressions. It is just too hard to use pattern matching for these tasks, as the tags can change regularly and have significant maintenance issues.
• BeautifulSoup is the go-to html parser for poor quality source. I have used this in the past and am pleased to hear that I was not too far off the money!
• Configuration of User Agent settings is discussed in detail, as well as other mechanisms that websites exploit to stop you from scraping content
• Good description of how to use the Live HTTP Headers add-on for Firefox.
• A thought-provoking discussion about APIs, and comments that suggest that their maintenance and support is woefully inadequate. I was interested to hear his views, as they imply that scraping may be the only alternative when you really need data that is highly inaccessible.

Other notes

The mechanise package features heavily in the examples for this presentation. The following link provides some good examples of how to use mechanise to automate forms:
http://wwwsearch.sourceforge.net/mechanize/forms.html

There was also some mention of how Javascript causes problems for web scrapers, although this problem can be overcome via the use of web-drivers such as Selenium (see http://pypi.python.org/pypi/selenium) and Watir. I have used safari-watir before, and from my experience it can perform many complex data gathering tasks with relative ease.

Please feel free to post your comments about your experiences with screen scraping, and other tools that you use to collect web data for R.

The R apply function – a tutorial with examples

Today I had one of those special moments that is uniquely associated with R. One of my colleagues was trying to solve what I term an 'Excel problem'. That is, one where the problem magically disappears once a programming language is employed. Put simply, the problem was to take a range, and randomly shift the elements of the list in order. For example, 12345 could become 34512 or 51234.

The list in question had forty-thousand elements, and this process needed to be repeated numerous times as part of a simulation. Try doing this in Excel and you will go insane: the shift function is doable but resource intensive. After ten minutes of waiting for your VBA script to run you will be begging for mercy or access to a supercomputer. However, in R the same can be achieved with the function:

translate<-function(x){
  if (length(x)!=1){
    r<-sample(1:(length(x)),1)
    x<-append(x[r:length(x)],x[1:r-1])
  }
  return(x)

}

My colleague ran this function against his results several thousand times and had the pleasure of seeing his results spit out in less than thirty seconds: problem solved. Ain't R grand.

More R magic courtesy of the apply function

The translate function above is not rocket science, but it does demonstrate how powerful a few lines of R can be. This is best exemplified by the incredible functionality offered by the apply function. However, I have noticed that this tool is often under-utilised by less experienced R users.

The usage from the R Documenation is as follows:

apply(X, MARGIN, FUN, ...)

where:

• X is an array or matrix;
• MARGIN is a variable that determines whether the function is applied over rows (MARGIN=1), columns (MARGIN=2), or both (MARGIN=c(1,2));
• FUN is the function to be applied.

In essence, the apply function allows us to make entry-by-entry changes to data frames and matrices. If MARGIN=1, the function accepts each row of X as a vector argument, and returns a vector of the results. Similarly, if MARGIN=2 the function acts on the  columns of X. Most impressively,  when MARGIN=c(1,2) the function is applied to every entry of X. As for the FUN argument, this can be anything from a standard R function, such as sum or mean, to a custom function like translate above.

An illustrative example
Consider the code below:

# Create the matrix
m<-matrix(c(seq(from=-98,to=100,by=2)),nrow=10,ncol=10)

# Return the product of each of the rows
apply(m,1,prod)

# Return the sum of each of the columns
apply(m,2,sum)

# Return a new matrix whose entries are those of 'm' modulo 10
apply(m,c(1,2),function(x) x%%10) 

In the last example, we apply a custom function to every entry of the matrix. Without this functionality, we would be at something of a disadvantage using R versus that old stalwart of the analyst: Excel. But with the apply function we can edit every entry of a data frame with a single line command. No autofilling, no wasted CPU cycles.

In the next edition of this blog, I will return to looking at R's plotting capabilities with a focus on the ggplot2 package. In the meantime, enjoy using the apply function and all it has to offer.

Multiple plots in R: lesson zero

Today, in one of my more productive days, I managed to create a sleek R script that plotted several histograms in a lattice, allowing for easy identification of the underlying trend. Although the majority of the time taken consisted of collecting the data and making various adjustments, it took a not inconsiderable amount of work to write the code.

As I was cursing the apply function – not for the last time I am sure – I suddenly realised the insane level of productivity that I have come to see as "par for the course". The level of computational analysis that can be conducted in a few hours with nothing more than a desktop PC,  a broadband connection, and copious amounts of caffeine is phenomenal. No longer can a lack of computing power or software be blamed for a lack of productivity growth: rate-limiting factors are now exclusively human.

Here at least is my contribution to the collective intelligence of biologicals. I have noticed that the most common reason that people avoid R is that they cannot rapidly make graphs that meet the high standards of their clients. In the next few editions of this blog we will build up a basic repertoire of plotting techniques, focusing on graphics that are sickeningly impressive.

Of course, Rome was not built in a day, and a thorough knowledge of R plotting cannot be built in one. Instead we will progress one layer at a time, adding additional levels of complexity and functionality. We start with a simple script that allows us to plot several graphs at the same time, each with a different value of a key variable. Here is the output:

Figure 1: Multiple plots using par

The code

# Clear all objects
rm(list=ls())

# Create a data set using random variables
df<-data.frame(x=rnorm(160),y=runif(160),a=sample(c(1,0,-1,10),
replace=TRUE,160))

df$z<-with(df,{3*a*y+x})

# Create a function that plots the value of "z" against the "y" value
plotM<-function(l){
  
  df.temp<-df[df$a==l,]
  plot(df.temp$y,df.temp$z,xlab="Y Value",ylab="Z Value",
  main=paste("Value of key variable: ",toString(l)))
  abline(lm(df.temp$z~df.temp$y),col="red")
  
}

# Create a grid to plot the different values of "a"
par(mfrow=c(2,2))

# Loop through each value of "a" and call the plotM function
for (i in c(1,0,-1,10)){
    plotM(i)
}

  
Walk-through
The first few lines of code create a data frame with 4 variables: x,y,z,a. Three of these variables are randomly generated, with the z variable dependent upon the other 3. Suppose that we are analysing this data set, unaware of the relationship between the x,y and z variables. A preliminary inspection shows that there are only 4 observed values of a. It seems sensible to plot z against y for each of these different a values.

The relevant code to create this plot starts with the function "plotM", standing for plot Multiple. This function takes an argument "l" that determines which value of the variable we are plotting. The line

df.temp<-df[df$a==l,]

filters the data frame to include only those rows where the variable a=l. The next line simply creates a standard R plot of y versus z. Finally, we use the abline function to plot a linear fit to highlight the trend. Too easy.

Now we come to the fun part. Using the par(mfrow=c(2,2)) we create a 2x2 grid in which to place the next four plots. Whatever plots we now create will be placed sequentially into this grid. Hence we can iterate over a vector containing the values of a,  calling plotM each time. This gives us our grid.

Comments
Okay, so the graph does not look like it came from NASA – or to be honest with NASA's ailing reputation maybe it does. But note that once we have created the plotM function, we only have to write 3 lines of code to make 4 separate charts. Moreover, the code would not increase even if we were plotting 100 charts.

Of course, we have not yet even drawn upon any of R's custom plotting packages. In the next edition of this blog we will look at how to use the ggplot2 package to add colour and a wide range of other features to our graphs.

Survival skills for today's analyst

I suffer a little from the age-old affliction of contrarianism. If a software package is used by the majority of the population, I assume it is flawed, highly limited, and its continued use will ultimately result in the downfall of the human race. Conversely, I am always extremely interested in a piece of software that has spread no further than the ivory tower in which it was first conceived.

The most longstanding example of this is my profound preference for the statistical computing language, R, over Microsoft Excel–a program in which I have begrudgingly developed an extremely high level of expertise. As every analyst knows, in the world of statistical software Excel is like McDonald's, Burger King, Pizza Hut, and KFC all rolled into one. It is so prepackaged and devoid of customization, yet so ubiquitous that we cannot do without it. Like the fast-food chains, we loathe Excel because it always produces the same graphs, the same simple statistical analyses. Yet when we find ourselves lost in a strange, unfriendly foreign country we go running back to the grid lines of excel to order a Big Mac. As soon as we enter the jungle, our survival skills are found wanting.

Yes my friends, like it or not, Excel is here to stay although not for lack of alternatives. The fact is that it is the user-friendly nature of this program that has been the key to its success. A friend of mine once put it thus: "Excel has allowed a generation of knowledge workers to survive without being able to program."

In truth, the driving force behind Excel's success is simple: Excel is easy. Oh I know that the die-hards will talk about how it is a superior visual tool, and that spreadsheets allow for increased transparency in financial models. But this argument falls flat on its face when we introduce macros to the equation; if spreadsheets are about transparency then why do we add VBA scripts that the user can neither see nor understand. And if we are happy to use scripts at some level, why on earth do we need to do everything else in a cumbersome visual environment.

Furthermore, the simple ends to which Excel users put their tools is demonstrated by the tiny fraction of users use the (admittedly limited) functionality afforded by VBA. That so many users can get by without loops, functions, or any notion of encapsulation is testament to the primitive uses to which Excel is put: it is just a big button calculator with an autofill feature. Surely there are more skills that we need to survive in the analytical jungle.

Finding an alternative
As I write this, I am sure that I have just alienated the entire community of so called "Power Excel Users". But I am sure that many engineers, economists, and scientists will agree that Excel is too limited to be the only quantitative tool that you have available in your office. The problem is finding software that you can successfully use in an office environment, and that is worth investing the time in learning.

The most important obstacle to overcome is the cost barrier. One of my friends, Rex, works for a major insurance group in their risk division. As far as I can tell, there are few organisations as willing to shell out money on analytical software as an insurance company. As a result, Rex regularly tells me about the wonderful software package that they just bought for $X million. These packages are highly customised and very user friendly (that's why they cost big dollars). The problem is, what happens when the company's systems change, or when you need to solve a new problem? Moreover, how does Rex do his job when he no longer has access to the software (ie, if he moves to another job). The cost of these highly customized packages means that they are not useful tools to acquire for your repertoire. As a rule of thumb, if it costs more than the latest version of Excel then assume that it is not portable: you cannot take it with you.

Enter R
Since I am a contrarian, I am sure that my advice should be taken with a grain (if not a barrel) of salt. However, I believe that there is now a viable alternative to Excel: R. R has been around for a long time, but it has taken a while to gain the following that it so rightfully deserves.

R is completely free and thus available at your fingertips wherever you go. No need to negotiate with the boss about breaking the budget for some fancy new piece of software. Download the binary, install it, and you are good to go. The advantage of this is not just that it is freely available, but that you can rely on it being available.

That just leaves its functionality, and my friends the good news is that R has functionality in spades. Take a quick look at its graphical features and you will see that almost any chart or graph you can dream of can be generated in R. In addition, the R community is continually adding new packages with new functions. In the last few years, the development of these packages has exploded in line with growth in the user base.

Transcending Excel and transitioning to R

Having used R for a reasonable amount of time, I find it hard to see why other analysts struggle day-in day-out with Excel. However, the great barrier to using R is that it is one step closer to all-out coding. Run through an interpreter, R seems strange and frightening to the non-programmer. If you have never learned a programming language, then chances are it will take you some time to shift to R.

Another issue is the need for other people to have the ability to review, check, and edit your work. Unless your boss is up to speed with R or is willing for your work to be checked by another R-literate colleague, you may have to stick with Excel for the moment.

There is, however, great scope for the analyst to grow their organisation into R over time. Whenever you are asked to do a self-contained piece of work independently, try doing it in R. I tend to go overboard and try to create advanced graphics that showcase R's capabilities. The majority of the time, people ask how I made the graph and are then keen to see what else R can do.

Into the jungle
As the old saying goes, "to the man that has only a hammer, every problem looks like a nail". At the moment, there are an awful lot of organisations that are filled with people who only have Excel and every problem sure looks like a spreadsheet.

I believe that analysts that fail to expand their toolkit tend to lose the ability to solve new problems. The generation of knowledge workers who are now in their 40s may have been lucky enough to survive on nothing more than their spreadsheet skills. However, as a twenty-something making my way in the business world, I cannot see how an analyst will be able to survive without some high-powered programming in their utility belt. R may not be enough on its own, but it seems like a good starting point.

Good luck in the analytical jungle.