quadDep.R

library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(ggplot2)

# Lê um shapefile com UFs do Brasil - http://www.usp.br/nereus/?dados=brasil
uf <- readOGR("../_temp/ufe.shp",encoding="UTF-8",stringsAsFactors = FALSE)
uf_centroides <- gCentroid(uf,byid = TRUE)
ufs = gSimplify(uf,tol=0.01)

# Lê uma lista com número de deputados por Estado
dep <- read.csv("../_temp/dep.csv",sep=";",encoding = "UTF-8")

# Calcula o número de células necessário para ter ao menos 513 deputados
# sobre a área do Brasil. A seguir calcula o tamanho da célula (lado - raiz 
# quadrada)
areaq <- sum(area(uf))/513 #Soma as áreas das UFs e divide por 513 deputados
ladoq <- sqrt(areaq) #Lado mínimo

# Calcula os parâmetros da grade
# Posição inicial (offset): usa o bounding box do shape das UFs, acrescenta metade do tamanho de célula
cc <- ceiling(bbox(uf)[,1]) + (ladoq/2)

# Dimensão da célula
cd <- ceiling(diff(t(bbox(uf)))/ladoq) 

# Cria a grade como topologia
grade <- GridTopology(cellcentre.offset=cc, cells.dim=cd, cellsize = c(ladoq,ladoq))

# Cria a grade como SpatialGrid, usando mesmo CRS das UFs
grade_sp <- SpatialGridDataFrame(grade,data.frame(ID=c(1:(cd[1]*cd[2]))),proj4string=CRS(proj4string(uf)))

# Cria uma grade de polígonos.
grade_poly <- as(grade_sp,"SpatialPolygonsDataFrame")

# Calcula os polígonos em intersecção.
qa <- intersect(ufs,grade_poly)

#Mantém apenas 520 polígonos
qa$area = area(qa)
qMa = qa[order(qa$area,decreasing = TRUE),] # Ordena por Maior área.
qa <- qMa[1:520,] 
grade_poly <- grade_poly[unlist(qa$ID),]

# Inicia um dataframe que vai conter a distância de cada centróide 
# da UF à coordenada central de cada célula da grade
dist <- data.frame(i = "",j = "", dist = 0.0,stringsAsFactors = FALSE)

# Passa para um dataframe as coordenadas dos centróides, facilitando iteração no loop
gradeC <- gCentroid(qa,byid = TRUE)

# Iteração para cada UF (centróides), calcula a distância para todos os 
# pontos da grade.
# (Substituir o uso do for por apply?)
for (i in 1:nrow(uf_centroides@coords))
{
	for(j in 1:nrow(gradeC@coords))
	{
		# d - distância calculada
		d <- sqrt( (uf_centroides[i,]$x - gradeC[j,]$x)**2 + (uf_centroides[i,]$y - gradeC[j,]$y)**2 )
		# Acrescenta na matriz de distâncias: o geocódigo da UF, o número 
		# da célula de grade e a distância.
		dist <- rbind(dist,c(as.character(uf[i,]$CD_GEOCODU),qa[j,]$ID,d))
	}
}

# Converte as distâncias para número.
dist$dist <- as.numeric(dist$dist)

# Cria um vetor que irá conter as células já usadas (já atribuídas a um estado)
cel_usadas <- as.vector(0)

# Cria uma matriz que conterá a célula e o estado a que estiver associada
matriz <- data.frame("","",stringsAsFactors = FALSE)
# Obs.: parece uma tarefa duplicada, mas tentei separar a matriz final 
# de um controle de células usadas.

# Faz a iteração para cada linha do quadro de deputados por UF.
# Para cada UF, irá calcular as N células mais próximas ao centróide,
# que não tiverem sido previamente usadas por outra UF.
for(w in 1:nrow(dep))
{
	# Atribui a uma variável totaldep o total de deputados na UF
	# que está sendo iterada.
	totaldep <- dep[w,2]
	
	# Separa somente os registros da tabela de distância da UF em questão.
	l <- dist[which(dist$i == as.character(dep[w,1])),]
	
	# Matriz de proximidade ordenada (distâncias crescentes).
	proximidade <- l[order(l$dist,decreasing = FALSE),]
	
	# Nova iteração, percorre a lista obtendo as células que não foram
	# previamente usadas por outra UF.
	# x é o contador de deputados atribuídos
	# n é o contador de células da grade
	x = 1
	n = 1
	while((x <= totaldep) & (n < 520))
	{
		if( proximidade[n,2] %in% cel_usadas )
		{
			# ... passa para a próxima célula.
			n <- n + 1
		}
		else { # Se não foi usado...
			# Acrescenta na matriz final a atribuição do estado ao identificador
			# da célula na grade.
			matriz <- rbind(matriz, c(as.character(dep[w,1]),proximidade[n,2]))
			
			# Também acrescenta no vetor de células usadas.
			cel_usadas <- c(unlist(cel_usadas),proximidade[n,2])
			
			# Próximo deputado.
			x <- x + 1
		}
	}
}


grade_poly$Q = 0
for(i in 1:nrow(matriz))
{
	uf = matriz[i,1]
	quad = matriz[i,2]
	grade_poly@data[which(grade_poly$ID == as.numeric(quad)),]$Q = uf
}

plot(grade_poly,col=ceiling(as.numeric(grade_poly$Q)/5))