par.ori <- par(no.readonly=TRUE) require(geoComp) require(compositions) ## carregando o arquivo de dados do geoComp data(pivo) ## criando as coordenadas da borda bor <- cbind(c(0,seq(0, 200, l=100), 0), c(0,sqrt(200^2-seq(0, 200, l=100)^2), 0)) ## visualizando #postscript("pivo.eps", horizontal=F,width=3.5,height=3.5) jpeg("pivo.jpg")#,mar=c(3.5,3.5,0.5,0.5),mgp=c(1.8,0.8,0),oma=c(10,10,10,10) par(mfrow=c(1,1),mar=c(3,3,0.5,0.5),mgp=c(1.8,0.8,0)) plot(bor,asp=1,ty="n",xlim=c(0,200),ylim=c(0,200),xlab="",ylab="") points(pivo[,1:2],asp=1) polygon(bor) dev.off() head(pivo) ## selecionando os componentes e as coordenadas dados <- pivo[,c(6,7,8,1,2)] head(dados) ## visualizando #postscript("hists.eps", height=240, width=240) postscript("histDig.eps", height=6, width=6) par(mfrow=c(2,2), mar=c(3.5, 3.5, 0.5, 0.5), mgp=c(1.8,0.8,0)) hist(dados[,1], main="", xlab="Areia",ylab="Frequencia",ylim=c(0,25)) hist(dados[,2], main="", xlab="Silte",ylab="Frequencia",ylim=c(0,25)) hist(dados[,3], main="", xlab="Argila",ylab="Frequencia",ylim=c(0,25)) #na classe acomp os dados sao analisados na geometria relativa comp <- acomp(dados[,1:3]) names(comp) <- c("Areia", "Silte", "Argila") #construindo o diagrama ternario com regioes de confianca de 2 e 4 sigma plot(comp) plot(mean(comp),add=T,pch=20,col='red') ellipses(mean(comp),var(comp),col='red',r=2) ellipses(mean(comp),var(comp),col='red',r=4) dev.off() ## Convertendo o arquivo para o formato geoComp ## ALR com argila como denominador source("as.geoComp.R") dados <- as.geoComp(dados) #str(dados) jpeg("corr.jpg")#,width=6,height=2.5 #, mar=c(3.5, 3.5, 0.5, 0.5), mgp=c(1.8,0.8,0) par(mfrow=c(1,3)) hist(dados[[1]]$Y1, main="",ylab="Frequencia", xlab="ln(Areia/Argila)",ylim=c(0,20)) hist(dados[[1]]$Y2, main="",ylab="Frequencia", xlab="ln(Silte/Argila)") plot(dados[[1]]$Y1, dados[[1]]$Y2, xlab="ln(Areia/Argila)",ylab="ln(Silte/Argila)") dev.off() ## Estimacao do modelo geoestatistico bivariado composicional ## reparametrizacao/vero concentrada/hessiano/metodo delta estima <- mec(dados) estima ## Fazendo a "cokrigagem"/predicao bivariada na escala normal/MVnorm #md.cov.ck <- cokrigagem(estima[[1]]$Estimativas, dim.gride=15, dados) require(geoR) require(statmod) gr <- pred_grid(bor, by=4) source("cokrigagem.R") md.cov.ck <- cokrigagem(estima[[1]]$Estimativas, loc=gr, dados.comp=dados) str(md.cov.ck) ## Voltando por quadratura de Gauss Hermite preditos.gh <- volta.quad(md.cov.ck,n.pontos=6,Variancia=FALSE) head(preditos.gh[[1]]) ## Voltando por Simulacao preditos.simu <- volta.cokri(md.cov.ck,num.simu=50,int.conf=0.90) head(preditos.simu[[1]]) ## comparando quadratura e simulacao par(mfrow=c(1,3), mar=c(3.5, 3.5, 0.5, 0.5), mgp=c(1.8,0.8,0)) plot(preditos.gh[[1]][,1], preditos.simu[[1]][,1], xlab="Gauss-Hermite",ylab="Simulacao", main="\nAreia", asp=1); abline(0,1) plot(preditos.gh[[1]][,2], preditos.simu[[1]][,2], xlab="Gauss-Hermite",ylab="Simulacao", main="\nSilte", asp=1); abline(0,1) plot(preditos.gh[[1]][,3], preditos.simu[[1]][,3], xlab="Gauss-Hermite",ylab="Simulacao", main="\nArgila", asp=1); abline(0,1) par(mfrow=c(2,3), mar=c(3.5, 3.5, 0.5, 0.5), mgp=c(1.8,0.8,0)) image(structure(list(predict=preditos.gh[[1]][,1]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(100, 200), y.leg=c(170, 180)) image(structure(list(predict=preditos.gh[[1]][,2]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(100, 200), y.leg=c(170, 180)) image(structure(list(predict=preditos.gh[[1]][,3]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(100, 200), y.leg=c(170, 180)) image(structure(list(predict=preditos.simu[[1]][,1]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(100, 200), y.leg=c(170, 180)) image(structure(list(predict=preditos.simu[[1]][,2]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(100, 200), y.leg=c(170, 180)) image(structure(list(predict=preditos.simu[[1]][,3]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(100, 200), y.leg=c(170, 180)) ######################################################### ###### Criando os mapas com base no .txt do cluster predi.gh <- read.table("predGH4s200p7.txt", header=T) predi.simu <- read.table("predsimu4s200p7.txt",header=T) jpeg("mapaghsimu.jpg") par(mfrow=c(2,3), mar=c(12.5, 4.5, 0.1, 0.5), mgp=c(2,0.8,0)) # por quadratura image(structure(list(predict=predi.gh[,1]), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(125, 200), y.leg=c(170, 180),ylim=c(0,220)) text(x=145,y=200,"(A) Areia - QG") image(structure(list(predict=predi.gh[,2]), class="kriging"), loc=gr,bor, col=terrain.colors(21),main="Silte", x.leg=c(125, 200), y.leg=c(170, 180),ylim=c(0,220)) text(x=145,y=200,"(B) Silte - QG") image(structure(list(predict=predi.gh[,3]), class="kriging"), loc=gr,bor, col=terrain.colors(21),main="Argila", x.leg=c(125, 200), y.leg=c(170, 180),ylim=c(0,220)) text(x=145,y=200,"(C) Argila - QG") # por simulacao image(structure(list(predict=predi.simu$Areia), class="kriging"), loc=gr,bor, col=terrain.colors(21), x.leg=c(125, 200), y.leg=c(170, 180),ylim=c(0,210)) text(x=170,y=200,"(D) Areia - Sim") image(structure(list(predict=predi.simu$Silte), class="kriging"), loc=gr,bor, col=terrain.colors(21),main="Silte", x.leg=c(125, 200), y.leg=c(170, 180),ylim=c(0,210)) text(x=170,y=200,"(E) Silte - Sim") image(structure(list(predict=predi.simu$Argila), class="kriging"), loc=gr,bor, col=terrain.colors(21),main="Argila", x.leg=c(125, 200), y.leg=c(170, 180),ylim=c(0,210)) text(x=170,y=200,"(F) Argila - Sim") dev.off() # comparando quadratura e simulacao jpeg("gh_x_simu.jpg") par(mfrow=c(1,3), mar=c(3.5, 3.5, 0.5, 0.5), mgp=c(1.8,0.8,0)) plot(predi.gh$X1, predi.simu$Areia, xlab="Gauss-Hermite",ylab="Simulacao", main="\nAreia", asp=1); abline(0,1) plot(predi.gh$X2, predi.simu$Silte, xlab="Gauss-Hermite",ylab="Simulacao", main="\nSilte", asp=1); abline(0,1) plot(predi.gh$X3, predi.simu$Argila, xlab="Gauss-Hermite",ylab="Simulacao", main="\nArgila", asp=1); abline(0,1) dev.off()