### Análise geoestátística ##### # Carregando pacotes require(geoR) require(maptools) require(MASS) ls() #ver as variáveis na memoria # Definindo diretorio setwd("c://geo//SP") # Importando dados sp_clima = read.table("PONTOS_SP_ALBERS.txt", head=T) # Verificação de pontos influentes #par(mfrow=c(1,2), mar=c(2,2,1,1)) boxplot(sp_clima$PPT, ylab="Chuva Anual [mm]", main="Chuva Anual SP") hist(sp_clima$PPT,freq=TRUE ,main="Chuva Anual SP", xlab="Chuva Anual [mm]") sp_clima[which(sp_clima$PPT < 1000 | sp_clima$PPT > 1700),]#ouliers limites (1000 - 1800) sp_clima1 = sp_clima[which(sp_clima$PPT > 1000 & sp_clima$PPT < 1700),]#removendo outliers #visualizando novos dados boxplot(sp_clima1$PPT, ylab="Chuva Anual [mm]", main="Chuva Anual SP") hist(sp_clima1$PPT, freq=TRUE, main="Chuva Anual SP", xlab="Chuva Anual [mm]") # Calcula as Estatísticas descritivas dos dados # média, mediana, Q1, Q3, min, máx # Análise não espacial head(sp_clima1$PPT);tail(sp_clima1$PPT);summary(sp_clima1$PPT);length(sp_clima1$PPT) var(sp_clima1$PPT) # Apresenta a Variância sd(sp_clima1$PPT) # Apresenta o Desvio padrão cv_ppt = sd(sp_clima1$PPT)*100/mean(sp_clima$PPT) cv_ppt # mostra o coeficiente de variação summary(sp_clima$ALT) var(sp_clima1$ALT) # Apresenta a Variância sd(sp_clima1$ALT) # Apresenta o Desvio padrão cv_alt = sd(sp_clima1$ALT)*100/mean(sp_clima$ALT) #covariavel cv_alt # mostra o coeficiente de variação # Criando objetos geodata ppt_sp = as.geodata(sp_clima1, coords.col=c(8, 7), data.col=5, covar.col=3)#Altitude como covariável #Transformando as coordenadas de metros -> kilometros (*1/1000) ppt_sp$coords = ppt_sp$coords/1000 # Importando shapefile com os limites dos dados sp =readShapeSpatial("SP_SAD69_ALBERS_MTS.shp", proj4string=CRS("+proj=UTM")) #Extraindo a borda, e convertendo coordenadas para Km bord = sp@polygons[[1]]@Polygons[[1]]@coords/1000 plot(bord, type='p', cex=.05) points(ppt_sp, pt.div ="quarti", cex.min=1, cex.max=1, xlab="Long (Km)", ylab="Lat (Km)") plot(ppt_sp$coords[,1],ppt_sp$coords[,2], xlab="Long [Km]", ylab="Lat [Km]") polygon(bord,cex=4) #sumário do objeto geodata resumo = summary (ppt_sp) resumo; resumo$coords.summary maiord = resumo$distances.summary[2]; maiord dist1 = maiord * 0.6; dist1 #distancia recomendada resumo$covariate.summary #plotando geodatas plot(ppt_sp$coords[,1], ppt_sp$data, xlab="Long (Km)", ylab="Chuva Anual (mm)") plot(ppt_sp$coords[,2], ppt_sp$data, xlab="Lat (Km)", ylab="Chuva Anual (mm)") plot(ppt_sp, low=T) plot(ppt_sp, trend='1st', low=T)# plot(ppt_sp, trend=~coords + ALT, low=T)# trend = ~ x1 + x2 + ALT plot(ppt_sp, trend='2nd', low=T)# trend = ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1*x2) # Ajustando os variogramas dist1 distv = 600 v1 = variog(ppt_sp, max.dist=distv, trend='1st'); plot(v1, col=1, xlab="Distância (Km)") v2 = variog(ppt_sp, max.dist=distv, trend=~coords + ALT); plot(v2, col=2, xlab="Distância (Km)") v3 = variog(ppt_sp, max.dist=distv, trend='2nd'); plot(v1, col=3, xlab="Distância (Km)") par(mfrow=c(1,1), mar=c(4,4,4,4)) v1.env = variog.mc.env(ppt_sp, obj.variog=v1);plot(v1, env=v1.env, main="M1", xlab="Distância (Km)") v2.env = variog.mc.env(ppt_sp, obj.variog=v2);plot(v2, env=v2.env, main="M2", xlab="Distância (Km)") v3.env = variog.mc.env(ppt_sp, obj.variog=v3);plot(v3, env=v3.env, main="M3", xlab="Distância (Km)") x11() ef1 = eyefit(v1) ef2 = eyefit(v2) ef3 = eyefit(v3) ef1; ef2; ef3 vf1 = variofit(v1, ini=c(8292,421), nug= 12162, cov.model="spherical")#será usado o modelo calibrado no eyefit vf2 = variofit(v2, ini=c(7883,498), nug= 11562, cov.model="spherical") vf3 = variofit(v3, ini=c(6752,483), nug= 11428, cov.model="spherical") vf1;vf2;vf3 print(vf1$value); print(vf2$value); print(vf3$value) #Semivariogramas e Modelos ajustados (Esférico) plot.new() plot(v1, col=1, xlab="Distância (Km)") lines.variomodel(m1, col=1) plot(v2, col=2, xlab="Distância (Km)") lines.variomodel(m2, col=1) plot(v1, col=3, xlab="Distância (Km)") lines.variomodel(m3, col=1) #Determinação da Max Verosimilhança e AIC m1=likfit(ppt_sp,ini=c(7626,384), cov.model="spherical", trend='1st') m2=likfit(ppt_sp,ini=c(6880,405), cov.model="spherical", trend=~coords + ALT) m3=likfit(ppt_sp,ini=c(5639,817), cov.model="spherical", trend='2nd') m1;m2;m3 summary(m1);summary(m2);summary(m3) print(m1$loglik); print(m2$loglik); print(m3$loglik) print(m1$AIC); print(m2$AIC); print(m3$AIC) # m2 <- menor AIC plot((variog4(ppt_sp, max.dist=350, trend=~coords + ALT))) v2.i = variog(ppt_sp, max.dist=distv, trend=~coords + ALT); plot(v2.i, col=2, main="Distância (Km)") v2.0 = variog(ppt_sp, max.dist=distv, trend=~coords + ALT, dir=0); plot(v2.0, col=2, xlab="Distância (Km)") v2.45 = variog(ppt_sp, max.dist=distv, trend=~coords + ALT, dir=pi/4); plot(v2.45, col=2, xlab="Distância (Km)") v2.90 = variog(ppt_sp, max.dist=distv, trend=~coords + ALT, dir=pi/2); plot(v2.90, col=2, xlab="Distância (Km)") v2.135 = variog(ppt_sp, max.dist=distv, trend=~coords + ALT, dir=3*(pi/4)); plot(v2.135, col=2, xlab="Distância (Km)") vf2.i = variofit(v2.i, ini=c(10000,422), nug=12000, cov.model="sph");vf2.i vf2.0 = variofit(v2.0, ini=c(10000,422), nug=12000, cov.model="sph");vf2.0 vf2.45 = variofit(v2.45, ini=c(12200,406), nug=11000, cov.model="sph");vf2.45 vf2.90 = variofit(v2.90, ini=c(7000, 544), nug=12000, cov.model="sph");vf2.90 vf2.135 = variofit(v2.135, ini=c(9000, 122), nug=12000, cov.model="sph", max.d=400);vf2.135 plot(v2) print(vf2.0$value); print(vf2.45$value); print(vf2.90$value); print(vf2.135$value); print(vf2.i$value) par(mfrow=c(2,2)) par(mfrow=c(1,1)) plot(v2.135) plot(v2) plot(v2.135) lines(vf2.135) m2.a =likfit(ppt_sp, ini=c(10000,405), cov.model="sph", fix.psiR=F, fix.psiA=F, trend=~coords+ALT) m2.0 =likfit(ppt_sp, ini=c(10000,405), cov.model="sph", fix.psiR=T, psiA=(0), trend=~coords+ALT) m2.45 =likfit(ppt_sp, ini=c(10000,405), cov.model="sph", fix.psiR=T, psiA=(pi/4), trend=~coords+ALT) m2.90 =likfit(ppt_sp, ini=c(10000,405), cov.model="sph", fix.psiR=T, psiA=(pi/2), trend=~coords+ALT) m2.135=likfit(ppt_sp, ini=c(10000,405), cov.model="sph", fix.psiR=T, psiA=(3*pi/4), trend=~coords+ALT) m2.a$AIC; m2.0$AIC; m2.45$AIC; m2.90$AIC; m2.135$AIC ### validacao cruzada ### xvm2 <- xvalid(ppt_sp, model=m2) summary(xvm2) write.table(cbind(xvm2, data.frame(kc[c(1,4)])),file="teste.txt") ## definindo um "grid"de predição # definiu uma malha de 51 pontos gr = expand.grid(seq(500,1600, len=100), seq(600,1400, len=100)) ## visualizando o grido onde serão feitas as predições plot.new() points(ppt_sp, xlab="Long [Km]", ylab="Lat [Km]") points(gr, pch=19, cex=0.25, col=2) polygon(bord) ppt_sp$borders = bord ## fazendo a predição espacial (krigagem) # informar quem são os dados, onde se vai predizer (malha) e usando um # modelo de dependencia espacial- se salvar mais de um modelo [[1]] kc1 = krige.conv(ppt_sp, loc=gr, krige=krige.control(obj.m=m1)) kc2 = krige.conv(ppt_sp, loc=gr, krige=krige.control(obj.m=m2)) kc3 = krige.conv(ppt_sp, loc=gr, krige=krige.control(obj.m=m3)) par(mfrow=c(3,1), mar=c(4,4,4,4)) image(kc1, col=terrain.colors(19),x.leg=c(-41.9,-41),y.leg=c(-8.0,-7.9)) image(kc2, col=terrain.colors(19),x.leg=c(-41.9,-41),y.leg=c(-8.0,-7.9)) image(kc3, col=terrain.colors(19),x.leg=c(-41.9,-41),y.leg=c(-8.0,-7.9)) #SP[,2] image(kc) ## visualizando os valores preditos na forma de um mapa image(kc,x.leg=c(-41.9,-41),y.leg=c(-8.0,-7.9),ylim=c(-8.0,-6.9)) image(kc1, col=terrain.colors(19),x.leg=c(-41.9,-41),y.leg=c(-8.0,-7.9)) contour(kc,add=T) points(PPT_SP,add=T) image(kc,col=gray(seq(1,0,l=21))) ## visualizando os erros padrão de predição (mapa de incerteza de predição) image(kc, val=sqrt(kc$krige.var), coords.data=ppt_sp$coords) # escala comum ZL=range(c(kc$pred)) par(mfrow=c(1,1), mar=c(4,4,4,4)) image(kc,zlim=ZL) ## outras opções de visualização contour(kc) plot.new() image(kc, col=terrain.colors(21)) points(ppt_sp, add=T) contour(kc, add=T, nlev=21)