Este blog é resultado de um trabalho solicitado pelo professor Luis Paulo Vieira Braga em 2014 durante o curso de Introdução à Geoestatística da Faculdade de Geologia da Universidade Federal do Rio de Janeiro (UFRJ).
O trabalho tem por objetivo fazer uma análise estatística dos dados de concentração de Cobalto no dataset Jura utilizando o software R e os conceitos básicos de estatística.
Os dados dos poços serão usados para estimar os parâmetros de um modelo esférico de variograma omnidirecional. Por fim, o modelo teórico será utilizado para gerar o produto final: Um mapa de concentração de Cobalto. Para isso usaremos o método de interpolação por Krigagem ordinária.
domingo, 14 de dezembro de 2014
O Dataset Jura
O dataset Jura consiste em um banco de dados de 359 poços perfurados na Suiça onde foram coletados dados de tipo de rocha, tipo de uso do solo e concentração de metais pesados por poço.
Esse conjunto de dados é de domínio público, estando incluído no pacote GSTAT.
Esse pacote deve ser instalado no Software R para se ter acesso aos dados. O dataset é carregado digitando-se no terminal do R:
data(jura)
Os dados estão estruturados em um dataframe. Nele constam uma amostra de treinamento (prediction.dat) e uma amostra de validação (validation.dat). A primeira, que vamos utilizar neste trabalho, consiste em uma tabela com os seguintes atributos: coordenada Xloc(km), coordenada Yloc(km), tipo de rocha (1 para Argoviano, 2 Kimmeridgiano, 3 Sequanianao, 4 Portlandiano e 5 Quartenário), uso do solo (1 Floresta, 2 Pasto, 3 Prado e 4 Lavoura) e teores de alguns metais pesados em ppm. A Tabela 1 mostra parte dos dados representativos da amostragem prediction.dat.
Esse conjunto de dados é de domínio público, estando incluído no pacote GSTAT.
Esse pacote deve ser instalado no Software R para se ter acesso aos dados. O dataset é carregado digitando-se no terminal do R:
data(jura)
Os dados estão estruturados em um dataframe. Nele constam uma amostra de treinamento (prediction.dat) e uma amostra de validação (validation.dat). A primeira, que vamos utilizar neste trabalho, consiste em uma tabela com os seguintes atributos: coordenada Xloc(km), coordenada Yloc(km), tipo de rocha (1 para Argoviano, 2 Kimmeridgiano, 3 Sequanianao, 4 Portlandiano e 5 Quartenário), uso do solo (1 Floresta, 2 Pasto, 3 Prado e 4 Lavoura) e teores de alguns metais pesados em ppm. A Tabela 1 mostra parte dos dados representativos da amostragem prediction.dat.
| Xloc | Yloc | Landuse | Rock | Cd | Co | Cr | Cu | Ni | Pb | Zn |
| 2.386 | 3.077 | 3 | 3 | 1.740 | 9.320 | 38.32 | 25.720 | 21.32 | 77.36 | 92.56 |
| 2.544 | 1.972 | 2 | 2 | 1.335 | 10.000 | 40.20 | 24.760 | 29.72 | 77.88 | 73.56 |
| 2.807 | 3.347 | 2 | 3 | 1.610 | 10.600 | 47.00 | 8.880 | 21.40 | 30.80 | 64.80 |
| 4.308 | 1.933 | 3 | 2 | 2.150 | 11.920 | 43.52 | 22.700 | 29.72 | 56.40 | 90.00 |
Tabela 1
Desenvolvimento
1- Mapa Base
Primeiro vamos gerar o mapa base do conjunto de poços vistos em planta. Para isso vamos utilizar os comandos abaixo no terminal do R:
data(jura)
plot(prediction.dat[,1],prediction.dat[,2],xlab="X(Km)",ylab="Y(Km)",main="Mapa base dos poços")
2-Sumário com boxplot
Para gerar o sumário da variável Cobalto, vamos utilizar os comandos abaixo:
dados_de_cobalto <- prediction.dat$Co
summary(dados_de_cobalto)
O resultado obtido foi:
Valor mínimo: 1.552
1º quartil: 6.520
Mediana: 9.760
Média: 9.303
3º quartil: 11.980
Valor máximo: 17.720
Os valores podem ser conferidos graficamente ao lado no gráfico 2.
Comandos para gerar o boxplot (gráfico 2):
data(jura)
dados_de_cobalto <- prediction.dat$Co
boxplot(dados_de_cobalto, range=1.5, ylab="Cobalto(ppm)",main="Boxplot Cobalto")
3-Histograma
Para gerar um histograma dos dados e fazer um controle de qualidade, vamos comparar o histograma com a curva de distribuição normal de mesma média e desvio padrão. Para isso, vamos corrigir os dados de frequência relativa do histograma para a soma das áreas das barras ser igual a 1. Foram usados os comandos abaixo:
data(jura)
dados_de_cobalto <- prediction.dat$Co
h<-hist(dados_de_cobalto, breaks=15)
xhist<-c(min(h$breaks),h$breaks)
yhist<-c(0,h$density,0)
xfit<-seq(1, 18, by=0.01)
yfit<-dnorm(xfit,mean=mean(dados_de_cobalto),sd=sd(dados_de_cobalto))
plot(xhist,yhist,type="s",ylim=c(0,max(yhist,yfit)),xlab="Cobalto (ppm)",ylab="Frequência relativa corrigida",main="Histograma de densidade de Cobalto")
lines(xfit,yfit,col="red")
4-Variograma Experimental
Utilizando o R vamos gerar um variograma omnidirecional experimental para a variável Co e em seguida fazer um ajuste ao modelo teórico esférico de variograma.
Os comandos para gerar o variograma são:
data(jura)
g <- gstat(id="Co", formula=Co~1, locations=~Xloc+Yloc, data=prediction.dat)
graf <- variogram(g)
plot(graf, main="Variograma omnidirecional experimental para Cobalto", xlab="Distância", ylab="Semivariância")
5-Estimando um modelo teórico de variograma
Com os dados do variograma experimental, vamos fazer um ajuste para obtenção dos parâmetros do modelo teórico de variograma. Os comandos utilizados para tal foram:
g <- gstat(id="Co", formula=Co~1, locations=~Xloc+Yloc, data=prediction.dat)
graf <- variogram(g)
f<-fit.variogram(graf,vgm(15,"Sph",2.14,0.5))
Os valores ajustados são guardados na variável f. Abaixo podemos ver o conteúdo desta variável.
model psill range
1 Nug 1.304965 0.000000
2 Sph 12.524985 1.183543
Com os valores ajustados, geramos o modelo teórico (linha azul) e comparamos com os pontos do variograma experimental. Para isso, usamos os comandos a seguir:
ff<-variogramLine(f,maxdist=2.14,n=500,min=1.0e-6)
plot(ff,col="blue", main="comparação (variograma experimental e modelo esférico)")
points(graf[,2],graf[,3], col="red")
6-Mapa de concentração de teores de Co
Finalmente, vamos usar os parâmetros do modelo teórico calculados na etapa anterior para fazer a Krigagem e gerar o mapa de concentração de teor de Cobalto. Os comandos podem ser vistos abaixo.
Primeiro vamos gerar o mapa base do conjunto de poços vistos em planta. Para isso vamos utilizar os comandos abaixo no terminal do R:
data(jura)
plot(prediction.dat[,1],prediction.dat[,2],xlab="X(Km)",ylab="Y(Km)",main="Mapa base dos poços")
 |
