2.1 Pacotes necessários

Agora deve-se instalar os pacotes úteis para fazer os mapas. O R vem com alguns pacotes automaticamente instalados, mas outros pacotes púbicos devem ser chamados para realizar operações específicas. Um deles é o pacote rgdal.

if (!require(rgdal)) {
    install.packages("rgdal")
    library(rgdal)
}

O pacote rgdal já fornecerá a primeira função prática para fazer o mapa. Usa-se aqui a função readOGR e indica-se o nome do arquivo shape a ser buscado, neste caso, na pasta MS e o arquivo de nome 50MUE250GC_SIR. Neste caso, a base de municípios (indicado por MU). A saída (plot) mostra o mapa e as características básicas (79 municípios e 2 campos - nome e código do município).

MS.mu <- readOGR("MS", "50MUE250GC_SIR")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\amrof\Documents\disciplinas\economia regional\material de aula\laboratorio R\aula_mapas\mapas_no_R\MS", layer: "50MUE250GC_SIR"
with 79 features
It has 91 fields

plot(MS.mu)

Para a mesorregião homogênea:

MS.me <- readOGR("MS", "50MEE250GC_SIR")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\amrof\Documents\disciplinas\economia regional\material de aula\laboratorio R\aula_mapas\mapas_no_R\MS", layer: "50MEE250GC_SIR"
with 4 features
It has 2 fields

plot(MS.me)

Para a microrregião homogênea:

MS.mi <- readOGR("MS", "50MIE250GC_SIR")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\amrof\Documents\disciplinas\economia regional\material de aula\laboratorio R\aula_mapas\mapas_no_R\MS", layer: "50MIE250GC_SIR"
with 11 features
It has 2 fields

plot(MS.mi)

Para a Unidade da Federação de Mato Grosso do Sul (MS):

MS.uf <- readOGR("MS", "50UFE250GC_SIR")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\amrof\Documents\disciplinas\economia regional\material de aula\laboratorio R\aula_mapas\mapas_no_R\MS", layer: "50UFE250GC_SIR"
with 1 features
It has 3 fields

plot(MS.uf)

Agora propõe-se colorir o mapa das microrregiões com cores aleatórias, com janelas perguntando o arquivo origem das informações. O usuário será chamado a escolher numa janela o arquivo desejado. Para tanto, será necessário o pacote maptools.

if (!require(maptools)) {
    install.packages("maptools")
    library(maptools)
}

# abrirá a tela para escolher o shapefile
mapa <- readOGR(file.choose())

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\amrof\Documents\disciplinas\economia regional\material de aula\laboratorio R\aula_mapas\mapas_no_R\50MUE250GC_SIR.shp", layer: "50MUE250GC_SIR"
with 79 features
It has 23 fields

# gerar dados aleatórios e ajuntá-los ao arquivo shape
data <- rnorm(nrow(mapa@data))
# definir os intervalos e as cores com comandos
intervalos <- c(-Inf, -2, -1, 0, 1, 2, Inf)
cortes <- cut(data, intervalos, include.lowest = TRUE)
niveis <- levels(cortes)
cores <- heat.colors(length(levels(cortes)))
levels(cortes) <- cores
# com a função plot podemos desenhar o conjunto de informação introduzido no
# objetos
plot(mapa, border = gray(0.9), lwd = 0.1, axes = FALSE, las = 1, col = as.character(cortes))

# Se quiser, poderá agregar ainda legenda, contornos etc.

3.1 Vizinhança Queen (Rainha) e Rook (Torre)

# fazendo vizinhos por fronteira comum tipo Queen - Rainha library(spdep)
sids_nbq <- poly2nb(sids)  # default = queen = rainha
sids_nbr <- poly2nb(sids, queen = FALSE)  # vizinhos do tipo rook = torre
coords <- coordinates(sids)  # coordenadas

# split.screen(c(1,2))
plot(sids, main = "MS 79: vizinhanca Queen")
plot(sids_nbq, coords, add = T, col = "red", lty = 2)  # centroide ligado ao centroide vizinho

# screen(2)
plot(sids, main = "MS 79: vizinhanca Rook")
plot(sids_nbr, coords, add = T, col = "red", lty = 2)

# close.screen(all = TRUE)

3.2 Vizinhança pela distância

# fazendo vizinhos pela distancia: os k mais proximos
IDs <- row.names(as(sids, "data.frame"))
# Neighbours list from knn object: knn2nb K nearest neighbours for spatial
# weights: knearneigh a partir das coords ele define o k=1 mais próximo como
# vizinho ou o k=2 (os dois mais próximos como vizinhos) ou o k=4 (os quatro
# mais próximos como vizinhos)
sids_kn1 <- knn2nb(knearneigh(coords, k = 1), row.names = IDs)
sids_kn2 <- knn2nb(knearneigh(coords, k = 2), row.names = IDs)
sids_kn4 <- knn2nb(knearneigh(coords, k = 4), row.names = IDs)

