Caro Eder, Já havia feito a retirada de tendencia com o atributo altura da planta, que está ligado com a produção.
O variograma gerado ficou diferente se usar polinomial de 1st, mas não sei como validar isso depois... Abraço Em 31 de outubro de 2013 13:15, Éder Comunello <[email protected]>escreveu: > Boa tarde, Helio! > > Vi seu email, mas vou precisar reservar um tempo pra tratar dele. Você > chegou a usar as glebas como covariável na modelagem? Para isso cada ponto > deveria ter uma var a mais identificando a gleba origem. > > Atte., > > > > Éder Comunello <c <[email protected]>[email protected]> > Dourados, MS - [22 16.5'S, 54 49'W] > > > Em 29 de outubro de 2013 22:25, Hélio Gallo Rocha < > [email protected]> escreveu: > >> Caro Eder, >> >> Só hoje vi seu post com CRM anexo. >> >> Muito obrigado pela revisão do CRM. >> >> Como também está a procura de uma maior entendimento da geoestatística, >> vamos seguindo. >> >> >> Vou numerar para separar as observações: >> >> 1. Pra entender melhor a área, realmente é café, o recorte esta meio >> estranho, mas na realidade são três talhões de 0.5 ha cada. >> Cada um tem uma característica bem diferenciada com relação a >> produção:alta, baixa e média. >> >> Segundo li, numa situação em que a média é flutuante, como parece ser o >> caso, teríamos de usar krigagem ordinária. >> >> >> A maioria das observações tinha notado, a linha: >> plot(dados, lowess=T, trend='1st') >> achei interessante, mostrando o resultado da retirada de tendencia. >> >> 2. no final : >> >> ### Lembra do resíduo que separei anteriormente? >> ### Se fosse usar... >> va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, >> estimator="classical") >> va.r.env <- variog.mc.env(dados.res, obj.var = va.r) >> plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), >> env=va.r.env, main="Residual") >> >> Foi possível gerar o semivariograma sem tendencia, com os resíduos, >> fazer o ajuste assentimento usando: eyefit(va.r) >> >> Dai fazer a validação, fazendo a validação ou não e gerar o mapa sem >> tendencia. >> >> 3. Quanto ao envelopamento, não sei se estou certo, mas ele simula >> valores máximos e minimos para dependencia, mas como é uma simulação, a >> cada vez que vc. roda, pode ter pontos com dependencia, nem que seje um, e >> as vezes não dá nenhum, mas de qualquer forma indica que o valor do alcance >> é pequeno em relação ao todo e baixa dependencia espacial, como nos meus >> dados. >> >> Esse item do envelopamento vale discussão, principalmente em situações >> que a variancia sobe além do limite superior e desce, isso significa algo >> importante, mas o que? >> >> 4. Quanto aos links do Landim já tinha visto, usando o surfer, também >> vale atenção. >> Fiz o recomendado por ele, mas dobrou os valores do mapa., vc. tentou? >> >> 5. Quanto ao ajuste por verosimilhança colocada pelo Prof. Paulo: >> Na primeira vez que havia rodado os dados achei muito estranho, pois o >> semivariograma subia indefinidamente e o semivariograma teórico por >> verossimilhança, usando likfit, baixou muito o patamar. >> Interessante se usar likfit com trend="1st", o semivariograma teórico >> baixa mais ainda, das dai não consigo krigar >> >> 6. Não sei se posso fazer o ajuste assentimento usando o semivariograma >> com tendencia mas usar os parametros sem tendencia. >> >> 7. Tenho de retirar a tendencia sempre? >> Tem casos que há uma patamar claro, mas se rodar o variog com trend, fica >> efeito pepita puro...e ai? >> >> 8. Outra coisa que não sei se estou entendendo: >> após a retirada da tendencia, fazer o mapa com o residuo, mas os valores >> não tem nada a ver com os reais, onde estou entendendo errado? >> >> Abraço, >> >> Hélio >> >> >> >> >> >> Em 26 de outubro de 2013 11:05, Eder Comunello [via R-br] < >> [email protected]> escreveu: >> >>> Senhores, bom dia! >>> >>> Agradeço novamente ao Prof. Paulo pela atenção dispensada. >>> >>> Hélio, quanto ao script que eu tinha ficado de olhar, vou colocar minhas >>> considerações abaixo. >>> >>> Procurei não mudar muito, mas acabei alterando um pouco. Os comentários >>> estão no corpo do script. >>> >>> 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. >>> >>> 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'. >>> >>> Mas como o foco principal é o uso do R/geoR, vou tentar me deter nisso. >>> >>> Use por sua própria conta e risco! :D >>> >>> >>> ### <BEGIN> >>> >>> require(geoR) #;require(tcltk2) >>> setwd("C:/LAB/RGIS/trend"); dir() >>> par.ori <- par(no.readonly = TRUE) >>> >>> url1 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt' >>> url2 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt >>> ' >>> >>> ### Como criar uma pasta pública no Dropbox... >>> browseURL('https://www.dropbox.com/help/16/pt_BR') >>> >>> ### 1. Carregando dados >>> ------------------------------------------------------ >>> dados <- read.geodata(url1, coords.col=1:2, data.col=3, sep="", header=T) >>> dados$borders <- read.table(url2, header=T) >>> summary(dados); length(dados$data) >>> ls(dados);ls.str(dados) >>> >>> ### 2. Visualizando dados >>> ---------------------------------------------------- >>> >>> ### Visão geral >>> par(mfrow=c(1,1)); par()$mfrow >>> points(dados) >>> points(dados, pt.div="quart") >>> points(dados, pt.div="quint") >>> ### tem pelos três sub-áreas bem marcadas... >>> >>> ### Verificando tendências com o plot.geodata >>> -------------------------------- >>> ### Antes de partir pros variogramas, use o plot.geodata para >>> pré-avaliar tendências >>> ### Bom para esclarecer opções de trend, pois notei certa confusão >>> (cte,1st,2nd) >>> ### Recomendo ler o item 'trend' na ajuda da função plot.geodata >>> ?plot.geodata >>> >>> plot(dados, lowess=T) >>> plot(dados, lowess=T, trend='cte') ### mesmo que anterior >>> ### linhas de comando anteriores são equivalentes: trend='cte' é >>> 'default' da geoR! >>> ### nos gráf. 2 e 3, plota-se as coord vs. dados. Observar as distorções >>> na linha 'lowess' >>> >>> plot(dados, lowess=T, trend=~coords) >>> plot(dados, lowess=T, trend='1st') >>> plot(dados, lowess=T, trend=~dados$coords[,1]+dados$coords[,2]) >>> ### as três linhas anteriores são equivalentes >>> ### o uso da polinomial de primeira ordem melhora o comportamento do >>> histograma >>> ### mas a linha de lowess continua com comportamento anômalo >>> >>> plot(dados, lowess=T, trend='2nd') >>> plot(dados, lowess=T, >>> trend=~I(coords)+I(coords^2)+I(coords[,1]*coords[,2])) >>> ### as duas linhas anteriores são equivalentes (polinomal de segunda >>> ordem)! >>> ### não aparenta agregar nenhuma melhoria >>> >>> ### Mais sobre modelos de tendência com trend.spatial >>> ------------------------ >>> ### Pra melhorar o entendimento dos modelos, dá pra dar uma olhada nas >>> matrizes >>> ?trend.spatial >>> head(unclass(trend.spatial('cte', dados))) ### constante >>> head(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, >>> usando coords >>> head(unclass(trend.spatial('1st', dados))) ### mesmo que anterior >>> head(unclass(trend.spatial('2nd', dados))) ### polinomal 2a ordem, >>> usando coords >>> head(unclass(trend.spatial(~coords[,1], dados))) ### em função de X >>> head(unclass(trend.spatial(~coords[,2], dados))) ### em função de Y >>> head(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords >>> + X**2 >>> head(unclass(trend.spatial(~I(coords)+I(coords^2), dados))) >>> head(unclass(trend.spatial(~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]), >>> dados))) ### mesmo que 2nd >>> head(unclass(trend.spatial(~I(coords)+I(coords[,1]^2), dados))) >>> head(unclass(trend.spatial(~I(coords^2), dados))) >>> head(unclass(trend.spatial(~coords[,1]*coords[,2], dados))) >>> head(unclass(trend.spatial(~I(coords[,1]*coords[,2]), dados))) >>> head(unclass(trend.spatial(~coords[,1]+I(coords[,1]^2)+coords[,2], >>> dados))) >>> >>> ### Separando o drift com base na matriz trend.spatial >>> ----------------------- >>> trend.matrix <- unclass(trend.spatial(trend='1st', geodata=dados)); >>> nrow(trend.matrix) >>> trend.lm <- lm(dados$data ~ trend.matrix) >>> trend.res <- lm(dados$data ~ trend.matrix)$residuals >>> trend.fit <- lm(dados$data ~ trend.matrix)$fitted.values >>> >>> ### 'Empacotando' os resíduos em um >>> geodata----------------------------------- >>> dados.res <- dados >>> names(trend.res) <- NULL; trend.res >>> dados.res$data <- trend.res; dados.res$data >>> str(dados.res) >>> ### Pode usar dados.res pra fazer a análise variográfica no lugar de >>> dados, >>> ### caso tenha problemas com krige.conv()!!! >>> >>> ### Embora o uso da polinomial de 1a ordem tenha representado alguma >>> melhoria, >>> ### eu não fiquei muito entusiamado por conta do 'lowess' apresentado >>> ### Talvez fosse o caso de testar covariáveis. >>> ### Uma fácil de obter e que pode ter relação é a altitude... >>> ### Mas vamos investigar o efeito nos variogramas >>> >>> ### 3. Distancia máxima e parâmetros para semivariograma >>> --------------------- >>> hmax <- summary(dados)[[3]][[2]]*.8 ; hmax ### 80% da máxima distância >>> uvec <- 12; pairs <- 20 ### pairs.min >>> >>> ### 4. Variogramas >>> ----------------------------------------------------------- >>> ### Não precisaria rodar todos, afinal como visto, tá repetindo >>> modelos!!! >>> ### v0:média constante; v1:1st; v2:2nd; v3:cte; v4:coords >>> v0 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical") >>> v1 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend="1st") >>> v2 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend="2nd") >>> v3 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend="cte") >>> v4 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend=~coords) >>> >>> ### envelopes >>> v0.env <- variog.mc.env(dados, obj.var = v0) >>> v1.env <- variog.mc.env(dados, obj.var = v1) >>> v2.env <- variog.mc.env(dados, obj.var = v2) >>> v3.env <- variog.mc.env(dados, obj.var = v3) >>> v4.env <- variog.mc.env(dados, obj.var = v4) >>> >>> ### 5. Visualizando variogramas >>> ---------------------------------------------- >>> par(par.ori) >>> par(mfrow=c(2,3)) >>> >>> points(dados, pt.div="quartiles",xlab="Leste", ylab="Norte") >>> >>> plot(v0, xlab="Dist. (m)", ylim=c(0,summary(v0$v)[6]*1.1), env=v0.env, >>> main="Isotrópico com tendência") >>> >>> 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)) >>> >>> plot(v1, xlab="Dist. (m)", ylim=c(0,summary(v1$v)[6]*1.1), env=v1.env, >>> main="Isotrópico trend 1st") >>> >>> 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)) >>> >>> plot(v2, xlab="Dist. (m)", ylim=c(0,summary(v2$v)[6]*1.1), env=v2.env, >>> main="Isotrópico trend 2nd") >>> >>> 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)) >>> >>> plot(v3, xlab="Dist. (m)", ylim=c(0,summary(v3$v)[6]*1.1), env=v3.env, >>> main="Isotrópico trend cte") >>> >>> 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)) >>> >>> plot(v4, xlab="Dist. (m)", ylim=c(0,summary(v4$v)[6]*1.1), env=v4.env, >>> main="Isotrópico trend coords") >>> >>> 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)) >>> >>> ### Analisando os variogramas >>> ------------------------------------------------ >>> ### Realmente o comportamento do histograma de trend='cte' parece conter >>> uma tendência, >>> ### mas os modelos de 'retirada' usados deixam os dados sem uma >>> estrutura de >>> ### variação. Dificilmente vai conseguir ajustar um bom modelo nesses >>> casos, pois >>> ### o comportamento é muito errático. É possível que O erro que você >>> estava obtendo >>> ### seja decorrente da fragilidade do modelo ajustado. >>> >>> ### Realmente não ficou bom! Mas e aí? O que dá pra fazer? Tenta ssim >>> mesmo? >>> ### Ainda tem que testar a anisotropia, que não tem no seu script >>> original. >>> ?variog4 >>> >>> va <- variog4(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical") >>> par(par.ori) >>> plot(va) >>> >>> ### Pronto! Apareceu anisotropia... >>> ### Tem duas sugestões de leitura abaixo pra dar uma ajuda... >>> browseURL(' >>> http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf') >>> browseURL(' >>> https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298') >>> >>> ### Você vai ter que estudar adequadamente a questão da anisotropia em >>> seus dados >>> ### Isso não te isenta do efeito combinado anisotropia + tendência >>> ### Só pra testar as possibilidades, vamos fazer o variograma na direção >>> 135 graus >>> >>> ### Vai ter que entrar a direções, sendo mais fácil usar PI/n >>> ### Tabela das direções (PI/n) para variogramas >>> denom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadores >>> data.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), >>> deg=pi/denom/k) >>> ### 135 é 3*45 ou 3*pi/4 >>> >>> va.135 <- variog(dados, uvec=uvec, pairs.min=pairs, dir=3*pi/4, >>> max.dist=hmax, estimator="classical") >>> plot(va.135) >>> ### Aqui parece ser mais fácil ajustar. Talvez tenha que ajustar o 'lag' >>> ### Avaliar a questão da distância entre amostras... >>> >>> ### Lembra do resíduo que separei anteriormente? >>> ### Se fosse usar... >>> va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical") >>> va.r.env <- variog.mc.env(dados.res, obj.var = va.r) >>> plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), >>> env=va.r.env, main="Residual") >>> >>> ### Vou parar por aqui, porque antes de prosseguir vai precisar ver >>> ### adequadamente essa questão da anisotropia. >>> ### No caso do efeito ser só da anisotropia, vai simplificar a krigagem >>> e você >>> ### não deve mais ter o problema que vinha tendo. >>> >>> ### <END> >>> >>> Éder Comunello <[hidden >>> email]<http://user/SendEmail.jtp?type=node&node=4660724&i=0>[hidden >>> email] <http://user/SendEmail.jtp?type=node&node=4660724&i=1>> >>> Dourados, MS - [22 16.5'S, 54 49'W] >>> >>> >>> >>> >>> _______________________________________________ >>> R-br mailing list >>> [hidden email] <http://user/SendEmail.jtp?type=node&node=4660724&i=2> >>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br >>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça >>> código mínimo reproduzível. >>> >>> ------------------------------ >>> If you reply to this email, your message will be added to the >>> discussion below: >>> >>> http://r-br.2285057.n4.nabble.com/Re-R-br-erro-na-validacao-cruzada-com-pacote-geoR-tp4660619p4660724.html >>> To unsubscribe from R-br, click >>> here<http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3357982&code=aGVsaW9nYWxsb3JvY2hhQGdtYWlsLmNvbXwzMzU3OTgyfC0xMzQ3NTkwMDY4> >>> . >>> NAML<http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml> >>> >> >> >> >> -- >> Hélio Gallo Rocha >> IFSULDEMINAS - Câmpus Muzambinho >> >> _______________________________________________ >> R-br mailing list >> [email protected] >> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br >> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça >> código mínimo reproduzível. >> > > -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho
_______________________________________________ R-br mailing list [email protected] https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