| Gráfico 1 |
2-Sumário com boxplot
 |
| Gráfico 2 |
Para gerar o sumário da variável Cobalto, vamos utilizar os comandos abaixo:
dados_de_cobalto <- prediction.dat$Co
summary(dados_de_cobalto)
O resultado obtido foi:
Valor mínimo: 1.552
1º quartil: 6.520
Mediana: 9.760
Média: 9.303
3º quartil: 11.980
Valor máximo: 17.720
Os valores podem ser conferidos graficamente ao lado no gráfico 2.
Comandos para gerar o boxplot (gráfico 2):
data(jura)
dados_de_cobalto <- prediction.dat$Co
boxplot(dados_de_cobalto, range=1.5, ylab="Cobalto(ppm)",main="Boxplot Cobalto")
3-Histograma
Para gerar um histograma dos dados e fazer um controle de qualidade, vamos comparar o histograma com a curva de distribuição normal de mesma média e desvio padrão. Para isso, vamos corrigir os dados de frequência relativa do histograma para a soma das áreas das barras ser igual a 1. Foram usados os comandos abaixo:
data(jura)
dados_de_cobalto <- prediction.dat$Co
h<-hist(dados_de_cobalto, breaks=15)
xhist<-c(min(h$breaks),h$breaks)
yhist<-c(0,h$density,0)
xfit<-seq(1, 18, by=0.01)
yfit<-dnorm(xfit,mean=mean(dados_de_cobalto),sd=sd(dados_de_cobalto))
plot(xhist,yhist,type="s",ylim=c(0,max(yhist,yfit)),xlab="Cobalto (ppm)",ylab="Frequência relativa corrigida",main="Histograma de densidade de Cobalto")
lines(xfit,yfit,col="red")
 |