# split.screen(c(1,2))
plot(sids, main = "MS 79: k=2-nearest vizinhos")
plot(sids_kn2, coords, add = T, col = "red", lty = 2)  # centroide ligado aos dois centroides vizinhos

# screen(2 )
plot(sids, main = "MS 79: k=4-nearest vizinhos")
plot(sids_kn4, coords, add = T, col = "green", lty = 2)
legend("topright", lty = 1, pch = 1, col = 1:3, c("fronteiras", "k=2", "k=4"), 
    cex = 0.7)

# close.screen(all = TRUE)

# Pode atribuir vizinhos por uma valor determinado de distância

dist <- unlist(nbdists(sids_kn1, coords))
summary(dist)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1304  0.3210  0.3876  0.4235  0.5106  0.9192 

max_k1 <- max(dist)

sids_kd1 <- dnearneigh(coords, d1 = 0, d2 = 0.75 * max_k1, row.names = IDs)
sids_kd2 <- dnearneigh(coords, d1 = 0, d2 = 1 * max_k1, row.names = IDs)
sids_kd3 <- dnearneigh(coords, d1 = 0, d2 = 1.5 * max_k1, row.names = IDs)
# Ou pela distância bruta
sids_ran1 <- dnearneigh(coords, d1 = 0, d2 = 0.9196, row.names = IDs)

# split.screen(c(1,3))
plot(sids, main = "MS 79: d2=0.75*max_k1")
plot(sids_kd1, coords, add = T, col = "red", lty = 2)  # centroide ligado aos dois centroides vizinhos
legend("topright", lty = 1, pch = 1, col = 1:2, c("fronteira", "d2=0.75*max_k1"), 
    cex = 0.7)

# screen(2)
plot(sids, main = "MS 79: d2=1*max_k1")
plot(sids_kd2, coords, add = T, col = "green", lty = 2)
legend("topright", lty = 1, pch = 1, col = c(1, 4), c("fronteira", "d2=1*max_k1"), 
    cex = 0.7)

# screen(3)
plot(sids, main = "MS 79: d2=1.5*max_k1")
plot(sids_kd3, coords, add = T, col = "blue", lty = 2)
legend("topright", lty = 1, pch = 1, col = c(1, 4), c("fronteira", "d2=1.5*max_k1"), 
    cex = 0.7)

# close.screen(all = TRUE)

3.3 Criando matrizes de ponderação (W)

# criando matrizes de ponderação (W)
sids_nbq_w <- nb2listw(sids_nbq)  # usei criterio rainha sids_nbq
sids_nbq_w  # pesos com padronização para linha

Characteristics of weights list object:
Neighbour list object:
Number of regions: 79 
Number of nonzero links: 400 
Percentage nonzero weights: 6.409229 
Average number of links: 5.063291 

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 79 6241 79 34.13411 328.4183

# pesos binários
sids_nbq_wb <- nb2listw(sids_nbq, style = "B")
sids_nbq_wb

Characteristics of weights list object:
Neighbour list object:
Number of regions: 79 
Number of nonzero links: 400 
Percentage nonzero weights: 6.409229 
Average number of links: 5.063291 

Weights style: B 
Weights constants summary:
   n   nn  S0  S1   S2
B 79 6241 400 800 9152

# weights constants are: n = number of observations, n number of zones,
# S0=global sum of weights, S1= S1 sum of weights, S2=S2 sum of weights

# para contornar problemas de regioes sem vizinhos
sids_nbq_w <- nb2listw(sids_nbq, zero.policy = T)
sids_nbq_w

Characteristics of weights list object:
Neighbour list object:
Number of regions: 79 
Number of nonzero links: 400 
Percentage nonzero weights: 6.409229 
Average number of links: 5.063291 

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 79 6241 79 34.13411 328.4183

4.1 Moran’s I ou I de Moran

# I de Moran - Moran's I W queen padronizado
moran.test(sids$R2015, listw = sids_nbq_w, alternative = "two.sided")

    Moran I test under randomisation

data:  sids$R2015  
weights: sids_nbq_w    

Moran I statistic standard deviate = 2.841, p-value = 0.004497
alternative hypothesis: two.sided
sample estimates:
Moran I statistic       Expectation          Variance 
      0.125560803      -0.012820513       0.002372496 

