## Analise bivariada de IMA (Azul) e ARGILA (Vermelho) # ## Limpando a workspace rm(list=ls(all=TRUE)) # # Leitura do arquivo de dados # IMA.csv <- read.table("/media/usbdisk/ITAMAR/Azul/IMA_Arg_azul.csv", head=T, dec=",", sep="",col.names=c("X","Y","IMA","ARG1","Arg2")) IMA.csv <- read.table("http://wiki.leg.ufpr.br/data/media/projetos/mobasa/artigo/ima_arg_azul.csv", head=T, sep="", dec=",",col.names=c("X","Y","IMA","ARG1","Arg2")) head(IMA.csv) # # arg <- read.table("/media/usbdisk/ITAMAR/Vermelho/arg_vermelho.csv", head=T, dec=",", sep="") arg <- read.table("http://wiki.leg.ufpr.br/data/media/projetos/mobasa/artigo/arg_vermelho.csv", head=T, dec=",", sep="") head(arg) # # Lendo o arquivo de bordas borda <- read.table("http://wiki.leg.ufpr.br/data/media/projetos/mobasa/artigo/contorno.txt", head=T, dec=".", sep="") head(borda) # ## Avaliando a correlacao entre as variaveis # ## Construcao do scatterplot jpeg("IMA_Arg_scatter.jpg", quality=100) par(mar=c(3,3,0.5,0.5),mgp=c(2,.8,0)) plot(IMA.csv$IMA,IMA.csv$ARG1,pch=20, xlab='IMA',ylab='Argila') # abline(lm(IMA.csv$IMA~IMA.csv$ARG1)) dev.off() # cor.test(IMA.csv$IMA,IMA.csv$ARG1) # # carregando o pacote geoR require(geoR) # # Gerando o arquivo geodata IMA <- as.geodata(IMA.csv,coords.col=1:2, data.col=3,covar.col=4) Arg <- as.geodata(arg,coords.col=1:2, data.col=3) names(IMA) summary(IMA) summary(Arg) ## # Examinando os dados plot(Arg,borders=borda ) plot(IMA,borders=borda) # Verificando a necessidade de transformação require(MASS) boxcox(Arg) boxcox(IMA) # Construindo variogramas para a Argila (vermelho) Arg.v <- variog(Arg, max.dist=4000,lambda=0) plot(Arg.v) # Construindo variogramas para IMA IMA.v <- variog(IMA, max.dist=4000) plot(IMA.v) # Ajustando um modelo exponencial por máxima verossimilhança # Sem co-variável # Sys.time() -> inicio IMA.lik <- likfit(IMA, ini=c(50,2007),cov.model="matern", kappa=0.5) # Sys.time() -> final # (duracao <- final-inicio) # # Usando a argila como co-variavel IMAeARG.lik <- likfit(IMA, ini=c(50,2007),cov.model="matern", kappa=0.5, trend=~ARG1) IMAeARG.lik # IMA.lik Arg.lik <- likfit(Arg, ini=c(.1,2007),cov.model="matern",kappa=0.5,lambda=0) Arg.lik # Construindo um grid de predição gr <- pred_grid(borda,by=100) # Estabelecendo os controles de krigagem kc.IMA <- krige.control(obj.model=IMA.lik) kc.Arg <- krige.control(obj.model=Arg.lik,lambda=0) # Krigando ......... Arg.k <- krige.conv(Arg,loc=gr,krige=kc.Arg) IMA.k <- krige.conv(IMA,loc=gr,krige=kc.IMA) IMA.k # Gerando a imagem com krigagem par(mfrow=c(1,2)) image(IMA.k,loc=gr,col=gray((0:32)/32), borders=borda) image(Arg.k,loc=gr,col=gray((0:32)/32), borders=borda) par(mfrow=c(1,1)) # ********************************************* # AQUI COMECA A NOVIDADE........ # ********************************************* # names(Arg) inicio <- Sys.time() fit12 <- likfitBGCCM(IMA,Arg) final <- Sys.time() (final-inicio) class(fit12) ## Co-Krigando em HMED pred1 <- predict(fit12, loc=gr) names(pred1) mean(pred1$krige.var) # Gerando a imagem jpeg("IMA.jpg", quality=100) par(mar=c(3,3,0.5,0.5),mgp=c(2,.8,0)) image(Arg.k,loc=gr,col=gray((0:32)/32), borders=borda) dev.off() head(pred1$pred) ## na mesma escala jpeg("IMA_Cok.jpg",height=480,width=960, quality=100) # par(mfrow=c(1,2),mar=c(3.5,3.5,0.5,0.5),mgp=c(2,.8,0),cex=.5 par(mfrow=c(1,2),mar=c(3,3,0.5,0.5),mgp=c(2,.8,0)) zl <- range(c(IMA.k$pred, pred1$pred)) image(IMA.k,loc=gr,col=gray(seq(1,0.1,l=21)), borders=borda, zlim=zl,xlab="",ylab="", x.leg=c(638000, 643000), y.leg=c(7054000,7054400)) text(634000, 7054500, "(A)", cex=1.5) image(pred1,loc=gr,col=gray(seq(1,0.1,l=21)), borders=borda, zlim=zl,xlab="",ylab="",, x.leg=c(638000, 643000), y.leg=c(7054000,7054400)) text(634000, 7054500, "(B)", cex=1.5) dev.off() ## Construindo o diagrama de dispersão para correlacao ## entre os valores preditos pela krigagem convencional (modelo de 18 pontos) ## e os valores preditos pela co-krigagem considerando a variável argila (555 pontos) ## (Análise do Paulo) jpeg("IMA_Cok_scater.jpg", quality=100) par(mar=c(3.5,3.5,0.5,0.5),mgp=c(2,.8,0)) plot(IMA.k$pred, pred1$pred, xlim=zl, ylim=zl, asp=1, pch=20, cex=0.25) abline(c(0,1)) dev.off() ## examinar futuramente # # zlep <- range(c(HMED.k[[2]], pred1[[2]])) # plot(sqrt(HMED.k[[2]]), sqrt(pred1[[2]]), xlim=zlep, ylim=zlep, asp=1, pch=20, cex=0.25) #plot(sqrt(HMED.k[[2]]), sqrt(pred1[[2]]), asp=1, pch=20, cex=0.25) # abline(c(0,1))