rm(list=ls())
library(Hmisc)
library(RColorBrewer)
library(circlize)
library(data.table)
library(stargazer)
library(maptools)
library(classInt)

setwd('~/research/gTiger/Article_citsci/open_code')
source('settings/parameters.r')
source('settings/paths.r')
source('settings/functions.r')

setwd(path_to_data)

APS_ods_dos_interprov = readRDS('APS_ods_dos_interprov.Rds')

col_df = data.table(home_prov = unique(APS_ods_dos_interprov$home_prov), home_prov_color = c(brewer.pal(8,'Dark2'), brewer.pal(8,'Accent'), brewer.pal(8,'Set1'), brewer.pal(8,'Set2'), brewer.pal(8,'Paired'),brewer.pal(8,'Set3')))

APS_ods_dos_interprov = merge(APS_ods_dos_interprov, col_df, by='home_prov', all.x=T, all.y=T)

d = APS_ods_dos_interprov[ , list(home_prov,work_prov, transport_risk_od)]

make_aps_circle_plot = function(){
transp = .4
chordDiagram(d, annotationTrack = "grid", preAllocateTracks = list(track.height = 0.1), col=d$home_prov_color, grid.col=col_df$home_prov_color, order=col_df$home_prov, transparency=transp, directional=0)
circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
xlim = get.cell.meta.data("xlim")
xplot = get.cell.meta.data("xplot")
ylim = get.cell.meta.data("ylim")
sector.name = get.cell.meta.data("sector.index")
if(abs(xplot[2] - xplot[1]) < 20) {
circos.text(cex=.5, mean(xlim), ylim[1], sector.name, facing = "clockwise",
niceFacing = TRUE, adj = c(0, 0.5))
} else {
circos.text(cex=1, mean(xlim), ylim[1], sector.name, facing = "inside",
niceFacing = TRUE, adj = c(0.5, 0))
}
}, bg.border = NA)
}

pdf(paste0(path_to_figures, 'commuting_risk.pdf'))
make_aps_circle_plot()
dev.off()

png(paste0(path_to_figures, 'commuting_risk.png'))
make_aps_circle_plot()
dev.off()

svg(paste0(path_to_figures, 'commuting_risk.svg'))
make_aps_circle_plot()
dev.off()

source_sink = readRDS('source_sink.Rds')

provs = readShapePoly(paste0(path_to_external_data, 'mapas_completo_municipal/spain_provinces_ind_2.shp'))
proj4string(provs) = "+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
provs = provs[which(!is.na(provs@data$PROV) & !provs@data$PROV %in% c(35, 38, 51, 52)),]
provs_utm = spTransform(provs, PROJ4_UTM_SPAIN)
provs_utm@data$PROV = as.integer(provs_utm@data$PROV)


make_source_sink_map = function(){
chloro = merge(provs_utm, source_sink, by.x='PROV', by.y='o', all=T)
v = chloro@data$net
describe(v)
v = v/1000
this_pal = c(brewer.pal(3, "Blues")[2:1], brewer.pal(4, "Reds"))
these_breaks = c(-4,-2, 0, 2, 4, 6, 8)
nclass = classIntervals(v,style="fixed", fixedBreaks=these_breaks)
colcode = findColours(nclass, this_pal)
plot(provs_utm, col=colcode)
legend('bottomright', legend=names(attr(colcode, 'table')), fill=this_pal, bty="n", inset=c(.01, .15), cex=.6)
}


pdf(paste0(path_to_figures, 'source_since_map.pdf'))
par(mar=c(0,0,0,0))
make_source_sink_map()
dev.off()


sst = readRDS('APS_source_sink_table.Rds')

stargazer(sst[order(-sst[,4, with=F])], type="html", summary=F, digits=0, rownames=F, out=paste0(path_to_figures, "tab3.html"))

names(sst) = c("p", "o", "i", "n")

sst = sst[order(n)]

sst$pn = factor(1:nrow(sst))

sst

ggplot() + geom_bar(data=sst, aes(x=pn, y=o), stat="identity", fill="#ffffff", colour="black") + geom_bar(data=sst, aes(x=pn, y=-i), stat="identity", fill="#ffffff", colour="black") + geom_bar(data=sst, aes(x=pn, y=n), stat="identity") + theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,vjust=.5, colour = "black")) + scale_x_discrete(name="",labels=sst$p) + ylab("Inter-Province Tiger Mosquito Transfers (N)") + geom_abline(intercept=0, slope=0, colour="red") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
ggsave(paste0(path_to_figures, 'source_sink_bars.png'), height=8, width=8, dpi=1000)
ggsave(paste0(path_to_figures, 'source_sink_bars.svg'), height=8, width=8, dpi=1000)