# W queen binario
moran.test(sids$R2015, listw = sids_nbq_wb)

    Moran I test under randomisation

data:  sids$R2015  
weights: sids_nbq_wb    

Moran I statistic standard deviate = 2.3103, p-value = 0.01043
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.094594746      -0.012820513       0.002161642 

# preparando o moran scatterplot estabelecendo o ponto de partida para o
# Monte Carlo
set.seed(123456)
moran.montecarlo <- moran.mc(sids$R2015, sids_nbq_w, 99)
moran.montecarlo$res  # pegando o resultado do moran.mc

  [1] -0.023187760 -0.087711892 -0.066936966  0.027595111 -0.066123207
  [6] -0.038615763  0.034472549 -0.040418876 -0.023306819 -0.030569311
 [11]  0.108895743  0.046451675 -0.055173679 -0.018647436  0.032714419
 [16]  0.002348992 -0.024514045 -0.023576562 -0.071374293 -0.037244252
 [21] -0.051332712  0.177769888  0.026641260 -0.001527644 -0.024850271
 [26] -0.056484251 -0.021525655 -0.020150044  0.099141912 -0.079754592
 [31]  0.025252484 -0.022534351 -0.017610231 -0.083100317 -0.086846333
 [36] -0.048108128 -0.043746993 -0.077652470 -0.045021486 -0.023993092
 [41] -0.089161254  0.009796032 -0.049253377 -0.017139153 -0.026388982
 [46]  0.008134913  0.003652464  0.103046780 -0.012272921  0.036892639
 [51] -0.002564535 -0.001484142 -0.033879675 -0.117919378 -0.057039339
 [56] -0.012269654 -0.004463029  0.053202754 -0.032920478  0.048630094
 [61]  0.006734393 -0.058086615 -0.035628564 -0.010837500  0.009986281
 [66] -0.049359985 -0.032831137  0.056050764 -0.049185023  0.044701742
 [71]  0.042390053 -0.051084447  0.030366496 -0.028047725  0.076447712
 [76] -0.019964308 -0.071385349  0.029141166  0.051914207 -0.015843474
 [81]  0.019687494 -0.014916004 -0.015857618 -0.032654524 -0.042402831
 [86] -0.002181341 -0.100689365 -0.032815736  0.018119389 -0.000229275
 [91] -0.051693677 -0.054190695 -0.041968578  0.001342042  0.017806155
 [96]  0.020815835  0.039102223 -0.046384865 -0.046352871  0.125560803

# o ultimo valor eh o valor da estatistica I-Moran Global obtido

moran.plot(sids$R2015, listw = sids_nbq_w, labels = as.character(sids$NM_MUNICIP), 
    pch = 13)

4.2 LISA ou Local Moran

oid <- order(sids$CD_GEOCMU)
set.seed(123456)
resI <- localmoran(sids$R2015, listw = sids_nbq_w, alternative = "two.sided")
# outra opcao (identica a linha anterior) mais direta
resI_2 <- localmoran(sids$R2015, nb2listw(sids_nbq, zero.policy = T), alternative = "two.sided")
printCoefmat(data.frame(resI, row.names = sids$NM_MUNICIP), check.names = FALSE)

                                  Ii        E.Ii      Var.Ii        Z.Ii
ÁGUA CLARA                0.31945794 -0.01282051  0.06694425  1.28423806
ALCINÓPOLIS              -0.00294433 -0.01282051  0.11642994  0.02894389
AMAMBAI                   0.01378164 -0.01282051  0.06694425  0.10281586
ANASTÁCIO                 0.10711882 -0.01282051  0.07794107  0.42961403
ANAURILÂNDIA              0.03147710 -0.01282051  0.15491882  0.11254552
ANGÉLICA                  0.18988554 -0.01282051  0.09333662  0.66349964
ANTÔNIO JOÃO              0.10489495 -0.01282051  0.23189656  0.24444794
AQUIDAUANA                0.18628321 -0.01282051  0.05869663  0.82181256
ARAL MOREIRA             -0.03185011 -0.01282051  0.11642994 -0.05576955
BANDEIRANTES             -0.00065335 -0.01282051  0.07794107  0.04358188
BATAGUASSU                0.00285524 -0.01282051  0.11642994  0.04594055
BATAYPORÃ                 0.03046449 -0.01282051  0.15491882  0.10997283
BELA VISTA                0.17117162 -0.01282051  0.11642994  0.53922114
BODOQUENA                 0.07388083 -0.01282051  0.11642994  0.25409346
BONITO                    0.06936015 -0.01282051  0.06694425  0.31762377
BRASILÂNDIA               0.39222621 -0.01282051  0.11642994  1.18706031
CAARAPÓ                  -0.00486582 -0.01282051  0.07794107  0.02849312
CAMAPUÃ                  -0.04184901 -0.01282051  0.06694425 -0.11219354
CAMPO GRANDE             -0.00234627 -0.01282051  0.07794107  0.03751799
CARACOL                   0.21542967 -0.01282051  0.15491882  0.57990793
                         Pr.z....0.
