############################################################## ### Script - Análise Resistencia do solo ### a penetração ### ### ### ### Naimara Vieira do Prado ### ############################################################## #Carregando o m?dulo geoR require(geoR) # Carrega o pacote estat?sitico require(stats) require(e1071) ############################################################## ### Inserindo o arquivo de dados ### ############################################################## # Lendo o arquivo de dados dados<-read.geodata("prova.txt",head=F,coords.col=2:3,data.col=4) dados ################################################################ ### An?lise explorat?ria de dados ### ################################################################ # Calcula as Estat?sticas descritivas dos dados # m?dia, mediana, Q1, Q3, min, m?x summary(dados$data) mean(dados$data) # Apresenta a M?dia var(dados$data) # Apresenta a Vari?ncia sd(dados$data) # Apresenta o Desvio padr?o max(dados$data) max(dist(dados$coords)) min(dados$data) min(dist(dados$coords)) skewness(dados$data) #Assimetria kurtosis(dados$data) #Curtose # Calcular o Coeficiente de varia??o # CV = (DP/Media)x100 CV <- sd(dados$data)*100/mean(dados$data) CV # mostra o coeficiente de varia??o #Apresenta o gr?fico Box-plot boxplot(dados$data,main='Gr?fico Boxplot') #Apresenta o histograma hist(dados$data, breaks=5, main='Histograma') #Apresenta v?rios gr?ficos dos dados plot(dados) ########################################################### ### An?lise explorat?ria espacial ### ########################################################### # Gr?fico Post-plot para o estudo de tend?ncia direcional #Post Plot points(dados,xlab="", ylab="",main='RSP',font.main = 3,pt.div="quartile",col=gray(seq(0.8,0,l=5))) legend(40600,38050,c("M?n-Q1","Q1-Q2","Q2-Q3","Q3-M?x"),fill=c("#FFFFFF","#B2B2B2","#4C4C4C","#000000")) ######################################################## ### Constru??o dos semivariogramas direcionais ### ######################################################## max(dist(dados$coords)) #m?xima dist?ncia entre pontos min(dist(dados$coords)) #m?nima dist?ncia entre pontos plot(variog4(dados, coords = dados$coords, data = dados$data, ylims=c(0,5), uvec =seq(1,668.5,l=10), estimator.type = "classical", max.dist=668.5, pairs.min = 30, direction = c(0, pi/4, pi/2, 3*pi/4), tolerance = pi/8)) ########################################################### ### Constru??o do semivariograma ### ########################################################### d.var <- variog(dados,uvec=seq(1,668.5,l=10),estimador.type="classical", pairs.min=30, direction="omnidirectional",tolerance=pi/8) plot(d.var, xlab="dist?ncia", ylab="semivari?ncia",main='RSP',font.main = 3) d.var #obtendo informa??es sobre o semivariograma experimental dist?ncia<-d.var$u semivari?ncia<-d.var$v pares<-d.var$n autocorr<-d.var$sd tabela<-cbind(dist?ncia,semivari?ncia,autocorr,pares) tabela ########################################################### ### Gr?fico de envelopes ### ########################################################### ##envelopes (intervalos de confian?a dentro do semivariograma - os pontos fora do limite indica depend?ncia espacial) #.Random.seed <- save.seed #save.seed <- .Random.seed set.seed(02003) dados.env<-variog.mc.env(dados, obj.v=d.var,nsim = 99) plot(d.var,env=dados.env,xlab="dist?ncia", ylab="semivari?ncia",main='RSP',font.main = 3) ########################################################### ### Ajuste de modelos espaciais ### ########################################################### #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #%%% AJUSTE OLS %%% #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ######### M?todo M?nimos Quadrados - OLS ################### Chutar a contribui??o (f1+f2) e f3 ############################ ## Modelo exponencial exp.ols<-variofit(d.var,ini=c(0.1,100),weights="equal", cov.model="exp") exp.ols plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma OLS - Exponencial') lines(exp.ols, col="blue") summary(exp.ols) ## Modelo esf?rico sph.ols<-variofit(d.var,ini=c(0.1,250),weights="equal", cov.model="sph") sph.ols plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma OLS - Esf?rico') lines(sph.ols, col="blue") summary(sph.ols) ## Modelo gaussiano gaus.ols<-variofit(d.var,ini=c(0.1,150),weights="equal", cov.model="gaus") gaus.ols plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma OLS - Gaussiano') lines(gaus.ols, col="blue") summary(gaus.ols) ## Modelo Fam?lia Mat?rn k=0.7 matern.ols<-variofit(d.var,ini=c(0.1,100),weights="equal", cov.model="matern", k=0.7) matern.ols plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma OLS - F.Matern k=0.7') lines(matern.ols, col="blue") summary(matern.ols) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% AJUSTE WLS1 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ######### M?todo M?nimos Quadrados Ponderados - WLS1 ## Modelo exponencial exp.wls1<-variofit(d.var,ini=c(0.1,250), cov.model="exp") exp.wls1 plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - WLS1 - Exponencial') lines(exp.wls1, col="blue") summary(exp.wls1) ## Modelo esf?rico sph.wls1<-variofit(d.var,ini=c(0.1,300), cov.model="sph") sph.wls1 plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - WLS1 - Esf?rico') lines(sph.wls1, col="blue") summary(sph.wls1) ## Modelo gaussiano gaus.wls1<-variofit(d.var,ini=c(0.1,250), cov.model="gaus") gaus.wls1 plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - WLS1 - Gaussiano') lines(gaus.ols, col="blue") summary(gaus.wls1) ## Modelo Fam?lia Mat?rn k=0.7 matern.wls1<-variofit(d.var,ini=c(0.1,150), cov.model="matern", k=0.7) matern.wls1 plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - WLS1 - F.Matern k=0.7') lines(matern.wls1, col="blue") summary(matern.wls1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% AJUSTE ML %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ######### M?todo M?xima Verossimilhan?a - ML ## Modelo exponencial exp.ml<-likfit(dados,ini=c(0.1,100), method="ml", cov.model="exp") exp.ml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - ML - Exponencial') lines(exp.ml, col="blue") summary(exp.ml) ## Modelo esf?rico sph.ml<-likfit(dados,ini=c(0.1,250), lik.method="ml", cov.model="sph") sph.ml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - ML - Esf?rico') lines(sph.ml, col="blue") summary(sph.ml) ## Modelo gaussiano gaus.ml<-likfit(dados,ini=c(0.1,150), lik.method="ml", cov.model="gaus") gaus.ml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - ML - Gaussiano') lines(gaus.ml, col="blue") summary(gaus.ml) ## Modelo Fam?lia Mat?rn k=0.7 matern.ml<-likfit(dados,ini=c(0.09,450), lik.method="ml", cov.model="matern", k=0.7) matern.ml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - ML - F.Matern k=0.7') lines(matern.ml, col="blue") summary(matern.ml) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% AJUSTE RML %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ######### M?todo M?xima Verossimilhan?a Restrita - RML ## Modelo exponencial exp.rml<-likfit(dados,ini=c(0.1,100), lik.method="rml", cov.model="exp") exp.rml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - RML - Exponencial') lines(exp.rml, col="blue") summary(exp.rml) ## Modelo esf?rico sph.rml<-likfit(dados,ini=c(0.1,250), lik.method="rml", cov.model="sph") sph.rml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - RML - Esf?rico') lines(sph.rml, col="blue") summary(sph.rml) ## Modelo gaussiano gaus.rml<-likfit(dados,ini=c(0.1,150), lik.method="RML", cov.model="gaus") gaus.rml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - RML - Gaussiano') lines(gaus.rml, col="blue") summary(gaus.rml) ## Modelo Fam?lia Mat?rn k=0.7 matern.rml<-likfit(dados,ini=c(0.09,100), lik.method="rml", cov.model="matern", k=0.7) matern.rml plot(d.var, xlab='Alcance', ylab='Semivari?ncia', main='Semivariograma - RML - F.Matern k=0.7') lines(matern.rml, col="blue") summary(matern.rml) ##################################################################### ### Escolha do melhor modelo ajustado ### ##################################################################### ##################################################################### ### Valida??o cruzada para modelos ajustados por OLS ### ##################################################################### # Valida??o cruzada e erro absoluto de exp.ols (vc11= Exponencial com OLS) vc11 = xvalid(dados,model=exp.ols) EA11=sum(abs(vc11$predicted-vc11$data)) EA11 summary(vc11) # Valida??o cruzada e erro absoluto de sph.ols (vc12= Esf?rico com OLS) vc12 = xvalid(dados,model=sph.ols) EA12=sum(abs(vc12$predicted-vc12$data)) EA12 summary(vc12) # Valida??o cruzada e erro absoluto de gaus.ols (vc13= Gaussiano com OLS) vc13 = xvalid(dados,model=gaus.ols) EA13=sum(abs(vc13$predicted-vc13$data)) EA13 summary(vc13) # Valida??o cruzada e erro absoluto de matern.ols (vc14= Mat?rn com OLS) vc14 = xvalid(dados,model=matern.ols) EA14=sum(abs(vc14$predicted-vc14$data)) EA14 summary(vc14) ##################################################################### ### Valida??o cruzada para modelos ajustados por WLS1 ### ##################################################################### # Valida??o cruzada e erro absoluto de exp.wls1 (vc21= Exponencial com WLS1) vc21 = xvalid(dados,model=exp.wls1) EA21=sum(abs(vc21$predicted-vc21$data)) EA21 summary(vc21) # Valida??o cruzada e erro absoluto de sph.wls1 (vc22= Esf?rico com WLS1) vc22 = xvalid(dados,model=sph.wls1) EA22=sum(abs(vc22$predicted-vc22$data)) EA22 summary(vc22) # Valida??o cruzada e erro absoluto de gaus.wls1 (vc23= Gaussiano com WLS1) vc23 = xvalid(dados,model=gaus.wls1) EA23=sum(abs(vc13$predicted-vc23$data)) EA23 summary(vc23) # Valida??o cruzada e erro absoluto de matern.wls1 (vc24= Mat?rn com WLS1) vc24 = xvalid(dados,model=matern.wls1) EA24=sum(abs(vc24$predicted-vc24$data)) EA24 summary(vc24) ##################################################################### ### Valida??o cruzada para modelos ajustados por ML ### ##################################################################### # Valida??o cruzada e erro absoluto de exp.ml (vc31= Exponencial com ML) vc31 = xvalid(dados,model=exp.ml) EA31=sum(abs(vc31$predicted-vc31$data)) EA31 summary(vc31) # Valida??o cruzada e erro absoluto de sph.ml (vc32= Esf?rico com ML) vc32 = xvalid(dados,model=sph.ml) EA32=sum(abs(vc32$predicted-vc32$data)) EA32 summary(vc32) # Valida??o cruzada e erro absoluto de gaus.ml (vc33= Gaussiano com ML) vc33 = xvalid(dados,model=gaus.ml) EA33=sum(abs(vc33$predicted-vc33$data)) EA33 summary(vc33) # Valida??o cruzada e erro absoluto de matern.ml (vc34= Mat?rn com ML) vc34 = xvalid(dados,model=matern.ml) EA34=sum(abs(vc34$predicted-vc34$data)) EA34 summary(vc34) ##################################################################### ### Valida??o cruzada para modelos ajustados por RML ### ##################################################################### # Valida??o cruzada e erro absoluto de exp.rml (vc41= Exponencial com RML) vc41 = xvalid(dados,model=exp.rml) EA41=sum(abs(vc41$predicted-vc41$data)) EA41 summary(vc41) # Valida??o cruzada e erro absoluto de sph.rml (vc42= Esf?rico com RML) vc42 = xvalid(dados,model=sph.rml) EA42=sum(abs(vc42$predicted-vc42$data)) EA42 summary(vc42) # Valida??o cruzada e erro absoluto de gaus.rml (vc43= Gaussiano com RML) vc43 = xvalid(dados,model=gaus.rml) EA43=sum(abs(vc43$predicted-vc43$data)) EA43 summary(vc43) # Valida??o cruzada e erro absoluto de matern.rml (vc44= Mat?rn com RML) vc44 = xvalid(dados,model=matern.rml) EA44=sum(abs(vc44$predicted-vc44$data)) EA44 summary(vc44) ##################################################################### ### Constru??o de mapa tem?tico - 5 classes iguais ### ##################################################################### ##################################################################### ### Adicionando bordas ### ##################################################################### # Abre o aquivo bordas.txt bor <- read.table("provab.txt", head=F) bor # apresenta os dados da borda plot(bor) polygon(bor) apply(bor,2,range) #Mostra o m?nimo e m?ximo das coordenadas gr<-expand.grid(x=seq(39636.47,41044.52,by=5), y=seq(37043.59,37713.30, by=5)) plot(gr) points(dados, pt.div="equal") #monta o grid de interpola??o require(splancs) gi<- polygrid(gr,bor=bor) points(gi, pch="+", col=2) #o novo grid considerando apenas a regi?o limitada gi KC=krige.control(obj=gaus.rml,lam=1) #krigagem d.k=krige.conv(dados, loc=gi, krige=KC) #Faz o mapa da vari?vel d por krigagem ordin?ria image(d.k, loc=gr, border=bor, col=gray(seq(1,0,l=5)),x.leg=c(39800,41000), y.leg=c(36800,36900),zlim=range(d.k$predict)) title(main = " ", font.main=4) text(240800,7237850,expression(bold(NA%up%N))) ##################################################################### ### Porcentagem de cada classe ### ##################################################################### # Carregando o pacote classInt require(classInt) # ?rea total em ha areatotal<- areapl(bor)/10000 areatotal valores<-d.k$predict # criando o vetor classIntervals(valores, 5, style="equal", intervalClosure="right") ##################################################################### ### Constru??o de mapa tem?tico - 2 classes ### ##################################################################### KC=krige.control(obj=gaus.rml,lam=1) #krigagem d.k=krige.conv(dados, loc=gi, krige=KC) #Faz o mapa da vari?vel dados por krigagem ordin?ria image(d.k, loc=gr, border=bor, col=gray(seq(1,0,l=2)),zlim=c(1.10,2.59)) ##################################################################### ### Porcentagem de cada classe ### ##################################################################### # Carregando o pacote classInt require(classInt) # ?rea total em ha areatotal<- areapl(bor)/10000 areatotal valores<-d.k$predict # criando o vetor classIntervals(valores,2, style = "fixed",fixedBreaks=c(1.10,2.59,5.0))