rm(list=ls())
library(foreign)
library(Hmisc)
library(data.table)
library(raster)
library(parallel)
library(sp)
library(RColorBrewer)
library(geosphere)
library(rgeos)
library(classInt)
library(GISTools)
library(Cairo)

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

setwd(path_to_data)

ncores=20

gadm_prov = readRDS(paste0(path_to_external_data, "ESP_adm2.rds"))
gadm_prov_utm = spTransform(gadm_prov, PROJ4_UTM_SPAIN)
# gadm_mun = readShapePoly(paste0(path_to_external_data, 'mapas_completo_municipal/esp_muni_00.shp'))

gadm_mun = readShapePoly(paste0(path_to_external_data, 'municipality_shapeflie/Munis_IGN/Munis_IGN.shp'))

proj4string(gadm_mun) = "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs"

gadm_mun@data$PROVMUN = as.character(gadm_mun@data$CODIGOINE)

# getting rid of one duplicate row - a small polygon in Alicante
gadm_mun = gadm_mun[which(gadm_mun$ID!=153),]

# proj4string(gadm_mun) = "+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
gadm_mun_utm = spTransform(gadm_mun, PROJ4_UTM_SPAIN)

spain_gpdf = readRDS("spain_grid.Rds")
spain_gpdf = spTransform(spain_gpdf, PROJ4_UTM_SPAIN)

muni_daily_probs = readRDS("muni_daily_probs.Rds")

this_pal = c(rgb(255, 255, 190, maxColorValue=255), brewer.pal(9, 'YlOrRd')[3:5], hsv(0,1,.42))

# plot(1:length(this_pal), axes = FALSE, frame.plot = TRUE,xlab = "", ylab = "", pch = 21, cex = 3, bg = this_pal)
# dev.off()

gadm_prov_utm = gadm_prov_utm[which(!gadm_prov_utm@data$NAME_1 %in% c('Islas Canarias', 'Ceuta y Melilla')),]

gadm_mun_utm$prov = substring(gadm_mun_utm$PROVMUN, 1,2)

muni_daily_probs$prov = substring(muni_daily_probs$PROVMUN, 1,2)

##############################
# Daily prob animation files #
##############################
make_alert_prob_muni_map_animation = function(this_prov, this_spdf, these_values, these_breaks, this_day){
	par(mar=c(1,1,3,1))
	nclass = classIntervals(these_values,style="fixed", fixedBreaks=these_breaks)
	colcode = findColours(nclass, this_pal)
	plot(this_prov, col="grey")
	plot(this_spdf, col=colcode, pch=16, cex=.1, border=NA, add=T)
	plot(this_prov,  add=TRUE, border="white", lwd=5)
	plot(this_prov,  add=TRUE, lwd=3)
	par(new=TRUE, fig = c(0.5, .9, .75, 1))
	the_months = c("Jan.", "Feb.", "Mar.", "Apr.", "May", "Jun.", "Jul.", "Aug", "Sep.", "Oct.", "Nov", "Dec")[order(12:1)]
	pie(rep(1/12, 12), labels=the_months, col="white", border=F)
	par(new=TRUE, fig = c(0.5, .9, .75, 1))
	these_cols = rep("dark grey", 365)
	these_cols[(365-(this_day)):(365-(this_day+14))] = "white"
	pie(rep(1/365, 365), labels="", col=these_cols, border=F)
}

i = min(muni_daily_probs$sea_days)
dir.create(paste0(path_to_untracked_output, '/daily_alert_prob'), showWarning=F)
mclapply(min(muni_daily_probs$sea_day):max(muni_daily_probs$sea_day), function(i){
	print(i)
	flush.console()
	merged_gadm_mun = merge(gadm_mun_utm, muni_daily_probs[sea_days==i], by=c("PROVMUN", "prov"), all=T)
	this_date = as.Date(i, origin="2015-01-01")
	this_data = merged_gadm_mun[which(!is.na(merged_gadm_mun@data$tiga_prob) & !is.na(merged_gadm_mun@data$prov) & merged_gadm_mun@data$prov != "00"),]
	# Cairo(paste0(path_to_untracked_output, '/daily_alert_prob/mosquito_alert_muni_daily_alert_prob_', formatC(i, width=3, flag="0"), "_", this_date, ".png"), type="png", units="px", width=1200, height=1200)
	png(paste0(path_to_untracked_output, '/daily_alert_prob/mosquito_alert_muni_daily_alert_prob_', formatC(i, width=3, flag="0"), "_", this_date, ".png"), type="cairo", units="px", width=1200, height=1200)
	make_alert_prob_muni_map_animation(this_prov=gadm_prov_utm, this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1), this_day=this_data$sea_days[1])
	dev.off()
}, mc.cores=ncores)

