[R-br] Diferenças entre sistemas operacionais (geoR)
Yury Duarte
yurynepomuceno em gmail.com
Qui Out 4 20:12:15 -03 2018
Boa noite Paulo, como vai?
Agradeço desde já pela atenção.
Tentei remover ao máximo as partes que não afetam na função e anexei o
script no arquivo texto que segue.
Att
Em qui, 4 de out de 2018 19:23, Paulo Justiniano Ribeiro Junior <
paulojus em ufpr.br> escreveu:
> Olá!
>
> krige.conv() só utiliza simulações caso seja indicado em output.options
> o número de simulações da distribuição preditiva.
> Neste caso são geradas simulações.
> As simulações só são geradas se:
> 1. especificado explicitamente
> 2. se for utilizada alguma transformação com lambda diferente de 0, 1/2
>
> O seu exemplo se enquadra em algum caso?
>
> A msg tb indica que voce colocou algum em micro.scale e seria bom saber o
> que deseja com isto.
>
> Se dados e códigos são grande Envie o comando completo da chamada da
> krige.conv() para vermos as opções
>
> P.J.
>
>
>
>
> ------------------------------
> *From: *"Yury Duarte" <yurynepomuceno em gmail.com>
> *To: *"Thiago V. dos Santos" <thi_veloso em yahoo.com.br>
> *Cc: *"a lista Brasileira oficial de discussão do programa R." <
> r-br em listas.c3sl.ufpr.br>, paulojus em ufpr.br, "p diggle" <
> p.diggle em lancaster.ac.uk>
> *Sent: *Thursday, 4 October, 2018 6:08:32 PM
> *Subject: *Re: [R-br] Diferenças entre sistemas operacionais (geoR)
>
> Não sei se existe um gerador de números aleatórios dentro da função, mas
> estou passando pra ela todos os parâmetros já pré definidos antes dela ser
> chamada.
> Vou colocar um set.sead() no início do código, conforme sugestão do César
> para então comparar os valores gerados nos diferentes sistemas operacionais
> e ver o que acontece.
>
> Em qui, 4 de out de 2018 17:45, Thiago V. dos Santos <
> thi_veloso em yahoo.com.br> escreveu:
>
>> Da maneira que eu vejo, a unica explicacao plausivel seria se a
>> krige.conv usasse um algoritmo gerador de numeros aleatorios. Caso voce
>> tenha experiencia com esmagamento de bugs, vc pode conferir isso no codigo
>> dela. Ou entao esperar o feedback dos autores.
>>
>> Greetings,
>> -- Thiago V. dos Santos
>>
>> Postdoctoral Research Fellow
>> Department of Climate and Space Science and Engineering
>> University of Michigan
>>
>>
>> On Thursday, October 4, 2018, 4:39:08 PM EDT, Yury Duarte <
>> yurynepomuceno em gmail.com> wrote:
>>
>>
>> Sim, Thiago.
>> Até o último passo os resultados são iguais. A mesma grade com o mesmo
>> número de pontos é gerada. As mesmas informações são passadas como
>> parâmetro para a função, mas os resultados não batem.
>> Nem mesmo em casos onde a função não quebra, os valores gerados no Ubuntu
>> não batem com os valores gerados em Windows. Isso está me intrigando muito.
>>
>> Em qui, 4 de out de 2018 17:23, Thiago V. dos Santos <
>> thi_veloso em yahoo.com.br> escreveu:
>>
>> Nao deveria ser assim.
>>
>> O seu codigo produz resultados identicos ate o momento em que voce roda a
>> krige.conv?
>>
>> Em caso positivo, acho que seria uma boa ideia reportar esse erro ao
>> Paulo (paulojus em ufpr.br) ou ao Peter (p.diggle em lancaster.ac.uk).
>>
>> Greetings,
>> -- Thiago V. dos Santos
>>
>> Postdoctoral Research Fellow
>> Department of Climate and Space Science and Engineering
>> University of Michigan
>>
>>
>> On Thursday, October 4, 2018, 3:42:35 PM EDT, Yury Duarte via R-br <
>> r-br em listas.c3sl.ufpr.br> wrote:
>>
>>
>> Boa tarde a todos os coleas!
>>
>> Fiz um código simples para interpolar valores que tenho em uma grade que
>> construo dentro do próprio código. Para isso, estou utilizando um modelo
>> esférico de interpolação por krigagem através da função krige.conv().
>> Estou achando muito estranho o fato de que, quando rodo o script em
>> ambiente Windows, não tenho problemas e meu script gera as saídas
>> exatamente como o esperado, entretanto, quando rodo o mesmo script (com o
>> mesmo banco de dados) em ambiente Linux (Ubuntu), a função krige.conv()
>> corrompe e o erro gerado é:
>> krige.control: micro.scale must be in the interval [0, nugget]
>>
>> Não encontrei nada em outros fóruns que pudesse me ajudar para explicar
>> isso.
>> Caso alguém tenha interesse em testar exatamente o mesmo conjunto de
>> dados que estou utilizando, informem que passo um link no google drive, já
>> que o código e os dados de entrada são grandes demais para serem anexados
>> no corpo do e-mail/nos anexos.
>>
>> Desde já, agradeço a atenção e a ajuda de todos!
>>
>> Yury Duarte
>> Engenheiro Agrônomo - ESALQ/USP
>> _______________________________________________
>> R-br mailing list
>> R-br em listas.c3sl.ufpr.br
>> 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.
>>
>>
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20181004/63935a33/attachment.html>
-------------- Próxima Parte ----------
#-----Inicio do Script (rodar funcoes antes)-----#
#-----Funcoes se encontram no script 'funcoes_NDVI'-----#
rodar_bibliotecas(TRUE)
raiz = 'C:\\Users\\Yury\\Desktop\\InCeres\\Projeto_Fertilidade\\Sprint\\02\\'
#-----Ler arquivos de entrada-----#
shapefile_contorno = readOGR(paste0(raiz, 'fields_6_areas_validas.shp'))
#imagem_satelite = stack(paste0(raiz, 'imagem.tif'))
imagem_satelite_04 = stack(paste0(raiz, 'IMG_04.tif'))
imagem_satelite_05 = stack(paste0(raiz, 'IMG_05.tif'))
#-----Transformar projecao e recortar .TIF-----#
shapefile_contorno = transformar_em_coordenadas_utm(shapefile_contorno)
datum_correto = proj4string(shapefile_contorno)
numero_talhoes = 1:(length(shapefile_contorno[,1]))
#i=3
for (i in 1:length(numero_talhoes)){
#-----Definir a separacao dos poligonos para clusterizar um de cada vez-----#
shapefile_contorno_do_talhao_original = shapefile_contorno[i,1]
shapefile_contorno_do_talhao = shapefile_contorno[i,1]
shapefile_contorno_do_talhao = buffer(shapefile_contorno_do_talhao, width = 60)
shapefile_contorno_para_recorte = reprojetar_para_recortar_utm(shapefile_contorno_do_talhao,imagem_satelite_04)
imagem_recortada_04 = crop(imagem_satelite_04, shapefile_contorno_para_recorte)
imagem_recortada_05 = crop(imagem_satelite_05, shapefile_contorno_para_recorte)
reprojetado_04 = projectRaster(imagem_recortada_04, crs = datum_correto)
reprojetado_05 = projectRaster(imagem_recortada_05, crs = datum_correto)
#-----calcular NDVI-----#
a = reprojetado_05+reprojetado_04
b = reprojetado_05-reprojetado_04
raster_ndvi = b/a
raster_produtividade = 0.1665*exp(8.9245*(raster_ndvi))
raster_extracao = raster_produtividade*0.63
q2 = list(raster_ndvi, raster_produtividade, raster_extracao)
tipo = c('NDVI', 'Produtividade', 'Extracao_K')
#j=1
for(j in 1:length(q2)){
raster_ndvi_recortado = raster::mask(q2[[j]], shapefile_contorno_do_talhao)
shapefile_pontos_ndvi = rasterToPoints(raster_ndvi_recortado, spatial = TRUE)
shapefile_pontos_ndvi = intersect(shapefile_pontos_ndvi, shapefile_contorno_do_talhao_original)
dataframe_ndvi = rasterToPoints(raster_ndvi_recortado, spatial = FALSE)
write.table(dataframe_ndvi,(paste0(raiz, i, '_', tipo[j], '_dataframe_ndvi.txt')), sep = ';', row.names = F)
#plot(shapefile_pontos_ndvi)
#-----Carregando os dados para interpolacao-----#
dados_ndvi = data.frame(dataframe_ndvi)
attach(dados_ndvi)
#names(dados_ndvi)
nome_analise = 'Krigagem_NDVI'
#-----Carrega dados no geoR (pacote de geoestatistica)-----#
geo_dados_ndvi = read.geodata(paste0(raiz, i, '_', 'dataframe_ndvi.txt'), sep = ';', header = TRUE, coords.col=1:2, data.col=3)
#-----Estatistica analitica-----#
#estb = basicStats(layer, ci=0.95)
#tks = ks.test(layer, 'pnorm', mean=mean(layer), sd=sd(layer))
#tks
#shapiro.test(layer)
#jarqueberaTest(layer)
#boxplot(layer)
#-----Geoestatistica - escolhendo um melhor ajuste de modelo-----#
#-----Distancia Max e Min-----#
long_max = max(x)
lat_max = max(y)
mdta = sqrt(mx^2+my^2)
mdta
mdta/3
#-----Gera os diferentes semivariogramas-----#
#-----maxima distancia e numero de logs-----#
mdt = 130
nlg = 6
semi_variograma_teorico = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg)
#semi_variograma_teorico1 = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg, direction = pi/8)
#semi_variograma_teorico2 = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg, direction = pi/4)
#semi_variograma_teorico3 = variog(geo_dados_ndvi, max.dist = mdt, uvec = nlg, direction = pi/2)
h = semi_variograma_teorico$u
v = semi_variograma_teorico$v
npar = semi_variograma_teorico$n
ltmx = (max(h)+0.4*max(h))
ltmy = (max(v)+0.4*max(v))
ic = (1.96*(sqrt(2*v)/sqrt(npar)))
is = abs(v+ic)
ii = abs(v-ic)
#-----Ajuste dos semivariogramas teoricos-----#
esferico = variofit(semi_variograma_teorico, cov.model='sph', max.dist=max(h), messages=F)
#exponencial = variofit(semi_variograma_teorico, cov.model='exponential', max.dist=max(h), messages=F)
#gaussiano = variofit(semi_variograma_teorico, cov.model='gaussian', max.dist=max(h),messages=F)
#-----Separando os parametros dos modelos-----#
#-----Esferico-----#
c0.esf = esferico$nugget
c.esf = esferico$cov.pars[1]
a.esf = esferico$cov.pars[2]
#-----Exponencial-----#
#c0.exp = exponencial$nugget
#c.exp = exponencial$cov.pars[1]
#a.exp = exponencial$cov.pars[2]
#-----Gaussiano-----#
#c0.gau = gaussiano$nugget
#c.gau = gaussiano$cov.pars[1]
#a.gau = gaussiano$cov.pars[2]
#-----Gera o grid para receber os valores interpolados-----#
tamanho_do_pixel = 20 # Tamanho do pixel para interpolação
range_x = as.integer(range(x))
range_y = as.integer(range(y))
grid.map = expand.grid(x=seq(from=range_x[1], to=range_x[2], by = (tamanho_do_pixel)),
y=seq(from=range_y[1], to=range_y[2], by = (tamanho_do_pixel)))
#-----Krigagem-----#
semi_variograma_ajustado = esferico
krg_01 = krige.conv(geo_dados_ndvi, locations = grid.map, krige = krige.control(obj.model = semi_variograma_ajustado))
length(grid.map[,1])
fim = data.frame(grid.map$x, grid.map$y, krg_01$predict)
coordenadas = data.frame(grid.map$x, grid.map$y)
coordenadas = coordinates(coordenadas)
dados = data.frame(fim$krg_01.predict)
#write.table(fim, paste0(raiz, i, '_', 'teste_grid.txt'), sep = ',', row.names = F)
shp_final = SpatialPointsDataFrame(coords = coordenadas, data = dados)
shp_final = intersect(shp_final, shapefile_contorno_do_talhao_original)
salvar = "C:\\Users\\Yury\\Desktop\\InCeres\\Projeto_Fertilidade\\Sprint\\02"
writeOGR(shp_final,
dsn = salvar,
layer = paste0('talhao_',i, '_', tipo[j], '_NDVI_do_CACETE'),
driver = 'ESRI Shapefile',
overwrite_layer = TRUE)
detach(dados_ndvi)
}
}
Mais detalhes sobre a lista de discussão R-br