options(width=90) ## carregando pacotes necessários require(geoR) require(MASS) ## importando os dados de um arquivo formato texto dados <- read.table("profund020a.txt", head=T) head(dados) ## convertendo os dados de fósforo para o formato geodata ## considerando tipo de solo como covariável P <- as.geodata(dados, coords.col = 1:2, data.col=9, covar.col=5) summary(P) ## gráfico de pontos na área, tamanhos proporcionais aos teores de fósforo points(P) ## definindo uma borda (arbitrária) na area clicando em pontos bor <- locator(type="p", pch=21) polygon(bor) P$borders <- with(bor, cbind(x,y)) points(P) ## uma visualizacao dos dados plot(P) ## examinando relação da variável resposta (P) com a covariável with(P, boxplot(data~covariate$solo, varwidth=T)) ## examinando a necessidade de transformação boxcox(P) abline(v=c(0, 0.25, .5, 1), col=4) with(P, boxcox(data ~ covariate$solo)) ## removendo os valores muitos altos P1 <- subset(P, data <=70) summary(P1) boxcox(P1) ## visualizando os pontos excluídos points(P) points(P$coords[P$data > 70, ],pch=19, col=2) ## variogramas dados originais plot(variog(P, max.dist=2000, lam=0.5)) ## variogramas dados excluindo valores altos (note a erraticidade!!!!) plot(variog(P1, lam=0.5)) plot(variog(P1, max.dist=2000, lam=0.5)) plot(variog(P1, max.dist=1500, lam=0.5)) plot(variog(P1, max.dist=500, lam=0.5)) plot(variog(P1, max.dist=1000, lam=0.5)) plot(variog(P1, max.dist=750, lam=0.5)) plot(variog(P1, uvec=seq(0,750, by=50), lam=0.5)) plot(variog(P1, uvec=seq(0,500, by=50), lam=0.5)) plot(variog(P1, max.dist=1000, lam=0.5)) ## variograma considerando covariável plot(variog(P1, max.dist=1000, lam=0.5, trend=~solo)) ## estimação de parâmetros (por máxima verossimilhança) ## valores iniciais cf sugerido visualmente pelo variograma ml <- likfit(P1, ini=c(10, 30), nug=3, lam=0.5) ml ## vendo o função do modelo ajustado por máxima verossimilhança( MV) plot(variog(P1, max.dist=1000, lam=0.5)) lines.variomodel(ml) plot(variog(P1, max.dist=250, lam=0.5)) lines.variomodel(ml) ## avaliando o ajunte de MV summary(ml) ### outros dados ... argila boxcox(asin(dados$argila/100)~1) boxcox(asin(dados$argila/100)~1, lam=seq(-4, 2, length=100)) with(dados, boxcox(log((argila/100)/(1-argila/100))+5~1)) hist(dados$argila, prob=T) hist(log(dados$argila), prob=T) lines(density(log(dados$argila))) rug(log(dados$argila)) dados$argila summary(dados$argila) table(dados$argila) ### outros dados... Magnésio head(dados) ## vendo as opções da função que converte dados para formato geodata args(as.geodata.default) help(as.geodata) ## convertendo e resumindo os dados Mg <- as.geodata(dados, data.col=13, covar.col=5) summary(Mg) ## atribuindo as mesmas bordar selecionadas anteriormente Mg$borders <- P$borders ## gráficos exploratórios plot(Mg, low=T) ## vendo os tipos de solo points(Mg) points(Mg$coords[Mg$cov$solo=="AQ",], col=4, pch=19) points(Mg$coords[Mg$cov$solo=="LVa11",], col=2, pch=19) points(Mg$coords[Mg$cov$solo=="LVa12",], col=5, pch=19) ## verificando se há relação com covariável with(dados, plot(Magnésio ~ solo)) ## verificando tranformação boxcox(Mg) plot(Mg, low=T, lam=0) ## variogramas plot(variog(Mg, max.dist=1500, lam=0)) plot(variog(Mg, max.dist=800, lam=0)) plot(variog(Mg, uvec=seq(0,700, l=10), lam=0)) plot(variog(Mg, uvec=seq(0,700, l=20), lam=0)) plot(variog(Mg, uvec=seq(0,900, l=15), lam=0)) ## estimação de parâmetros via MV Mg.ml <- likfit(Mg, ini=c(0.15, 200), nug=0.1, lam=0) Mg.ml ## resumo do modelo ajustado summary(Mg.ml) ## definindo um grid para rpedição espacial GR <- pred_grid(Mg$borders, by=30) points(Mg) points(GR, pch=".") ## predição espacial: krigagem Mg.kr <- krige.conv(Mg, loc=GR, bor=Mg$bor, krige=krige.control(obj.model=Mg.ml)) image(Mg.kr) ## algumas opções de visualização image(Mg.kr, col=gray(seq(1, 0, len=100))) image(Mg.kr, col=gray(seq(1, 0, len=10))) image(Mg.kr, col=gray(seq(1, 0, len=5))) image(Mg.kr, col=gray(seq(1, 0, len=3)))