ccaas = c(cat="09",and="01",bal="04", mur="14", val="10" ,ara= "02" )
for(j in 1:length(ccaas)){
	dir.create(paste0(path_to_untracked_output, '/daily_alert_prob_', names(ccaas)[j]), showWarning=F)
}

make_alert_prob_muni_map_animation = function(this_prov, this_mun, this_spdf, these_values, these_breaks, this_day, clock_left=TRUE){
	par(mar=c(1,1,1,1))
	nclass = classIntervals(these_values,style="fixed", fixedBreaks=these_breaks)
	colcode = findColours(nclass, this_pal)
	plot(this_prov, border=NA, col="grey")
	plot(this_spdf, col=colcode, pch=16, cex=.1, border=NA, add=TRUE)
	plot(this_mun, lwd=.5, add=TRUE)
	plot(this_prov,  add=TRUE, border="white", lwd=5)
	plot(this_prov,  add=TRUE, lwd=3)
	if(clock_left){
		par(new=TRUE, fig = c(0, .4, .75, 1))
	} else{
		par(new=TRUE, fig = c(.5, .9, .75, 1))		
	}
	the_months = c("Jan.", "Feb.", "Mar.", "Apr.", "May", "Jun.", "Jul.", "Aug", "Sep.", "Oct.", "Nov", "Dec")[order(12:1)]
	pie(rep(1/12, 12), labels=the_months, col="white", border=F)
	if(clock_left){
		par(new=TRUE, fig = c(0, .4, .75, 1))
	} else{
		par(new=TRUE, fig = c(.5, .9, .75, 1))		
	}
	these_cols = rep("dark grey", 365)
	these_cols[(365-(this_day)):(365-(this_day+14))] = "white"
	pie(rep(1/365, 365), labels="", col=these_cols, border=F)
}

catprovs = c(bar = "08", tar = "43", gir="17", lle="25")

mclapply(min(muni_daily_probs$sea_days):max(muni_daily_probs$sea_days), function(i){
	print(i)
	flush.console()
	merged_gadm_mun = merge(gadm_mun_utm, muni_daily_probs[sea_days==i], by=c("PROVMUN", "prov"), all=T)
	this_date = as.Date(i, origin="2015-01-01")
	this_ccaa= "09"
	this_ccaa_gadm = 6
	this_data = merged_gadm_mun[which(!is.na(merged_gadm_mun@data$tiga_prob) & merged_gadm_mun@data$prov %in% catprovs),] 
		this_prov = gadm_prov_utm[which(gadm_prov_utm$ID_1==this_ccaa_gadm),]
		this_mun = gadm_mun_utm[which(gadm_mun_utm@data$prov %in% catprovs),]
	#	Cairo(paste0(path_to_untracked_output, '/daily_alert_prob_', names(ccaas)[j], '/tigatrapp_muni_daily_alert_prob_', formatC(i, width=3, flag="0"), "_", this_date, ".png"), type="png", units="px", width=1200, height=1200)
		png(paste0(path_to_untracked_output, '/daily_alert_prob_', names(ccaas)[j], '/tigatrapp_muni_daily_alert_prob_', formatC(i, width=3, flag="0"), "_", this_date, ".png"), type="cairo", units="px", width=1200, height=1200)
		make_alert_prob_muni_map_animation(this_prov, this_mun, this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1), this_day=this_data$sea_days[1], clock_left=FALSE)
		dev.off()
}, mc.cores=ncores)

# convert -delay 10 *_%03d_*.png animation.gif

# convert -delay 10 untracked/daily_alert_prob/*_%03d_* daily_alert_prob/animation.gif
# convert -delay 10 untracked/daily_alert_prob_cat/*.png daily_alert_prob_cat/animation.gif
# convert -delay 10 untracked/daily_alert_prob_val/*.png daily_alert_prob_val/animation.gif
# convert -delay 10 untracked/daily_alert_prob_mur/*.png daily_alert_prob_mur/animation.gif
# convert -delay 10 untracked/daily_alert_prob_bal/*.png daily_alert_prob_bal/animation.gif
# convert -delay 10 untracked/daily_alert_prob_and/*.png daily_alert_prob_and/animation.gif
# convert -delay 10 untracked/daily_alert_prob_ara/*.png daily_alert_prob_ara/animation.gif