| Gráfico 3 |
4-Variograma Experimental
Utilizando o R vamos gerar um variograma omnidirecional experimental para a variável Co e em seguida fazer um ajuste ao modelo teórico esférico de variograma.
Os comandos para gerar o variograma são:
data(jura)
g <- gstat(id="Co", formula=Co~1, locations=~Xloc+Yloc, data=prediction.dat)
graf <- variogram(g)
plot(graf, main="Variograma omnidirecional experimental para Cobalto", xlab="Distância", ylab="Semivariância")
 |
| Gráfico 4 |
5-Estimando um modelo teórico de variograma
Com os dados do variograma experimental, vamos fazer um ajuste para obtenção dos parâmetros do modelo teórico de variograma. Os comandos utilizados para tal foram:
g <- gstat(id="Co", formula=Co~1, locations=~Xloc+Yloc, data=prediction.dat)
graf <- variogram(g)
f<-fit.variogram(graf,vgm(15,"Sph",2.14,0.5))
Os valores ajustados são guardados na variável f. Abaixo podemos ver o conteúdo desta variável.
model psill range
1 Nug 1.304965 0.000000
2 Sph 12.524985 1.183543
Com os valores ajustados, geramos o modelo teórico (linha azul) e comparamos com os pontos do variograma experimental. Para isso, usamos os comandos a seguir:
ff<-variogramLine(f,maxdist=2.14,n=500,min=1.0e-6)
plot(ff,col="blue", main="comparação (variograma experimental e modelo esférico)")
points(graf[,2],graf[,3], col="red")
 |
| Gráfico 5 |
6-Mapa de concentração de teores de Co
Finalmente, vamos usar os parâmetros do modelo teórico calculados na etapa anterior para fazer a Krigagem e gerar o mapa de concentração de teor de Cobalto. Os comandos podem ser vistos abaixo.
data(jura)
coordinates(prediction.dat)=~Xloc+Yloc
gridded(jura.grid) = ~Xloc+Yloc
m <- vgm(12.524985, "Sph", 1.183543, 1.304965)
x <- krige(Co~1, prediction.dat, jura.grid, model=m)
spplot(x[1], col.regions=bpy.colors(40),main = "Mapa de concentração de Cobalto(ppm)", xlab="Xloc(0 a 6Km)", ylab= "Yloc(0 a 6Km)",xlim=c(0,6),ylim=c(0,6))
Conclusão
A partir do banco de dados "jura" foi possível aplicar as técnicas estatísticas e estimar um modelo teórico de variograma que se ajustou muito bem ao variograma experimental. Ao utilizar do modelo e da krigagem, podemos gerar mapas com dados mais próximos da realidade. Visto que, na Krigagem, em uma posição conhecida, o valor interpolado coincide com o valo medido.
No mapa de concentração de Cobalto, podemos identificar claramente que a continuidade das concentrações tendem a se alinhar com as direções NE-SW. Há uma faixa de altas concentrações que atravessa a região segundo a direção NE-SW. Além disso, na região NW temos um expressivo baixo teor de Cobalto. Em uma porção mais limitada à SE temos uma região de alta concentração.
As maiores frequência de concentração identificadas no histograma estão no intervalo de 3 a 5ppm (destacado em cores frias no mapa de concentração) e concentrações de 9 a 13ppm (destacado em cores quentes no mapa de concentração de Co).
O variograma experimental ajustado ao modelo teórico nos informa que, para a variavel Cobalto, os poços distantes até 1km de distância tem boa correlação. Analisando o mapa base, podemos medir a distância entre poços e encontrar valores em torno de 500m ou menos, o que nos dá uma boa continuidade dos dados nesta malha de amostragem.
No mapa de concentração de Cobalto, podemos identificar claramente que a continuidade das concentrações tendem a se alinhar com as direções NE-SW. Há uma faixa de altas concentrações que atravessa a região segundo a direção NE-SW. Além disso, na região NW temos um expressivo baixo teor de Cobalto. Em uma porção mais limitada à SE temos uma região de alta concentração.
As maiores frequência de concentração identificadas no histograma estão no intervalo de 3 a 5ppm (destacado em cores frias no mapa de concentração) e concentrações de 9 a 13ppm (destacado em cores quentes no mapa de concentração de Co).
O variograma experimental ajustado ao modelo teórico nos informa que, para a variavel Cobalto, os poços distantes até 1km de distância tem boa correlação. Analisando o mapa base, podemos medir a distância entre poços e encontrar valores em torno de 500m ou menos, o que nos dá uma boa continuidade dos dados nesta malha de amostragem.
Assinar:
Comentários (Atom)