Curso Geoestatistica

Prévia do material em texto

SERVIÇO PÚBLICO FEDERAL 
UNIVERSIDADE FEDERAL DO PARÁ 
INSTITUTO DE GEOCIENCIAS 
CURSO DE ESPECIALIZAÇÃO EM GEOMINAS 
 
Prof. Joaquim Carlos Barbosa Queiroz 
 
Unidade I: Introdução, Análise exploratória de dados, Descrição 
univariada, Descrição espacial bi-variada. Aplicações 
 
1. Introdução 
A estatística é uma ferramenta bastante utilizada em diversos 
campos da ciência na análise de dados. Dependendo da área em que é 
aplicada, os métodos estatísticos precisam ser considerados segundo as 
características com que os dados se apresentam. Muitas ferramentas 
estatísticas são úteis no desenvolvimento da compreensão qualitativa de 
fenômenos naturais; outras podem ser usadas no desenvolvimento de 
respostas quantitativas para questões específicas. Entretanto, muitos 
métodos estatísticos não utilizam a informação espacial, que é um ponto 
importante em conjuntos de dados utilizados nas ciências da terra. A 
Geoestatística oferece um modo de descrever e quantificar a 
continuidade espacial que é uma característica essencial de muitos 
fenômenos naturais. O formalismo matemático usado na geoestatística 
não é trivial e, mesmo evitando tal formalismo a apresentação não é 
simples. Além disso, sabe-se que a análise de dados de dados das 
ciências da terra podem ser frustrantes e carregados de dificuldades. O 
objetivo deste curso é, portanto, fazer uma introdução aos métodos 
geoestatísticos, sem aprofundamento do formalismo matemático, com 
exemplos práticos e orientação no uso de programas computacionais 
geoestatísticos. 
i 
2. Análise exploratória de dados 
 
A organização e apresentação é um ponto importante na análise 
das características de um conjunto de dados. Em geral se inicia com 
uma descrição univariada, que consiste basicamente em apresentar os 
dados em forma de tabelas de freqüências, histogramas, Box-plot, qq- 
plot, no caso de variáveis quantitativas e gráficos de colunas ou barras e 
pizza, no caso de variáveis qualitativas; e na representação dos dados 
em termos de estatísticas resumo, representadas por medidas de 
localização (média, mediana e moda), medidas de dispersão (amplitude, 
variância, desvio-padrão) e medidas de forma (assimetria e curtose). 
Quando se realiza uma analise em um conjunto de dados 
multivariados em geral se faz um estudo de possíveis relações e 
dependências entre as variáveis. Neste caso, utiliza-se ferramentas de 
analise bi-variada que incluem diagramas de dispersão, covariância, 
correlação e semivariograma. 
 
3. Descrição univariada 
3.1. Conceitos básicos e Notação utilizada: 
 
1. Variável Aleatória (VA) – é uma variável que pode apresentar uma 
série de resultados de acordo com alguma distribuição de probabilidades 
[freqüências] (pi, i = 1,...,N). As N probabilidades de ocorrências devem 
satisfazer as condições (1) pi ≥ 0, para todo i = 1,...,N e (2) ∑
N 
p = 1 . 
i=0 
 
2. Z(u): variável aleatória relacionada a alguma localização no espaço, 
em que u é o vetor de coordenadas da localização. 
3. Função de distribuição acumulada (FDA): caracteriza a distribuição de 
probabilidades de uma VA Z(u) e é definida como 
F(u;z) = Prob{ Z(u) _ z }. (1) 
A FDA apresenta as seguintes propriedades: 
 
a) F(u;z) é não decrescente; 
 
b) F(u;z) ϵ [0,1]; 
 
c) F( – ∞ ) = 0 e F(∞ ) = 1. 
 
4. A função densidade de probabilidade (FDP) é a derivada da FDA, se 
esta for diferenciável, ou seja: 
 
F (u; z + dz) 
_ 
F (u; z) 
f (u; z) = F ´(u; z) = l im 
dz→0 dz 
(2) 
5. Valor esperado – é a média ponderada dos n possíveis resultados, 
onde cada resultado é ponderado por sua probabilidade de ocorrência: 
– no caso discreto: 
 
 
N 1 N 
E(Z ) = m = ∑p z = ∑z 
i i N i i=1 i=1 
(3) 
– no caso contínuo, sob condições de existência das integrais: 
 
 
∞ ∞ 
E(Z ) = m = ∫zdF(u; z) = ∫zdf (u; z)dz ∞ ∞ (4) 
onde F(u;z) e f(u;z) são a FDA e FDP, respectivamente. 
 
6. Variância da VA Z – é definida como o quadrado do desvio esperado 
de Z em relação à sua média 
 
Var(Z) = σ 
2 
= E{(Z 
_ 
m)
2 
} = E[Z}
_ 
m
2
 
N 
= ∑p (z m)2 caso discreto 
i i 
i=1 
 
∞ 
= ∫(z m)2 f (u; z)dz caso contínuo _ ∞ 
 
 
(5) 
A raiz quadrada da variância, , é chamada de desvio-padrão, e 
sua razão em relação à média, σ/m, para variáveis não negativas, é 
denominada coeficiente de variação (ou desvio-padrão relativo) o qual é 
adimensional. Coeficientes variação maiores do que 1 podem indicar a 
presença de valores altos (outliers) na distribuição. 
4. Descrição espacial bi-variada 
4.1. Funções de Covariância, Correlação e Variograma 
 
Em ciências da Terra é freqüentemente importante conhecer o 
padrão de dependência de uma variável X em relação a outra Y, por 
exemplo, quando se deseja saber a relação entre pares de 
concentrações de metais medidas na mesma localização. 
A distribuição conjunta de resultados de um par de variáveis 
aleatórias X e Y é caracterizada pela FDA conjunta (ou bivariada) 
definida como: 
 
 (6) 
 
estimada, na prática, pela proporção de pares de dados conjuntamente 
abaixo dos respectivos valores (valores de corte) x e y. Isto pode ser 
mostrado no diagrama de dispersão (Figura 1) onde cada par de dado 
(xi,yi) é plotado como um ponto. 
O grau de dependência entre as duas variáveis X e Y pode ser 
caracterizado pela dispersão em torno da linha de 45o no diagrama de 
dispersão. A dependência perfeita (X = Y) corresponde a todos os pares 
experimentais (xi,yi), i = 1,..., N plotados sobre a linha de 45o. 
 
 
Figura 1 – Diagrama de dispersão para os pares (xi , yi) 
 
O momento de inércia do diagrama de dispersão em torno da linha 
de 45o – chamado de “semivariograma” do conjunto de pares (xi , yi) – 
é definido como a metade da média das diferenças quadráticas dos 
teores entre as coordenadas de cada par, ou seja: 
 
 (7) 
Quanto maior o valor do semivariograma, maior a dispersão e 
menos relacionadas são as duas variáveis X e Y. 
A covariância centrada (na média), ou simplesmente covariância, 
é definida por: 
 
Cov{XY}=σXY = E{XY} = E{[X – mX][Y– mY]} 
1 N 
estimada na prática por ∑(x _ m )( y _ m ) 
N i X i Y i=1 
 
(8) 
A covariância padronizada (adimensional) entre duas VA’s X e Y é 
conhecida como coeficiente de correlação linear, ou seja: 
 
(9) 
 
A relação experimental entre o semivariograma e a covariância 
pode ser obtida pelo desenvolvimento da equação (7), cujo resultado é 
mostrado a seguir: 
 
 
 (10) 
Para variáveis padronizadas, cujas estatísticas são mX’ = mY’ = 0 
e σ 2 = σ 2 = 1 , teremos 
X ' Y ' 
 
 
(11) 
 
As equações (10) e (11) mostram que um aumento no 
semivariograma é acompanhado de um decréscimo na covariância e 
correlação. Quanto maior for a dispersão dos pares (xi,yi) em torno da 
linha 45o no diagrama de dispersão, maior será o semivariograma e 
menor a covariância e o coeficiente de correlação _XY, o que caracteriza 
o semivariograma como uma medida de variabilidade (ou 
dissimilaridade) e a covariância e correlação como medidas de 
similaridade. 
5. Programas de Geoestatística 
 
A Tabela 1 apresenta as principais características de alguns 
programas de geoestatística mais comunmente utilizados. 
 
Tabela 1. Lista programas de geoestatística 
 
 
5.1 Uso do programa R 
 
. Módulos GeoR, gstat : Fornecem funções para análise de dados 
geoestatísticos usando o programa R. 
 
Observações: 
Os comandos no programa R são apresentados em vermelho 
As saídas no programa R são apresentadas em azul 
O caractere # representa comentários 
 
Símbolos ou comandos importantes 
o Sair do programa q() 
o Salva o trabalho realizado save.image() 
o Lista todos os objetos da área de trabalho atual ls() 
o Remove o objeto x rm(x) 
o Dado ausente (data missing) NA 
o Mostra todos os pacotes instalados library() 
o Carregar (p.ex.) o pacote nlme library(nlme) 
5.1. Instruções básicas para execuçãodo programa R. 
Uso do programa R (Módulos GeoR, gstat) 
Fornecem funções para análise de dados geoestatísticos usando o programa R. 
 Os comandos no programa R são apresentados em vermelho
 As saídas no programa R são apresentadas em azul
 O caractere # representa comentários, representados em preto
 Símbolos ou comandos importantes
