rm(list=ls(all=TRUE)) par.ori <- par(no.readonly=T) # ## Entrada de dados da Camada 2 # qc2 <- read.table("Quimicos_C2.txt", head=T, dec=",") # ## Entrada de dados das bordas da área e sub-áreas # borda <- read.table("Borda.txt") head(borda) a1 <- read.table("Area1.txt") a2 <- read.table("Area2.txt") a3 <- read.table("Area3.txt") # ## Inspecionando os dados # head(qc2) names(qc2) # ## Convertendo a variavel Regiao do tipo numerica para tipo fator # qc2$Regiao <- as.factor(qc2$Regiao) # summary(qc2) # ## Obtendo os valores minimos e maximos para as coordenadas # apply(qc2[,1:2], 2 , range) # ## Ajustando os valores das coordenadas # qc2$X <- (qc2$X-688000) qc2$Y <- (qc2$Y-7100000) qc2$X qc2$Y # require(geoR) # ## Convertendo os dados para o formato geodata # pha <- as.geodata(qc2, coords.col=1:2, data.col=6,covar.col=13:14) pha # par(mfrow=c(1,1)) points(pha) polygon(borda) polygon(a1) polygon(a2) # pha$borders <- borda # names(pha) # ## Excluindo outliers # plot(pha) phaog1 <- subset(pha, pha$coords[,1] != 571) plot(phaog1) pha1 <- subset(phaog1,phaog1$coords[,1] != 529) # ## Analises Exploratorias plot(pha1) plot(pha1,low=T) # ## Correlacao entre a variavel sem os outliers e as covariaveis profundidade e regiao # par(mfrow=c(1,1)) plot(pha1$covar$Prof,pha1$data) lines(lowess(pha1$data~pha1$covar$Prof)) # boxplot(pha1$data~pha1$covar$Regiao, notch=T,varwidth=T) plot(pha1, trend=~Regiao) # ## Investigando a necessidade de transformacao nos dados # require(MASS) # ## Resolvemos adotar lambda=1, sem transformacao # boxcox(pha1,lam=seq(-4,4,l=100)) boxcox(pha1,trend=~Regiao,lam=seq(-4,4,l=100)) boxcox(pha1,trend=~Prof,lam=seq(-4,4,l=100)) # ### Iniciando a modelagem ### # par(par.ori) par(mfrow=c(1,2)) points(pha1) points(pha1,cex.min=0.1,cex.max=2) # ## Delimitando as Regioes # polygon(a1) polygon(a2) # ## Explorando a dependencia espacial atraves do variograma # summary(pha1) v0 <- variog(pha1, max.dist=250) v1 <- variog(pha1,max.dis=250,trend=~Regiao) v2 <- variog(pha1,max.dis=250,trend=~Prof) par(mfrow=c(1,3)) plot(v0) plot(v1) plot(v2) # # Garantindo a mesma escala para y em ambos os graficos # vmax <- max(c(v0$v,v1$v,v2$v)) plot(v0,ylim=c(0,vmax)) plot(v1,ylim=c(0,vmax)) plot(v2,ylim=c(0,vmax)) # ## Escolhendo o modelo # # Sem as covariaveis # pha1.msc <- list() inis <- cbind(cbind(seq(0.02,0.07,l=6),seq(50,150,l=6))) pha1.msc$pha11 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=0.5) pha1.msc$pha12 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=1.5) pha1.msc$pha13 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=2.5) pha1.msc$pha14 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=3.5) pha1.msc$pha15 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=4.5) # # Com a covariavel Regiao # pha1.mccr <- list() inis <- cbind(cbind(seq(0.02,0.07,l=6),seq(50,150,l=6))) pha1.mccr$pha11 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=0.5,trend=~Regiao) pha1.mccr$pha12 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=1.5,trend=~Regiao) pha1.mccr$pha13 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=2.5,trend=~Regiao) pha1.mccr$pha14 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=3.5,trend=~Regiao) pha1.mccr$pha15 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=4.5,trend=~Regiao) # # Com a covariavel Profundidade # pha1.mccp <- list() inis <- cbind(cbind(seq(0.02,0.07,l=6),seq(50,150,l=6))) pha1.mccp$pha11 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=0.5,trend=~Prof) pha1.mccp$pha12 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=1.5,trend=~Prof) pha1.mccp$pha13 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=2.5,trend=~Prof) pha1.mccp$pha14 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=3.5,trend=~Prof) pha1.mccp$pha15 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=4.5,trend=~Prof) # # Com as covariaveis REgiao e Profundidade pha1.mccrp <- list() inis <- cbind(cbind(seq(0.02,0.07,l=6),seq(50,150,l=6))) pha1.mccrp$pha11 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=0.5,trend=~Regiao+Prof) pha1.mccrp$pha12 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=1.5,trend=~Regiao+Prof) pha1.mccrp$pha13 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=2.5,trend=~Regiao+Prof) pha1.mccrp$pha14 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=3.5,trend=~Regiao+Prof) pha1.mccrp$pha15 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=4.5,trend=~Regiao+Prof) # ## Modelo escolhido e com covariavel regiao e kappa=1.5: # sapply(pha1.msc,logLik) sapply(pha1.mccr,logLik) sapply(pha1.mccp,logLik) sapply(pha1.mccrp,logLik) # ## Parametros estimados # inis <- cbind(cbind(seq(0.02,0.07,l=6),seq(50,150,l=6))) pha1.mccr$pha12 <- likfit(pha1,ini=inis,nug=seq(0.01,0.04,l=5),cov.model='matern',kappa=1.5,trend=~Regiao) pha1.mccr$pha12 # ## Iniciando a interpolacao # par(par.ori) par(mfrow=c(1,1)) points(pha1) # # Definindo o grid de predicao # gr <- pred_grid(borda,by=10) dim(gr) points(gr,pch=19, cex=0.25) gr0 <- polygrid(gr,border=borda,bound=T) dim(gr0) points(gr0,pch=19,cex=0.25,col=1) # 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 table(ind.reg) ind.reg <- as.factor(ind.reg) # set.seed(512) KC <- krige.control(trend.d=~Regiao,trend.l=~ind.reg,obj.model=pha1.mccr$pha12) # # Especificando 3 mapas de simulacao # OC <- output.control(n.pred=3) # # Interpolando # pha1 pha1.kc <- krige.conv(pha,loc=gr,krige=KC,output=OC) # ## Construindo os mapas # par(mfrow=c(2,2),mar=c(1,2,1,1)) image(pha1.kc,col=gray(seq(0.95,0.1,l=51))) image(pha1.kc,val=pha1.kc$simul[,1],col=gray(seq(0.95,0.1,l=51))) image(pha1.kc,val=pha1.kc$simul[,2],col=gray(seq(0.95,0.1,l=51))) image(pha1.kc,val=pha1.kc$simul[,3],col=gray(seq(0.95,0.1,l=51))) par(mfrow=c(1,1)) #