Create Flow Map in R using ggplot2

Flow Map in R

## This gist shows how to create Flow Maps in R using ggplot2. 
## source: This is based on different bits of code from other with amazing R skills: 

@ceng_l          : http://web.stanford.edu/~cengel/cgi-bin/anthrospace/great-circles-on-a-recentered-worldmap-in-ggplot
@3wen            : http://egallic.fr/maps-with-r/
@spatialanalysis : http://spatialanalysis.co.uk/2012/06/mapping-worlds-biggest-airlines/
@freakonometrics : http://freakonometrics.hypotheses.org/48184


# Libraries
  library(maps)
  library(geosphere)
  library(dplyr)
  library(ggplot2)
  library(rworldmap)
  library(plyr)
  library(data.table)
  library(ggthemes)


# Get World map
  worldMap <- getMap()
  mapworld_df <- fortify( worldMap )


# Read data on airports and flights
  airports <- read.csv("http://www.stanford.edu/~cengel/cgi-bin/anthrospace/wp-content/uploads/2012/03/airports.csv", as.is=TRUE, header=TRUE)
  flights <- read.csv("http://www.stanford.edu/~cengel/cgi-bin/anthrospace/wp-content/uploads/2012/03/PEK-openflights-export-2012-03-19.csv", as.is=TRUE, header=TRUE)

# get airport locations
  airport_locations <- airports[, c("IATA","longitude", "latitude")]

# aggregate number of flights (frequency of flights per pair)
  flights.ag <- ddply(flights, c("From","To"), function(x) count(x$To))


# Link airport lat  long to origin and destination
  OD <- left_join(flights.ag, airport_locations, by=c("From"="IATA") )
  OD <- left_join(OD, airport_locations, by=c("To"="IATA") )
  OD$id <-as.character(c(1:nrow(OD))) #create and id for each pair



##### Two Simple Maps ##### 

# 1. Using straight lines
  ggplot() + 
    geom_polygon(data= mapworld_df, aes(long,lat, group=group), fill="gray30") +
    geom_segment(data = OD, aes(x = longitude.x, y = latitude.x, xend = longitude.y, yend = latitude.y, color=freq),
                 arrow = arrow(length = unit(0.01, "npc"))) +
    scale_colour_distiller(palette="Reds", name="Frequency", guide = "colorbar") +
    coord_equal()


# 2. Using Curved Lines
  ggplot() + 
    geom_polygon(data= mapworld_df, aes(long,lat, group=group), fill="gray30") +
    geom_curve(data = OD, aes(x = longitude.x, y = latitude.x, xend = longitude.y, yend = latitude.y, color=freq),
               curvature = -0.2, arrow = arrow(length = unit(0.01, "npc"))) +
    scale_colour_distiller(palette="Reds", name="Frequency", guide = "colorbar") +
    coord_equal()



##### A more professional map ####  
# Using shortest route between airports considering the spherical curvature of the planet

# get location of Origin and destinations airports
  setDT(OD) # set OD as a data.table for faster data manipulation
  beijing.loc <- OD[ From== "PEK", .(longitude.x, latitude.x)][1] # Origin
  dest.loc <- OD[ , .(longitude.y, latitude.y)] # Destinations

# calculate routes between Beijing (origin) and other airports (destinations)
  routes <- gcIntermediate(beijing.loc, dest.loc, 100, breakAtDateLine=FALSE, addStartEnd=TRUE, sp=TRUE)
  class(routes) # SpatialLines object


# Convert a SpatialLines object into SpatialLinesDataFrame, so we can fortify and use it in ggplot
  # create empty data frate  
  ids <- data.frame()
  # fill data frame with IDs for each line
  for (i in (1:length(routes))) {         
    id <- data.frame(routes@lines[[i]]@ID)
    ids <- rbind(ids, id)  }

  colnames(ids)[1] <- "ID" # rename ID column

# convert SpatialLines into SpatialLinesDataFrame using IDs as the data frame
  routes <- SpatialLinesDataFrame(routes, data = ids, match.ID = T)

# Fortify routes (convert to data frame)  +++  join attributes
  routes_df <- fortify(routes, region= "ID") # convert into something ggplot can plot
  gcircles <- left_join(routes_df, OD, by= ("id"))
  head(gcircles)

### Recenter ####

center <- 115 # positive values only - US centered view is 260

# shift coordinates to recenter great circles
  gcircles$long.recenter <-  ifelse(gcircles$long  < center - 180 , gcircles$long + 360, gcircles$long) 

# shift coordinates to recenter worldmap
  worldmap <- map_data ("world")
  worldmap$long.recenter <-  ifelse(worldmap$long  < center - 180 , worldmap$long + 360, worldmap$long)

### Function to regroup split lines and polygons
# takes dataframe, column with long and unique group variable, returns df with added column named group.regroup
  RegroupElements <- function(df, longcol, idcol){  
    g <- rep(1, length(df[,longcol]))
    if (diff(range(df[,longcol])) > 300) {          # check if longitude within group differs more than 300 deg, ie if element was split
      d <- df[,longcol] > mean(range(df[,longcol])) # we use the mean to help us separate the extreme values
      g[!d] <- 1     # some marker for parts that stay in place (we cheat here a little, as we do not take into account concave polygons)
      g[d] <- 2      # parts that are moved
    }
    g <-  paste(df[, idcol], g, sep=".") # attach to id to create unique group variable for the dataset
    df$group.regroup <- g
    df
  }

