Outubro de 2017

Pacotes Utilizados

require(gstat)
require(sp)

O que é GEOESTATÍSTICA

Análise gráfica

require(RgoogleMaps)
require(googleVis)

Dados Meuse

  • O conjunto de dados meuse fornecido pelo pacote sp

  • Composto por quatro metais pesados medidos no solo superior em uma planície de inundação ao longo do rio Meuse

  • Possui outras covariáveis

  • Os sedimentos poluídos é transportado pelo rio, e na maior parte depositado perto da margem do rio, e em áreas com baixa altitude.

– Publicação: Burrough, P.A., R.A. McDonnell, 1998. Principles of Geographical Information Systems, 2nd Edition. Oxford University Press.

Introdução

  • A mineração histórica de metais causou a dispersão generalizada de chumbo, zinco, cobre e cádmio no solo aluvial

  • Os poluentes podem restringir o uso da terra nessas áreas

  • Assim, são necessários mapas detalhados que identifiquem zonas com altas concentrações

  • Nosso objetivo específico será gerar um mapa de um metal pesado (zinco) no solo usando observações pontuais.

Lendo os dados

data(meuse)
class(meuse)
## [1] "data.frame"
names(meuse)
##  [1] "x"       "y"       "cadmium" "copper"  "lead"    "zinc"    "elev"   
##  [8] "dist"    "om"      "ffreq"   "soil"    "lime"    "landuse" "dist.m"
  • Veja que no help(meuse) tem todas as informações a respeito dos dados

Os dados

str(meuse)
## 'data.frame':    155 obs. of  14 variables:
##  $ x      : num  181072 181025 181165 181298 181307 ...
##  $ y      : num  333611 333558 333537 333484 333330 ...
##  $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##  $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
##  $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
##  $ zinc   : num  1022 1141 640 257 269 ...
##  $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
##  $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
##  $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##  $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##  $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##  $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##  $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...

Estudando a variável zinc

zinco <- meuse$zinc
hist(zinco, breaks=20, prob=T, col="purple", ylim=c(0,0.004))
lines(density(zinco), lwd=3, col="gray20")

Estudando a variável zinc

boxplot(zinco, col="purple", horizontal = TRUE)

gstat

  • Neste curso será utilizado o pacote gstat para análise geoestatística.

  • No Brasil, o mais conhecido é o seu "concorrente": geoR

No gstat é possível:

  • Variogramas

  • Ajustes de modelos

  • Krigagem

  • CoKrigagem

  • Entre outras funcionalidade

sp

- O pacote gstat trabalha em conjunto com a biblioteca sp

- O sp fornece classes e métodos para dados espaciais

- Os dados meuse faz parte da sua biblioteca

- A função coordinates: Definir coordenadas espaciais para criar um objeto Espacial

Especificando que são nossas coordenadas

coordinates(meuse) = ~x + y

Mudança da classe dos dados

class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Análise Gráfica Descritiva

bubble(meuse, "zinc")

Alterando os quartis

zq <- quantile(meuse$zinc, seq(0,1,0.2))
zq
##     0%    20%    40%    60%    80%   100% 
##  113.0  186.8  246.4  439.6  737.2 1839.0
bubble(meuse, "zinc", key.entries = zq)

Análise Gráfica Descritiva

spplot(meuse['zinc'])

Análise Gráfica Descritiva

spplot(meuse['zinc'], scales=list(draw=T))

Análise Gráfica Descritiva

spplot(meuse['zinc'], scales=list(draw=T), key.space="right")

Análise Gráfica Descritiva

spplot(meuse['zinc'], scales=list(draw=T), key.space="right",
       colorkey=T)

Estudando Tendência

  • Em alguns casos a variável (atributo) possui uma relação com as próprias coordenadas
par(mfrow=c(1,2))
plot(zinc~x, data=meuse)
plot(zinc~y, data=meuse)

Log(meuse)

par(mfrow=c(1,2))
plot(log(zinc)~x, data=meuse)
plot(log(zinc)~y, data=meuse)

O Semivariograma

lzn.vgm <- variogram(log(zinc)~1, meuse)
lzn.vgm
##     np       dist     gamma dir.hor dir.ver   id
## 1   57   79.29244 0.1234479       0       0 var1
## 2  299  163.97367 0.2162185       0       0 var1
## 3  419  267.36483 0.3027859       0       0 var1
## 4  457  372.73542 0.4121448       0       0 var1
## 5  547  478.47670 0.4634128       0       0 var1
## 6  533  585.34058 0.5646933       0       0 var1
## 7  574  693.14526 0.5689683       0       0 var1
## 8  564  796.18365 0.6186769       0       0 var1
## 9  589  903.14650 0.6471479       0       0 var1
## 10 543 1011.29177 0.6915705       0       0 var1
## 11 500 1117.86235 0.7033984       0       0 var1
## 12 477 1221.32810 0.6038770       0       0 var1
## 13 452 1329.16407 0.6517158       0       0 var1
## 14 457 1437.25620 0.5665318       0       0 var1
## 15 415 1543.20248 0.5748227       0       0 var1

O Semivariograma

plot(lzn.vgm)

Ajustando o modelo

