Monday, August 3, 2009

Maps in R

The maptools package in R allows plotting data in a map --- besides keeping your work in only one environment, it saves you the like 10 minutes it takes for ArcGIS just to open, and the 10 minutes that it takes to render the map each time you modify something.

I have here a little example plotting percapita gdp growth in the world (data from the World Band Development Indicators Online). The final map looks like this (the code is below):




## read gdp growth in 2007 originally taken from the World Bank Development Indicators On-line:
dat<-read.table("http://web.ics.purdue.edu/~nvillori/blog/gdpgrowth2007.txt",header=TRUE)
## Load package maptools --- wrld_simpl is a polygon map of the world included in this package.
require(maptools)
data(wrld_simpl)
## Index ctry.names by position in the map.
dat$pos<-match(dat$Country.Code,unique(wrld_simpl$ISO3))
dat<-dat[order(dat$pos),]
## Classes & Palettes: Use packages RColorBrewer for cool palettes and classInt for easy sectioning of the values to be mapped:
### Variable to be plotted
plotvar <- dat$YR2007
### Palettes and intervals:
require(RColorBrewer)
nclr<-9
plotclr <- brewer.pal(nclr,"YlOrRd")
require(classInt)
class <- classIntervals(plotvar, nclr, style="pretty")
colcode <- findColours(class, plotclr)
### Define margins and other graphic parameters:
par(mar=c(1,1,1,1), lheight = .8)
### Plot an "empty" map with borders in grey --- aspect ratio is modified to distort image a bit:
plot(wrld_simpl,axes=TRUE,border="darkgrey",xlim=c(-110,150),bg="lightcyan", ylim=c(-45,60),asp=1.5,main="GDP Growth in 2007")
### Plot the gdp growth rates:
plot(wrld_simpl[wrld_simpl$ISO3 %in% dat$Country.Code,],col=colcode,axes=TRUE,add=TRUE)
### Add country codes to the map:
text(wrld_simpl$LON[wrld_simpl$ISO3 %in% dat$Country.Code],wrld_simpl$LAT[wrld_simpl$ISO3 %in% dat$Country.Code],wrld_simpl$ISO3[wrld_simpl$ISO3 %in% dat$Country.Code],cex=.4,col="black")
### Add a legend:
legend(title="GDP growth",-45, -15, legend=names(attr(colcode, "table")), fill=attr(colcode, "palette"), cex=.6, bty="y",bg="lightgrey")


2 comments:

  1. What about putting the code in the comments like we said??

    ReplyDelete
  2. Holy rules of the blog! I was not aware of it... in any case, if the post is about code, it seems better to have it there. I have added a little scroll-down box in which you can put the code and thus not take all the space. I'll put a post about the HTML code soon so we can use it.

    ReplyDelete