rm(list=ls(all=TRUE)) par.ori <- par(no.readonly=T) # ## Entrada de dados da Camada_1 # Quimicos<-read.table("Quimicos_C1.txt", header=TRUE) # ## Entrada de dados das bordas da area e sub-areas # borda<-read.table("Borda.txt", head=F) head(borda) A1<-read.table("Area1.txt", head=F) A2<-read.table("Area2.txt", head=F) A3<-read.table("Area3.txt", head=F) # ## Inspecionando os dados # head(Quimicos) names(Quimicos) # ## Convertendo a variavel regiao do tipo numerica para tipo fator # Quimicos$regiao <- as.factor(Quimicos$regiao) # summary(Quimicos) # require(geoR) # ## Convertendo os dados para o formato geodata # Mg<-as.geodata(Quimicos,coords.col=1:2,data.col=5,covar.col=12) Mg # par(mfrow=c(1,1)) points(Mg) polygon(borda) polygon(A1) polygon(A2) # Mg$borders <- borda # names(Mg) # ## Analise Exploratoria dos dados # summary(Mg) plot(Mg) plot(Mg,low=T) # ## Correlacao entre a variavel e a covariavel regiao # boxplot(Mg$data~Mg$covar$regiao, notch=T, varwidth=T) plot(Mg, trend=~regiao) # ## Investigando a necessidade de transformacao nos dados # require(MASS) # ## Resolvemos adotar lambda=1 # par(mfrow=c(1,2)) boxcox(Mg,lam = seq(0,2.5, l=100))-> BC1 boxcox(Mg, trend=~regiao, lam = seq(0,2.5, l=100))-> BC2 par(mfrow=c(1,1)) # plot(Mg, trend=~regiao) # ### Iniciando a modelagem # par(mar=c(3.5,3,0,0.5),mgp=c(2,0.8,0)) points(Mg, lambda=1,cex.min=.3,cex.max=1, axes=F, ann=F,col=gray(seq(0.9,0,l=length(Mg$data)))) polygon(A1) polygon(A2) polygon(A3) par(mar=c(3.5,3.5,0.5,0.5),mgp=c(2,0.8,0)) # par(mar=c(3.5,3,0,0.5),mgp=c(2,0.8,0)) points(Mg, trend=~regiao, lambda=1,cex.min=.3,cex.max=1, axes=F, ann=F,col=gray(seq(0.9,0,l=length(Mg$data)))) polygon(A1) polygon(A2) polygon(A3) par(mar=c(3.5,3.5,0.5,0.5),mgp=c(2,0.8,0)) # ## Explorando a dependencia espacial atraves do variograma # Mg.v<-variog(Mg,max.dist=350,lambda=1, trend=~regiao) plot(Mg.v, pch=19,cex=0.8) # ## Escolhendo o modelo # # Sem a covariavel # Mg.ml <- list() Mg.ml$Mg.1<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=0.5, lambda=1) Mg.ml$Mg.2<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=1.5, lambda=1) Mg.ml$Mg.3<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=2.5, lambda=1) Mg.ml$Mg.4<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=3.5, lambda=1) Mg.ml$Mg.5<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=4.5, lambda=1) # # Com a covariavel regiao # Mg.ml1 <- list() Mg.ml1$Mg.1<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=0.5, lambda=1, trend=~regiao) Mg.ml1$Mg.2<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=1.5, lambda=1, trend=~regiao) Mg.ml1$Mg.3<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=2.5, lambda=1, trend=~regiao) Mg.ml1$Mg.4<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=3.5, lambda=1, trend=~regiao) Mg.ml1$Mg.5<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=4.5, lambda=1, trend=~regiao) # ## Modelo escolhido e com a covariavel e kappa=1.5 # sapply(Mg.ml, logLik) sapply(Mg.ml1, logLik) # ## Parametros estimados # Mg.ml1$Mg.2<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=1.5, lambda=1, trend=~regiao) Mg.ml1$Mg.2 # ## Iniciando a interpolacao # # Definindo o grid de predicao # Mg$borders <- borda gr <- pred_grid(Mg$borders, by=10) gr0<- polygrid(gr, borders=Mg$borders, bound=T) ind.reg<-numeric(nrow(gr0)) ind.reg[.geoR_inout(gr0,A1)]<-1 ind.reg[.geoR_inout(gr0,A2)]<-2 ind.reg[.geoR_inout(gr0,A3)]<-3 ind.reg<-as.factor(ind.reg) # set.seed(512) KC1<-krige.control(trend.d=~regiao, trend.l=~ind.reg, obj.model=Mg.ml1$Mg.2) # ## Especificando e mapas de simulacao # OC=output.control(n.pred=3) # ## Interpolando # pred<-krige.conv(Mg, loc=gr, krige=KC1,output=OC) # ## Construindo os mapas # par(mfrow=c(2,2), mar=c(1,2,1,1)) image(pred, col=gray(seq(0.95,0.1,l=51))) image(pred, val=pred$simul[,1],col=gray(seq(0.95,0.1,l=51))) image(pred, val=pred$simul[,2],col=gray(seq(0.95,0.1,l=51))) image(pred, val=pred$simul[,3],col=gray(seq(0.95,0.1,l=51))) par(mfrow=c(1,1))