head(vgm())
##   short                               long
## 1   Nug                       Nug (nugget)
## 2   Exp                  Exp (exponential)
## 3   Sph                    Sph (spherical)
## 4   Gau                     Gau (gaussian)
## 5   Exc Exclass (Exponential class/stable)
## 6   Mat                       Mat (Matern)

Modelo Esférico

Sph.fit = fit.variogram(lzn.vgm, 
                        model = vgm(psill = 0.5, model = "Sph",
                                    range = 900, nugget = 0.2))
Sph.fit
##   model      psill    range
## 1   Nug 0.05066017   0.0000
## 2   Sph 0.59060556 897.0044

Modelo Esférico

plot(lzn.vgm, Sph.fit)

Modelo Exponencial

Exp.fit = fit.variogram(lzn.vgm, 
                        model = vgm(psill = 0.5, model = "Exp",
                                    range = 900, nugget = 0.2))
Exp.fit
##   model     psill    range
## 1   Nug 0.0000000   0.0000
## 2   Exp 0.7186599 449.7668

Modelo Exponencial

plot(lzn.vgm, Exp.fit)

O grid

data("meuse.grid")
gridded(meuse.grid) = ~x+y
summary(meuse.grid)
## Object of class SpatialPixelsDataFrame
## Coordinates:
##      min    max
## x 178440 181560
## y 329600 333760
## Is projected: NA 
## proj4string : [NA]
## Number of points: 3103
## Grid attributes:
##   cellcentre.offset cellsize cells.dim
## x            178460       40        78
## y            329620       40       104
## Data attributes:
##      part.a           part.b            dist        soil     ffreq   
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   1:1665   1: 779  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.1193   2:1084   2:1335  
##  Median :0.0000   Median :1.0000   Median :0.2715   3: 354   3: 989  
##  Mean   :0.3986   Mean   :0.6014   Mean   :0.2971                    
##  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.4402                    
##  Max.   :1.0000   Max.   :1.0000   Max.   :0.9926

Krigagem Ordinária

Ord.kriged = krige(log(zinc)~1,
                   meuse,
                   meuse.grid, 
                   model = Exp.fit)
## [using ordinary kriging]

Krigagem Ordinária

spplot(Ord.kriged['var1.pred'], scales=list(draw=T), key.space="right",
       colorkey=T,
       main="Krigagem Ordinária")

Krigagem Simples

Spl.kriged = krige(log(zinc)~1,
                   meuse,
                   meuse.grid, 
                   model = Exp.fit,
                   beta = 7)
## [using simple kriging]

Krigagem Simples

spplot(Spl.kriged['var1.pred'], scales=list(draw=T), key.space="right",
       colorkey=T,
       main="Krigagem Simples")

Krigagem Universal

Krigagem Universal

  • Modelando a tendência
plot(log(zinc)~sqrt(dist), data=meuse)

Krigagem Universal

  • Ajustando um modelo
mod <- lm(log(zinc)~sqrt(dist), data=meuse)
summary(mod)
## 
## Call:
## lm(formula = log(zinc) ~ sqrt(dist), data = meuse)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04624 -0.29060 -0.01869  0.26445  1.59685 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.99438    0.07593   92.12   <2e-16 ***
## sqrt(dist)  -2.54920    0.15498  -16.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4353 on 153 degrees of freedom
## Multiple R-squared:  0.6388, Adjusted R-squared:  0.6364 
## F-statistic: 270.6 on 1 and 153 DF,  p-value: < 2.2e-16

Krigagem Universal

  • Gráfico do ajuste
plot(log(zinc)~sqrt(dist), data=meuse, cex=0.8, pch=19)
curve(coef(mod)[1]+coef(mod)[2]*x, col='blue', add=T, lwd=2)

Krigagem Universal

  • Variograma Resisual
meuse$res <- residuals(mod)
vt <- variogram(res~1, meuse)
plot(vt)

Krigagem Universal

  • Ajustando um modelo
Exp.fit = fit.variogram(vt, 
                        model = vgm(psill = 0.15, model = "Exp",
                                    range = 900, nugget = 0.05))
Exp.fit
##   model     psill    range
## 1   Nug 0.0571200   0.0000
## 2   Exp 0.1764149 340.3013

Krigagem Universal

  • Ajustando um modelo
plot(vt, Exp.fit)

Krigagem Universal

  • Valores a serem preditos
new <- data.frame(meuse.grid)
names(new)
## [1] "part.a" "part.b" "dist"   "soil"   "ffreq"  "x"      "y"
new$pred <- predict(mod, new)

Krigagem Universal

  • Krigagem nos resíduos
lz.ukRes <- krige(res~1, meuse, meuse.grid, Exp.fit)
## [using ordinary kriging]
res_kg <- lz.ukRes$var1.pred

Krigagem Universal

  • Adicionando a tendencia
lz.ukRes$pred_kg <- new$pred + res_kg

Krigagem Universal

spplot(lz.ukRes['pred_kg'], scales=list(draw=T), key.space="right",
       colorkey=T,
       main="Krigagem Universal")

Krigagem Universal: direto no gstat

  • Krigagem gstat
lz.ukGstat <- krige(log(zinc)~sqrt(dist), meuse, meuse.grid, Exp.fit)
## [using universal kriging]

Krigagem Universal: direto no gstat

  • Krigagem gstat
pred_kg_gstat <- lz.ukGstat$var1.pred
hist(pred_kg_gstat-lz.ukRes$pred_kg)