Plotting Great Circle Paths

2019-02-07 14:32发布

问题:

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?

回答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