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.

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)


To leave a comment for the author, please follow the link and comment on his blog: Omnia sunt Communia! » R-english.

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...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , , , ,

Comments are closed.