### script de comandos em R
### abrir em R e executar cada linha com 
### <control> + <R>

### Instalar pacotes 'rgdal' e 'spdep' e tambem 
### os pacotes que estes usam: 'sp', etc...

### instalando pacote 'rgdal' e suas dependencias
install.packages("rgdal", dep=TRUE)

### instalando pacote 'spdep' e suas dependencias
install.packages("spdep", dep=TRUE)

### Objetivos:
### 1 - ler mapa .MIF
### 2 - ler mapa 'shapefile'
### 3 - ler dados de municipios
### 4 - vizualisar os dados no mapa
### 5 - testar autocorrelacao espacial (I Moran)

#### mudando diretorio de trabalho (onde estao os arquivos)
setwd("C:/Documents and Settings/Elias/Meus documentos/mapas")

#####################################################################
### 1 - ler mapa .MIF
### carregando pacote rgdal
### install.packages("rgdal") ### instala
require(rgdal) ### carrega

### importando mapa em formato .MIF
### arquivo obtido do site do geominas
map1 <- readOGR("MG_MUN96.MIF", "MG_MUN96")

### verifiando a classe do objeto map1
class(map1)

### um SpatialPolygonsDataFrame tem os elementos 'data' e 'polygons'
### inspecionando o 'data' (tabela de atributos associada)
dim(map1@data)
head(map1@data)

### plotando o mapa
plot(map1)


#####################################################################
### 2 - ler mapa 'shapefile'

### carregando o pacote spdep
### neste pacote ha funcoes para ler shapefiles 
### e para analises estatisticas espaciais
require(spdep)

### importando shapefile
### arquivo obtido do site do geominas
map2 <- readShapePoly("MG_MUN96")

### verifiando a classe do objeto map2
class(map2)

### inspecionando o 'data' (tabela de atributos associada)
dim(map2@data)
head(map2@data)

### plotando o mapa
plot(map2)

#####################################################################
### 3 - ler dados de municipios

### lendo os dados de percentual de residentes em domicilios 
### com rede geral de agua
### arquivo gerado com dados obtidos do site do datasus
dad <- read.csv2("redeGeralMG.csv")
dim(dad)
head(dad)

### sumario dos dados
summary(dad)
### variancia da variavel redeGeral
var(dad$redeGeral)
### desvio padrao da variavel redeGeral
sd(dad$redeGeral)
### histograma da variavel redeGeral
hist(dad$redeGeral)

#####################################################################
### 4 - vizualisar os dados no mapa
### colocar os dados na mesma ordem dos polygons

### extraindo o codigo dos municipios dos dados
cod.mun.dat <- as.numeric(substr(as.character(dad$Mun), 3, 7))
head(cod.mun.dat)

### vamos colocar os dados na ordem dos polygons
ord <- sapply(cod.mun.dat, function(x)
  which(is.element(map2@data$CODMUN,x)))
head(ord)

### ordenando os dados
dad.ord <- dad[ord, ]

### vizualisando os dados no mapa
cores <- c("brown", "red3", "salmon2", "salmon1", "yellow")
summary(dad$rede)
lims <- c(0, 20, 40, 60, 80, 100)
classe <- findInterval(dad.ord$rede, lims)
table(classe)

par(mar=c(0,0,2,0))
plot(map2, col=cores[classe])
title("Percentual residentes em domicilios com rede de įgua")
legend(c(-51, -48), c(-14, -18), 
  leglabs(lims, "Menor que", "Maior que"), 
  bty='n', fill=cores)

#####################################################################
### 5 - testar autocorrelacao espacial (I Moran)

### encontrando quem e' vizinho de quem (demora um pouco)
nb.mi.mg <- poly2nb(map2)
nb.mi.mg
nb.mi.mg[[1]]
nb.mi.mg[[2]]

### ponderacao
nbw <- nb2listw(nb.mi.mg)
names(nbw)
nbw$nei[[1]]
nbw$wei[[1]]

nbw$nei[[2]]
nbw$wei[[2]]

### teste de Moran
args(moran.mc)
imoran <- moran.mc(dad.ord$rede, nbw, nsim=999)
imoran

### vizualizando os resultados
par(mar=c(4,4,2,2))
hist(imoran$res, xlab='Indice', main='', col=gray(.5), border=gray(.7))
arrows(imoran$stat,-2,imoran$stat,10,lwd=2,col=2,leng=.1,code=1)
segments(imoran$stat, 3, 0.4, 120, lty=2)
text(.4, 150, paste("I Moran =", format(imoran$stat,dig=4)))
text(.4, 130, paste("valor-p =", format(imoran$p.val, dig=4)))