## importando dados file.show("dadosIntro.txt") #file.show("dadosintro.txt") dados <- read.table("http://www.leg.ufpr.br/~paulojus/geo2012/dadosIntro.txt", head=T) dados ## instalando o pacote geoR (se necessário) e suas dependencias #install.packages("geoR", dep=TRUE) ## carregando o pacote geoR require(geoR) dG <- as.geodata(dados, coords.col=c(1,2), data.col=3) class(dG) ## vendo os conjuntos de dados disponíveis no pacote #data(package="geoR") x11() ## carregando o conjunto de dados dG ####data(dG) ## visualizando os dados points(dG) args(points.geodata) points.geodata(dG,pt.div="quint") points.geodata(dG,pt.div="quint", cex.min=1, cex.max=1) ##points.geodata(coords=dados[,1:2], data=dados[,3]) ## um resumo dos dados summary(dG) summary(dG, lambda=0) dG ## outra visualização dos dados x11() plot(dG) plot(dG, lambda=0) plot(dG, low=T) args(plot.geodata) points(dG) points(dG, lambda=0) points(dG, pt.div="quint", cex.min=1, cex.max=1) x11() ?points.geodata points(dG, pt.div="quint", cex.min=1, cex.max=1,permute=T) ## obtendo o variograma v <- variog(dG) plot(v) vC <- variog(dG, option="cloud") vC <- variog(dG, max.dist=100, option="cloud") plot(vC) v <- variog(dG, max.dist=100) plot(v) v <- variog(dG, max.dist=80) plot(v) v <- variog(dG, uvec=seq(0, 80, by=10)) plot(v) v <- variog(dG, uvec=seq(0, 80, by=5)) plot(v) vS <- variog(dG, max.dist=100, option="smooth") plot(vS, ylim=c(0, 650)) v.mc <- variog.mc.env(dG, obj.variog=v) plot(v, env=v.mc) ## escolhendo funcao de variograma plot(v) ef <- eyefit(v) ef vf <- variofit(v, ini=ef) vf ## estimando a Media ONES <- rep(1, nrow(dados)) D <- as.matrix(dist(dG$coords)) names(vf) vf$cov.pars SIGMA <- vf$nugget * diag(nrow(dados)) + vf$cov.pars[1] * exp(-D/vf$cov.pars[2]) dim(SIGMA) iS.ONES <- solve(SIGMA, ONES) (mu.gls <- solve(crossprod(iS.ONES, ONES), crossprod(iS.ONES, dG$data))) ## estimativa por máxima verossimilhança ## função de (Log)verossimilhança (escrita e forma ingenua) ltheta <- function(pars, gd, printpars=FALSE){ ## pars = mu, sigma, phi, tau n <- nrow(gd$coords) D <- as.matrix(dist(gd$coords, upper=T, diag=T)) SIGMA <- pars[4]^2 * diag(n) + pars[2]^2 * exp(-D/pars[3]) logdetS <- determinant(SIGMA, log=TRUE)$modulus res <- gd$data - pars[1] Q <- crossprod(res, solve(SIGMA, res)) logVero <- -0.5*(n*log(2*pi) + logdetS + Q) if(printpars) print(c(pars, logVero)) return(logVero) } lthetaConc <- function(pars, gd, printpars=FALSE){ ## pars : phi e nu = tau/sigma n <- nrow(gd$coords) D <- as.matrix(dist(gd$coords, upper=T, diag=T)) SIGMA <- pars[1]^1 * diag(n) + exp(-D/pars[1]) logdetS <- determinant(SIGMA, log=TRUE)$modulus mu <- res <- gd$data - pars[1] Q <- crossprod(res, solve(SIGMA, res)) logVero <- -0.5*(n*log(2*pi) + logdetS + Q) if(printpars) print(c(pars, logVero)) return(logVero) } names(vf) (parsVF <- with(vf, c(mu.gls, sqrt(cov.pars[1]), cov.pars[2], sqrt(nugget)))) ltheta(parsVF, gd = dG) ltheta(c(50, 15, 20, 0), gd = dG) args(optim) parsML <- optim(c(50, 15, 20, 0), ltheta, gd = dG, printpars=T, method="L-BFGS-B", lower=c(-Inf, 0, 0, 0), control=list(fnscale=-1)) parsML$par^c(1,2,1,2) vf ## usano a fuincao da geoR ml <- likfit(dG, ini=c(300, 20)) ml ## definindo um "grid" de predição gr <- expand.grid(seq(0,100, len=51), seq(0,100, len=51)) ## visualizando o grido onde seráo feitas as predições points(dG) points(gr, pch=19, cex=0.25, col=2) ## fazendo a predição espacial (krigagem) kc1 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[1]])) kc2 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[2]])) kc3 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[3]])) kc.ml <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ml)) names(kc.ml) ## visualizando os valores preditos na forma de um mapa par(mfrow=c(2,2), mar=c(2,2,1,1)) image(kc1) image(kc1, col=terrain.colors(21)) image(kc1, col=gray(seq(1,0, l=21))) image(kc1, col=gray(seq(0,1, l=4))) points(dG, add=T) par(mfrow=c(1,1), mar=c(2,2,1,1)) image(kc.ml) ## visualizando os erros padrão de predição (mapa de incerteza de predição) image(kc1, val=sqrt(kc1$krige.var), coords.data=dG$coords) image(kc2, val=sqrt(kc1$krige.var), coords.data=dG$coords) image(kc.ml, val=sqrt(kc.ml$krige.var)) ## escala comum ZL <- range(c(kc1$pred, kc2$pred, kc3$pred)) par(mfrow=c(1,3), mar=c(2,2,1,1)) image(kc1, zlim=ZL) image(kc2, zlim=ZL) image(kc3, zlim=ZL) ## outras opções de visualização contour(kc1) image(kc1, col=terrain.colors(21)) points(dG, add=T) contour(kc1, add=T, nlev=21) persp(kc1) persp(kc1, theta=20, phi=35) ## ## exportando predições ## ## exportando somente um vetor dos atributos (formato texto) write(kc1$pred, file="kc1.txt", ncol=1) file.show("kc1.txt") ## tecle c para terminar a exibição do arquivo ## exportando somente os atributos como matriz (formato texto) write(kc1$pred, ncol=51, file="kc2.txt") ## exportando atributos e locações de predição names(kc1) write.table(cbind(gr, data.frame(kc1[c(1,2)])),file="kc3.txt", row.names = FALSE) ## ou... write.csv(cbind(gr, data.frame(kc1[c(1,2)])),file="kc4.csv", row.names = FALSE) ## ou... write.csv2(cbind(gr, data.frame(kc1[c(1,2)])),file="kc5.csv", row.names = FALSE) ## simulacoes condicionais args(krige.conv) args(output.control) kc.ml <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ml), output = output.control(n.pred=1000)) names(kc.ml) dim(kc.ml$simulations) ## a krigagem (pos simulacao) correponde a predSim <- apply(kc.ml$simulations, 1, mean) plot(kc.ml$pred, predSim) abline(0,1) summary(kc.ml$pred - predSim) ## mapa de prob de estar acima de 50 p50 <- apply(kc.ml$simulations, 1, function(x) mean(x>50)) image(kc.ml, val=p50) contour(kc.ml, val=p50) image(kc.ml, val=p50) ## suponha agora que definimos como regiao de alto risco ## as quye possuiram prob > 0.7 de estar acima de 50 p50.alto <- ifelse(p50 < 0.7, 0, 1) image(kc.ml, val=p50.alto, col=c(1,0)) ## agora definindo 5 niveis de risco em fc da prob image(kc.ml, val=p50, breaks=seq(0, 1, by=0.2), col=c("blue","green","yellow","red","black")) bor <- locator(type="l") polygon(bor, col=4) ## completar aqui a krigagem dentrop do polygono plot(dG) v <- variog(dG, max.dist=80) plot(v) args(variog) v0 <- variog(dG, max.dist=80, direction=0) v45 <- variog(dG, max.dist=80, direction=pi/4) v90 <- variog(dG, max.dist=80, direction=pi/2) v135 <- variog(dG, max.dist=80, direction=3*pi/4) lines(v0) lines(v45, col=2) lines(v90, col=3) lines(v135,col=4) v4 <- variog4(dG, max.dist=100) plot(v4) plot(v4, omn=T) ml.iso <- likfit(dG, ini=c(500, 20)) ml.aniso <- likfit(dG, ini=c(500, 20), fix.psiR=F, fix.psiA=F) ml.anisoR <- likfit(dG, ini=c(500, 20), fix.psiR=F, fix.psiA=T, psiA=0) logLik(ml.iso) logLik(ml.anisoR) logLik(ml.aniso) ## Inf. Bayesiana summary(dG) args(krige.bayes) args(model.control) MC <- model.control() args(prior.control) PC <- prior.control(phi.discrete=seq(0, 40, len=41)) args(output.control) OC <- output.control(n.post=5000) kb <- krige.bayes(dG, model=MC, prior=PC, output=OC) names(kb) plot(kb) ## outra priorti PC <- prior.control(phi.discrete=seq(0, 40, len=41), phi.prior="rec") kb <- krige.bayes(dG, model=MC, prior=PC, output=OC) plot(kb) ## priori tb tem tau^2_Rel PC <- prior.control(phi.discrete=seq(0, 50, len=51), phi.prior="rec", tausq.rel.discrete=seq(0, 1.5, length=16), tausq.rel.prior="rec") kb <- krige.bayes(dG, model=MC, prior=PC, output=OC) par(mfrow=c(1,2)) plot(kb) kb <- krige.bayes(dG, loc=gr, model=MC, prior=PC, output=OC) names(kb) names(kb$post) names(kb$pred) image(kb) p50bayes <- apply(kb$pred$simul, 1, function(x) mean(x>50)) image(kb, values = p50bayes) ## comparando as predicoes image(kc.ml) image(kb) plot(kc.ml, kb$predictive$pred)