[igraph] Maps and Networks

[Lorenzo Isella]

Dear All,
Apologies for posting this to two different mailing lists, but I really  
need some cross expertise.
I am trying to do the following: plot a geographical map colored according  
to a scalar value (choropleth map) and to overimpose a network on the map.
To fix the ideas, you can imagine that I am plotting the GDP in different  
countries (choropleth map) and the flow of payments among them (here is  
where the network kicks in).
Please consider the following snippet

##################################################################################à
library(igraph)
library(maps)


t <- map("italy", plot=FALSE)
pos1 <- c(16,40.3)
pos2 <- c(12,44)
pos3 <- c(13,42)

el <- matrix( c("foo", "bar", "bar", "foobar"), nc=2, byrow=TRUE)
g <- graph.edgelist(el)

g <- as.undirected(g)

l <- matrix(ncol=2, nrow=3)

l[1,] <- pos1
l[2,] <- pos2
l[3,] <- pos3

pdf("italy.pdf")
plot(t$x, t$y, "l", xaxt='n', yaxt='n',ann=FALSE,  axes=F)

autocurve.edges (g, start = 0.5)

plot(g, add=TRUE, rescale=FALSE, layout=l,
     vertex.label.dist=0.5,vertex.label=NA, vertex.size=15
     )


dev.off()

##############################################################################

Colors aside, it gives the idea of what I am looking for.
At the end of the email, you find a rather longish script, based on  
ggplot2, that plots a map of Europe with several countries colored  
according to a scalar.
You need to download the world map from

http://www.mappinghacks.com/data/TM_WORLD_BORDERS_SIMPL-0.2.zip


At this stage, I have two problems
(1) I cannot manage to combine the previous network with the ggplot object  
I generate and plot at the end of the script.
(2) I would like to specify the nodes' position by passing their latitude  
and longitude instead of some "coordinates" relative to the size of the  
plot.

Sorry for the long post and any help is appreciated.
Cheers

Lorenzo


##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################

library(maptools)
library(reshape)
library(ggplot2)
## require(gpclib)
library(RColorBrewer)
library(grid)



my_ggplot_theme <- function(legend_coord){  theme( panel.background =  
element_rect(fill="gray",
                          colour = "black", size = 0.5, linetype = 1),
                panel.grid.minor = element_blank(),
                axis.ticks = element_line(colour = "black", size=1),
                axis.ticks.length = unit(0.15, "cm"),
                strip.background = element_rect(colour = 'blue',
                           fill = 'white', size = 1, linetype=1),
                strip.text.x = element_text(colour = 'red', angle = 0,
                                     size = 12, hjust = 0.5,
                                     vjust = 0.5, face = 'bold'),
                axis.title.x = element_text(size = 20),
                axis.title.y = element_text(size = 20, angle=90, vjust=.4),
                axis.text.x = element_text(size=15, colour="black",  
vjust=1),
                axis.text.y = element_text(size=15, colour="black",  
hjust=1),
                legend.text = element_text(size = 14, vjust=0.4),
                legend.title = element_text(size = 24, hjust=0),
                legend.position = legend_coord,
                legend.title = element_blank(),
                legend.background=element_rect(color=NA, fill=NA),
                legend.key = element_rect(colour = NA, fill=NA) )
                                                 }






mono_scale_discrete <- function(col_name, n_cols){

res <- brewer.pal(n_cols, col_name)

return(res)

}


double_scale_discrete <- function(col_name1, col_name2,n_cols_tot,  
n_cols1){

res1 <- brewer.pal(n_cols_tot, col_name1)

n_cols2 <- n_cols_tot-n_cols1

res2 <- brewer.pal(n_cols_tot, col_name2)

res <- c(res1[n_cols_tot:(n_cols_tot-n_cols1+1)],
         rev(res2[n_cols_tot:(n_cols_tot-n_cols2+1)]))

return(res)

}

double_scale_discrete2 <- function(col_name1, col_name2,n_cols_tot,  
n_cols1){

res1 <- rev(brewer.pal(n_cols1, col_name1))

n_cols2 <- n_cols_tot-n_cols1

res2 <- brewer.pal(n_cols2, col_name2)

## res <- c(res1[n_cols_tot:(n_cols_tot-n_cols1+1)],
##          res2[(n_cols2+1):n_cols_tot])


res <- c(res1, res2)

return(res)

}



col_gen <- function(n) {
  black <- "#000000"
  if (n <= 9) {
    c(black,brewer.pal(n-1, "Set2"))
  } else {
    c(black,hcl(h=seq(0,(n-2)/(n-1),
                  length=n-1)*360,c=100,l=65,fixup=TRUE))
  }
}

##############################################

#Read and set the world map
worldmap <- readShapeSpatial("TM_WORLD_BORDERS_SIMPL-0.3.shp")


worldmap <- fortify(worldmap, region="ISO2")


################################################


country_list <- c("EU27", "BE", "BG","CZ", "DK", "DE", "EE", "IE", "GR",
                   "ES", "FR", "IT", "CY", "LV", "LT", "LU", "HU", "MT",
                   "NL", "AT","PL", "PT", "RO", "SI", "SK", "FI","SE","GB"  
)



set.seed(1234)

t <- runif(length(country_list), 40, 280)

prc_ppp_ind <- cbind(as.data.frame(country_list),
                     as.data.frame(t))

names(prc_ppp_ind) <- c("geo", "value")


prc_ppp_ind$geo <- as.character(prc_ppp_ind$geo)

len <- length(worldmap$id)

worldmap$value <- rep(NA, len)

for (i in seq(length(country_list))){

country <- country_list[i]

sel <- which(worldmap$id==country)

sel2 <- which(prc_ppp_ind$geo==country)

my_val <- prc_ppp_ind$value[sel2]



worldmap$value[sel] <- rep(my_val,length(sel)) #regions
#are at 3 levels, and I need to make sure that I am giving the
#same value to every country as a whole (not just different values
# to different areas of the same country.)

}



############################################################
## Do the plotting

#NB: remove the constraints on latitude and longitude to plot
# the whole world!!!

latlimits <- c(30, 75)
longlimits <- c(-15, 35)

break_seq <- c(40,  60,  80, 100, 120,140,240, 260, 280)

n_cols <- length(break_seq)-1




my.cols <- double_scale_discrete("Reds", "Blues",n_cols, 3)


gpl <- ggplot(worldmap, aes(long,lat,group=group)) +
geom_polygon(aes(fill=cut(value,break_seq))) +

    geom_polygon(data = worldmap, aes(long,lat),
               fill= NA,
               color = "black",
               size=0.1) + # white borders
    scale_fill_manual("PPP (%)",values = my.cols,na.value =  "grey85")+

labs(title = expression(paste("Fictitious Data")))+
       theme(plot.title = element_text(size = 20))+
    coord_map("ortho", orientation=c(41.9, 12.5, 0),
              xlim = longlimits, ylim = latlimits) +
labs(x=" ", y=" ")+
theme(axis.ticks = element_blank(), axis.text.x = element_blank(),  
axis.text.y = element_blank())

pdf("map_colored_EU.pdf"## ,width=5,height=5
    )
plot(gpl)
dev.off()