### Function to close regrouped polygons
# takes dataframe, checks if 1st and last longitude value are the same, if not, inserts first as last and reassigns order variable
  ClosePolygons <- function(df, longcol, ordercol){
    if (df[1,longcol] != df[nrow(df),longcol]) {
      tmp <- df[1,]
      df <- rbind(df,tmp)
    }
    o <- c(1: nrow(df))  # rassign the order variable
    df[,ordercol] <- o
    df
  }

# now regroup
  gcircles.rg <- ddply(gcircles, .(id), RegroupElements, "long.recenter", "id")
  worldmap.rg <- ddply(worldmap, .(group), RegroupElements, "long.recenter", "group")

# close polys
  worldmap.cp <- ddply(worldmap.rg, .(group.regroup), ClosePolygons, "long.recenter", "order")  # use the new grouping var



# Flat map
  ggplot() +
    geom_polygon(data=worldmap.cp, aes(long.recenter,lat,group=group.regroup), size = 0.2, fill="#f9f9f9", color = "grey65") +
    geom_line(data= gcircles.rg, aes(long.recenter,lat,group=group.regroup, color=freq), size=0.4, alpha= 0.5) +
    scale_colour_distiller(palette="Reds", name="Frequency", guide = "colorbar") +
    theme_map()+
    ylim(-60, 90) +
    coord_equal()


# Spherical Map
  ggplot() +
    geom_polygon(data=worldmap.cp, aes(long.recenter,lat,group=group.regroup), size = 0.2, fill="#f9f9f9", color = "grey65") +
    geom_line(data= gcircles.rg, aes(long.recenter,lat,group=group.regroup, color=freq), size=0.4, alpha= 0.5) +
    scale_colour_distiller(palette="Reds", name="Frequency", guide = "colorbar") +
    # Spherical element
    scale_y_continuous(breaks = (-2:2) * 30) +
    scale_x_continuous(breaks = (-4:4) * 45) +
    coord_map("ortho", orientation=c(61, 90, 0)) 


# Any ideas on how to color the oceans ? :)

As a crude workaround, this works to create set of 'ocean' polygons to plot behind the world map (excuse the clunky code!):

Create set of 'ocean' polygons of a certain longitudinal width:

width = 5 # Set string width in degrees
n = 360/width
ocean = data.frame(long = numeric(n*5), lat = numeric(n*5), group = numeric(n*5), order = integer(n*5), region = character(n*5), stringsAsFactors = FALSE)

ocean$group = rep(1:n, each = 5)
ocean$order = rep(1:5, n)
ocean$region = paste("Ocean",ocean$group, sep = "")
ocean$long = c(-180, -180, -180+width, -180+width, -180) + (ocean$group-1)*width
ocean$lat = rep(c(-89.9, 89.9, 89.9, -89.9, -89.9), n) 

Recentre and regroup this set of polygons as with the world map and it paints the map blue in n strips (making the strips narrow, say 5 degrees) avoids overlap issues at the 'back' of the spherical view...

ocean$long.recenter <- ifelse(ocean$long < center - 180 , ocean$long + 360, ocean$long) ocean.rg <- ddply(ocean, .(group), RegroupElements, "long.recenter", "group") ocean.cp <- ddply(ocean.rg, .(group.regroup), ClosePolygons, "long.recenter", "order")

And add to the ggplot code:

ggplot() + geom_polygon(data=ocean.cp, aes(long.recenter,lat,group=group.regroup), size = 0.2, fill="blue", colour = "blue") + geom_polygon(data=worldmap.cp, aes(long.recenter,lat,group=group.regroup), size = 0.2, fill="#f9f9f9", color = "grey65") + geom_line(data= gcircles.rg, aes(long.recenter,lat,group=group.regroup, color=freq), size=0.4, alpha= 0.5) + scale_colour_distiller(palette="Reds", name="Frequency", guide = "colorbar") + scale_y_continuous(breaks = (-2:2) * 30) + scale_x_continuous(breaks = (-4:4) * 45) + coord_map("ortho", orientation=c(61, 90, 0))

And thanks for the excellent script (and to those whose work you drew on)!

D

blue <- rgb(190/255, 232/255, 255/255)

then simply add at the end of your ggplot()

theme(panel.background = element_rect(fill = blue))

for example:

ggplot() +
  geom_polygon(data=worldmap.cp, aes(long.recenter,lat,group=group.regroup), size = 0.2, fill="#f9f9f9", color = "grey65") +
  geom_line(data= gcircles.rg, aes(long.recenter,lat,group=group.regroup, color=freq), size=0.4, alpha= 0.5) +
  scale_colour_distiller(palette="Reds", name="Frequency", guide = "colorbar") +
  theme_map()+
  ylim(-60, 90) +
  coord_equal()+
  theme(panel.background = element_rect(fill = blue))