###

rm(list = ls())
rm(list=lsf.str())

library(rgeos)
library(stringr)
library(reshape)
library(maptools)
library(ggplot2)
library(SmarterPoland)
library(sp)
library(grid)
library(rgdal)

# choose your working directory for shapefile download
setwd("C:/Users/Milena/Documents/R")
current.dir <- getwd()

#download the shapefile
download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip", 
               destfile="TM_WORLD_BORDERS-0.3.zip")

#unzip to SpatialPolygonsDataFrame
unzip("TM_WORLD_BORDERS-0.3.zip")

world <- readOGR(dsn = current.dir, layer = "TM_WORLD_BORDERS-0.3")

# extract from the world shapefile only the EU countries
c.eu <- c("AUT", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN",
           "FRA", "DEU", "GRC", "HUN", "IRL", "ITA", "LVA", "LTU", "LUX", 
           "MLT", "NLD", "POL", "PRT", "ROU", "SVK", "SVN", "ESP", "SWE", 
           "GBR")

vec.data.eu <- rep(0, times=length(c.eu))

for (i in 1:length(c.eu) ) 
{vec.data.eu[i] = which(world@data[, "ISO3"]==c.eu[i])}

world.eu1 <- world[vec.data.eu,]
world.eu <- world.eu1

# country centroids (middle points in each EU country)
LAT  <- c(47.33, 50.83, 43.00, 45.17, 35.00, 49.75, 56.00, 
          59.00, 64.00, 46.00, 51.50, 39.00, 47.00, 53.00, 
          42.83, 57.00, 56.00, 49.75, 35.92, 52.50, 52.00, 
          39.50, 46.00, 48.67, 46.25, 40.00, 62.00, 54.00)

LONG <- c(13.33,  4.00, 25.00, 15.50, 33.00, 15.00, 10.00,
          26.00, 26.00,  2.00, 10.50, 22.00, 20.00, -8.00, 
          12.83, 25.00, 24.00,  6.17, 14.43,  5.75, 20.00, 
          -8.00, 25.00, 19.50, 15.17, -4.00, 15.00, -4.00)

cent <- cbind(LAT,LONG)
rownames(cent) <- c.eu

# choose two countries and display them in two different colors (color A and color B)
# I create a vector of ones length equal to the number of countries (28)

world.eu <- world.eu1
a <- rep(1, length(c.eu))

# and overwrite the value for the 1st selected country with value 2 
# and the 2nd selected country with 3 

cor1 <<- 15 #Italy
cor2 <<- 21 #Poland

a[cor1] <- 2
a[cor2] <- 3

# combine the vector of levels with country names
# and call the factor column "score"

b <- cbind(c.eu, a)
dataframe2 <- data.frame(b)  
colnames(dataframe2) <-c("Country.Code", "score")

# merge score with the spatial points data frame 
# and pass it to ggplot for visualization

matched.indices.eu <- match(world.eu@data[, "ISO3"], dataframe2[, "Country.Code"])
world.eu@data <- data.frame(world.eu@data, dataframe2[matched.indices.eu, ])
world.f.eu <- fortify(world.eu, region = "ISO3")
world.m.eu <- merge(world.f.eu, world.eu@data, by.x = "id", by.y = "Country.Code")

# draw curved arrows from country A to country B in color A (an arrow head pointing at country B)
# and from country B to country A in color B (an arrow head pointing at country A)
# the size of the arrow should correspond with the size of trade flow 
# from country A to B and from B to A 
# (that is however not the object of that question and could be done later)

# since the ggplot will be encapsulated in a shiny application 
# code should draw correct arrows in every country A and country B combination
# and should execute really fast. 

# I made a trial with geom_curve but the arrows either did not appear at all 
# (even though no error message was displayed)
# or if appear that would be after a long wait

# I found function curveGrob (which worked quite fast on my machine)
# and annotation_custom to pass it to ggplot
# I thought that it would be a solution. 
# my problem is similar I also want to visualize export between countries

# all I had to change was the curvature of the first arrow to the -0.3
# and assign the arrow head of to the beginning not the end of the arrow

myCurve<-curveGrob(0, 0, 1, 1, default.units = "npc",
                   curvature = -0.3, angle = 60, ncp = 20, shape = 1,
                   square = FALSE, squareShape = 1,
                   inflect = FALSE, arrow = arrow(length = unit(0.15, "inches"), type="closed"), open = TRUE,
                   debug = FALSE,
                   name = NULL, gp = gpar(col="blue", lwd=10, lineend="round", fill="blue"), vp = NULL)

myCurve2<-curveGrob(0, 0, 1, 1, default.units = "npc",
                    curvature = 0.3, angle = 60, ncp = 20, shape = 1,
                    square = FALSE, squareShape = 1,
                    inflect = FALSE, arrow = arrow(length = unit(0.15, "inches"), type="closed", ends="first"), 
                    open = TRUE, debug = FALSE,
                    name = NULL, gp = gpar(col="magenta", lwd=2, lineend="round", fill="magenta"), vp = NULL)

ggplot(world.m.eu, aes(long, lat, group = group))+
geom_polygon(aes(fill = world.m.eu$score),show.legend=FALSE)+ 
geom_polygon(data = world.m.eu, aes(long,lat), 
             fill="NA", color = "white",  size=0.01) +
      coord_cartesian(xlim = c(-17, 37), ylim = c(34, 72))+ 
      scale_fill_manual(values = c("lightblue", "blue", "magenta")) +
      annotation_custom(grob=myCurve, xmin=cent[cor1,2], xmax=cent[cor2,2],
                      ymin=cent[cor1,1], ymax=cent[cor2,1]) +
      annotation_custom(grob=myCurve2,xmin=cent[cor2,2], xmax=cent[cor1,2],
                        ymin=cent[cor2,1], ymax=cent[cor1,1])
    