o Sair do programa q() 
o Salva o trabalho realizado save.image() 
o Lista todos os objetos da área de trabalho atual ls() 
o Remove o objeto x rm(x) 
o Remove os objetos x e y rm(x,y) 
o Dado ausente (data missing) NA 
o Mostra todos os pacotes instalados library() 
o Carregar (p.ex.) o pacote nlme library(nlme) 
ou require(nlme) (neste caso, se o pacote já 
estiver instalado o comando não é executado) 
o Para a execução de uma linha de comando deve-se pressionar a tecla 
“Enter” 
Instalar o programa R.3.5.1 (disponibilizado na pasta ou na internet: 
https://cran.r-project.org/bin/windows/base/old/3.5.1/ ) 
Depois de instalado usa-se o ícone na área de trabalho 
 Procedimento básico para instalar um pacote
o Escolher espelho CRAN (É necessário estar conectado à internet) 
https://cran.r-project.org/bin/windows/base/old/3.5.1/%20)
 
 
o Instalar (carregar) o pacote desejado 
 
 
 
 
 
 Procedimento comum a todos os ensaios
># Carregar os pacotes necessários (Programa R versão 3.1.1) 
>library(geoR) 
>library(gstat) 
>library(lattice) 
>library(sp) 
 
#Definir a área de trabalho: pasta do computador onde estao 
os arquivos que serão utilizados 
setwd("C:/PPGCA/Geominas/CursoR/Ensaios") 
 
6. Ensaios práticos 
Para a realização de ensaios práticos será utilizado o banco de 
dados jura disponível em Goovaerts, 1997 (em anexo). Esse banco 
consiste de 9 variáveis, sendo 7 variáveis contínuas (concentrações de 
Cd, Cu, Pb, Co, Cr, Ni e Zn e 2 variáveis categóricas (Uso da terra – 
florestas(1), pastos(2), prados(3) e lavoura(4) – e tipos de rocha – 
Argoviano(1), Kimeridgiano(2), Sequaniano(3), Portlandiano(4) e 
Quaternário(5)). 
 
6.1. Aplicações: 
Banco de dados Jura (Goovaerts, 1997): 
• Contem 9 variáveis: 
• Concentrações (em ppm) dos metais de Cd, Cu, Pb, Co, Cr, Ni e Zn 
(Contínuas) 
• Uso da terra e tipo de rocha (Categóricas) 
TABELA 1. Estatísticas descritivas e Tolerância máxima 
 
Estatísticas Cd Co Cr Cu Ni Pb Zn 
Média 1.31 9.3 35.07 23.73 19.73 53.92 75.08 
Mediana 1.07 9.76 34.84 17.6 20.56 46.4 73.56 
Desvio padrão 0.92 3.58 10.96 20.71 8.23 29.79 29.02 
Variância amostral 0.84 12.79 120.07 429.01 67.78 887.57 842.12 
Curtose 2.58 -0.8 0.02 12.12 0.24 11.6 2.16 
Assimetria 1.51 -0.18 0.29 2.88 0.16 2.91 1.03 
Intervalo 4.99 16.17 58.88 162.44 49 210.6 194.12 
Mínimo 0.14 1.55 8.72 3.96 4.2 18.96 25.2 
Máximo 5.13 17.72 67.6 166.4 53.2 229.56 219.32 
Tolerancia máxima 0.8 25 75 50 50 50 200 
Contagem 259 259 259 259 259 259 259 
Ensaio 1 (Script: Ensaio_1): 
 Descrição univariada e bi-variada.
 Gráficos das variáveis quantitativas (histograma, boxplot, qqplot):
 Gráficos das variáveis qualitativas
 Construção do mapa de localizações das amostras
a) Variáveis selecionadas: Cadmium e Niquel (quantitativas) e uso da 
terra (qualitativa) 
 
PROGRAMA R 
# Procedimento comum a todos os ensaios 
# Carregar os pacotes necessários (Programa R versão 3.2.2) 
>library(geoR) 
>library(gstat) 
>library(lattice) 
>library(sp) 
 
>#Definir a área de trabalho 
>setwd("C:/PPGCA/Geominas/CursoR/Ensaios") 
 
#Carregar podem ser lidos no formato ASCII (.txt) 
>Geo <- read.table("jura.txt", head=T) #banco de dados 
 
>names(Geo) #mostrar o nome das variáveis do banco Geo 
 
[1] "Landuse" "Rock" "Cd" "Co" "Cr" "Cu" "Ni" 
[8] "Pb" "Zn" 
>summary(Geo) #Resumo do arquivo Geo 
# Selecionar a variavel Cadmium 
>Cd <- as.geodata(Geo, data.col = 3) 
 
# Descrição Univariada 
 
#Principais estatísticas descritivas da variável Cadmium 
>summary(Cd$data) 
 
Min. 1st Qu. Median Mean 3rd Qu. Max. 
0.1350 0.6375 1.0700 1.3090 1.7150 5.1290 
 
#Histograma 
>hist(Cd$data, breaks = 20, main = "Cadmio", xlab = "Cd", 
ylab = "Freqüência", col = 4) 
 
 
#Box PLot 
>boxplot(Cd$data, col = "gray", main = "Cadmio", ylab = 
"CdVal") 
 
 
 
 
#QQ Plot 
>qqnorm(Cd$data, main = "Cadmio", xlab = "Cd", ylab = 
"CdVal") 24 
 
 
 
 
# Descrição Bi-variada 
# Selecionar a variavel Niquel 
>Ni <- as.geodata(Geo, data.col = 9) 
>summary(Ni$data) #Principais estatísticas descritivas do Ni 
Min. 1st Qu. Median Mean 3rd Qu. Max. 
4.20 13.80 20.56 19.73 25.42 53.20 
 
#Diagrama de dispersao 
>plot(Ni$data,Cd$data, pch="*") 
 
 
 
#Covariancia 
>cov(Ni$data,Cd$data) 
[1] 3.672182 
#Correlação 
>cor(Ni$data,Cd$data) 
[1] 0.4873752 
 
#Teste de correlação 
>cor.test(Ni$data,Cd$data,method=c("pearson")) 
Pearson's product-moment correlation 
data: Ni$data and Cd$data 
t = 8.9479, df = 257, p-value < 2.2e-16 
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval: 
0.3885697 0.5750997 
sample estimates: 
cor 
0.4873752 
 
#Conclusão: A correlação entre o Niquel e o Cadmio é 
significatia, pois o p-value, p < 0.001. 
 
>#Mapa de pontos 
>points(Cd, xlab = "W-E", main = "Cadmio", ylab = "N-S", 
cex.min = 1.5,pt.divide="quart",x.leg=c(0.2,1.8),y.leg=c(4.3, 
5.7)) 
 
 
Variável selecionada (Uso da terra (Landuse) : categórica 
Categorias : florestas(1), pastos(2), prados(3) e lavoura(4) 
>#Variavel qualitativa Landuse 
># Selecionar a variavel Landuse 
>#carregar arquivo de dados categoricos (uso da terra) 
 
>land <- read.table("land.txt", head=T) 
>land[1:5,] # mostrar os dados na area de trabalho 
X Y Landuse Flore Pasto Prado Lavoura 
1 2.699 1.199 1 1 0 0 0 
2 4.232 1.588 1 1 0 0 0 
3 4.459 2.624 1 1 0 0 0 
4 2.916 4.916 1 1 0 0 0 
5 3.649 2.674 1 1 0 0 0 
>summary(land) 
 
X Y Landuse Flore 
Min. :0.626 Min. :0.580 Min. :1.000 Min. :0.0000 
1st Qu.:2.282 1st Qu.:1.487 1st Qu.:2.000 1st Qu.:0.0000 
Median :3.043 Median :2.581 Median :3.000 Median :0.0000 
Mean :2.980 Mean :2.665 Mean :2.548 Mean :0.1274 
3rd Qu.:3.665 3rd Qu.:3.752 3rd Qu.:3.000 3rd Qu.:0.0000 
Max. :4.920 
Pasto 
Min. :0.0000 
1st Qu.:0.0000 
Median :0.0000 
Mean :0.2162 
3rd Qu.:0.0000 
Max. :1.0000 
Max. :5.690 
Prado 
Min. :0.0000 
1st Qu.:0.0000 
Median :1.0000 
Mean :0.6371 
3rd Qu.:1.0000 
Max. :1.0000 
Max. :4.000 
Lavoura 
Min. :0.0000 
1st Qu.:0.0000 
Median :0.0000 
Mean :0.0193 
3rd Qu.:0.0000 
Max. :1.0000 
Max. :1.0000 
 
># Figuras (Grafico de pizza e coluna) 
 
>attach(land) #vincula os nomes das variáveis ao caminho de busca 
>tland <- table(Landuse) #cria tabela da variavel Landuse 
>tland <- as.data.frame(tland) #tabela com varios tipos de dados 
>tland$Porc <- round(tland$Freq/nrow(land)*100, dig = 2) # cria variavel 
percentagem 
>tland$Landuse <-c("Floresta", "Pasto", "Prado", "Lavoura") #rótulos na 
variavel Landuse 
 
>#Grafico de Pizza 
>pie(tland$Freq,labels=paste(tland$Porc, "%"), col =c(4,5,6,7 )) 
>legend("topleft", legend=tland$Landuse, fill = c(4,5,6,7), 
horiz = FALSE, bty = "n", cex=1) 
 
 
 
 
 
 
 
 
 
 
 
 
>#Gráfico de Colunas 
>barplot(tland$Freq,labels=paste(tland$Porc,"%"),col=c(4,5,6, 
7)) 
>legend("topleft",legend=tland$Landuse,fill=c(4,5,6,7),horiz= 
FALSE, bty = "n", cex=1) 
 
Exercício 1: Escolha uma variável quantitativa do banco de dados jura.txt. 
a) Faça o mapa de localização das amostras e gráficos das estatísticas descritivas 
(histograma, boxplot e qqplot para a variavel quantitativa e tabelas de frequencias 
b) Faça o grafico de barras e gráfico de pizza da variável categórica (tipos de rocha) 
do banco Jura. 
 
Siga os mesmos procedimentos do Ensaio_1 
 
 
 
Unidade II: Inferência estatística, Análise variográfica, Modelo linear de 
regionalização, Modelo linear de co-regionalização 
1. Inferência estatística 
O principal objetivo da inferência estatística, neste caso, é a caracterizaçãodo grau de dependência espacial entre duas variáveis aleatórias Z(u) e Z(u+h) 
separados de um vetor h. 
As VA’s X e Y podem representar: 
(i) duas diferentes propriedades medidas em uma mesma localização, por exemplo, 
concentração (ou teor) de arsênio e manganês no solo; 
(ii) uma mesma propriedade medida em duas diferentes localizações do espaço, 
característica das variáveis regionalizadas, por ex., concentração de arsênio nas 
localizações x e x + h separadas por um vetor h, sendo: X = Z(x), Y = Z( x + h); 
(iii) duas diferentes propriedades medidas em duas diferentes localizações, por 
exemplo, concentração de arsênio na localização x e concentração de manganês 
separada por um vetor h, sendo: X = ZAs(x), Y = ZMn(x + h). 
Em todos os casos, o semivariograma γXY ou a correlação ρXY medirão o grau 
de variabilidade/similaridade entre as duas VA’s X e Y. 
O segundo caso (ii) é de particular interesse em problemas de interpolação 
espacial, onde se deseja inferir (mapear) uma determinada área para uma 
determinada propriedade, Z(u), u ϵ área A, que apresenta continuidade espacial, a 
partir de uma amostra n de Z(u). A Combinação de todos os n(h) pares de dados de 
Z(u) sobre a mesma área/zona/camada/população, com tais pares separados 
aproximadamente pelo mesmo vetor h (em comprimento e direção), permite 
estimar o semivariograma característico (ou experimental) da variabilidade espacial 
na área A: 
 
 
 
(2.1) 
2. Análise variográfica 
Para o modelamento do semivariograma, realizado depois de construído o 
semivariograma experimental, é necessário que a hipótese de estacionariedade 
seja considerada. Esta hipótese estabelece, em resumo, que os dois primeiros 
momentos (média e variância) da diferença [Z(u) – Z(u+h)] são independentes da 
localização u e função somente do vetor h. O segundo momento dessa diferença 
corresponde ao semivariograma, ou seja, 
2γ(h) = E{[Z(u) – Z(u+h)]2} é independente de u ϵ A (2.2) 
O desenvolvimento da equação acima conduz a 
γ (h) = Co – C(h) (2.3) 
 
onde, 
Var{Z(u)} = Var{Z(u+h)} = σ2 = Co para todo u ϵ A 
Cov{Z(u),Z(u+h)} = C(h) para todo u ϵ A 
A relação (2.3) é utilizada então para a determinação do modelo 
 
semivariográfico. A variância Co é chamada na geoestatística como patamar (ou 
sill). O semivariograma é a ferramenta preferida para inferência estatística porque 
oferece algumas vantagens sobre a covariância, entre as quais (CHRISTAKOS, 1984): 
(1) seu cálculo empírico está sujeito a erros menores; 
(2) oferece melhor caracterização da variabilidade espacial 
i j 
(3) exige a chamada hipótese de estacionaridade intrínseca, ou seja, que Z(u) é 
uma função aleatória com incrementos Z(u+h) – Z(u) estacionários, mas não 
necessariamente ela própria estacionária. 
CHRISTAKOS (1984) mostra ainda que, em geral, para uma função contínua 
escolhida como semivariograma, é necessário e suficiente satisfazer a propriedade 
de positivo definida que garante que a variância de certas combinações lineares 
n 
descritas por Y = ∑λi Z (u i ) , sejam positivas. O modelo semivariográfico γ(h) deve 
i=1 
 
ser condicionalmente negativo definido, isto é, 
 
n n n 
∑∑λ λ γ(h) ≥0 quando ∑λi = 0 ∀n 
i=1 j =1 i=1 
(2.4) 
 
Outra condição a ser satisfeita pelo semivariograma é: 
lim γ(h) / | h |2 = 0 
|h|→∞ 
 
(2.5) 
 
Uma vez que não é fácil verificar essas condições diretamente 
(CHRISTAKOS, 1984), a solução prática é usar combinações lineares de modelos 
básicos que são válidos, isto é, permissíveis. Os modelos básicos, mais utilizados 
em geoestatística serão apresentados. 
 
2.1- Variografia Experimental 
O gráfico do semivariograma experimental, γ(h) versus h, dado pela 
equação (2.1), mostrará uma série de pontos discretos de γ(h) correspondendo a 
cada valor de h, e para o qual, uma função contínua (modelo) deve ser ajustada. 
Na construção do semivariograma experimental recomenda-se utilizar um 
mínimo número de pares por lag em torno de 20 a 30 e a distância máxima para h 
deve ser metade da separação máxima 
Pode-se construir o semivariograma experimental em qualquer direção. 
Apesar do afastamento entre os pontos ser o único parâmetro definido para se 
calcular o semivariograma, a escolha do vetor de afastamento entre os pares 
afetará na direção seguida pelo semivariograma. Esta direção poderá ser 
horizontal, vertical ou diagonais. Para o caso tri-dimensional, o semivariograma é 
construído, incluindo uma terceira dimensão. A Figura 2.1 mostra 
esquematicamente a construção do semivariograma experimental no caso de 
amostras regularmente espaçadas. Cada par de pontos (γ(h), h) corresponde à 
soma das diferenças ao quadrado, dividido por 2*N(h), média, dos valores da 
variável Z em espaçamentos (lags), h, diferentes. 
 
 
 
 
 
 
 
Figura 2.1. Amostras regularmente espaçadas 
A Figura 2.2 mostra um exemplo do cálculo de um ponto do 
senmivariograma em duas direções, N-S e S-W, para a distancia de 25 m. 
 
* 
N -S 
= 
1 [(10_11 )2 +(12_11 )2 +(11_11 )2 +(10_10 
18 
)2 +(10_10 )2 +(11_13 )2 + 
+ (13_13 )2 +(13_11 )2 + (11_12 )2 ] = 0,6111 
 
 
* 
E -W 
= 
1 [(12 _11 )2 +(13_12 
18 
)2 +(11_10 )2 +(10_11 )2 +(11_11 )2 +(11_11 )2 + 
+ (12_10 )2 +(10_14 )2 + (14 _13 )2 ]= 1,4444 
γ 
γ 
 
 
Figura 2.2. Exemplo de dados para amostras regularmente espaçadas 
 
 
Porém, na prática são raros os experimentos em que as amostras são 
regularmente espaçadas. DEUTSCH E JOURNEL (1997, p. 48-49), apresentam um 
método de seleção de pares para cálculo do semivariograma em amostras 
irregularmente espaçadas, em até três dimensões. Algo parecido com as 
características físicas de um feixe de luz (Figura 2.4), onde alguns critérios devem 
ser observados e alguns parâmetros definidos pelo pesquisador, de acordo com o 
objeto e a área de estudo, tais como: 
 Azimute ou ângulo de rotação horizontal, que é o ângulo em relação ao eixo 
norte no sentido horário; 
 Inclinação ou ângulo de rotação vertical, que é o ângulo em relação a 
elevação no sentido horário; 
 Tolerância angular horizontal, que é a metade da “janela” de abrangência 
horizontal; 
 Tolerância angular vertical, que é a metade da “janela” de abrangência 
vertical (a tolerância angular horizontal e a tolerância angular vertical 
representam os raios de uma elipse ou circulo, caso sejam iguais, que 
determina a forma do foco do feixe de luz); 
 Afastamento (lag) entre os pares de pontos; 
 A tolerância dos afastamento, que define os limites para cada distância, 
geralmente é a metade do afastamento, para que não haja interseções; 
 A distância máxima; 
 Banda horizontal, é o raio delimitador horizontal para a tolerância angular 
horizontal; 
 Banda vertical, é o raio delimitador vertical para a tolerância angular vertical 
(a banda horizontal e vertical fixam o foco do _feixe de luz_ até a distância 
máxima); 
Por convenção os ângulos mais utilizados para definir o azimute e a inclinação, 
caso haja, são 0º, 45º, 90º e 135º.. 
 
Figura 2.4. Modelo de seleção de pares para cálculo do semivariograma em 
amostras irregularmente espaçadas em três dimensões. 
 
A expectativa é que o semivariograma mantenha o mesmo padrão para 
qualquer que seja ângulo escolhido, sofrendo influência apenas do afastamento 
entre os pares de pontos. Muitas vezes, na prática, isto não ocorre. Quando o 
semivariograma sofre influência da direção dos pares, diz-se que o processos é 
anisotrópico. A anisotropia é uma característica física encontrada em diversos tipos 
de materiais, substâncias ou fenômenos da natureza. Logo, não é raro a ocorrência 
deste fenômenos quando se analisa variáveis ambientais e geológicas. Quando o 
semivariograma não tem influencia da direção dos pares, o processo é isotrópico. 
Por vezes, faz-se necessário analisar a variabilidade espacial de v.a 
categóricas ou categorizadas. Nestes casos usa-se o semivariogramaindicador, que 
consiste em calcular o semivariograma experimental tradicional com as seguintes 
funções indicadoras : 
 
 
 
(2.6) 
quando z(ui) é uma v.a. categórica, onde zcat é o código ou índice da categoria de 
interesse, e 
 
 
 
(2.7) 
para v.a. quantitativas, onde zcut é uma valor de corte para v.a. quantitativa Z(u). 
Portanto, o semivariograma indicador é definido como: 
 
 
 
 
(2.8) 
e mede a proporção de pares de pontos separados pelo vetor h, em que ocorreu o 
evento de interesse, no caso z(u) = zcat ou z(u) ≤ zcut. 
Teoricamente, o semivariograma é zero se a distância entre o pares de 
pontos for zero, mas na prática, geralmente ocorre uma descontinuidade na 
origem do semivariograma, chamado de Efeito pepita. Ou seja, a medida que h 
tende a zero o semivariograma assumirá um valor positivo Co. Isto se deve a erros 
amostrais ou a presença de variabilidade em curtas escalas, onde amostras 
separadas em pequenas distância pode apresentar diferenças significantes. (ISAAKS 
& SRIVASTAVA, 1990, p. 143). 
Outro conceito importante, ao se analisar o comportamento do 
semivariograma em função de um determinado intervalo de distâncias h's, é que a 
partir de uma determinada distância ϕ os pares de pontos z(u) e z(u + h) deixam de 
influenciar-se mutuamente (ρ(h) = 0) e a variabilidade permanece constante. Esta 
distância a é denominada de Alcance. Já ao valor que o semivariograma assume 
quando está em função de uma distância igual ou superior ao alcance, a, é 
denominado de Patamar. A Figura 2.5 mostra os parâmetros citados na estrutura 
do semivariograma ajustado. 
Figura 2.5: Parâmetros do semivariograma 
 
 
2.2- Modelamento variográfico 
Para se realizar predições espaciais por meio de krigagem ou simulação 
estocástica, precisa-se substituir o semivariograma experimental por um modelo 
de semivariograma aceitável. Parte da razão para isto é que o algoritmo de 
krigagem terá acesso à valores de semivariância para distâncias diferentes das 
utilizados no semivariograma experimental. 
2.2.1 Modelagem dos Semivariogramas e da Anisotropia 
Os modelos mais utilizados em geoestatística são apresentados a seguir: 
 Modelo Esférico: É o modelo mais utilizado na geoestatística, tanto para uma, 
como para duas ou três dimensões. Tem um comportamento linear para 
distâncias próximas a origem, mas converge rapidamente em distâncias mais 
 
longas, quando h a e é considerado um modelos transitivo, pois atinge o 
patamar. Definido como: 
 
(2.9) 
 
 
 
 Modelo Exponencial 
 
 
 
(2.10) 
 Modelo Gaussiano 
 
 
 
 
(2.11) 
A Figura 2.6 representa o comportamentos dos modelos teóricos esférico, 
exponencial e gaussiano, submetidos aos mesmos parâmetros, efeito pepita ( Co), 
patamar parcial (C) e alcance (a). Observa-se que o modelo esférico realmente 
atinge o valor especificado do patamar no alcance especificado. Já o modelo 
exponencial e o modelo gaussiano, atingem o patamar assintoticamente. 
Figura 2.6: Comportamento dos modelos teóricos esférico, exponencial e 
gaussiano, submetidos aos mesmos parâmetros. 
Em alguns casos, o semivariograma apresenta padrões ou escalas 
diferentes dependendo da distância entre os pares. Então, para se obter um 
melhor ajuste, faz-se necessário combinar dois ou mais modelos teóricos, 
chamados modelos aninhados. 
 Modelo Potência: É um modelo não transitivo (Figura 2.7), definido como: 
 (2.12) 
 
onde 0 < w < 2 e c é o coeficiente de declividade. 
 
 
Figura 2.7. Comportamento do modelo teórico potência, com w = 0,5, w = 1,0 e 
w = 1,5, e c = 1. 
 
A anisotropia pode ser modelada para minimizar o erro de mensuração da 
variabilidade espacial e, consequentemente, os erros nas predições. Primeiro deve- 
se detectar a existência da mesma, construindo-se o semivariograma para 
diferentes direções (Figura 2.8). Isto vale tanto para a rotação horizontal, através 
do azimute, quanto para a rotação vertical, através da inclinação. Depois se verifica 
se há diferenças significativas no alcance e no patamar dos semivariogramas para 
cada eixo. Existem três tipos de anisotropia: (a) geométrica, em que o alcance dos 
eixos na direção de anisotropia são diferentes, enquanto o patamar permanece 
constante (Figura 2.8 a); (b) zonal, neste caso, o alcance dos eixos na direção de 
anisotropia permanecem iguais enquanto os patamares são distintos (Figura 2.8 b); 
(c) anisotropia combinada, nesta caso, tanto o alcance quanto o patamar do 
semivariograma mudam na direção do eixo de anisotropia (Figura 2.8 c). 
No caso tridimensional, em que há indícios de anisotropia, pode-se aplicar 
uma transformação nas coordenadas originais, antes de calcular o semivariograma 
experimental. Os eixos das coordenadas transformadas devem ficar alinhados com 
as direções de anisotropia. Isto permite um calcular normalmente o 
semivariograma experimental unidirecional ao longo dos eixos da anisotropia 
(ISAAKS; SRIVASTAVA, 1990, p. 181). 
 
 
 
 
 
 
 
Figura 2.8: (a) Modelo semivariográfico padrão de anisotropia geométrica. (b) 
Modelo semivariográfico padrão de anisotropia zonal. (c) e (d) Modelos 
semivariográficos padrões de anisotropia combinada 
l 
L L 
3. Modelo Linear de Regionalização 
Primeiro são estabelecidas as condições (estacionariedade e modelo 
semivariográfico condicionalmente negativo definido) que qualquer modelo de 
covariância ou semivariograma devem satisfazer e, então o modelo linear de 
regionalização é introduzido. 
Procedimento: 
1. Ajuste da forma do semivariograma experimental ou função de covariância : dois 
ou mais modelos básicos, g(h) ou c(h) devem ser combinados (aninhados). 
2. Entretanto, nem todas as combinações de semivariogramas permissíveis ou 
modelos de covariância resultam em uma função de covariância ou 
semivariograma permissivel. 
3. O modo mais fácil de construir um modelo permissivo consiste primeiro em 
construir uma função aleatória. A função de covariância ou semivariograma da 
função aleatória é então, por definição, permissível. 
4. O MLR constrói a função aleatória Z(u) como uma combinação linear de (L + 1) 
funções aleatórias independentes Yl(u) cada uma com média zero e função básica 
de covariância cl (h) 
 
 
Z (u)  
l 0 
a
l
Y 
l 
(u)  m 
onde 
E[Z (u)]  m 
E[Y
l 
(u)]  0 
Cov[Y 
l 
(u).Y 
l ' 
(u  h)]  c (h) 
 0 
 
 
 
se l  l' 
caso contrário 
 
independência mútua (duas a duas) das (L + 1) funções aleatórias Yl(u) 
A decomposição da equação acima requer que a função de covariância da 
função aleatória Z(u) seja uma combinação linear das (L + 1) funções de 
covariâncias básicas : C(h)  CovZ (u), Z (u  h)

C(h)  al al ' CovY l (u),Y l ' (u  h)
l 0 l '0 
L 
l 
C(h)  b c ( h) 
l 
2 
 (h)  b g ( h) l 
L L 
C(h)  a l alc ( h) 
l 0 l '0 
L 
l 
l 
 
com 
 
b
l 
 (a
l 
)
2 
 0 
l '0 
 
bl : sill (patamar) positivo do modelo de covariância básico 
Um desenvolvimento similar é feito em termos de semivariograma : 
Para gl(h), semivariograma da FA, Y
l 
(u) , em que o semivariograma cruzado entre 
duas FA’s diferentes, Yl (u) e Y l ' (u) é igual a zero: 
E[Y l (u)  Y l (u  h)].[Y l ' (u)  Y l ' (u  h)] c (h) se l  l' 
 0 caso contrário 
 
O modelo de semivariograma  (h) é então expresso como uma combinação 
 
linear positiva dos modelos básicos de semivariograma gl(h) 
 (h)  1 EZ (u), Z (u  h)2 

L 
l 
com b
l 
 (a
l 
)
2 
 0 
l 0 
 
 
Prática do modelamento 
1. Semivariograma experimental ou valores de covariância 
2. Semivariogramas ou modelos de covariância permissíveis 
3. Utilização de informação secundária, tais como, conhecimento físico da área 
Processo de modelamento: baseado em decisões do usuário: 
1. Escolha de modelos isotrópicos ou anisotrópicos 
– Usar no mínimo três direções 
– Usar mapas de semivariogramas 
2. Qual o numero (L + 1) e tipos de modelos semivariográficos basicos usar– Evitar o uso de semivariogramas se suporpondo. O objetivo é capturar as 
principais características espaciais do atributo. Quando modelos aninhados 
diferentes fornecem ajustes similares, selecionar o mais simples. Em geral, modelos 
mais complicados não conduzem a estimativas mais precisas. 
3. Quais os parâmetros (sill, alcance, anisotropia,.) de cada modelo básico usar 
– Efeito pepita (nugget): determinado pela extrapolação do comportamento dos 
primeiros valores do semivariograma em relação ao eixo vertical. 
– Sill (Patamar): contribuição de estruturas não aleatórias. Uma pratica não 
confiável é forçar a soma dos (L+1) sills a se igualar à variância amostral 
– Alcance: para semivariogramas experimentais bem comportados é facil de inferir 
em modelos de transição. 
 
4. Modelo linear de co-regionalização 
Usado em situações em que se tem variáveis correlacionadas. Por exemplo, em 
ciências do solo, pode-se ter a condutividade hidráulica e/ou retenção de água 
(difíceis e caras para se medir) com teores de partículas na camada superficial do 
solo (fáceis de medir). 
 
4.1. Semivariogramas cruzados 
Sejam {Z1(u1i), i=1,..., N1} e {Z2(u2j), j=1,..., N2} realizações particulares das 
funções aleatórias Z1(u1i) e Z2(u2j). Em geral, Z2 é subamostrada, em relação a Z1, 
i.e.,, N1 > N2. Então Z2 deverá ser estimado para os locais não amostrados. 
Os semivariogramas diretos e cruzados de Z1(u1i) e Z2(u2j) são dados por 
 
1 2 
 11 (h) = 
2 
E {Z1 (u1 j + h)  Z1 (u1 j )} 
 
1 2 
 22 (h) = 
2 
E {Z 2 (u2 j + h)  Z 2 (u2 j )} 
 12 (h) = 
1 
21 (h) = 
2 
E {[Z1 (u1 j + h)  Z1 (u1 j )][Z 2 (u2 j + h)  Z 2 (u2 j )} 
k ij k 
ij 
2 
(u 
2 j 
Estimativa do Semivariograma cruzado 
 
(h) 
1 
 
 
2N (h) 
N

( h )
Z
 
i1 
 
 h)  Z )Z  h)  Z )


 Semivariograma : caso particular do semivariograma cruzado, quando as duas 
