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

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

setwd(path_to_data)

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'))
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")
muni_summer_probs = readRDS("muni_summer_probs.Rds")
muni_sept_probs = readRDS("muni_sept_probs.Rds")



jetpal = c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")

YlOrRd.colors = colorRampPalette(c(brewer.pal(8, "YlOrRd")[1], brewer.pal(8, "YlOrRd")[8]))
jet.colors <- colorRampPalette(jetpal)

this_pal = jet.colors(10)[c(1, 6:10)]

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

make_alert_prob_muni_map = function(this_spdf, these_values, these_breaks){
	par(mar=c(1,1,1,1))
nclass = classIntervals(these_values,style="fixed", fixedBreaks=these_breaks)
colcode = findColours(nclass, this_pal)
plot(gadm_prov_utm, col="grey")
plot(this_spdf, col=colcode, pch=16, cex=.1, border=NA, add=T)
legend('bottomright', legend=names(attr(colcode, 'table')), fill=this_pal, bty="n", inset=c(.01, .15))
}


########################
# plotting summer tiga #
########################

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

merged_gadm_mun = merge(gadm_mun_utm, muni_summer_probs, by="PROVMUN", all=T)
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"),] 
describe(this_data@data$tiga_prob)

pdf(paste0(path_to_figures, 'mosquito_alert_muni_summer_risk_daily_scale.pdf'))
	make_alert_prob_muni_map(this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1))
dev.off()

png(paste0(path_to_figures, 'mosquito_alert_muni_summer_risk_daily_scale.png'))
	make_alert_prob_muni_map(this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1))
dev.off()

#########
# FIG 3 #
#########
pdf(paste0(path_to_figures, 'fig3.pdf'))
	make_alert_prob_muni_map(this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1))
dev.off()
png(paste0(path_to_figures, 'fig3.png'))
	make_alert_prob_muni_map(this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1))
dev.off()

######################
# plotting sept tiga #
######################
muni_sept_probs$prov = substring(muni_sept_probs$PROVMUN, 1,2)

merged_gadm_mun = merge(gadm_mun_utm, muni_sept_probs, by="PROVMUN", all=T)
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"),] 

describe(this_data@data$tiga_prob)

pdf(paste0(path_to_figures, 'mosquito_alert_muni_sept_risk_daily_scale.pdf'))
	make_alert_prob_muni_map(this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1))
dev.off()

png(paste0(path_to_figures, 'mosquito_alert_muni_sept_risk_daily_scale.png'))
	make_alert_prob_muni_map(this_spdf=this_data, these_values=this_data@data$tiga_prob, these_breaks = c(seq(0, .4, .1), 1))
dev.off()