# however what worked for the combination Italy - Poland 
# and perhaps a few others 
# does not work for example for Austria - Cyprus
# or France - Greece. 
# it seems that regardless the assignment of xmin and xmax
# or ymin and ymax
# it looks like the ys are switched 
# for the lines to be drawed upwards


### Austria - Cyprus

world.eu <- world.eu1
a <- rep(1, length(c.eu))

cor1 <<- 1 #Austria
cor2 <<- 5 #Cyprus

a[cor1] <- 2
a[cor2] <- 3

b <- cbind(c.eu, a)
dataframe2 <- data.frame(b)  
colnames(dataframe2) <-c("Country.Code", "score")

matched.indices.eu <- match(world.eu@data[, "ISO3"], dataframe2[, "Country.Code"])
world.eu@data <- data.frame(world.eu@data, dataframe2[matched.indices.eu, ])
world.f.eu <- fortify(world.eu, region = "ISO3")
world.m.eu <- merge(world.f.eu, world.eu@data, by.x = "id", by.y = "Country.Code")

myCurve<-curveGrob(0, 0, 1, 1, default.units = "npc",
                   curvature = -0.3, angle = 60, ncp = 20, shape = 1,
                   square = FALSE, squareShape = 1,
                   inflect = FALSE, arrow = arrow(length = unit(0.15, "inches"), type="closed"), open = TRUE,
                   debug = FALSE,
                   name = NULL, gp = gpar(col="blue", lwd=10, lineend="round", fill="blue"), vp = NULL)

myCurve2<-curveGrob(0, 0, 1, 1, default.units = "npc",
                    curvature = 0.3, angle = 60, ncp = 20, shape = 1,
                    square = FALSE, squareShape = 1,
                    inflect = FALSE, arrow = arrow(length = unit(0.15, "inches"), type="closed", ends="first"), 
                    open = TRUE, debug = FALSE,
                    name = NULL, gp = gpar(col="magenta", lwd=2, lineend="round", fill="magenta"), vp = NULL)

ggplot(world.m.eu, aes(long, lat, group = group))+
  geom_polygon(aes(fill = world.m.eu$score),show.legend=FALSE)+ 
  geom_polygon(data = world.m.eu, aes(long,lat), 
               fill="NA", color = "white",  size=0.01) +
  coord_cartesian(xlim = c(-17, 37), ylim = c(34, 72))+ 
  scale_fill_manual(values = c("lightblue", "blue", "magenta")) +
  annotation_custom(grob=myCurve, xmin=cent[cor1,2], xmax=cent[cor2,2],
                    ymin=cent[cor1,1], ymax=cent[cor2,1]) +
  annotation_custom(grob=myCurve2,xmin=cent[cor2,2], xmax=cent[cor1,2],
                    ymin=cent[cor2,1], ymax=cent[cor1,1])


### France - Greece

world.eu <- world.eu1
a <- rep(1, length(c.eu))

cor1 <<-10 #France 
cor2 <<-12 #Greece

a[cor1] <- 2
a[cor2] <- 3

b <- cbind(c.eu, a)
dataframe2 <- data.frame(b)  
colnames(dataframe2) <-c("Country.Code", "score")

matched.indices.eu <- match(world.eu@data[, "ISO3"], dataframe2[, "Country.Code"])
world.eu@data <- data.frame(world.eu@data, dataframe2[matched.indices.eu, ])
world.f.eu <- fortify(world.eu, region = "ISO3")
world.m.eu <- merge(world.f.eu, world.eu@data, by.x = "id", by.y = "Country.Code")

myCurve<-curveGrob(0, 0, 1, 1, default.units = "npc",
                   curvature = -0.3, angle = 60, ncp = 20, shape = 1,
                   square = FALSE, squareShape = 1,
                   inflect = FALSE, arrow = arrow(length = unit(0.15, "inches"), type="closed"), open = TRUE,
                   debug = FALSE,
                   name = NULL, gp = gpar(col="blue", lwd=10, lineend="round", fill="blue"), vp = NULL)

myCurve2<-curveGrob(0, 0, 1, 1, default.units = "npc",
                    curvature = 0.3, angle = 60, ncp = 20, shape = 1,
                    square = FALSE, squareShape = 1,
                    inflect = FALSE, arrow = arrow(length = unit(0.15, "inches"), type="closed", ends="first"), 
                    open = TRUE, debug = FALSE,
                    name = NULL, gp = gpar(col="magenta", lwd=2, lineend="round", fill="magenta"), vp = NULL)

ggplot(world.m.eu, aes(long, lat, group = group))+
  geom_polygon(aes(fill = world.m.eu$score),show.legend=FALSE)+ 
  geom_polygon(data = world.m.eu, aes(long,lat), 
               fill="NA", color = "white",  size=0.01) +
  coord_cartesian(xlim = c(-17, 37), ylim = c(34, 72))+ 
  scale_fill_manual(values = c("lightblue", "blue", "magenta")) +
  annotation_custom(grob=myCurve, xmin=cent[cor1,2], xmax=cent[cor2,2],
                    ymin=cent[cor1,1], ymax=cent[cor2,1]) +
  annotation_custom(grob=myCurve2,xmin=cent[cor2,2], xmax=cent[cor1,2],
                    ymin=cent[cor2,1], ymax=cent[cor1,1])

###