variáveis são idênticas 
 Semivariograma cruzado : calculado usando as informações existentes para 
posições geográficas coincidentes 
Características 
– não é obvio que o valor do semivariograma cruzado para h=0, deva ser nulo 
– alcance : representa o final ou a distância máxima de dependência espacial entre 
as variáveis 
– patamar : se existir, deve aproximar-se da covariância entre as duas variáveis 
– Os modelos utilizados para o semivariograma cruzado são os mesmos já 
requeridos para o semivariograma 
4.2. Modelamento da co-regionalização 
Covariância multivariada ou inferência do semivariograma fornece um 
conjunto de matrizes Ĉ(hk )  Ĉij (hk ) ou ̂ (h )  ̂ (h )
Número de modelos de semivariograma (ou covariância) diretos e cruzados 
é dados por, 
Nv (Nv 1) 
2 
A Maior dificuldade é que os modelos não podem ser construídos 
independentemente um do outro e existem condições que qualquer matriz de 
modelos semivariográficos ou covariância devem obedecer. Entretanto, o 
programa R tem uma função (fit.lmc) que ajusta os variogramas. 
S 
Notação matricial: (h)  Bs
s0 
s (h) 
Bs = [bs ] : matriz de co-regionalização 

12 1 
(u 
1 j 1 
(u 
1 j 2 
(u 
2 j 
ii 
b 
ij 
s 
Condições suficientes para que o chamado modelo linear de co-regionalização 
(MLC) seja permissível: 
1. Funções s(h) : modelos semivariográficos permissíveis ; 
2. Matriz de co-regionalização é positivo-definida. Isso implica que os coeficientes 
devem satisfazer as seguintes restrições: 
a) bs  0 
b) jj  0 
c) | bs |


