## carregando a geoR require(geoR) ## vendo os conjuntos de dados disponíveis na geoR data(package="geoR") ## carregando o conjunto de dados ksat data(Ksat) ## obtendo informações sobre os dados # help(ksat) ## carregando o conjunto de dados ca20 data(ca20) ## vendo informações sobre os dados # help(ca20) ## resumo dos dados summary(Ksat) ## vendo os elementos do objeto Ksat names(Ksat) ## visualizando os dados points(Ksat, bor=borders) # outra visualização dos dados plot(Ksat, bor=borders) # os dados são fortemente assimétricos e além disto a transformação # log() deve ser razoável dada a natureza da variável # Novo plot com dados transformados plot(Ksat, bor=borders, lambda=0) ## carregando o conjunto de dados parana data(parana) ## vendo informações sobre os dados # help(parana) names(parana) summary(parana) points(parana, bor=borders) plot(parana, bor=borders) ## os dados exibem uma clara tendência com as coordenadas ## reexaminando os dados, descontando o efeito da tendência plot(parana, bor=borders, trend="1st") ## dados de cálcio do solo data(ca20) # help(ca20) summary(ca20) names(ca20) points(ca20, bor=borders) polygon(ca20$reg1) polygon(ca20$reg2) plot(ca20, bor=borders) ## investigando relação com covariáveis with(ca20, plot(covariate$altitude, data)) with(ca20, plot(covariate$area, data)) ## refazendo plot removend efeito das áreas plot(ca20, bor=borders, trend=~area) ## Buscando a melhor transformação para os dados require(MASS) boxcox(Ksat) boxcox(Ksat, lambda=seq(-1,1,lam=101)) boxcox(parana, trend="1st") boxcox(ca20, trend=~area) data(SIC) boxcox(sic.100) boxcox(sic.all, lambda=seq(0, 1, l=101)) plot(Ksat, bor=borders) plot(Ksat, bor=borders, lam=0) plot(parana, bor=borders, trend="1st") plot(parana, bor=borders, trend="1st", lam=0.5) plot(sic.all, bor=sic.borders, lambda=0.5) ## explorando a presença de dependência espacial ## compare os gráficos dos dados originais com os obtidos ## após permutar os dados nas diferentes localizações par(mfrow=c(1,2), mar=c(3,3,.5, .5), mgp=c(2,1,0)) points(Ksat, lam=0, bor=borders, cex.max=1, cex.min=1, pt.div="quint") points(Ksat, lam=0, bor=borders, cex.max=1, cex.min=1, pt.div="quint", permute=T) par(mfrow=c(1,1)) par(mfrow=c(1,2)) points(parana, lam=0.5, bor=borders, trend="1st", cex.min=1,cex.max=1, pt.div="quint") points(parana, lam=0.5, bor=borders, trend="1st", cex.min=1,cex.max=1, pt.div="quint", permute=T) par(mfrow=c(1,1)) par(mfrow=c(1,2)) points(ca20, bor=borders, trend=~area, cex.min=1, cex.max=1, pt.div="quint") points(ca20, bor=borders, trend=~area, cex.min=1, cex.max=1, pt.div="quint", permute=T) par(mfrow=c(1,1)) ## variogramas Ksat.v <- variog(Ksat, max.dist=8.5, uvec=seq(0,8,by=1.25), lam=0) plot(Ksat.v) parana.v <- variog(parana, trend="1st", lam=0.5, max.dist=500) plot(parana.v) ca20.v <- variog(ca20, trend=~area, max.dist=450, op="cloud") plot(ca20.v) ca20.v <- variog(ca20, trend=~area, max.dist=450) plot(ca20.v) ## ajusta variogramas "a olho" com a função eyefit ca20.ef <- eyefit(ca20.v) ca20.ef Ksat.ef <- eyefit(Ksat.v) Ksat.ef ## Interpolando valores ## definindo um grid de pontos cobrindo a área grid0 <- expand.grid(seq(0, 22.5, by=0.25) , seq(0, 13, by=0.25)) ## visualizando o grid points(Ksat, bor=borders) points(grid0, pch="+", col="green") # selecionando os pontos dentro das bordas da area grid1 <- polygrid(grid0, borders=Ksat$borders) points(grid1, pch="+", col="blue") ## procedendo a interpolação Ksat.in <- krige.conv(Ksat, loc=grid0, borders=Ksat$borders, krige=krige.control(obj=Ksat.ef)) # visualizando resultados image(Ksat.in, col=gray(seq(1, 0.1, l=21))) image(Ksat.in, val=sqrt(Ksat.in$krige.var),col=gray(seq(1, 0.1, l=21))) save.image() ## estimativa por ajuste de mínimos quadrados Ksat.vf <- variofit(Ksat.v, ini=c(1.7, 3), nug=0.5) ## ou pode usar o eyefit como valor inicial... Ksat.vf <- variofit(Ksat.v, ini=Ksat.ef) ## estimativa de parametros pela verossimilhanca Ksat.ml <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0) ## ou pode usar o eyefit como valor inicial... Ksat.ml <- likfit(Ksat, ini=Ksat.ef, lam=0) ## comparando modelos Ksat.ml1 <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0, cov.model="sph") Ksat.ml2 <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0, kappa=1) Ksat.ml3 <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0, kappa=1.5) Ksat.ml4 <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0, kappa=2.5) Ksat.ml5 <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0, kappa=3.5) Ksat.ml6 <- likfit(Ksat, ini=c(1.7, 3), nug=0.5, lam=0, kappa=4.5) options(digits=7) Ksat.ml$loglik Ksat.ml2$loglik Ksat.ml3$loglik Ksat.ml4$loglik Ksat.ml5$loglik Ksat.ml6$loglik kps <- c(0.5, 1, 1.5, 2.5, 3.5, 4.5) llk <- c(Ksat.ml$loglik,Ksat.ml2$loglik,Ksat.ml3$loglik,Ksat.ml4$loglik,Ksat.ml5$loglik, Ksat.ml6$loglik) plot(kps, llk) ## selecionando modelo para ca20 ca20.ml1 <- likfit(ca20, ini=c(60, 100), nug=30) ca20.ml2 <- likfit(ca20, ini=c(60, 100), nug=30, trend=~area) ca20.ml3 <- likfit(ca20, ini=c(60, 100), nug=30, trend=~area) ca20.ml4 <- likfit(ca20, ini=c(60, 100), nug=30, trend=~area+altitude) ca20.ml5 <- likfit(ca20, ini=c(60, 100), nug=30, trend=~area+ca20$coords[,2]) ca20.ml6 <- likfit(ca20, ini=c(60, 100), nug=30, trend=~area+altitude+ca20$coords[,2]) ca20.fits <- rbind( c(ca20.ml1$npars,ca20.ml2$npars,ca20.ml3$npars,ca20.ml4$npars,ca20.ml5$npars,ca20.ml6$npars), c(ca20.ml1$loglik,ca20.ml2$loglik,ca20.ml3$loglik,ca20.ml4$loglik,ca20.ml5$loglik,ca20.ml6$loglik), c(ca20.ml1$AIC,ca20.ml2$AIC,ca20.ml3$AIC,ca20.ml4$AIC,ca20.ml5$AIC,ca20.ml6$AIC), c(ca20.ml1$BIC,ca20.ml2$BIC,ca20.ml3$BIC,ca20.ml4$BIC,ca20.ml5$BIC,ca20.ml6$BIC)) colnames(ca20.fits) <- paste("modelo", 1:6) rownames(ca20.fits) <- c("np", "ll", "AIC", "BIC") ca20.fits ## mais opções na interpolação espacial args(krige.conv) args(output.control) OC <- output.control(n.post=500, quantile=c(0.1, 0.9), thres=2) Ksat.in <- krige.conv(Ksat, loc=grid0, borders=Ksat$borders, krige=krige.control(obj=Ksat.ef), output=OC) names(Ksat.in) dim(Ksat.in$simulations) dim(Ksat.in$quantile) par(mfcol=c(2,3), mar=c(2,2,0,0), mgp=c(1.8, .8, 0)) image(Ksat.in, col=gray(seq(1, 0, l=21)), x.leg=c(5,20), y.leg=c(-4,-2.5)) text(12, 15, "Mapa de E[Y_0|y]") image(Ksat.in, val=sqrt(Ksat.in$krige.var),col=gray(seq(1, 0, l=21)), x.leg=c(5,20), y.leg=c(-4,-2.5)) text(12, 15, "Mapa de SD[Y_0|y]") image(Ksat.in, val=Ksat.in$simulations[,23],col=gray(seq(1, 0, l=21)), x.leg=c(5,20), y.leg=c(-4,-2.5)) text(12, 15, "Mapa de uma simulação de [Y_0|y]") image(Ksat.in, val=Ksat.in$prob,col=gray(seq(1, 0.2, l=21)), x.leg=c(5,20), y.leg=c(-4,-2.5)) text(12, 15, "Mapa de probabilidade P[Y_0 < 2]") image(Ksat.in, val=Ksat.in$quant[,1],col=gray(seq(1, 0, l=21)), x.leg=c(5,20), y.leg=c(-4,-2.5)) text(12, 15, "Mapa de do quantil 0.10") image(Ksat.in, val=Ksat.in$quant[,2],col=gray(seq(1, 0, l=21)), x.leg=c(5,20), y.leg=c(-4,-2.5)) text(12, 15, "Mapa de do quantil 0.90")