Great circles, raster, sp and lattice

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.

First, let’s load the packages.

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)
flights <- read.csv("http://datasets.flowingdata.com/tuts/maparcs/flights.csv", header=TRUE, as.is=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,
                          sp=TRUE, addStartEnd=TRUE)
  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)
                      )

Finally, the boundaries. You have to download and unzip this file.

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)