ÁGUA CLARA                   0.1991
ALCINÓPOLIS                  0.9769
AMAMBAI                      0.9181
ANASTÁCIO                    0.6675
ANAURILÂNDIA                 0.9104
ANGÉLICA                     0.5070
ANTÔNIO JOÃO                 0.8069
AQUIDAUANA                   0.4112
ARAL MOREIRA                 0.9555
BANDEIRANTES                 0.9652
BATAGUASSU                   0.9634
BATAYPORÃ                    0.9124
BELA VISTA                   0.5897
BODOQUENA                    0.7994
BONITO                       0.7508
BRASILÂNDIA                  0.2352
CAARAPÓ                      0.9773
CAMAPUÃ                      0.9107
CAMPO GRANDE                 0.9701
CARACOL                      0.5620
 [ reached getOption("max.print") -- omitted 59 rows ]

# a funcao retorna valores de Ii, E.Ii, Var.Ii, z.Ii e Pr.z.Ii Ii
# estatistica local moran E.Ii Valor esperado da estatistica local moran
# (media) Var.Ii variancia da estatistica local moran Z.Ii desvio padrao da
# estatistica local moran Pr() p-value da estatistica local moran posso
# enviar para um arquivo do excel

# library(xlsx)
dados.LISA.CSV <- write.csv(cbind(sids, resI), "dados.LISA.csv", col.names = TRUE, 
    row.names = FALSE)


hist(resI[, 5])

mean(resI[, 1])

[1] 0.1255608

4.3 Mapa do Local Moran (LISA)

# codigo de Roger Bivand em
# http://r-sig-geo.2731867.n2.nabble.com/Local-Moran-p-values-in-Geoda-Python-and-R-td7587887.html
nat.data <- data.frame(sids)
attach(nat.data)
# coords<-coordinates(sids) # coordenadas, já feito anteriormente
# IDs<-row.names(as(sids, 'data.frame')) # IDs, já feito anteriormente
# sids_nbq_w<-nb2listw(sids_nbq, zero.policy=T) # matriz W queen já feito
# anteriormente
list_w <- sids_nbq_w
# LISA values
lisa <- localmoran(R2015, list_w, zero.policy = T, alternative = "two.sided")
# centers the variable of interest around its mean
cDV <- R2015 - mean(R2015)
# centers the local Moran's around the mean
mI <- lisa[, 1]
C_mI <- mI - mean(mI)  # MAS N?O QUEREMOS CENTRAR! Apenas o sinal importa
quadrant <- vector(mode = "numeric", length = nrow(lisa))
quadrant[cDV > 0 & mI > 0] <- 1
quadrant[cDV < 0 & mI > 0] <- 2
quadrant[cDV > 0 & mI < 0] <- 3
quadrant[cDV < 0 & mI < 0] <- 4
# set a statistical significance level for the local Moran's
signif <- 0.05
# places non-significant Moran's in the category '5'
quadrant[lisa[, 5] > signif] <- 5
# map
png(file = "Lisa_map_R.png", width = 680, res = 100)
colors <- c("red", "blue", "lightpink", "skyblue2", rgb(0.95, 0.95, 0.95))
par(mar = c(0, 0, 1, 0))  # sets margin parameters for plot space 
plot(sids, border = "grey", col = colors[quadrant], main = "LISA Cluster Map, 2015 Produtividade Relativa (R)")
legend("bottomright", legend = c("alto-alto", "baixo-baixo", "alto-baixo", "baixo-alto"), 
    fill = colors, bty = "n", cex = 0.7, y.intersp = 1, x.intersp = 1)
dev.off()

png 
  2 

######################################################################################## na tela
colors <- c("red", "blue", "lightpink", "skyblue2", rgb(0.95, 0.95, 0.95))
par(mar = c(0, 0, 1, 0))  # sets margin parameters for plot space 
plot(sids, border = "grey", col = colors[quadrant], main = "LISA Cluster Map, 2015 Produtividade Relativa (R)")
legend("bottomright", legend = c("alto-alto", "baixo-baixo", "alto-baixo", "baixo-alto"), 
    fill = colors, bty = "n", cex = 1, y.intersp = 1, x.intersp = 1)