Ensaio 2 (Script: Ensaio_2): 
 Análise variográfica 
 Modelo linear de regionalização 
 Modelo linear de co-regionalização 
a) Variáveis selecionadas: Niquel e Cobalto (quantitativas) 
PROGRAMA R 
# Carregar os pacotes necessários (Programa R 3.5.1) 
library(geoR) 
library(gstat) 
library(sp) 
 
#Definir a área de trabalho 
setwd("D:/PPGCA/Hidrogeo") 
#os dados de entrada podem ser lidos no formato .txt 
Geo <- read.table("jura.txt", head=T) #banco de dados 
names(Geo) #mostrar o nome das variáveis do banco Geo 
Geo[1:5,] #mostrar 5 primerias linhas do arquivo Geo 
summary(Geo) #Resumo do arquivo Geo 
 
### Selecionar as variaveis Niquel(coluna 9) e 
### Cobalto(coluna 6) 
 
Ni <- as.geodata(Geo, coords.col = 1:2, data.col = 9) 
Ni #mostrar a variavel Niquel area de trabalho 
 
Co <- as.geodata(Geo,coords.col = 1:2, data.col = 6) 
Co #mostrar a variavel Cobalto na area de trabalho 
b s .b
s 
ii jj 
#### 1. Analise Variográfica ### 
################################ 
### 1.Verificação da estacionaridade (ausencia de tendencia) 
# Grafico do Niquel na direção x (W-E) 
plot(Ni$coords[, 1], Ni$data, xlab = "W-E", ylab = "Niquel", 
pch = 20, cex = 1.5, main = "Niquel") 
lines(lowess(Ni$data ~ Ni$coords[, 1])) 
 
