###################################### ## instalando os pacotes (se necessário) ###################################### install.packages(c("maptools", "sp", "spdep", "classInt", "RColorBrewer"), dep=T) ##################################### ## Carregando pacotes necessários ##################################### require(maptools) ## funções para importação/expostação e manipulação de mapas e dados geográficos gpclibPermit() ## função para habilitar licença de uso require(sp) ## (classes) para representação de dados espaciais no R require(spdep) ## funções análises de dados de áreas require(classInt) ## rotinas para faclitar a divisão de dados em classes por vários critérios require(RColorBrewer) ## usada aqui para criar palhetas de cores nas visualizações em mapas par.ori <- par(no.readonly=TRUE) ## ## Casos de AIDS no RJ ## arquivos necessários (copiar da página): ## rj_mun_region_pol.shp, rj_mun_region_pol.shx, rj_mun_region_pol.dbf, casosaids.txt ## importando mapa do RJ rj <- readShapePoly("rj_mun_region_pol.shp", IDvar="GEOCODIGO") rj$IDvar plot(rj) slotNames(rj) ## vendo os dados/atributos importados com os mapas dim(rj) names(rj) ## ou... dim(rj@data) names(rj@data) head(rj@data) rj@data[,1] str(rj@data) ## importanto arquivo com dados dos casos casos <- read.csv2("casosaids.txt") dim(casos) names(casos) str(casos) casos[,1] ## pareando os dados dos dois arquivos para poder "juntá-los" rj$GEOCODIGO <- as.character(rj$GEOCODIGO) casos$GEOCODIGO_STRING <- as.character(casos$GEOCODIGO_STRING) cbind(rj$GEOCODIGO, casos$GEOCODIGO_STRING) all.equal(rj$GEOCODIGO, casos$GEOCODIGO_STRING) ## note que estes dados já estão na mesma ordem!!! all.equal(rj$GEOCODIGO, casos$GEOCODIGO_STRING) ## mesmo assim, atítuolo de ilustração, vamos mostrar os comandos necessarios para combinar os dados ## assegurando a ordem correta ind <- match(rj$GEOCODIGO, casos$GEOCODIGO_STRING) ind casos <- casos[ind,] ## para "juntar" os dados é necessário colocar os nomes das colunas row.names(casos) <- rj$GEOCODIGO rj <- spCbind(rj, casos) ## ## Bairros de Curitiba ## ## lendo dados tipo shapefiles ctba <- readShapePoly("bairros.shp", IDvar="CODE") plot(ctba) head(ctba@data) ctba@data$NOME <- as.character(ctba@data$NOME) ## mudança de codificação de caracteres (só use se necessário, ## depende do sistema operacional e do enconding do sistema) ## Veja os nomes dos unicípios e se os acentos aparecem corretamente #Encoding(ctba@data$NOME) #Encoding(ctba@data$NOME) <- "latin1" #ctba@data$NOME <- enc2native(ctba@data$NOME) #head(ctba@data) tb <- read.csv("tabe.csv", enc="latin1") head(tb) ## ops nao importou corretamente assim! ## vamos importar de outra forma: tb <- read.table("tabe.csv", sep=";", head=T, enc="latin1") head(tb) str(tb) ## e agora ordenando os dados corretamente, para compatibilizar a ordem ## dos minicípios no shape e na tabela de atributos ind <- match(ctba@data$CODE, tb$CODE) ind tb <- tb[ind,] row.names(tb) <- ctba$CODE ctba <- spCbind(ctba, tb) names(ctba) head(ctba@data) ## Visualizando um mapa do atributo: Segurança #INT <- classIntervals(ctba$Segurança, n=5, style="quantile") INT <- classIntervals(ctba$Segurança, style="fixed", fixedBreaks=c(0, 20, 40, 60, 80, 100)) CORES.5 <- c(rev(brewer.pal(3, "Blues")), brewer.pal(3, "Reds"))[-1] ## tirando uma cor no final COL <- findColours(INT, CORES.5) plot(ctba, col=COL) title("Segurança") TB <- attr(COL, "table") legtext <- paste(names(TB)) #, " (", TB, ")", sep="") legend("bottomright", fill=attr(COL, "palette"), legend=legtext, bty="n") ## calculando vizinhanças (segundo um critério de vizinhança -- note que há outros!!!) ctba.nb1 <- poly2nb(ctba) args(poly2nb) ## OBS: outros tipos de vizinhanças podem ser calculados: ## - mudando argumentos de poly2nb ## - usando outras funções de vizinhança tal como: dnearneigh() and knearneigh() ## Moran global moran.test(ctba$Segurança, listw=nb2listw(ctba.nb1)) args(moran.test) moran.mc(ctba$Segurança, listw=nb2listw(ctba.nb1), nsim=99) args(moran.mc) ## Uma alternativa ao Moran: estatística de Geary geary.test(ctba$Segurança, listw=nb2listw(ctba.nb1)) ## Moran na presença de covariáveis (substituri 1 por covariaveis abaixo) ## faz Moran em resíduos lm.morantest(lm(ctba$Segurança ~1), listw=nb2listw(ctba.nb1)) lm.morantest.exact(lm(ctba$Segurança ~1), listw=nb2listw(ctba.nb1)) ## Moran local(um tipo de LISA - Local indicator of spatial association) ctba.mloc <- localmoran(ctba$Segurança, listw=nb2listw(ctba.nb1, style="C")) ctba.mloc ## Fazendo um mapa dos índices de Moran local head(ctba.mloc) range(ctba.mloc[,1]) INT <- classIntervals(ctba.mloc[,1], style="fixed", fixedBreaks=seq(-1, 2.5, by=0.5)) CORES.7 <- c(rev(brewer.pal(3, "Blues")), brewer.pal(4, "Reds")) COL <- findColours(INT, CORES.7) plot(ctba, col=COL) title("I de Moran Local para Segurança") TB <- attr(COL, "table") legtext <- paste(names(TB)) #, " (", TB, ")", sep="") legend("bottomright", fill=attr(COL, "palette"), legend=legtext, bty="n") #, cex=0.9, y.inter=0.8) ## uma visualização com outro tipo de fatiamento (por quantis) INT <- classIntervals(ctba$Segurança, n=7, style="quantile", dataPrecision=2) #INT <- classIntervals(ctba$Segurança, n=5, style="sd") CORES.7 <- c(rev(brewer.pal(3, "Blues")), brewer.pal(4, "Reds")) COL <- findColours(INT, CORES.7) plot(ctba, col=COL) title("I de Moran Local para Segurança") TB <- attr(COL, "table") legtext <- paste(names(TB)) #, " (", TB, ")", sep="") legend("bottomright", fill=attr(COL, "palette"), legend=legtext, bty="n") #, cex=0.9, y.inter=0.8) ## Poderia-se fazer tb um mapa dos índices de Moran local padronizados! ## Mapa das probabilidades (significâncias do I de Moral local) head(ctba.mloc) INT <- classIntervals(ctba.mloc[,5], style="fixed", fixedBreaks=c(0, 0.001, 0.01, 0.05, 0.10, 0.20, 1)) CORES.6 <- c(rev(brewer.pal(5, "Reds")), brewer.pal(5, "Blues")) COL <- findColours(INT, CORES.6) plot(ctba, col=COL) title("p-valores do I de Moran Local para Segurança") TB <- attr(COL, "table") legtext <- paste(names(TB)) #, " (", TB, ")", sep="") legend("bottomright", fill=attr(COL, "palette"), legend=legtext, bty="n", cex=0.9, y.inter=0.8) ## LISA Map (+/+) , (-,-), (+,-), (-,+) ## calculando a médias dos atributos locais (padronizados) com os dos vizinhos head(ctba.mloc) ## montando matrix W de vizinhança ctba.nb1.mat <- nb2mat(ctba.nb1) ## Segurança Padronizada Segurança_SD <- scale(ctba$Segurança) ## Cálculo da m;édia da Seguranca Padronizada nos visinhos Segurança_W <- ctba.nb1.mat %*% Segurança_SD ## Diagrama de espalhamento de Moran plot(Segurança_SD, Segurança_W) abline(v=0, h=0) ## Montando o mapa do espalhamento de Moran Q <- ifelse((Segurança_SD>=0 & Segurança_W >=0), 1, 0) Q[(Segurança_SD<0 & Segurança_W < 0)] <- 2 Q[(Segurança_SD>=0 & Segurança_W < 0)] <- 3 Q[(Segurança_SD<0 & Segurança_W >= 0)] <- 4 table(Q) CORES.4 <- c("blue", "green" , "red", "yellow") plot(ctba, col=CORES.4[Q]) title("Mapa LISA") legend("bottomright", c("Q(+/+)", "Q(-/-)", "Q(+/-)", "Q(-/+)"), fill=CORES.4) ## Fazer análises para outras variáveis destes dados!!!table ## OBS: ## Para dados não contínuos (e normais) modificações e alternativas são usadas!!!