Análise geoestatística do elemento Dy no estado do Piauí

O disprósio (Dy) é um elemento químico metálico de símbolo Dy e de número atómico igual a 66, com massa atómica 162,5 u. À temperatura ambiente, o disprósio encontra-se no estado sólido. Faz parte do grupo das terras raras.

Mapa do estado e a presença de Dy em seu território

Visualizando o mapa:

Gráfico com a inclusão de diferentes variáveis

Primeiramente, sem nenhum efeito de tendência

plot(geosolo,low=T)

## Componente de tendência na coordenada X

plot(geosolo,low=T, trend=~X1)

Componente de tendência na coordenada Y

plot(geosolo,low=T, trend=~Y1)

Componente de tendência na coordenada X e Y

plot(geosolo,low=T, trend=~X1+Y1)

Vimos acima que quando não aplicamos efeitos de tendência os dados nas coordenadas X e Y parecem ter um menor efeito de tendência. No quarto gráfico vemos a densidade dos resíduos, que possui uma grande concentração em torno do 0 e alguns valores positivos muito altos.

fazendo um gráfico de envelope

Se há ao menos um ponto fora dos intervalos do envelope simulado, então existe evidência de dependência espacial. Vemos abaixo que há uma dependência espacial

v1 = variog(geosolo, max.dist=400, uvec=seq(0,400,l=50))
## variog: computing omnidirectional variogram
env.var = variog.mc.env(geosolo, obj.v=v1, nsim=50) 
## variog.env: generating 50 simulations by permutating data values
## variog.env: computing the empirical variogram for the 50 simulations
## variog.env: computing the envelops
plot(v1, env=env.var)

Variograma para diferentes modelos

v1 = variog(geosolo, max.dist=400) 
## variog: computing omnidirectional variogram
plot(v1, xlim=c(0,333))

v1 = variog(geosolo, max.dist=400, trend=~X1) 
## variog: computing omnidirectional variogram
plot(v1, xlim=c(0,300))

v1 = variog(geosolo, max.dist=400, trend=~Y1) 
## variog: computing omnidirectional variogram
plot(v1, xlim=c(0,400))

v1 = variog(geosolo, max.dist=400, trend=~X1+Y1) 
## variog: computing omnidirectional variogram
plot(v1, xlim=c(0,400))

Usando o comando eyefit para visualizar valores iniciais para “sigmasq”, “phi” e “tausq”, temos que: “sigmasq” = 48.73, “phi”=128.2 e “tausq”=6.09

faremos então os primeiros modelos para o semivariograma

ml1  <- likfit(geosolo, ini=c(48.73,128.2), nug=6.09, cov.model = "matern")
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml2  <- likfit(geosolo, ini=c(48.73,128.2), nug=6.09, cov.model = "spherical")
## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

Para tendência em X, temos: “sigmasq”=33.85, “phi”=310.45 e “tausq”=4.23

ml3  <- likfit(geosolo, ini=c(33.85,310.45), nug=4.23, trend=~X1, cov.model = "matern")
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml4  <- likfit(geosolo, ini=c(33.85,310.45), nug=4.23, trend=~X1, cov.model = "spherical")
## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

Para tendência em Y, temos: “sigmasq”=48.73, “phi”=128.2 e “tausq”=6.09

ml5  <- likfit(geosolo, ini=c(48.73,128.2), nug=6.09, trend=~Y1, cov.model = "matern")
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml6  <- likfit(geosolo, ini=c(48.73,128.2), nug=6.09, trend=~Y1, cov.model = "spherical")
## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

Para tendência e X e Y, temos: “sigmasq”=32.63,“phi”=128.2 e “tausq”=4.08

ml7  <- likfit(geosolo, ini=c(32.63,128.2), nug=4.08, trend=~X1+Y1, cov.model="matern")
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml8  <- likfit(geosolo, ini=c(32.63,128.2), nug=4.08, trend=~X1+Y1, cov.model="spherical")
## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

Escolhendo o melhor modelo pelo AIC

ml1[17]
## $AIC
## [1] 1483.672
ml2[17]
## $AIC
## [1] 1483.124
ml3[17]
## $AIC
## [1] 1482.126
ml4[17]
## $AIC
## [1] 1482.548
ml5[17]
## $AIC
## [1] 1485.558
ml6[17]
## $AIC
## [1] 1485.058
ml7[17]
## $AIC
## [1] 1480.685
ml8[17]
## $AIC
## [1] 1479.097
#Modelo 8 com menor AIC

O melhor modelo é o oitavo, pois possui o menor AIC

plot(variog(geosolo, max.dist = 400,nugget.tolerance = 20))
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
lines.variomodel(ml8,col="black")
title("Ajuste do variograma ")

Predição espacial (krigagem)

GR <- pred_grid(geosolo$borders, by=6)
points(geosolo)
points(GR, pch=".",cex=1.5)

geosolo.kr <- krige.conv(geosolo, loc=GR, bor=geosolo$bor, krige=krige.control(obj.model=ml8))
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(main="Mapa dos quantis",geosolo.kr,col=terrain.colors(15),x.leg=c(500, 850), y.leg=c(0,50), xlim=c(-100,950))

explorando mais os resultados:

mapa de probabilidades: P(S(x)|y > 0.9), maior do que a mediana = 0.9

GR <- pred_grid(geosolo$borders, by=15)
OC <- output.control(thres=0.9,quantile=c(0.1, 0.9))
geosolo.kc <- krige.conv(geosolo, loc=GR, krige=krige.control(obj.model=ml8), out=OC)
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: Kriging performed using global neighbourhood
#names(geosolo.kc)

#dim(geosolo.kc$simula)

image(main="Mapa de Probabilidade para elemento Dy, considerando sua mediana=0.9",geosolo.kc, loc=GR, bor=geosolo$borders, col=terrain.colors(21), 
      val=1-geosolo.kc$prob,x.leg=c(500, 1000), y.leg=c(0,50))

Vemos que na região sudeste do estado de Piauí é onde tem maior incidência e maior probabilidade de encontrarmos o elemento Dy nos solos do estado, levando em consideração que a quantidade mediana de Dy no solo do Piauí é de 0.9.