# 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,
[email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */[[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([email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */[[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: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Tags: , , , , , ,