ANÁLISE DESCRITIVA ESPACIAL COM O
PACOTE geobr (PARTE II)
1ª ESCOLA DE VERÃO
PROF. GUILHERME AUGUSTO VELOSO
UFF - 20, 22 E 24 DE MARÇO DE 2023
Os dados de área são associados a levantamentos populacionais, como censos e estatísticas de saude e que originalmente se referem a indivíduos localizados em pontos específicos do espaço. Por razões, por exemplo, de confidencialidade, estes dados são agregados em unidades de análise, usualmente delimitadas por polígonos fechados (setores censitários, zonas de endereçamento postal, municípios).
A forma usual de apresentação de dados agregados por áreas é o uso de
mapas coloridos com o padrão espacial do
fenômeno. O objetivo dessa segunda parte do curso é a
construção e manipulação de mapas coropléticos no
R. Além disso, essa parte tem como objetivo a
introdução de conceitos importantes da análise
descritiva e exploratória de dados espaciais.
Precisaremos dos pacotes:
geobr: Permite o acesso aos shapefiles de diferentes unidades geográficas brasileiras.ggplot2: construção e customização de mapas coropléticos.RColorBrewer: Paleta de Cores.openxlsx: leitura de arquivos no formato.xlsx.ggiraph: permite a construção de mapas dinâmicos.sp: classes e métodos para dados espaciais.spdep: construção das matrizes de vizinhançadplyr: manipulação de bases de dados
Dessa forma, segue o código para a instalação dos pacotes:
if(!require(geobr)){ install.packages("geobr"); require(geobr)}
if(!require(ggplot2)){ install.packages("ggplot2"); require(ggplot2)}
if(!require(RColorBrewer)){ install.packages("RColorBrewer"); require(RColorBrewer)}
if(!require(openxlsx)){ install.packages("openxlsx"); require(openxlsx)}
if(!require(ggiraph)){ install.packages("ggiraph"); require(ggiraph)}
if(!require(sp)){ install.packages("sp"); require(sp)}
if(!require(spdep)){ install.packages("spdep"); require(spdep)}1 Exemplo Prático
O exemplo que utilizaremos para a visualização e exploração de dados espaciais vem do estudo GBD (Global Burden of Disease). Para mais detalhes sobre o estudo GBD, basta consultar os sites GBD e GBD Compare
Nesta parte do curso, utilizaremos a base de dados contendo as
taxas de mortalidade padronizadas por idade associadas a
acidentes envolvendo motociclistas nos municípios de Minas Gerais e do
Rio de Janeiro durante o triênio 2009-2011. Os dados podem ser
baixados nos links DADOMG
e DADORJ.
Basta ir em arquivo e fazer download.
2 Visualização dos Dados de Área
O primeiro passo na utilização dos dados de área é a escolha da unidade geográfica a ser explorada. Essa escolha está diretamente associada à disponibilidade dos dados que serão fornecidos e/ou coletados e, portanto, dependerá do fenômeno a ser estudado.
As técnicas de análise descritiva exploratória aplicadas a dados espaciais de área são essenciais ao desenvolvimento das etapas da modelagem estatística espacial, em geral sensível ao tipo de distribuição, à presença de valores extremos.
As técnicas empregadas são, em geral, adaptações das ferramentas usuais. Assim, se na investigação de valores extremos se utiliza ferramentas gráficas como histogramas ou boxplots, na análise espacial é importante também investigar outliers não só no conjunto dos dados mas também em relação aos vizinhos.
O uso de diferentes pontos de corte da variável induz a visualização de diferentes aspectos. Basicamente, podem ser utilizados três formas distintas de pontos de corte:
intervalos iguais: Os valores máximo e mínimo são divididos pelo número de classes. Se a variável tem um distribuição muito concentrada de um lado, a maior parte das áreas será alocada a poucas cores.
percentis: O uso de percentis para definação de classes obriga a alocação dos polígonos em quantidades iguais pelas cores. Isto pode mascarar diferenças significativas em valores extremos e pode dificultar a identificação de áreas críticas.
desvios padrões: A distribuição da variável é apresentada em gradações de cores diferentes para valores acima e abaixo da média. Faz a suposição da normalidade da distribuição da variável. Esta hipótese é pouco realista no caso de variáveis em países de grande desigualdade social como o Brasil.
Em resumo, é parte importante da análise exploratória experimentar diferentes pontos de corte da variável na visualização dos mapas.
As diferentes técnicas de visualização estão ilustradas no exemplo a seguir, em que mostramos a distribuição espacial do indicador que mede a proporção de recém-nascidos que nasce em boas condições de saúde (Índice de APGAR) para os bairros do Rio de Janeiro, no ano de 1994. Foram geradas duas visualizações, ambas com 5 pontos de corte e 5 cores. Na primeira figura, utilizou-se quintis e, na próxima, cinco classes de igual tamanho.
Como a distribuição da variável não é simétrica, quando se divide em classes de amplitudes iguais as de valores mais baixos (ou piores), assinaladas em vermelho ficam reduzidas a poucas áreas, enquanto que na divisão em quintis, por definição, um quinto das áreas ficará em cada classe.
Outra questão interessante é a comparação de mapas. Supondo a distribuição espacial de um indicador em diferentes anos: como visualizar a evolução temporal? Certamente os pontos de corte da variável nos diferentes períodos devem ser os mesmos. Observe na figura a evolução temporal da mortalidade por homicídios para os triênios 79-81 e 90-92, no Estado do Rio de Janeiro.
A apresentação dos quintis da distribuição conjunta dos indicadores permite visualizar bem o espalhamento desta “doença”.
Para a customização dos mapas nesse curso,
utilizaremos a tabela de cores oriunda do pacote
RColorBrewer, disponível no R. O pacote
RColorBrewer é uma ferramenta para gerenciar
cores no R, oferecendo várias paletas de
cores, como pode ser visto abaixo:
Existem 3 tipos de paletas:
As paletas sequenciais são adequadas para dados ordenados que progridem de baixo para cima. As etapas de claridade dominam a aparência desses esquemas, com cores claras para valores de dados baixos e cores escuras para valores de dados altos.
As paletas qualitativas não implicam em diferenças de magnitude entre as classes de legenda, e os matizes são usados para criar as diferenças visuais primárias entre as classes.
As paletas divergentes colocam igual ênfase nos valores críticos intermediários e extremos em ambas as extremidades do intervalo de dados. A classe crítica ou quebra no meio da legenda é enfatizada com cores claras e os extremos baixo e alto são enfatizados com cores escuras que possuem matizes contrastantes.
Para ajudar a categorizar os pontos de corte, precisaremos da seguinte função age.cat:
## Função utilizada para categorizar os pontos de corte
age.cat <- function(x, breaks,
sep = "-", above.char = "+",min,max) {
labs <- c(paste(c(min,breaks)[-(length(breaks)+1)],
breaks-0.01,
sep = sep),
paste(breaks[length(breaks)], max, sep = "-"))
cut(x, breaks = c(min,breaks,Inf),
right = FALSE, labels = labs)
}A construção do mapa coroplético considerando percentis como pontos de corte é feito a partir do seguinte código
## Leitura dos dados
Dados_MG <- read.xlsx("DadosMG.xlsx") ## Leitura do Shape dos municípios Mineiros
MG<-read_municipality(code_muni=31,year=2018,simplified=T)## Pontos de corte
breaks <- round(quantile(Dados_MG$taxa,probs = c(0.05,0.20,0.40,0.60,0.80,0.95)),2)
## Escolha da paleta de cores
mycolors <- rev(brewer.pal(7, "RdBu"))
## Categorizando os pontos de corte
Dados_MG$Categoria <- age.cat(round(Dados_MG$taxa,2),breaks,min=min(round(Dados_MG$taxa,2)),max=max(round(Dados_MG$taxa,2)))
## Merge do arquivo contendo as taxas com o shapefile
Dados_MG <- merge(MG,Dados_MG,by="code_muni",all.x=T)## Construção do mapa
Mapa_MG <- ggplot() + geom_sf(data=Dados_MG,aes(fill=Categoria),colour="black") +
theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
Mapa_MG
Os comandos a seguir permitem a construção do mapa dinâmico:
## O que será mostrado ao encostar em cada município
Dados_MG$tooltip <- c(paste0(Dados_MG$name_muni, "\n ", round(Dados_MG$taxa,2)))
## Construção do mapa
gg <- ggplot(Dados_MG) + geom_sf(aes(fill=Categoria)) +
geom_sf_interactive(aes(fill=Categoria, tooltip = tooltip, data_id = code_muni),color="black",size=0.1)+ theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de Mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
x <- girafe(ggobj = gg)
print(x)2.0.1 Atividade
Reproduza o exemplo anterior para os municípios do estado do Rio de Janeiro.
3 Centroides
Uma importante ferramenta capaz de resumir cada polígono fechado em uma única localização geográfica é a utilização dos centroides. Eles são frequentemente definidos como o centro de gravidade dos polígonos e têm muitas funções importantes na análise espacial.
Os centroides são frequentemente usados para “representar” o polígono, por exemplo, em cálculos de distância e ao atribuir valores de algum fenômento de interesse para um ponto (por exemplo, renda per capita, incidência de doenças, composição do solo). O efeito dessa atribuição é assumir que a variável de interesse está bem aproximada atribuindo valores a um único ponto.
No R, para os dados de Minas gerais do GBD, os
centroides de todos os municípios podem ser calculados com os
seguintes comandos:
## Obtenção dos centroides
sf_cent_MG <- st_centroid(Dados_MG)
## Mapa dos centroides
ggplot() +
geom_sf(data = Dados_MG, fill = 'white') +
geom_sf(data = sf_cent_MG, color = 'red') + theme_bw() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
3.1 Atividade
Obter os centroides para os municípios do estado do Rio de Janeiro e representar em um mapa.
3.2 Gráfico de médias e medianas
Os gráficos de médias e medianas segundo linhas e colunas permitem investigar a dependência espacial entre vizinhos em duas direções: Norte-Sul e Leste-Oeste. Para construir estes gráficos, utiliza-se as coordenadas dos centróides das áreas, aproximando-as para um espaçamento regular de forma a montar uma matriz.
Calcula-se então as médias e as medianas do fenômeno de interesse ao longo das linhas (eixo Leste-Oeste) e colunas (eixo Norte-Sul) desta matriz. Esta técnica permite identificar a flutuação das medidas ao longo de duas direções, sugerindo a presença de valores discrepantes quando a diferença entre estas é grande, e a tendência ao longo de uma direção quando os valores variam suavemente.
Na figura abaixo, apresenta-se o resultado desta técnica aplicada a dois indicadores socioeconômicos do censo 1991 - renda média do chefe da família e proporção de chefes de família com escolaridade igual ou superior ao segundo grau - para setores censitários da Ilha do Governador, no Rio de Janeiro.
Essa ilha é composta por 225 setores censitários, cujos centróides estão assinalados no primeiro quadro da figura: observe que nas extremidades do “mapa” a quantidade de pontos é muito pequena, e, consequentemente, qualquer medida nesta área será pouco robusta.
No eixo Norte-Sul (colunas), pode-se observar que a renda média do chefe da família apresenta tendência variável, bem menor no centro da região. À mesma coisa acontece para escolaridade, embora com maior flutuação.
No eixo Leste-Oeste (linhas), também parece haver algum deslocamento para valores mais altos no sentido leste, mas o descolamento de médias (x) e medianas (o) sugere a presença de valores extremos dos indicadores.
4 Matrizes de Vizinhança
Para estimar a variabilidade espacial de dados de área, uma ferramenta básica é a matriz de proximidade espacial, também chamada matriz de vizinhança. Dado um conjunto de \(n\) áreas \((A_1,A_2, ..., A_n)\), construímos a matriz W(n x n), onde cada um dos elementos \(w\), representa uma medida de proximidade entre \(A_i\), e \(A_j\), Esta medida de proximidade pode ser calculada a partir de um dos seguintes critérios:
\(w_{ij}=1\) se o centroide de \(A_i\) está a uma determinada distância de \(A_j\); caso contrário \(w_{ij}=0\)
\(w_{ij}=1\) se \(A_i\) compartilha um lado comum com \(A_j\), caso contrário, \(w_{ij}=0\). Nesse caso, a contiguidade pode ser estabelecida de três maneiras:
- \(w_{ij}=l_{ij}/l_i\), onde \(l_{ij}\) é o comprimento da fronteira entre \(A_i\) e \(A_j\) e \(l_{i}\) é o perímetro de \(A_i\).
Como a matriz de proximidade é utilizada em cálculos de indicadores na fase de análise exploratória, é muito útil normalizar suas linhas, para que a soma dos pesos de cada linha seja igual a 1.
Isto simplifica vários cálculos de índices de autocorrelação espacial. A figura abaixo ilustra um exemplo simples de matriz de proximidade espacial, em que os valores dos elementos da matriz refletem o critério de adjacência e foram normalizados.
A idéia da matriz de proximidade espacial pode ser generalizada para vizinhos de maior ordem (vizinhos dos vizinhos). Com critério análogo ao adotado para a matriz de vizinhança de primeira ordem, pode-se construir as matrizes \(W^1\), \(W^2\), …, \(W^k\). Por exemplo, na figura acima, as áreas A e C são vizinhas na matriz de proximidade espacial de ordem 2.
A seguir, seguem os códigos para a construção da matriz de vizinhança para os municípios de Minas Gerais (considerando a contiguidade Queen e Rook) além do mapa representando as conexões:
## Mudando de arquivo tipo SF para SP
Dados_MG_sp <- sf:::as_Spatial(Dados_MG$geometry)
## Calculando os vizinhos pela contiguidade Queen
Queen_MG<-spdep::poly2nb(Dados_MG_sp, queen=T)
## Construindo a matriz de vizinhança Queen
WQueen_MG<-matrix(0,853,853)
for(i in 1:853)
WQueen_MG[i,Queen_MG[[i]]]=1
## Calculando os vizinhos pela contiguidade Rook
Rook_MG<-spdep::poly2nb(Dados_MG_sp, queen=F)
## Construindo a matriz de vizinhança Rook
WRook_MG<-matrix(0,853,853)
for(i in 1:853)
WRook_MG[i,Rook_MG[[i]]]=1
## Plotando as duas estruturas de vizinhança lado a lado
par(mfrow=c(1,2))
sp::plot(Dados_MG_sp, border="lightgrey",main="QUEEN")
sp::plot(Queen_MG, coordinates(Dados_MG_sp), pch = 19, cex = 0.6, add = TRUE, col= "blue")
sp::plot(Dados_MG_sp, border="lightgrey",main="ROOK")
sp::plot(Rook_MG, coordinates(Dados_MG_sp), pch = 19, cex = 0.6, add = TRUE, col= "blue")
4.1 Atividade
Construa a matriz de vizinhança e represente as conexões para os municípios do estado do Rio de Janeiro, considerando as contiguidades Queen e Rook.
5 Média Móvel Espacial
Uma forma simples de explorar a variação da tendência espacial dos dados é calcular a média dos valores dos vizinhos. Isto reduz a variabilidade espacial , pois a operação tende a produzir uma superfície com menor flutuação que os dados originais. A média móvel \(\hat{\mu_i}\) associada ao atributo z, relativo à i-ésima área, pode ser calculada a partir dos elementos \(w\), da matriz normalizada de proximidade espacial \(W\), tomando-se simplesmente a média dos vizinhos:
A figura abaixo ilustra o uso do estimador de média móvel para o percentual de idosos (mais de 70 anos) para os 96 distritos da cidade de São Paulo.
Estes dados são indicadores da grande disparidade social da cidade, com uma grande variação entre o centro (onde a proporção de idosos chega a 8%) com a periferia (onde há várias regiões com menos de 1%).
O valor máximo do percentual de idosos é de 8,2% e o mínimo de 0,8%, com um desvio padrão de aproximadamente 2%. Com a média local, há um alisamento: o valor mínimo é de 1% e o máximo é reduzido a 6,8%. Pode-se notar, ao comparar os dois mapas da figura, que a média móvel local fornece uma visão das grandes tendências do fenômeno em estudo e no caso do percentual de idosos, mostra um forte gradiente centro-periferia.
Nas próximas linhas, tem-se os comandos necessários para fazer a suavização dos dados de mortalidade por acidentes envolvendo motociclistas em Minas Gerais, a partir da média móvel com a contiguidade Queen:
## Criando a matriz de vizinhança Queen normalizada
NQueen_MG<-WQueen_MG/apply(WQueen_MG,1,sum)
## Extraindo o vetor de taxas
taxa<-Dados_MG$taxa
## Médias Móveis
Movel_MG<-NQueen_MG%*%taxa
## Inserindo as médias móveis no shapefile
Dados_MG$Movel<-Movel_MG
breaks <- round(quantile(c(Dados_MG$Movel),probs = c(0.05,0.20,0.40,0.60,0.80,0.95)),2)
Dados_MG$Categoria2 <- age.cat(round(Dados_MG$Movel,2),breaks,min=min(round(Dados_MG$Movel,2)),max=max(round(Dados_MG$Movel,2)))
Mapa_Movel_MG<-ggplot() + geom_sf(data=Dados_MG,aes(fill=Categoria2),colour="black") + theme_bw() +
scale_fill_manual(values = mycolors,name="Estimador de média móvel \n da taxa de Mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
Mapa_MG
Mapa_Movel_MG
5.1 Atividade
Reproduza o exemplo anterior para os dados dos municípios do Rio de Janeiro
6 Tópicos Adicionais
Nesta seção, serão apresentados alguns tópicos adicionais:
- Sobreposição de shapefiles: útil na combinação de mais que uma shapefile para construir mapas mais completos.
- Comparação de Mapas: como plotar dois mapas avaliados em tempos distintos e colocar uma mesma legenda.
6.1 Sobreposição de shapefiles
Considere o exemplo de uma região que não esteja
definida no geobr, como por exemplo a região dos
municípios da Bacia Hidrográfica do Rio Paraopeba, em
Minas Gerais.
O rompimento da Barragem I de contenção de rejeitos da Mina de Córrego do Feijão, no município de Brumadinho, Minas Gerais, em 2019, operada pela Companhia Vale, provocou um dos mais graves desastres do mundo relacionados a barragens de mineração. Considerado o maior acidente de trabalho do Brasil, vitimou fatalmente 272 pessoas, entre as quais 250 funcionários diretos e terceirizados da Vale.
O desastre ocasionou o lançamento de, pelo menos, 12 milhões de metros cúbicos de rejeitos no solo e no Rio Paraopeba, que atingiram mais de 160 km de extensão.
Esses rejeitos resultaram em danos ambientais à vegetação e à fauna, bem como no lançamento de metais pesados na água, como manganês, alumínio, ferro, arsênio, restringindo seu uso e impedindo o abastecimento da água para a região metropolitana de Belo Horizonte.
Nesse curso de impactos ambientais, além de Brumadinho, foram considerados atingidos outros 25 municípios do estado, que correspondem a uma população de cerca de 1,1 milhão de habitantes.
Os comandos a seguir ajudam a como localizar essa região em
Minas Gerais, considerando os shapefiles disponíveis no
geobr:
## Instalando o pacote dplyr
if(!require(dplyr)){ install.packages("dplyr"); require(dplyr)}
## ShapeFile do estado de Minas Gerais
Estado_MG<-read_state(code_state=31,year=2018,simplified=T,showProgress=F)
## ShapeFile dos municípios do estado de Minas Gerais
Municipio_MG<-read_municipality(code_muni=31,year=2018,simplified=T,showProgress=F)
## Municípios da Bacia Hidrográfica do Rio Paraopeba
Municipios_Bacia <- c(3100203, 3106705, 3107000, 3109006, 3109907, 3120904, 3124104,
3125705, 3126000, 3126406, 3130101, 3136652, 3139706, 3140159,
3140704, 3143500, 3146404, 3146909, 3147105, 3147402, 3149606,
3152006, 3161700, 3162922, 3163102, 3169356)
## Criando a shapefile dos municípios da Bacia Hidrográfica do Rio Paraopeba
Bacia <- Municipio_MG %>% filter(code_muni %in% Municipios_Bacia)
## Construção do mapa
ggplot() + geom_sf(data=Estado_MG,color="black")+
geom_sf(data=Bacia,fill="#2D3E50", color="#FEBF57")+
theme_bw() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
Vamos voltar ao exemplo dos dados de taxas de mortalidade padronizadas por idade considerando acidentes envolvendo motociclistas nos municípios do estado de Minas Gerais em 2009-2011. Agora, vamos sobrepor a esse mapa as regiões intermediárias de Minas Gerais, com os seguintes comandos:
## Regiões Intermediárias de Minas Gerais
Intermediate_MG<-read_intermediate_region(code_intermediate = 31,year=2020,simplified=T,showProgress = F)
## Construção do Mapa
Mapa_MG_Int <- ggplot() + geom_sf(data=Dados_MG,aes(fill=Categoria)) +
geom_sf(data=Intermediate_MG,fill=NA,colour="black",lwd=1)+
theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
Mapa_MG_Int
6.2 Comparação de Mapas
Os comandos a seguir, permitem a construção de dois mapas avaliados nos triênios (2000-2002) e (2009-2011) para os municípios do estado do Rio de Janeiro, considerando a mortalidade por acidentes envolvendo motociclistas.
## Leitura dos dados
Dados_RJ_T1_T2 <- read.xlsx("DadosRJ_T1_T2.xlsx") ## Leitura do Shape dos municípios fluminenses
RJ<-read_municipality(code_muni=33,year=2018,simplified=T)## Pontos de corte
breaks <- round(quantile(Dados_RJ_T1_T2$taxa,probs = c(0.05,0.20,0.40,0.60,0.80,0.95)),2)
## Escolha da paleta de cores
mycolors <- rev(brewer.pal(7, "RdBu"))
## Categorizando os pontos de corte
Dados_RJ_T1_T2$Categoria <- age.cat(round(Dados_RJ_T1_T2$taxa,2),breaks,min=min(round(Dados_RJ_T1_T2$taxa,2)),max=max(round(Dados_RJ_T1_T2$taxa,2)))
## Merge do arquivo contendo as taxas com o shapefile
Dados_RJ_T1_T2 <- merge(RJ,Dados_RJ_T1_T2,by="code_muni",all.x=T)## Construção do mapa (função facet faz o plot duplo)
Mapa_RJ_T1_T2 <- ggplot() + geom_sf(data=Dados_RJ_T1_T2,aes(fill=Categoria),colour="black") +
theme_bw() + facet_grid(~trienio)+
scale_fill_manual(values = mycolors,name="Taxa de mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
Mapa_RJ_T1_T2
7 Referências
Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004 (ISBN: 85-7383-260-6).