<div dir="ltr">Senhores, bom dia!<div><br></div><div>Agradeço novamente ao Prof. Paulo pela atenção dispensada.</div><div><br></div><div>Hélio, quanto ao script que eu tinha ficado de olhar, vou colocar minhas considerações abaixo.</div>
<div><br></div><div>Procurei não mudar muito, mas acabei alterando um pouco. Os comentários estão no corpo do script.</div><div><br></div><div>Lembrando sempre que também estou aprendendo sobre o assunto, então pode ter muita coisa errada ou mal interpretada aí. Mas acho que pode te ajudar.</div>
<div><br></div><div>Outra coisa é que não tenho muito ideia do que os dados são, como foram obtidos e as condições da área de onde vieram, o que pra mim é fundamental pra análise. Pelos indícios nos dados (localização, estrutura), vou arriscar que são dados da produção de café ou citros e que o terreno está numa pendente da direita pra esquerda, com maior produção na parte baixa. O que chamou a atenção é o formato da área que é bem particular e não parece muito 'natural'.</div>
<div><br></div><div>Mas como o foco principal é o uso do R/geoR, vou tentar me deter nisso.</div><div><br></div><div>Use por sua própria conta e risco! :D</div><div><br></div><div><br></div><div><div><font face="courier new, monospace">### <BEGIN></font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">require(geoR) #;require(tcltk2)</font></div><div><font face="courier new, monospace">setwd("C:/LAB/RGIS/trend"); dir()</font></div>
<div><font face="courier new, monospace">par.ori <- par(no.readonly = TRUE)</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">url1 <- '<a href="https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt">https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt</a>'</font></div>
<div><font face="courier new, monospace">url2 <- '<a href="https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt">https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt</a>'</font></div><div>
<font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Como criar uma pasta pública no Dropbox...</font></div><div><font face="courier new, monospace">browseURL('<a href="https://www.dropbox.com/help/16/pt_BR">https://www.dropbox.com/help/16/pt_BR</a>')</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### 1. Carregando dados ------------------------------------------------------</font></div><div><font face="courier new, monospace">dados <- read.geodata(url1, coords.col=1:2, data.col=3, sep="", header=T)</font></div>
<div><font face="courier new, monospace">dados$borders <- read.table(url2, header=T)</font></div><div><font face="courier new, monospace">summary(dados); length(dados$data)</font></div><div><font face="courier new, monospace">ls(dados);ls.str(dados)</font></div>
<div><font face="courier new, monospace"> </font></div><div><font face="courier new, monospace">### 2. Visualizando dados ----------------------------------------------------</font></div><div><font face="courier new, monospace"><br>
</font></div><div><font face="courier new, monospace">### Visão geral</font></div><div><font face="courier new, monospace">par(mfrow=c(1,1)); par()$mfrow</font></div><div><font face="courier new, monospace">points(dados)</font></div>
<div><font face="courier new, monospace">points(dados, pt.div="quart")</font></div><div><font face="courier new, monospace">points(dados, pt.div="quint")</font></div><div><font face="courier new, monospace">### tem pelos três sub-áreas bem marcadas...</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Verificando tendências com o plot.geodata --------------------------------</font></div><div><font face="courier new, monospace">### Antes de partir pros variogramas, use o plot.geodata para pré-avaliar tendências</font></div>
<div><font face="courier new, monospace">### Bom para esclarecer opções de trend, pois notei certa confusão (cte,1st,2nd)</font></div><div><font face="courier new, monospace">### Recomendo ler o item 'trend' na ajuda da função plot.geodata</font></div>
<div><font face="courier new, monospace">?plot.geodata </font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">plot(dados, lowess=T)</font></div><div><font face="courier new, monospace">plot(dados, lowess=T, trend='cte') ### mesmo que anterior</font></div>
<div><font face="courier new, monospace">### linhas de comando anteriores são equivalentes: trend='cte' é 'default' da geoR!</font></div><div><font face="courier new, monospace">### nos gráf. 2 e 3, plota-se as coord vs. dados. Observar as distorções na linha 'lowess'</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">plot(dados, lowess=T, trend=~coords)</font></div><div><font face="courier new, monospace">plot(dados, lowess=T, trend='1st') </font></div>
<div><font face="courier new, monospace">plot(dados, lowess=T, trend=~dados$coords[,1]+dados$coords[,2])</font></div><div><font face="courier new, monospace">### as três linhas anteriores são equivalentes</font></div><div>
<font face="courier new, monospace">### o uso da polinomial de primeira ordem melhora o comportamento do histograma</font></div><div><font face="courier new, monospace">### mas a linha de lowess continua com comportamento anômalo</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">plot(dados, lowess=T, trend='2nd') </font></div><div><font face="courier new, monospace">plot(dados, lowess=T, trend=~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]))</font></div>
<div><font face="courier new, monospace">### as duas linhas anteriores são equivalentes (polinomal de segunda ordem)!</font></div><div><font face="courier new, monospace">### não aparenta agregar nenhuma melhoria</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Mais sobre modelos de tendência com trend.spatial ------------------------</font></div><div><font face="courier new, monospace">### Pra melhorar o entendimento dos modelos, dá pra dar uma olhada nas matrizes</font></div>
<div><font face="courier new, monospace">?trend.spatial</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial('cte', dados)))   ### constante</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, usando coords</font></div>
<div><font face="courier new, monospace">head(unclass(trend.spatial('1st', dados)))   ### mesmo que anterior</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial('2nd', dados)))   ### polinomal 2a ordem, usando coords</font></div>
<div><font face="courier new, monospace">head(unclass(trend.spatial(~coords[,1], dados))) ### em função de X</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial(~coords[,2], dados))) ### em função de Y</font></div>
<div><font face="courier new, monospace">head(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords + X**2</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial(~I(coords)+I(coords^2), dados)))</font></div>
<div><font face="courier new, monospace">head(unclass(trend.spatial(~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]), dados))) ### mesmo que 2nd</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial(~I(coords)+I(coords[,1]^2), dados)))</font></div>
<div><font face="courier new, monospace">head(unclass(trend.spatial(~I(coords^2), dados)))</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial(~coords[,1]*coords[,2], dados)))</font></div><div>
<font face="courier new, monospace">head(unclass(trend.spatial(~I(coords[,1]*coords[,2]), dados)))</font></div><div><font face="courier new, monospace">head(unclass(trend.spatial(~coords[,1]+I(coords[,1]^2)+coords[,2], dados)))</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Separando o drift com base na matriz trend.spatial -----------------------</font></div><div><font face="courier new, monospace">trend.matrix <- unclass(trend.spatial(trend='1st', geodata=dados)); nrow(trend.matrix)</font></div>
<div><font face="courier new, monospace">trend.lm     <- lm(dados$data ~ trend.matrix)</font></div><div><font face="courier new, monospace">trend.res    <- lm(dados$data ~ trend.matrix)$residuals</font></div><div><font face="courier new, monospace">trend.fit    <- lm(dados$data ~ trend.matrix)$fitted.values</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### 'Empacotando' os resíduos em um geodata-----------------------------------</font></div><div><font face="courier new, monospace">dados.res <- dados</font></div>
<div><font face="courier new, monospace">names(trend.res) <- NULL; trend.res</font></div><div><font face="courier new, monospace">dados.res$data <- trend.res; dados.res$data</font></div><div><font face="courier new, monospace">str(dados.res)</font></div>
<div><font face="courier new, monospace">### Pode usar dados.res pra fazer a análise variográfica no lugar de dados, </font></div><div><font face="courier new, monospace">### caso tenha problemas com krige.conv()!!!</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Embora o uso da polinomial de 1a ordem tenha representado alguma melhoria,</font></div><div><font face="courier new, monospace">### eu não fiquei muito entusiamado por conta do 'lowess' apresentado</font></div>
<div><font face="courier new, monospace">### Talvez fosse o caso de testar covariáveis. </font></div><div><font face="courier new, monospace">### Uma fácil de obter e que pode ter relação é a altitude...</font></div><div>
<font face="courier new, monospace">### Mas vamos investigar o efeito nos variogramas</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### 3. Distancia máxima e parâmetros para semivariograma ---------------------</font></div>
<div><font face="courier new, monospace">hmax <- summary(dados)[[3]][[2]]*.8 ; hmax ### 80% da máxima distância</font></div><div><font face="courier new, monospace">uvec <- 12; pairs <- 20 ### pairs.min</font></div>
<div><font face="courier new, monospace">           </font></div><div><font face="courier new, monospace">### 4. Variogramas -----------------------------------------------------------</font></div><div><font face="courier new, monospace">### Não precisaria rodar todos, afinal como visto, tá repetindo modelos!!!</font></div>
<div><font face="courier new, monospace">### v0:média constante; v1:1st; v2:2nd; v3:cte; v4:coords</font></div><div><font face="courier new, monospace">v0 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical")</font></div>
<div><font face="courier new, monospace">v1 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="1st")</font></div><div><font face="courier new, monospace">v2 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="2nd")</font></div>
<div><font face="courier new, monospace">v3 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="cte")</font></div><div><font face="courier new, monospace">v4 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend=~coords)</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### envelopes</font></div><div><font face="courier new, monospace">v0.env <- variog.mc.env(dados, obj.var = v0)</font></div>
<div><font face="courier new, monospace">v1.env <- variog.mc.env(dados, obj.var = v1)</font></div><div><font face="courier new, monospace">v2.env <- variog.mc.env(dados, obj.var = v2)</font></div><div><font face="courier new, monospace">v3.env <- variog.mc.env(dados, obj.var = v3)</font></div>
<div><font face="courier new, monospace">v4.env <- variog.mc.env(dados, obj.var = v4)</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### 5. Visualizando variogramas ----------------------------------------------</font></div>
<div><font face="courier new, monospace">par(par.ori)</font></div><div><font face="courier new, monospace">par(mfrow=c(2,3))</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">points(dados, pt.div="quartiles",xlab="Leste", ylab="Norte")</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">plot(v0, xlab="Dist. (m)", ylim=c(0,summary(v0$v)[6]*1.1), env=v0.env, main="Isotrópico com tendência")</font></div>
<div><font face="courier new, monospace">text(v0$u,v0$v,round(v0$n,1),pos=1);text(v0$u,v0$v,round(v0$u,1),pos=3);lines(v0$u,v0$v);abline(h=mean(v0$v))</font></div><div><font face="courier new, monospace"><br></font></div>
<div><font face="courier new, monospace">plot(v1, xlab="Dist. (m)", ylim=c(0,summary(v1$v)[6]*1.1), env=v1.env, main="Isotrópico trend 1st")</font></div><div><font face="courier new, monospace">text(v1$u,v1$v,round(v1$n,1),pos=1);text(v1$u,v1$v,round(v1$u,1),pos=3);lines(v1$u,v1$v);abline(h=mean(v1$v))</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">plot(v2, xlab="Dist. (m)", ylim=c(0,summary(v2$v)[6]*1.1), env=v2.env, main="Isotrópico trend 2nd")</font></div>
<div><font face="courier new, monospace">text(v2$u,v2$v,round(v2$n,1),pos=1);text(v2$u,v2$v,round(v2$u,1),pos=3);lines(v2$u,v2$v);abline(h=mean(v2$v))</font></div><div><font face="courier new, monospace"><br></font></div>
<div><font face="courier new, monospace">plot(v3, xlab="Dist. (m)", ylim=c(0,summary(v3$v)[6]*1.1), env=v3.env, main="Isotrópico trend cte")</font></div><div><font face="courier new, monospace">text(v3$u,v3$v,round(v3$n,1),pos=1);text(v3$u,v3$v,round(v3$u,1),pos=3);lines(v3$u,v3$v);abline(h=mean(v3$v))</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">plot(v4, xlab="Dist. (m)", ylim=c(0,summary(v4$v)[6]*1.1), env=v4.env, main="Isotrópico trend coords")</font></div>
<div><font face="courier new, monospace">text(v4$u,v4$v,round(v4$n,1),pos=1);text(v4$u,v4$v,round(v4$u,1),pos=3);lines(v4$u,v4$v);abline(h=mean(v4$v))</font></div><div><font face="courier new, monospace"><br></font></div>
<div><font face="courier new, monospace">### Analisando os variogramas ------------------------------------------------</font></div><div><font face="courier new, monospace">### Realmente o comportamento do histograma de trend='cte' parece conter uma tendência,</font></div>
<div><font face="courier new, monospace">### mas os modelos de 'retirada' usados deixam os dados sem uma estrutura de </font></div><div><font face="courier new, monospace">### variação. Dificilmente vai conseguir ajustar um bom modelo nesses casos, pois </font></div>
<div><font face="courier new, monospace">### o comportamento é muito errático. É possível que O erro que você estava obtendo </font></div><div><font face="courier new, monospace">### seja decorrente da fragilidade do modelo ajustado.</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Realmente não ficou bom! Mas e aí? O que dá pra fazer? Tenta ssim mesmo?</font></div><div><font face="courier new, monospace">### Ainda tem que testar a anisotropia, que não tem no seu script original.</font></div>
<div><font face="courier new, monospace">?variog4</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">va <- variog4(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical")</font></div>
<div><font face="courier new, monospace">par(par.ori)</font></div><div><font face="courier new, monospace">plot(va)</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Pronto! Apareceu anisotropia...</font></div>
<div><font face="courier new, monospace">### Tem duas sugestões de leitura abaixo pra dar uma ajuda...</font></div><div><font face="courier new, monospace">browseURL('<a href="http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf">http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf</a>')</font></div>
<div><font face="courier new, monospace">browseURL('<a href="https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298">https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298</a>')</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Você vai ter que estudar adequadamente a questão da anisotropia em seus dados</font></div><div><font face="courier new, monospace">### Isso não te isenta do efeito combinado anisotropia + tendência</font></div>
<div><font face="courier new, monospace">### Só pra testar as possibilidades, vamos fazer o variograma na direção 135 graus</font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Vai ter que entrar a direções, sendo mais fácil usar PI/n</font></div>
<div><font face="courier new, monospace">### Tabela das direções (PI/n) para variogramas</font></div><div><font face="courier new, monospace">denom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadores</font></div>
<div><font face="courier new, monospace">data.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), deg=pi/denom/k)</font></div><div><font face="courier new, monospace">### 135 é 3*45 ou 3*pi/4</font></div><div>
<font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">va.135 <- variog(dados, uvec=uvec, pairs.min=pairs, dir=3*pi/4, max.dist=hmax, estimator="classical")</font></div>
<div><font face="courier new, monospace">plot(va.135)</font></div><div><font face="courier new, monospace">### Aqui parece ser mais fácil ajustar. Talvez tenha que ajustar o 'lag'</font></div><div><font face="courier new, monospace">### Avaliar a questão da distância entre amostras...</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### Lembra do resíduo que separei anteriormente?</font></div><div><font face="courier new, monospace">### Se fosse usar...</font></div>
<div><font face="courier new, monospace">va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical")</font></div><div><font face="courier new, monospace">va.r.env <- variog.mc.env(dados.res, obj.var = va.r)</font></div>
<div><font face="courier new, monospace">plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), env=va.r.env, main="Residual")</font></div><div><font face="courier new, monospace"><br></font></div>
<div><font face="courier new, monospace">### Vou parar por aqui, porque antes de prosseguir vai precisar ver </font></div><div><font face="courier new, monospace">### adequadamente essa questão da anisotropia.</font></div>
<div><font face="courier new, monospace">### No caso do efeito ser só da anisotropia, vai simplificar a krigagem e você</font></div><div><font face="courier new, monospace">### não deve mais ter o problema que vinha tendo.</font></div>
<div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">### <END></font></div></div><div><br></div><div class="gmail_extra"><div><div dir="ltr">Éder Comunello <<a href="mailto:comunello.eder@gmail.com" target="_blank">c</a><a href="mailto:omunello.eder@gmail.com" target="_blank">omunello.eder@gmail.com</a>> <br>
Dourados, MS - [22 16.5'S, 54 49'W]<br></div></div>
<br><br><div class="gmail_quote"><br></div></div></div>