# Great circles, raster, sp and lattice

(This article was first published on Omnia sunt Communia! » R-english, and kindly contributed to R-bloggers)

Recently I found a post at FlowingData with a detailed tutorial to map connections with great circles with R. The tutorial of FlowingData is excellent, but I feel more comfortable with the sp classes and methods, and with the lattice and latticeExtra packages. Besides, I want to use the free spatial data available from the DIVA-GIS project, from the developer of the raster and the geosphere packages.

Here is what I have got.

library(lattice)
library(latticeExtra)
library(maps)
library(geosphere)
library(sp)
library(maptools)
library(raster)


Now it’s time to get the data. First, airports and flights:

airports <- read.csv("http://datasets.flowingdata.com/tuts/maparcs/airports.csv", header=TRUE)


With this information, following the code from the FlowingData post, let’s define a list of SpatialLines. The makeLines function defines a SpatialLines for each connection, and stores the number of flights in the ID slot. The linesAA is a list with the result of makeLines over fsub with the apply function. Previously, the fsub data.frame has been ordered by cnt, as the FlowingData post teaches.

makeLines <- function(x){
air1 <- airports[airports$iata == x[2],]##start air2 <- airports[airports$iata == x[3],]##end
inter <- gcIntermediate(c(air1[1,]$long, air1[1,]$lat),
c(air2[1,]$long, air2[1,]$lat),
n=100,
inter@lines[[1]]@ID <- as.character(x[4])##cnt
inter
}

fsub <- flights[flights$airline == "AA",] fsub <- fsub[order(fsub$cnt),]
linesAA <- apply(fsub, 1, makeLines)


Each component of the list can be plot with the sp.lines function. The colIndex function assigns a color from a palette to the values of cnt stored in the ID slot.

colIndex <- function(x, cnt, palette)palette[match(x, cnt)]
makePlot <- function(x, cnt, palette){
idx <- as.numeric(x@lines[[1]]@ID)
sp.lines(x, col.line=colIndex(idx, cnt, palette))
}


I define the palette for the list of lines with a combination of colorRampPalette, brewer.pal and level.colors, following the example of the help file of this last function.

cnt <- fsub\$cnt
breaks <- do.breaks(range(cnt), 30)
palLines <- colorRampPalette(brewer.pal('Greens', n=9))
colors <-level.colors(cnt,
at = breaks,
col.regions = palLines)


Let’s get the elevation data. You have to download and unzip this file (you can use the getData function from the raster package, but I am having problems with it.). The next code builds a sp object with this information:

elevUSA <- raster('USA_alt/USA1_msk_alt.grd')
prj <- CRS("+proj=longlat +datum=WGS84")
projection(elevUSA) <- prj
elevUSAsp <- as(sampleRegular(elevUSA, 5e5, asRaster=TRUE), 'SpatialGridDataFrame')


The altitude map is constructed with the spplot function over elevUSAsp. I use the combination of colorRampPalette and brewer.pal to get the palette.

palMap=colorRampPalette(brewer.pal('Greys', n=9))
altitudeMap <- spplot(elevUSAsp,
col.regions=palMap,
colorkey=list(space='bottom'),
scales=list(draw=TRUE)
)


mapUSA <- readShapeLines('USA_adm/USA_adm0.shp', proj4string=prj)


Now it’s time for joining all together. I use the layer function from latticeExtra to draw the boundaries (with sp.lines) and the list of lines (with lapply and makePlot).

p <- altitudeMap +
layer(sp.lines(mapUSA, col.line='black')) +
layer(lapply(linesAA, makePlot, cnt, colors))


Here it is! (click on the image for a high resolution PDF file)

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...