# Grafico do Niquel na direção y (N-S) 
plot(Ni$coords[, 2], Ni$data, xlab = "S-N", ylab = "Niquel", 
pch = 20, cex = 1.5, main = "Niquel") 
lines(lowess(Ni$data ~ Ni$coords[, 2])) 
 
# Grafico do Cobalto na direção x (W-E) 
plot(Co$coords[, 1], Co$data, xlab = "W-E", ylab = "Cobalto", 
pch = 20, cex = 1.5, main = "Cobalto") 
lines(lowess(Co$data ~ Co$coords[, 1])) 
 
# Grafico do Cobalto na direção y (N-S) 
plot(Co$coords[, 2], Co$data, xlab = "S-N", ylab = "Cobalto", 
pch = 20, cex = 1.5, main = "Cobalto") 
lines(lowess(Co$data ~ Co$coords[, 2])) 
 
#plotando os 4 graficos juntos 
par(mfrow=c(2,2)) # Dividir a tela em 2 linha e 2 colunas 
plot(Ni$coords[, 1], Ni$data, xlab = "W-E", ylab = "Niquel", 
pch = 20, cex = 1.5, main = "Niquel") 
lines(lowess(Ni$data ~ Ni$coords[, 1])) 
 
plot(Ni$coords[, 2], Ni$data, xlab = "S-N", ylab = "Niquel", 
pch = 20, cex = 1.5, main = "Niquel") 
lines(lowess(Ni$data ~ Ni$coords[, 2])) 
 
plot(Co$coords[, 1], Co$data, xlab = "W-E", ylab = "Cobalto", 
pch = 20, cex = 1.5, main = "Cobalto") 
lines(lowess(Co$data ~ Co$coords[, 1])) 
 
plot(Co$coords[, 2], Co$data, xlab = "S-N", ylab = "Cobalto", 
pch = 20, cex = 1.5, main = "Cobalto") 
lines(lowess(Co$data ~ Co$coords[, 2])) 
 
 
 
#### !!! Dados estacionarios #### 
 
### 2. Verificação de anisotropia 
## 2.1. Semivariogramas experimentais direcionais 
#Obs: depois do comando em parenteses é o nome da subrotina 
# Comando variog4 (geoR): computa o semivariograma em 4 
direções para avaliação da anisotropia 
# uvec: vetor de parametros: min: 0 max: 2 (km); 
# length: distancia maxima=20; Portanto teremos 20 pontos 
cujos passos serao 2/20 = 0.1 km 
# direction: define as direções em um angulo (tolerance) de 
45 graus (degrees). 
 
#Variavel Niquel 
v4Ni <- variog4(Ni, uvec = seq(0, 2, length = 20), 
direction = c(0, 45, 90, 135), tolerance = 45, 
unit.angle = "degrees") 
v4Ni 
par(mfrow=c(1,1)) 
plot(v4Ni) #plotar o semivariograma experimental 
nas 4 direções 
title("Niquel") # titulo do grafico 
 
#Variavel Cobalto 
v4Co <- variog4(Co, uvec = seq(0, 2, length = 20), direction 
= c(0, 45, 90, 135), tolerance = 45, 
unit.angle = "degrees") 
v4Co 
plot(v4Co) 
title("Cobalto") 
 
## 2.2 Semivariograma de superficie: verificar a anisotropia 
# Comando trellis.par.set (lattice): padrao de cores 
# Comando gstat(gstat): gerar objeto para construir o 
# variograma experimental 
#NULL : novo objeto 
#id: nome da variavel 
#formula: equação da variavel dependente 
#data: nome do arquivo de dados 
#Comando variogram(gstat): semivariograma experimental na 
# subrotina gstat 
#vgni: objeto gerado no gstat 
#cuttof=2 : maxima distancia ou separaçao em que os pontos 
#estimados sao incluidos 
#width=.2 : distancia do passo 
#cex: tamanho dos caracteres graficos 
#map=TRUE: opção para o semivariograma de superficie 
#main: titulo do grafico 
 
#Variavel Niquel 
coordinates(Geo) = ~ x + y; #definição das coordenadas 
trellis.par.set(sp.theme()) 
vgni = gstat(NULL, id="Niquel", formula= I(Ni)~1, data=Geo) 
plot(variogram(vgni, cutoff=2, width=.2, map=TRUE),cex=1.5, 
main = "(cross) semivariance maps: Niquel") 
# Sem indícios de anisotropia!!!!! 
#Variavel Cobaltotrellis.par.set(sp.theme()) 
vgco = gstat(NULL, "Cobalto", I(Co)~1, Geo) 
plot(variogram(vgco, cutoff=2, width=.2, map=TRUE),cex=1.5, 
main = "(cross) semivariance maps: Cobalto") 
 
 
 
### OBS. Indicios de anisotropia na variavel Cobalto !!! 
 
#### 2. Modelo de Regionalização ##### 
######################################## 
# Analise semivariografica: Niquel 
variogNi <- variog(Ni, uvec = seq(0, 2, length = 20), option = "bin") 
variogNi 
#Semivariograma unidirecional experimental 
plot(variogNi, xlab = "Distância", ylab = 
"Semivariância",type="l") 
title("Semivariograma Experimental: Niquel") 
#Semivariograma unidirecional e modelo ajustado 
plot(variogNi, xlab = "Distância", ylab = "Semivariância") 
lines.variomodel(cov.model = "spherical" , cov.pars = 
rbind(c(0,0.2),c(71,1.3)), nugget = 11, max.dist = 2, lwd 
= 2, col = "red") 
text(1,0.1, "MODELO: 11 + Sph(h/0,2) + 71Sph (h/1,3)") 
title("Semivariograma Experimental e modelo ajustado:Niquel") 
 
 
 
# Analise semivariografica: Cobalto : anisotropia !!!! 
#Semivariogramas experimentais nas duas direções 
coordinates(Geo) = ~ x + y 
vgmCo <- variogram(Co~1, Geo,coords.col = 1:2, 
alpha=c(67.5,157.5)) 
plot(vgmCo,type="b", col="black",main="Cobalto") 
# Modelo anisotropico ajustado 
model.Co <- vgm(11.5,"Sph",1.5,2.0,anis=c(67.5,0.67)) 
plot(vgmCo,model=model.Co,type="p", 
lty=1,col=c("black"),main="Cobalto") 
text(0.5,"MODELO:2+11.5Sph(h/1.5)+11.5 Sph(h/1,0) dir 67.5") 
#### 3. Modelo de Co-regionalização ##### 
########################################### 
 
#Cálculo da correlação entre as variaveis Niquel e Cobalto 
Cor <- cor.test(Ni$data,Co$data) #correlações 
Cor 
0.7507369 
## Modelos variograficos ajustados pelo gstat 
# Niquel 
Ni <- gstat(id = "Ni", formula = Ni ~ 1, data = Geo) 
vNi<- variogram(Ni, cutoff = 2, width = 2/20)#experimental Ni 
plot(vNi,type="o",main = "Variograma experimental:Niquel") 
mNi <-vgm(0, "Sph", 0.2, add.to = 
vgm(71, "Sph", 1.3, nugget = 11)) #Modelo ajustado 1 
plot(vNi, mNi,main = "Modelo ajustado: Niquel") 
 
#modelo ajustado pelo gstat 
fNi <- fit.variogram(vNi, mNi, fit.ranges = FALSE) 
plot(vNi, fNi,main = "Modelo ajustado gstat: Niquel") 
text(1,0.1, "MODELO: 11 + 9,0 Sph(h/0,12) + 64 Sph (h/1,4)") 
fNi 
model psill range 
1 Nug 11.20039 0.0 
2 Sph 70.03977 1.3 
3 Sph 0.00000 0.2 
# Cobalto 
Co <- gstat(id = "Co", formula = Co ~ 1, data = Geo) 
vCo <- variogram(Co, cutoff = 2, width = 
2/20,alpha=c(67.5,157.5)) # experimental Co 
plot(vCo,type="o",main = "Variograma experimental:Cobalto") 
 
mCo <- vgm(11.5,"Sph",1.5,2.0,anis=c(67.5,0.67)) #Modelo 
ajustado 1 
plot(vCo, mCo,main = "Modelo ajustado: Cobalto") 
 
#modelo ajustado pelo gstat 
fCo <- fit.variogram(vCo, mCo, fit.ranges = FALSE) 
plot(vCo, fCo,main = "Modelo ajustado gstat: Cobalto") 
fCo 
model psill range ang1 anis1 
1 Nug 1.14819 0.0 0.0 1.00 
2 Sph 12.13941 1.5 67.5 0.67 
 
 
 
