Plotting Great Circle Paths

Question

I'm trying to plot some path/connection based maps but am unable to figure out how.

I see a lot of possibilities for one point based metrics (crime hotspots in London, etc. with googleVis, ggmap, etc.) but I can't find too many examples of two points based metrics (immigrations between cities, train routes, etc.) There is an example with the package geosphere, but it seems to not be available for R 3.0.2.

Ideally I would like something like this D3 example and would also like to customise the thickness, colour, etc. of the lines and the circles.

PS: I don't suppose rCharts can do this just yet, right?

Larimore

geosphere is available for R 3.0.2, if you are having troubles installing I would focus on that part first. In R geosphere would be my first port of call for this question. – Dortch

You're right. It had something to do with the package sp. I thought dependencies get updated automatically. Thanks for the heads up. – Larimore

Answer 1

Here's something I started last year and never got around to polishing properly but hopefully it should answer your question of joining points on a map with great circle lines, with the flexibility to customise lines, circles etc.

I use the (my) rworldmap package for mapping, WDI for world bank data and geosphere for the great circle lines. The aim was to plot aid flows from all donor countries to all recipient countries (one plot per donor). Below is an example plot, and below that the code. Hope it helps. Would be nice to find time to pick it up again! Andy

library(rworldmap)
library(WDI) # WORLD BANK INDICATORS

## lines of either type may obscure more than they add
##**choose line option here
addLines <- 'gc' #'none''straight' 'gc'
if ( addLines == 'gc' ) library(geosphere)

# setting background colours
oceanCol = rgb(7,0,30,maxColorValue=255) 
landCol = oceanCol 

#produces a list of indicator IDs and names as a matrix
indicatorList <- WDIsearch('aid flows')

#setting up a world map shaped plot window
#*beware this is windows specific
mapDevice('windows',width=10,height=4.5)


year <- 2000
#for(indNum in 1:2)
for(indNum in 1:nrow(indicatorList))
{
  indID <- indicatorList[indNum][1]
  donorISO3 <- substr(indID,start=8,stop=10)

  dFdonor <- WDI(indicator=indID,start=year,end=year)
  #divide by 10^6 for million dollars
  dFdonor[indID] <- dFdonor[indID] * 1/10^6

  sPDFdonor <- joinCountryData2Map(dFdonor,nameJoinColumn='country',joinCode='NAME')
  #take out Antarctica
  sPDFdonor <- sPDFdonor[-which(row.names(sPDFdonor)=='Antarctica'),]

  legendTitle=paste("aid flow from",donorISO3,year,"(millions US$)") 
  mapBubbles(sPDFdonor, nameZSize=indID, plotZeroVals=FALSE, legendHoriz=TRUE, legendPos="bottom", fill=FALSE, legendTitle=legendTitle, oceanCol=oceanCol, landCol=landCol,borderCol=rgb(50,50,50,maxColorValue=255),lwd=0.5,lwdSymbols=1)
  #removed because not working , main=paste('donor', donorISO3,year)

  #now can I plot lines from the centroid of the donor to the centroids of the recipients
  xDonor <- sPDFdonor$LON[ which(sPDFdonor$ISO3==donorISO3) ]
  yDonor <- sPDFdonor$LAT[ which(sPDFdonor$ISO3==donorISO3) ] 
  xRecips <- sPDFdonor$LON[ which(sPDFdonor[[indID]] > 0) ]
  yRecips <- sPDFdonor$LAT[ which(sPDFdonor[[indID]] > 0) ]
  amountRecips <- sPDFdonor[[indID]][ which(sPDFdonor[[indID]] > 0) ]


  ## straight lines
  if ( addLines == 'straight' )
  {
    for(line in 1:length(xRecips))
    {  
       #col <- 'blue'
       #i could modify the colour of the lines by the size of the donation
       #col=rgb(1,1,1,alpha=amountRecips[line]/max(amountRecips))
       #moving up lower values
       col=rgb(1,1,0,alpha=sqrt(amountRecips[line])/sqrt(max(amountRecips)))
       lines(x=c(xDonor,xRecips[line]),y=c(yDonor,yRecips[line]),col=col, lty="dotted", lwd=0.5)   #lty = "dashed", "dotted", "dotdash", "longdash", lwd some devices support <1
    }
  }

  ## great circle lines
  ## don't work well when donor not centred in the map
  ## also the loop fails at CEC & TOT because not ISO3 codes
  if ( addLines == 'gc' & donorISO3 != "CEC" & donorISO3 != "TOT" )
  {  
    for(line in 1:length(xRecips))
    {
      #gC <- gcIntermediate(c(xDonor,yDonor),c(xRecips[line],yRecips[line]), n=50, breakAtDateLine=TRUE)
      #30/10/13 lines command failed with Error in xy.coords(x, y) : 
      #'x' is a list, but does not have components 'x' and 'y'
      #adding sp=TRUE solved
      gC <- gcIntermediate(c(xDonor,yDonor),c(xRecips[line],yRecips[line]), n=50, breakAtDateLine=TRUE, sp=TRUE)

      #i could modify the colour of the lines by the size of the donation
      #col=rgb(1,1,1,alpha=amountRecips[line]/max(amountRecips))
      #moving up lower values
      col=rgb(1,1,0,alpha=sqrt(amountRecips[line])/sqrt(max(amountRecips)))

      lines(gC,col=col,lwd=0.5)
    }
  }  

  #adding coasts in blue looks nice but may distract
  data(coastsCoarse)
  plot(coastsCoarse,add=TRUE,col='blue')

  #repeating mapBubbles with add=T to go on top of the lines
  mapBubbles(sPDFdonor, nameZSize=indID, plotZeroVals=FALSE, fill=FALSE, addLegend=FALSE, add=TRUE, ,lwd=2)
  #removed because not working : , main=paste('donor', donorISO3,year)

  #looking at adding country labels
  text(xRecips,yRecips,sPDFdonor$NAME[ which(sPDFdonor[[indID]] > 0) ],col=rgb(1,1,1,alpha=0.3),cex=0.6,pos=4) #pos=4 right (1=b,2=l,3=ab)

  #add a title 
  nameDonor <- sPDFdonor$NAME[ which(sPDFdonor$ISO3==donorISO3) ]
  mtext(paste("Aid flow from",nameDonor,year), cex = 1.8, line=-0.8)

  #savePlot(paste("C:\\rProjects\\aidVisCompetition2012\\Rplots\\greatCircles\\wdiAidFlowLinesDonor",donorISO3,year,sep=''),type='png')
  #savePlot(paste("C:\\rProjects\\aidVisCompetition2012\\Rplots\\greatCircles\\wdiAidFlowLinesDonor",donorISO3,year,sep=''),type='pdf')

} #end of indNum loop

Protonema

This looks nice, good stuff. I'm not sure how to check but is the underlying map image vectorised? geosphere seems to not encourage zooming in/out. Also, since it looks like the package isn't maintained actively, I shall only upvote and not accept the answer. I hope the 15 points encourage you to take it up again :) – Larimore

Thanks! Yes the rworldmap map is vector (an sp SpatialPolygonsDataFrame). You can tell this by saving the plot as a pdf (uncomment the penultimate line of the code and edit the path), then zoom in and you'll see lines not pixels. – Protonema