require(geoR) require(MASS) Quimicos<-read.table("Quimicos_C1.txt", header=TRUE) borda<-read.table("Borda.txt", head=F) A1<-read.table("Area1.txt", head=F) A2<-read.table("Area2.txt", head=F) A3<-read.table("Area3.txt", head=F) Quimicos$regiao <- as.factor(Quimicos$regiao) Alum<-as.geodata(Quimicos,coords.col=1:2,data.col=3,covar.col=12) Alum$borders <- borda Al4<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=3.5, lambda=0, trend=~regiao) gr <- pred_grid(Alum$borders, by=10) gr0<- polygrid(gr, borders=Alum$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=Al4) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(Alum, loc=gr, krige=KC1,output=OC) pred1 <- pred$simulations[,1] pred1 Ca<-as.geodata(Quimicos,coords.col=1:2,data.col=4,covar.col=12) Ca$borders <- borda Ca2<-likfit(Ca, ini=c(1.2,75),cov.model="matern", kappa=1.5, lambda=1, trend=~regiao) gr <- pred_grid(Ca$borders, by=10) gr0<- polygrid(gr, borders=Ca$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=Ca2) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(Ca, loc=gr, krige=KC1,output=OC) pred2 <- pred$simulations[,1] pred2 Mg<-as.geodata(Quimicos,coords.col=1:2,data.col=5,covar.col=12) Mg$borders <- borda Mg2<-likfit(Mg, ini=c(1.2,75),cov.model="matern", kappa=1.5, lambda=1, trend=~regiao) 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=Mg2) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(Mg, loc=gr, krige=KC1,output=OC) pred3 <- pred$simulations[,1] pred3 PhA<-as.geodata(Quimicos,coords.col=1:2,data.col=6,covar.col=12) PhA$borders <- borda PhA1 <- subset(PhA,PhA$coords[,1] !=571) PhA5<-likfit(PhA1, ini=c(0.0006,50),cov.model="matern", kappa=4.5, lambda=0, trend=~regiao) gr <- pred_grid(PhA$borders, by=10) gr0<- polygrid(gr, borders=PhA$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=PhA5) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(PhA, loc=gr, krige=KC1,output=OC) pred4 <- pred$simulations[,1] pred4 PhK<-as.geodata(Quimicos,coords.col=1:2,data.col=7,covar.col=12) PhK$borders <- borda PhK5<-likfit(PhK, ini=c(0.03,100),cov.model="matern", kappa=4.5, lambda=1, trend=~regiao) gr <- pred_grid(PhK$borders, by=10) gr0<- polygrid(gr, borders=PhK$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=PhK5) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(PhK, loc=gr, krige=KC1,output=OC) pred5 <- pred$simulations[,1] pred5 K<-as.geodata(Quimicos,coords.col=1:2,data.col=8,covar.col=12) K$borders <- borda K1<-likfit(K, ini=c(0.05,40),cov.model="matern", kappa=0.5, lambda=0, trend=~regiao) gr <- pred_grid(K$borders, by=10) gr0<- polygrid(gr, borders=K$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=K1) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(K, loc=gr, krige=KC1,output=OC) pred6 <- pred$simulations[,1] pred6 P<-as.geodata(Quimicos,coords.col=1:2,data.col=9,covar.col=12) P$borders <- borda P2<-likfit(P, ini=c(0.2,100),cov.model="matern", kappa=1.5, lambda=0, trend=~regiao) gr <- pred_grid(P$borders, by=10) gr0<- polygrid(gr, borders=P$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=P2) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(P, loc=gr, krige=KC1,output=OC) pred7 <- pred$simulations[,1] pred7 C<-as.geodata(Quimicos,coords.col=1:2,data.col=10,covar.col=12) C$borders <- borda C1 <- subset(C,C$coords[,1] != 775) Car1<-likfit(C1, ini=c(1,20),cov.model="matern", kappa=0.5, lambda=1, trend=~regiao) gr <- pred_grid(C$borders, by=10) gr0<- polygrid(gr, borders=C$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=Car1) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(C, loc=gr, krige=KC1,output=OC) pred8 <- pred$simulations[,1] pred8 H<-as.geodata(Quimicos,coords.col=1:2,data.col=11,covar.col=12) H$borders <- borda H5<-likfit(H, ini=c(0.015,40),cov.model="matern", kappa=4.5, lambda=0, trend=~regiao) gr <- pred_grid(H$borders, by=10) gr0<- polygrid(gr, borders=H$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=H5) OC=output.control(n.pred=3) # define 3 simulacoes pred<-krige.conv(H, loc=gr, krige=KC1,output=OC) pred9 <- pred$simulations[,1] pred9 # Unindo as 9 colunas de predicao predS <-as.table(round(cbind(pred1,pred2,pred3,pred4,pred5,pred6,pred7,pred8,pred9),4)) predS # Resumo rapido das nove colunas de predicao head(predS) # guardando estes resultados em um arquivo txt para posterior analise write.table(predS,file='predS.txt')