#Variogramas cruzados 
gCoNi <- gstat(id = "Co", formula = Co ~ 1, data=Geo) 
gCoNi <- gstat(gCoNi, "Ni", Ni ~ 1, Geo) 
vCoNi <- variogram(gCoNi, cutoff = 2, width = 2/20) 
plot(vCoNi) #variogramas simples e cruzado 
 
 
 
# variograma Cruzado CoNi 
vCoNi2 <- vCoNi[1:20,] #Co-variograma experimental CoNi 
plot(vCoNi2,type="o",main = "Co-variograma 
experimental:CoNi") 
#modelo ajustado pelo gstat 
fCoNi <- fit.variogram(vCoNi2, 
model = vgm(0, "Sph", 0.2, add.to = 
vgm(28, "Sph", 1.3, nugget = 2)), 
fit.ranges = FALSE) #modelo ajustado CoNi 
plot(vCoNi2, fCoNi,col="red",lw=2,main = "Co-variograma 
ajustado:CdNi") 
fCoNi 
 
 
model psill range 
1 Nug 2.735748 0.0 
2 Sph 26.649184 1.3 
3 Sph -3.395847 0.2 
 
#Modelo Linear de Coregionalização 
gCoNi <- gstat(gCoNi, id = "Co", model = fCo) 
gCoNi <- gstat(gCoNi, id = "Ni", model = fNi) 
gCoNi <- gstat(gCoNi, id = c("Co","Ni"), model = fCoNi) 
 
CoNi.fit = fit.lmc(vCoNi, gCoNi, 
model = vgm(0, "Sph", 0.2, add.to = 
vgm(28,"Sph", 1.3, nugget = 2))) 
CoNi.fit 
data: 
Co : formula = Co`~`1 ; data dim = 259 x 10 
Ni : formula = Ni`~`1 ; data dim = 259 x 10 
variograms: 
model psill range 
Co[1] Nug 1.094104 0.0 
Co[2] Sph 12.738695 1.3 
Co[3] Sph 1.958586 0.2 
Ni[1] Nug 11.200386 0.0 
Ni[2] Sph 70.039773 1.3 
Ni[3] Sph 1.693443 0.2 
Co.Ni[1] Nug 2.735748 0.0 
Co.Ni[2] Sph 26.649184 1.3 
Co.Ni[3] Sph -1.821196 0.2 
 
 
Modelo de Co-regionalização: 
 
  Co Co (h)  Co Ni (h)  
 
 1.094 2.735 g  1.958 (h) 
1.821  h  
 
 12.74 
 
 
26.65   h 

 
 (h)  (h) 


 
2.735 11.20 
 0
 
 
1.821 1.693 Sph 
0.2 

 
26.65 70.04 
Sph
1.3 


 Ni Co Ni Ni           



Exercício 2: Escolher uma variável quantitativa (Cu, Pb, Co, Cr e Zn ) do banco 
jura.txt e realizar a análise variográfica, seguindo os mesmos procedimentos do 
Ensaio _2 
Unidade III: Métodos de estimação. krigagem simples, ordinária e 
indicativa. 
1. Métodos de estimação 
 
Para estimativas de valores em localizações não amostradas, existem diversos 
algoritmos de interpolação, tais como: krigagem, inverso do quadrado da distância, 
curvatura mínima, funções de base radial, triangulação, etc. Todos esses métodos 
estimam o valor em um determinado local, como uma soma ponderada de valores 
de dados em locais adjacentes. Quase todos atribuem pesos de acordo com as 
funções que dão um peso diminuindo com o aumento da distância de separação. A 
krigagem atribui pesos de acordo com os valores de covariância espacial. Em 
particular: 
 Se os dados são bastante densos e uniformemente distribuídas ao longo da 
área de estudo, em geral se tem boas estimativas, independentemente do 
algoritmo de interpolação 
 Se os dados estão agrupados (em clusters) e com grandes lacunas entre si, 
pode-se obter estimativas não confiáveis, independentemente do algoritmo de 
interpolação. 
 Quase todos os algoritmos de interpolação subestimam os altos valores e 
superestimam os baixos; isto é inerente ao cálculo da média e se um algoritmo de 
interpolação não considera a média não seria considerado razoável 
 
Algumas vantagens da krigagem: 
 Ajuda para compensar os efeitos do agrupamento de dados, atribuindo a 
pontos individuais dentro de um cluster menos peso do que pontos isolados de 
dados (ou, tratando aglomerados mais como pontos únicos) 
 Dá estimativa do erro de estimação (variância de krigagem), juntamente com a 
estimativa da variável, Z. 
 A disponibilidade de erro de estimação fornece base para a simulação 
estocástica de possíveis realizações de Z(u). 
 Considera a anisotropia do meio. 
2. Krigagem 
2.1. Krigagem Simples 
To dos os estimadores de krigagem são variantes do estimador de 
regressão linear básico Z*(u) definida como (GOOVAERTS, 1997): 
 3.1 
onde λα(u) é o peso atribuído a cada valor observado Z(uα) localizado dentro de 
determinada vizinhança W(u) centrado em u. Os n(u) pesos são escolhidos de 
forma a minimizar a estimação ou erro da variância sob a condição de não- 
tendenciosidade do estimador. Os pesos de krigagem, λα(u), são derivados da 
função de covariância ou semivariograma As diferenças entre os vários tipos de 
krigagem residem no modelo considerado para a tendência m(u) na expressão 
(3.1). 
A krigagem simples (SK) considera a média m(u) conhecida e constante 
por toda a área de estudo, enquanto a krigagem ordinária (KO) limita o domínio de 
estacionaridade da média à vizinhança local W(u) e, ao contrário da KS, a média é 
desconhecida. 
2.2. Krigagem Ordinária (OK) 
 
A KO considera a variação local da média limitada ao domínio de 
estacionaridade da média a uma vizinhança local W(u) centrada sobre a localização 
u a ser estimada. O estimador da krigagem ordinária pode ser escrito como 
 
3.2 
 
 
Neste caso, a média local m(u’) é constante, mas desconhecido. 
A minimização da variância do erro, sujeita à restrição da soma dos pesos 
igual a 1, permitea determinação dos pesos λα(u) a partir do seguinte sistema de 
equações chamado de sistema de krigagem ordinária (equações normais com 
restrições): 
 
3.3 
 
 
onde C(uβ – uα) e C(u – uα) são, respectivamente, a covariância entre os pontos uβ 
e uα e u e uα e μ(u) é o parâmetro de Lagrange associado com a restrição 
 
. 
Substituindo a covariância por sua expressão C(h) = Co – γ(h), o sistema (3.3) e a 
variância podem ser escritos em função do modelo semivariográfico γ(h). 
2.3. Krigagem Indicativa (IK) 
A krigagem indicativa é utilizada para o mapeamento de variáveis nas 
seguintes situações: 
 Mapeamento de risco de ocorrência de uma variável contínua acima de 
determinado ponto de corte (cutoff) 
 Mapeamento de variáveis categóricas. Ex. uso da terra, tipos de solos,... 
 Predição da variação espacial da probabilidade da variável estar dentro de um 
determinado critério (mapas probabilisticos). 
São utilizadas variáveis indicativas para definição dos valores limites 
(threshold), definidas como: 
 
 
quando z(ui) é uma v.a. categórica, onde zcat é o código ou índice da categoria de 
interesse, e 
para v.a. quantitativas, onde zcut é uma valor de corte para v.a. quantitativa Z(u). 
No modelamento da variabilidade espacial são utilizados variogramas 
indicativos definidos como, 
 
 
Os variogramas indicativos medem a proporção de pares de pontos 
separados pelo vetor h, em que ocorreu o evento de interesse, no caso z(u) = zcat 
ou z(u) ≤ zcut. 
Teoricamente, o semivariograma é zero se a distância entre o pares de 
pontos for zero, mas na prática, geralmente ocorre uma descontinuidade na 
origem do semivariograma, chamado de Efeito pepita. Ou seja, a medida que h 
tende a zero o semivariograma assumirá um valor positivo Co. Isto se deve a erros 
amostrais ou a presença de variabilidade em curtas escalas, onde amostras 
separadas em pequenas distância pode apresentar diferenças significantes. (ISAAKS 
& SRIVASTAVA, 1990, p. 143). 
No caso de variáveis quantitativas o semivariograma indicativo mede o 
quanto dois valores z separados por um vetor h estão em lados opostos do valor de 
corte zk. Em outras palavras, mede a frequencia de transição entre duas classe de 
valores z em função de h. Quanto maior o semivariograma indicativo, menos 
conectados no espaço são os pequenos ou grandes valores. 
Para variáveis categóricas o semivariograma indicativo mede a freqüência 
com que duas localizações separadas por um vetor h pertencem a diferentes 
categorias. Quanto menor o valor do semivariograma, melhor a conectividade 
espacial da categoria. 
Aproximação indicativa no modelamento da incerteza local 
 
• A função é modelada para um conjunto de K valores limites zk, que 
discretizam a faixa de variação de z 
 
 
 
 A resolução da fdac é então aumentada pela interpolação dentro de cada classe 
(zk;zk+1] e extrapolação dos dois valores extremos, z1 e zk 
 A estimação geoestatística não-paramétrica dos valores da fdac é baseada na 
 
interpretação da probabilidade condicional como a esperança 
condicional de uma variável aleatória indicativa: 
 
onde 
 
 Os valores da fdac podem ser obtidos por krigagem ordinária de 
 
PROCEDIMENTO (O programa R realiza essas tarefas automaticamente  ) 
1. Fazer a codificação de cada observação em um vetor de K valores indicativos. O 
numero de classes deve ser escolhido de modo que as K+1 classes tenham 
aproximadamente frequencias iguais. Recomenda-se usar mais de 5 classes e 
menos de 15. 
α 
∑λ α 
α 
2. Calcular os semivariogramas experimentais para cada ponto de corte zk e, em 
seguida, fazer o modelamento dos mesmos. 
3. Em cada localização não amostrada u, 
– Estimar os K valores da fdac usando, por exemplo, a krigagem ordinária 
indicativa, i.e, 
– Interpolar ou extrapolar os valores da fdac para construir um modelo continuo 
para a fdac, que permita recuperar a probabilidade de qualquer ponto z. 
 
3. Krigagem multivariada: co-krigagem 
• A krigagem, é um caso particular do método co-krigagem 
• Requisitos : 
– dependência espacial para cada uma das variáveis Z1 e Z2, 
– dependência espacial entre Z1 e Z2 
No caso de uma simples variável secundária (Y), o estimador de cokrigagem 
ordinária é escrito como: 
 
* 
COK 
n1 
(u) = ∑λ 
1 
α1 =1 
 
(u)Z (u 
1 
n2 
) + ∑λ' 
2 
α2 =1 
 
(u) Y (u' ) 
2 
 
n 
n2 
condições de não-tendenciosidade 
1 
(u) = 1 
1 
∑λ' = 1 
2 
α2 =1 
α=1 
 
 
Sistema de cokrigagem para determinação dos n pesos 
 
λα1 e λ 2 
 
Composto de (n1 + n2 + 2) equações com 2 parâmetros de Lagrange μ1 e μ2 
associados às duas condições de não-tendenciosidade: 
∑
n1  (u)C (u - u )  ∑
n2 
 C (u - u )    C (u - u ) para   1,..., n 
  1  Z    '1  ' ZY  '  1 Z 1 
∑n1  (u)C (u - u )  ∑
n2
  C (u - u )    C (u - u ) para  '  1,..., n 
  1 
 n 
ZY   '  '1  ' Y  '  ' 2 ZY  ' 2 
∑1   1 
  1 

∑
n ' 
  1 
  '  ' 
Z α α α 
α 
 
 
 
 
 
Este sistema de cokrigagem apresenta uma e somente uma solução se: 
(i) não há dois dados totalmente redundantes, isto é: 
u  u se:    = 1,..., n. 
u’  u’ se: ’  ’ = 1,..., n’. 
 
CZ (h) CZY (h) 
(ii) a matriz de funções de covariância 
CZY (h) CY (h) 
é positiva-definida. 
 
 
Ensaio 3 (Script: Unidade 3): 
 krigagem ordinária isotrópica para o Niquel e anisotrópica para o Cobalto 
 Krigagem indicativa 
o mapas probabilistiscos e temáticos para o uso da terra 
o mapas de riscos para o cádmio (valor máximo permitido ou ponto de 
corte: 0.8 ppm) 
a) Variáveis selecionadas: Niquel, Cobalto e Cadmio (quantitativas) 
b) landuse (qualitativas) 
PROGRAMA R 
1. krigagem ordinária: 
1.1 isotrópica para o Niquel 
 
#Definir a área de trabalho 
setwd("D:/PPGCA/Hidrogeo") 
#dados de entrada podem ser no formato ASCII (.txt) 
#Carregar dados 
jura.pred <- read.table("jura.txt", head=T) 
jura.pred[1:5,] 
## lendo arquivo com bordas da área 
bor <- read.table("jurabor.txt", head=T) 
bor[1:5,] # mostrar 5 linhas do aqrquivo 
x y 
1 2.1497 4.3352 
2 2.1151 4.2197 
3 2.1151 4.0927 
4 2.1728 3.9888 
5 2.1843 3.8849 
 
 
####### Krigagem Ordinária usando o geoR ####### 
### Analise semivariografica: Niquel 
 
jura.pred<-read.table("jura.txt",head=T) 
Ni<-as.geodata(jura.pred,coords.col = 1:2, data.col = 9) 
Ni #mostrar a variavel Niquel area de trabalho 
varNi <- variog(Ni, uvec = seq(0, 2, length = 20), 
option = "bin") 
varNi 
#Semivariograma unidirecional experimental 
#par(mfrow = c(1,2)) # set up the graphics 
plot(varNi, xlab = "Distância", ylab = "Semivariância", 
type="l") 
title("Semivariograma Experimental: Niquel") 
#Semivariograma unidirecional e modelo ajustado 
plot(varNi, xlab = "Distância", ylab = "Semivariância") 
lines.variomodel(cov.model = "spherical" , cov.pars = 
rbind(c(9,0.12),c(64,1.4)), nugget = 11, max.dist = 
2, lwd = 2, col = "red") 
text(1,0.1,"MODELO:11+9,0 Sph(h/0,12) + 64 Sph (h/1,4)") 
title("Semivariograma Experimental e modelo ajustado: 
Niquel") 
 
 
 
###Mapa de Krigagem do Niquel 
locNi <- expand.grid(seq(0,5.5,l=100), seq(0,6,l=120)) 
kcNi<-ksline(Ni,locations=locNi,cov.model= "spherical", 
cov.pars=rbind(c(9,0.12),c(64,1.4)),nugget = 11) 
contour(kcNi,filled=TRUE,bor=bor,col=terrain.colors(20), 
nlevels = 13) 
title("Krigagem ordinaria Niquel") 
 
# definindo uma palleta de cores 
require(maptools) 
bluepal<-colorRampPalette(c=terrain.colors(20)) 
brks <- c(0,5,10,15,20,25,30,35,40,45) 
cols <- bluepal(length(brks) - 1) 
names(kcNi) 
image(kcNi, border=bor, loc=locNi, val=kcNi$predict, col 
= terrain.colors(20)) 
legend(-0.5,6, fill = cols, legend = leglabs(brks,"<", 
">="),bty = "n", cex = 0.8,col=terrain.colors(20)) 
title("Krigagem Ordinaria Ni") 
 
 
1.2 anisotrópica para o Cobalto 
 
###### 1.2 Krigagem ordinária com anisotropia ##### 
Co <-as.geodata(jura.pred,coords.col=1:2,data.col = 6) 
#Modelo semivariografico do cobalto (Unidade 2): 
#(h)= 2 + 11.5 Sph[(h/1.5)2 + 11.5 Sph(h/1.0)2 dir 67.5 
###Mapa de Krigagem do Cobalto 
locCo <- expand.grid(seq(0,5.5,l=100), seq(0,6,l=120)) 
kcCo <- ksline(Co, locations=locCo, cov.model= " 
spherical", cov.pars = rbind(c(11.5,1.5)),nugget = 2.0, 
aniso.pars = c(67.5,1.5)) 
contour(kcCo, filled=TRUE,borders=bor,col=terrain. 
colors(30)) 
title("Krig.Ord.:Cobalto(anisotropia") 
 
 
 
 
2. Krigagem indicativa 
3.1 Mapas probabilistiscos e temáticos para o uso da terra 
###### 3. Krigagem Indicativa ######### 
###### 3.1 Mapas de proabilidades 
###### para gerar os mapas temáticos ######### 
#carregar arquivo de dados categoricos (uso da terra) 
Landuse <- read.table("land.txt", head=T) 
Landuse[1:5,] # mostrar os dados na area de trabalho 
X Y Landuse Flore Pasto Prado Lavoura 
1 2.699 1.199 1 1 0 0 0 
2 4.232 1.588 1 1 0 0 0 
3 4.459 2.624 1 1 0 0 0 
4 2.916 4.916 1 1 0 0 0 
5 3.649 2.674 1 1 0 0 0 
 
summary(Landuse) 
###Semivariograma indicativo experimental:Floresta 
Flore <-as.geodata(Landuse,coords.col=1:2,data.col = 4) 
 
par(mfrow = c(1,2)) # set up the graphics 
v.flor <- variog(Flore, uvec = seq(0, 2, length = 20), 
option = "bin") 
plot(v.flor, xlab = "Distância", ylab = "Semivariância", 
main = "Semivariograma Floresta", 
type = "b", col = "black") 
 
#Semivariograma unidirecional e modelo ajustado 
plot(v.flor, xlab = "Distância", ylab = "Semivariância") 
lines.variomodel(cov.model = "spherical" , cov.pars = 
rbind(c(0.04,0.5),c(0.08,1.4)), nugget = 0.01, max.dist 
= 2, lwd = 2, col = "red") 
text(1,0.01, "MODELO: 0.01 + 0.04 Sph(h/0,5) + 0.08 Sph 
(h/1,4)") 
title("Semivariograma Experimental e modelo ajustado") 
# Ajuste usando mínimos quadrados ponderados 
par(mfrow = c(1,1)) 
locin <- matrix(c(0.04,0.08,0.5,1.4), ncol=2) 
flor.ols <- variofit(v.flor, ini=locin, nugget=0.01, 
min="optim") 
plot(v.flor) 
lines(flor.ols, lty=1, col="red") 
title("Semivariograma ajustado Floresta") 
 
 
 
 
###Semivariograma experimental da variavel Pasto 
Pasto <-as.geodata(Landuse,coords.col =1:2,data.col = 5) 
par(mfrow = c(1,2)) # set up the graphics 
v.pas <- variog(Pasto, uvec = seq(0, 2, length = 20), 
option = "bin") 
plot(v.pas, xlab = "Distância", ylab = "Semivariância", 
main="Semivariograma Pasto", type = "b", col = "black") 
 
#Semivariograma unidirecional e modelo ajustado 
plot(v.pas, xlab = "Distância", ylab = "Semivariância") 
lines.variomodel(cov.model = "spherical" , cov.pars = 
rbind(c(0.17,0.4)), nugget = 0.01, max.dist = 2, lwd = 
2, col = "red") 
text(1,0.01, "MODELO: 0.01 + 0.17 Sph(h/0,4)") 
title("Semivariograma Experimental e modelo ajustado") 
 
# Ajuste usando mínimos quadrados ponderados 
par(mfrow = c(1,1)) # set up the graphics 
pas.ols <- variofit(v.pas, ini=c(0.17,0.4), nugget=0.01, 
min="optim") 
plot(v.pas) 
lines(pas.ols, lty=1, col="red") 
title("Semivariograma ajustado Pasto") 
 
 
 
 
###Semivariograma experimental da variavel Prado 
Prado <-as.geodata(Landuse,coords.col=1:2, data.col = 6) 
par(mfrow = c(1,2)) 
v.pra <- variog(Prado, uvec = seq(0, 2, length = 20), 
option = "bin") 
plot(v.pra, xlab = "Distância", ylab = "Semivariância", 
main = "Semivariograma Prado", type = "b", col = 
"black") 
#Semivariograma unidirecional e modelo ajustado 
plot(v.pra, xlab = "Distância", ylab = "Semivariância") 
lines.variomodel(cov.model = "spherical" , cov.pars = 
rbind(c(0.23,0.4)), nugget = 0.01, max.dist = 2, lwd = 
2, col = "red") 
text(1,0.01, "MODELO: 0.01 + 0.23 Sph(h/0,4)") 
title("Semivariograma Experimental e modelo ajustado") 
## Ajuste usando mínimos quadrados ponderados 
par(mfrow = c(1,1)) # set up the graphics 
pra.ols <- variofit(v.pra, ini=c(0.23,0.4), 
nugget=0.02,min="optim") 
plot(v.pra) 
lines(pra.ols, lty=1, col="red") 
title("Semivariograma ajustado Prado") 
 
 
 
 
###Semivariograma experimental da variavel Lavoura 
Lavoura<-as.geodata(Landuse,coords.col=1:2,data.col = 7) 
par(mfrow = c(1,2)) 
v.lav <- variog(Lavoura, uvec = seq(0, 2, length = 20), 
option = "bin") 
plot(v.lav, xlab = "Distância", ylab = "Semivariância", 
main = "Semivariograma Lavoura", type = "b", col = 
"black") 
#Semivariograma unidirecional e modelo ajustado 
plot(v.lav, xlab = "Distância", ylab = "Semivariância") 
lines.variomodel(cov.model = "spherical" , cov.pars = 
rbind(c(0.020,0.4)), nugget = 0.001, max.dist = 2, lwd = 
2, col = "red") 
text(1,0.001, "MODELO: 0.001 + 0.020 Sph(h/0,4)") 
title("Semivariograma Experimental e modelo ajustado") 
 
# Ajuste usando mínimos quadrados ponderados 
par(mfrow = c(1,1)) # set up the graphics 
lav.ols <- variofit(v.lav, ini=c(0.020,0.4), 
nugget=0.00,min="optim") 
plot(v.lav) 
lines(lav.ols, lty=1, col="red") 
title("Semivariograma ajustado Lavoura") 
 
 
 
 
 
##Krigagem indicativa para os mapas probabilisticos 
#Definir grid (malha) 
x <- seq(0,5.5,l=100) 
y <- seq(0, 6, l=120) 
locn <- expand.grid(x=x,y=y) 
 
kc.flor <- krige.conv(Flore, loc=locn, 
krige=krige.control(obj=flor.ols)) 
summary(kc.flor) 
 
kc.pas <- krige.conv(Pasto, loc=locn, 
krige=krige.control(obj=pas.ols)) 
summary(kc.pas) 
 
kc.pra <- krige.conv(Prado, loc=locn, 
krige=krige.control(obj=pra.ols)) 
summary(kc.pra) 
 
kc.lav <- krige.conv(Lavoura, loc=locn, 
krige=krige.control(obj=lav.ols)) 
summary(kc.lav) 
# Mapas probabilisticos para cada categoria 
par(mfrow = c(2,2)) 
flor<-kc.flor$predict 
image(kc.flor, border=bor, loc=locn, 
val=flor,col=rainbow(9)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(2.5,5.5),flor,col=rainbow(9), vert=TRUE) 
title("Floresta") 
 
pas<-kc.pas$predict 
image(kc.pas, border=bor, loc=locn, 
val=pas,col=rainbow(9)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(2.5,5.5),pas,col=rainbow(9), vert=TRUE) 
title("Pasto") 
 
pra<-kc.pra$predict 
image(kc.pra, border=bor, loc=locn, 
val=pra,col=rainbow(9)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(2.5,5.5),pra,col=rainbow(9), vert=TRUE) 
title("Prado") 
 
lavo<-kc.lav$predict 
image(kc.lav, 
loc=locn,border=bor,val=lavo,col=rainbow(9)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(2.5,5.5),val=lavo,col=rainbow(8), vert=TRUE) 
title("Lavoura") 
 
 
 
 
### Construçao de mapas temáticos com uso de Mapas de 
### probabilidade 
#Construção das tabelas com valores preditos 
prob <- cbind(kc.flor$predict,kc.pas$predict, 
kc.pra$predict,kc.lav$predict) 
colnames(prob) <- c("Floresta", "Pasto", "Prado", 
"Lavoura") #Renomeando as colunas de cada categoria 
max <- max.col(prob) #cat.com valor máximo de probabilidade 
mapa <- cbind(locn,max) #Tabela com localização e valor máximo 
colnames(mapa) <- c("X","Y","Prob") #Renomeando as colunas 
summary(mapa) 
 
class(mapa) <- "kriging" # classe das variaveis 
mapa$predict <- mapa$Prob # redefinindo o valor "Prob" 
mapa$coords <- cbind(mapa$X,mapa$Y) #coordenadas 
 
#Mapa de probabilidades por classes 
par(mfrow = c(1,1)) 
par(cex = 1.1) #definfindo o tamanho dos caracteres 
image(mapa, locations = locn, borders=bor, xlab="W-E", 
ylab="S-N", col = rainbow(4)) 
legend(-0.5, 5.5,c("Floresta", "Pasto", 
"Prado","Lavoura"), fill = rainbow(4)) 
 
 
 
 
 
 
3.2 Mapa de riscos para o Cádmio (valor máximo permitido (Zc : 0.8 ppm)) 
###### 3.2 Mapas de riscos 
###### Variável selecionada: Cadmio 
###### Valor máximo permitido (Zc = 0.8 ppm) 
### Selecionar a variavel Cadmio 
jura.pred <- read.table("jura.txt", head=T) 
Cd<-as.geodata(jura.pred,coords.col = 1:2, data.col = 5) 
# Analise variografica 
par(mfrow = c(1,2)) # set up the graphics 
variogCd <- variog(Cd, uvec = seq(0, 2, length = 20), 
option = "bin") 
plot(variogCd,type="l",main="Semivariograma 
experimental") 
#Semivariograma unidirecional e modelo ajustado (geoR) 
plot(variogCd, xlab = "Distância", ylab = 
"Semivariância",main="Modelo ajustado") 
lines.variomodel(cov.model = "spherical" , cov.pars= 
rbind(c(0.35,0.2),c(0.26,1.3)), nugget = 0.25, 
max.dist = 2, lwd = 2, col = "red") 
text(1,0.1, "MODELO: 0,25 + 0,35 Sph(h/0,2) + 0,26 Sph 
(h/1,3)") 
## Ajuste usando mínimos quadrados ponderados 
par(mfrow = c(1,1)) # set up the graphics 
locin <- matrix(c(0.35,0.26,0.2,1.3), ncol=2) 
Cd.ols <- variofit(variogCd, ini=locin, nugget=0.25, 
min="optim") 
Cd.ols 
plot(variogCd) 
lines(Cd.ols, lty=1, col="blue", lwd=3) 
 
 
 
 
 
##Execução com um grid (50x60:menos tempo de computação 
##Krigagem indicativa 
#Definir grid (malha) 
x <- seq(0,5.5,l=50) 
y <- seq(0, 6, l=60) 
locCdi <- expand.grid(x=x,y=y) 
summary(locCdi) 
## Krigagem Indicativa: predição de probabilidades e 
quantis 
## Probabilidade de Cd ser maior que 0.8 ppm em cada um 
dos pontos 
OC <- output.control(thres=0.8, 
quan=c(0.34,0.50,0.60,0.80,0.90)) 
Cd.kc <- krige.conv(Cd, loc=locCdi, 
krige=krige.control(obj=Cd.ols),output=OC) 
 
# examinando o objeto retornado pela função 
names(Cd.kc) 
Zp<-Cd.kc$pred # valores preditos 
Zv<-Cd.kc$krige.var # variancia 
Zprob=Cd.kc$prob # probabilidades estimados do valor ser 
menor que 0.8 nos pontos selecionados 
Zc=1-Cd.kc$prob # probabilidades estimados do valor ser 
maior que 0.8 nos pontos selecionados 
Ztype=Cd.kc$mean # E-type (médias) 
Zq=Cd.kc$quant # quantis 
 
#Salvando os resultados em uma tabela 
ki <- cbind(Zp,Zv,Zc,Zprob,Ztype,Cd.kc$quant) 
 
kind<-cbind(locCdi,ki) 
names(kind) 
write.table(kind,"tkind0.dat") 
 
###Mapas 
#Mapa de probabilidades Zc>0.8 
image(Cd.kc, border=jbor, loc=locCdi, 
val=Zc,col=terrain.colors(30)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(3,6),Zc,col=terrain.colors(30), vert=TRUE) 
title("Cadmio: Prob Zc>0.8") 
 
 
 
 
 
# Mapa dos valores médios (E-type) 
image(Cd.kc, border=bor, loc=locCdi, 
val=Cd.kc$Ztype,col=terrain.colors(30)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(3,6),Cd.kc$Ztype,col=terrain.colors(30), 
vert=TRUE) 
title("Cadmio: E-type") 
 
 
 
 
## Mapeando quantis 
names(Cd.kc$quant) 
# Mapa dos valores medianos (quantil-0.5) 
Zmed=Cd.kc$q50 
image(Cd.kc, border=bor, loc=locCdi, 
val=Zmed,col=terrain.colors(30)) 
legend.krige(x.leg=c(-0.4,-0.05), 
y.leg=c(3,6),Zmed,col=terrain.colors(30), vert=TRUE) 
title("Cadmio: Mediana") 
 
 
 
 
 
Exercício 3: Escolher uma ou duas variáveis quantitativas (Cu, Pb, Co, Cr e Zn ) do 
banco jura.txt 
a) Construir o mapa da variavel com uso da krigagem ordinaria 
b) Com uso da krigagem indicativa, construir o mapa temático para a variável tipo e 
o mapa de risco para a variável escolhida em a). As tolerancias máximas permitidas 
se encontram na Tabela 1 (Unidade 1) 
Use os mesmos procedimentos do Ensaio _3 
	Uso do programa R (Módulos GeoR, gstat)
	 Símbolos ou comandos importantes
	 Procedimento básico para instalar um pacote
	1. Inferência estatística
	2. Análise variográfica
	2.1- Variografia Experimental
	2.2- Modelamento variográfico
	2.2.1 Modelagem dos Semivariogramas e da Anisotropia
	3. Modelo Linear de Regionalização
	4. Modelo linear de co-regionalização
	4.1. Semivariogramas cruzados
	Características
	4.2. Modelamento da co-regionalização
	Ensaio 2 (Script: Ensaio_2):
	PROGRAMA R
	1. Métodos de estimação
	Algumas vantagens da krigagem:
	2. Krigagem
	2.2. Krigagem Ordinária (OK)
	3.2
	2.3. Krigagem Indicativa (IK)
	Aproximação indicativa no modelamento da incerteza local
	3. Krigagem multivariada: co-krigagem
	Ensaio 3 (Script: Unidade 3):
	PROGRAMA R (1)
	class(mapa) <- "kriging" # classe das variaveis