######################################################################## ## Taxa de infiltração básica de água no solo para 3 diferentes solos ## ## Jefferson Vieira José jfvieira@hotmail.com.br ## ######################################################################## rm(list=ls()) ## carregando os pacotes require(foreign) require(geoR) require(MASS) ## carregando o conjunto de dados dados <- read.table ('C://Documents and Settings/Administrador/Desktop/vib.txt',header=TRUE,sep="",dec=".") # Diferentes Classes de solo Z==1,2,3 TIBO1 <- as.geodata(subset(dados, Z==1), coords.col=c(2,3), data.col = 10) TIBO2 <- as.geodata(subset(dados, Z==2), coords.col=c(2,3), data.col = 10) TIBO3 <- as.geodata(subset(dados, Z==3), coords.col=c(2,3), data.col = 10) ############################################################### ## Análise para Z==1 CTI_Nitossolo vermelho ## ############################################################### summary(TIBO1$data) plot(TIBO1, low=T) boxplot(TIBO1$data,low=T) x11() par(mfrow=c(1,2),tcl=-0.2) points(TIBO1, pt.div="quint", cex.min=0.5, cex.max=1.5,main="Dados originais") points(TIBO1, pt.div="quint", cex.min=0.5, cex.max=1.5, permute=T, main="Sem padrão espacial") boxcox(TIBO1) # inclui o zero -> transformação log. ## Variograma... determinar os pontos inciais para o likfi summary(TIBO1) plot(variog(TIBO1, max.dist=14)) plot(variog(TIBO1, max.dist=11)) plot(variog(TIBO1, max.dist=7)) v <- variog(TIBO1, max.dist=11,lambda=0) # escolha subjetiva. preferi este. x11() ef <- eyefit(v) # exponencial. ## Ajuste dos modelo fit <- likfit(TIBO1, ini=ef,cov.model="exp",lambda=0) fit ## definindo um "grid" de predição gr <- expand.grid(seq(1,11, len=50), seq(1,11, len=50)) ## visualizando o grid onde serão feitas as predições points(TIBO1) points(gr, pch=19, cex=0.1, col=2) ## fazendo a predição espacial (krigagem) kc <- krige.conv(TIBO1, loc=gr, krige = krige.control(obj.model=fit)) ## visualizando os valores preditos na forma de um mapa image(kc, col=terrain.colors(21)) points(TIBO1, add=T) ## visualizando os erros padrão de predição (mapa de incerteza de predição) image(kc, val=sqrt(kc$krige.var)) ## outras opções de visualização contour(kc) image(kc, col=terrain.colors(21),useRaster=T) leg(-5,max(TIBO1$coords[,2]),zlim=range(kc$predict),col=terrain.colors(21) ) points(TIBO1, add=T) contour(kc, add=T, nlev=29) persp(kc) persp(kc, theta=20, phi=35) ## mapas de probabilidades OC <- output.control(quantile=c(0.1, 0.5, 0.9)) kc <- krige.conv(TIBO1, loc=gr, krige=krige.control(obj.model=fit), output=OC) names(kc) par(mfrow=c(1,2)) PROBS <- apply(kc$simulation, 2, function(x) mean(x < 12)) hist(PROBS, prob=T, main="Distribuição preditiva na da proporção da área com P[S|y < 10]", panel.first=lines(density(PROBS))) P1 <- apply(kc$simulations, 1, function(x) mean(x < 12)) image(kc, val=P1,useRaster=T, col=gray(seq(1,0,l=21)),,main="Localizações na área com P[S|y < 10]") leg(-5,max(TIBO1$coords[,2]),zlim=range(PROBS),col=gray(seq(1,0,l=21))) ############################################################### ## Análise para Z==2 FEI_Latossolo Vermelho ## ############################################################### summary(TIBO2$data) plot(TIBO2, low=T) boxplot(TIBO2$data,low=T) x11() par(mfrow=c(1,2),tcl=-0.2) points(TIBO2, pt.div="quint", cex.min=0.5, cex.max=1.5,main="Dados originais") points(TIBO2, pt.div="quint", cex.min=0.5, cex.max=1.5, permute=T, main="Sem padrão espacial") boxcox(TIBO2) # inclui o zero -> transformação log. ## Variograma... determinar os pontos inciais para o likfi summary(TIBO2) plot(variog(TIBO2, max.dist=17)) plot(variog(TIBO2, max.dist=13)) plot(variog(TIBO2, max.dist=12)) v <- variog(TIBO2, max.dist=10,) # escolha subjetiva. preferi este. x11() ef <- eyefit(v) # exponencial. ## Ajuste dos modelo fit <- likfit(TIBO2, ini=ef,cov.model="exp",lambda=0) fit ## definindo um "grid" de predição gr <- expand.grid(seq(1,11, len=50), seq(1,11, len=50)) ## visualizando o grid onde serão feitas as predições points(TIBO2) points(gr, pch=19, cex=0.1, col=2) ## fazendo a predição espacial (krigagem) kc <- krige.conv(TIBO2, loc=gr, krige = krige.control(obj.model=fit)) ## visualizando os valores preditos na forma de um mapa image(kc, col=terrain.colors(21)) points(TIBO2,lambda=0, add=T) ## visualizando os erros padrão de predição (mapa de incerteza de predição) image(kc, val=sqrt(kc$krige.var)) ## outras opções de visualização contour(kc) image(kc, col=terrain.colors(21),useRaster=T) leg(-5,max(TIBO2$coords[,2]),zlim=range(kc$predict),col=terrain.colors(21) ) points(TIBO2, add=T) contour(kc, add=T, nlev=29) persp(kc) persp(kc, theta=20, phi=35) ## mapas de probabilidades OC <- output.control(quantile=c(0.1, 0.5, 0.9) kc <- krige.conv(TIBO2, loc=gr, krige=krige.control(obj.model=fit), output=OC) names(kc) par(mfrow=c(1,2)) PROBS <- apply(kc$simulation, 2, function(x) mean(x > 70)) hist(PROBS, prob=T, main="Distribuição preditiva \n da proporção da área com P[S|y < 10]", panel.first=lines(density(PROBS))) P1 <- apply(kc$simulations, 1, function(x) mean(x > 70)) image(kc, val=P1,useRaster=T, col=gray(seq(1,0,l=21)),,main="Localizações na área com P[S|y < 10]") leg(-5,max(TIBO2$coords[,2]),zlim=range(PROBS),col=gray(seq(1,0,l=21))) ############################################################### ## Análise para Z==3 IAPAR_Latossolo Vermelho ## ############################################################### summary(TIBO3$data) plot(TIBO3, low=T) boxplot(TIBO3$data,low=T) x11() par(mfrow=c(1,2),tcl=-0.2) points(TIBO3, pt.div="quint", cex.min=0.5, cex.max=1.5,main="Dados originais") points(TIBO3, pt.div="quint", cex.min=0.5, cex.max=1.5, permute=T, main="Sem padrão espacial") boxcox(TIBO3) # inclui o zero -> transformação log. ## Variograma... determinar os pontos inciais para o likfi summary(TIBO3) plot(variog(TIBO3, max.dist=14)) plot(variog(TIBO3, max.dist=13)) plot(variog(TIBO3, max.dist=7)) v <- variog(TIBO3, max.dist=13,lambda=0) # escolha subjetiva. preferi este. x11() ef <- eyefit(v) # exponencial. ## Ajuste dos modelo fit <- likfit(TIBO3, ini=ef,cov.model="exp",lambda=0) fit ## definindo um "grid" de predição gr <- expand.grid(seq(1,11, len=50), seq(1,11, len=50)) ## visualizando o grid onde serão feitas as predições points(TIBO3) points(gr, pch=19, cex=0.1, col=2) ## fazendo a predição espacial (krigagem) kc <- krige.conv(TIBO3, loc=gr, krige = krige.control(obj.model=fit)) ## visualizando os valores preditos na forma de um mapa image(kc, col=terrain.colors(21)) points(TIBO3, add=T) ## visualizando os erros padrão de predição (mapa de incerteza de predição) image(kc, val=sqrt(kc$krige.var)) ## outras opções de visualização contour(kc) image(kc, col=terrain.colors(21),useRaster=T) leg(-5,max(TIBO3$coords[,2]),zlim=range(kc$predict),col=terrain.colors(21) ) points(TIBO2, add=T) contour(kc, add=T, nlev=20) persp(kc) persp(kc, theta=20, phi=35) ## mapas de probabilidades OC <- output.control(quantile=c(0.1, 0.5, 0.9)) kc <- krige.conv(TIBO3, loc=gr, krige=krige.control(obj.model=fit), output=OC) names(kc) par(mfrow=c(1,2)) PROBS <- apply(kc$simulation, 2, function(x) mean(x < 6.9)) hist(PROBS, prob=T, main="Distribuição preditiva na da proporção da área com P[S|y < 10]", panel.first=lines(density(PROBS))) P1 <- apply(kc$simulations, 1, function(x) mean(x < 6.9)) image(kc, val=P1,useRaster=T, col=gray(seq(1,0,l=21)),,main="Localizações na área com P[S|y < 10]") leg(-5,max(TIBO1$coords[,2]),zlim=range(PROBS),col=gray(seq(1,0,l=21)))