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

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
# ncores = 1

# first using gadm data
gadm_mun = readRDS(paste0(path_to_external_data, "ESP_adm4.rds"))
gadm_mun_utm = spTransform(gadm_mun, PROJ4_UTM_SPAIN)

spain_gpdf = readRDS("spain_grid.Rds")

spain_gpdf = spTransform(spain_gpdf, PROJ4_UTM_SPAIN)

tigacell_muni_fracs_list = mclapply(1:length(gadm_mun_utm), function(i){
	print(paste0(i, ' of ', length(gadm_mun_utm)))
	flush.console()
	result = NULL
	test = intersect(spain_gpdf, gadm_mun_utm[i,])
	if(! is.null(test) && length(test)>0){
	total_area = sum(unlist(lapply(test@polygons, function(x) x@Polygons[[1]]@area)))	
	result = data.table(NAME_2=test@data$NAME_2,ID_2=test@data$ID_2, NAME_4=test@data$NAME_4,ID_4=test@data$ID_4, TigacellID=as.character(test@data$TigacellID), frac=unlist(lapply(1:length(test), function(j) test@polygons[[j]]@Polygons[[1]]@area/total_area)))
	}
	return(result)
}, mc.cores=ncores)

tigacell_muni_fracs = rbindlist(tigacell_muni_fracs_list)

saveRDS(tigacell_muni_fracs, file="tigacell_muni_fracs.Rds")

# now INE polys and codes
library(maptools)
options(stringsAsFactors = FALSE)

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)

gadm_mun_utm$PROVMUN = as.character(gadm_mun_utm$PROVMUN)

gadm_mun_utm = gadm_mun_utm[which(gadm_mun_utm$PROVMUN!="00000" & ! gadm_mun_utm$COM %in% c("05", "18", "19")),] # excluding holes, canaries, ceuta, and melilla

PROVMUNs = unique(gadm_mun_utm$PROVMUN)

tigacell_muni_fracs_list = mclapply(1:length(PROVMUNs), function(i){
	print(paste0(i, ' of ', length(PROVMUNs)))
	flush.console()
	result = NULL
	this_PROVMUN = PROVMUNs[i]

	test = intersect(spain_gpdf, gadm_mun_utm[which(gadm_mun_utm$PROVMUN==this_PROVMUN),])
	if(! is.null(test) && length(test)>0){
	total_area = sum(unlist(lapply(test@polygons, function(x) x@Polygons[[1]]@area)))	
	result = data.table(MUNICIPIO0=test@data$MUNICIPIO0,PROVMUN=test@data$PROVMUN, COM=test@data$COM, TigacellID=as.character(test@data$TigacellID), frac=unlist(lapply(1:length(test), function(j) test@polygons[[j]]@Polygons[[1]]@area/total_area)))
	}
	return(result)
}, mc.cores=ncores)

tigacell_muni_fracs = rbindlist(tigacell_muni_fracs_list)

saveRDS(tigacell_muni_fracs, file="tigacell_muni_fracs_ine.Rds")

