Predição lenta - Interpolação Kriging
Bom dia pessoal, Espero que possam me ajudar. Entro em contato para pedir uma ajuda em uma interpolação que estou tentando fazer. Já usei o mesmos comandos abaixo para realizar interpolações de variáveis de um lago por ordinary kriging com 240 pontos, sem problemas. Agora estou tentando interpolar a batimetria do lago, para tanto criei pontos no limite do lago com valor zero, cerca de 1200 pontos. Usando a mesma metodologia, a predição “predict( )" para o grid passa de 3 horas e não tive paciência para esperar mais, visto que as outras vezes que usei a predição para o grid demorou cerca de 2 minutos para 240 pontos interpolados. Imagino que a quantidade de pontos seja o problema. Alguém teria alguma sugestão de caminho que tivesse menos custo computacional, pois acho que a quantidade de pontos seja o problema. Agradeço desde já a atenção. Segue abaixo os comandos utilizados:
library(geoR) library(raster) library(gstat) library(rgdal) interpol_pcsand<-data_interpol[,c(1,2,6)]
coordinates(interpol_pcsand)<- ~long+lat
#setting the CRS for the geospatial data, same as the mask and grid
projection(interpol_pcsand)<-'+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs '
#search for duplicated coords
dup.coords(interpol_pcsand$coords)
#Creating a gstat object
g_pcsand <- gstat(id="pcsand", formula=pcsand ~ 1, data=interpol_pcsand)
# variogram with eye ball fitting
v.eye_sand <- eyefit(variog(as.geodata(interpol_pcsand["pcsand"]), max.dist = 20000))
#saving a readable file for gstat
ve.fit <- as.vgm.variomodel(v.eye[[1]])
#updating the gstat object
g_pcsand <- gstat(g_pcsand, id="pcsand", model=ve.fit_pcsand )
#predicting to new data # grid4 diemensions - #features : 1630464 #extent : 468333.3, 499983.3, 6634933, 6681223 (xmin, xmax, ymin, ymax) #coord. ref. : +proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs
p_pcsand <- predict(g_pcsand, model=ve.fit_pcsand , newdata=grid4)
#write a raster with the interpolated matrix rasterpcsand_pred <- rasterFromXYZ(as.data.frame(p_pcsand)[, c("x", "y", "data.pred")]) #Create raster tif file
writeRaster(raster, filename="raster.tif", overwrite=TRUE)
_________________________________ Thiago Cesar Lima Silveira PhD candidate - Zoology PUCRS - Brazil Visitor PhD. student - School of Ocean Sciences - Bangor University Building Westbury Mouth Room Nautilus 327 Menai Bridge - Isle of Anglesey - UK Telephone.: 07843 115244 e-mail: thiago.cesar@acad.pucrs.br CV: http://lattes.cnpq.br/5960267776845701
Ir de 240 para 1200 pontos de preodicao nao deve causar toda esta diferenca no tempo de predicao se voce estiver simplesmente obtendo media e variancia da preditiva. Mesmo simulando de preditiva eu nao esperaria este tempo todo nao tenho familiaridade com alguns objetos que voce está usando mas talvez valha a pena isoldar dados em formatos mais simples e/ou rodar rotinas alternativas de predicao em outros pacotes On Tue, 17 Jun 2014, Thiago Cesar Lima Silveira wrote:
Bom dia pessoal,Espero que possam me ajudar.
Entro em contato para pedir uma ajuda em uma interpolação que estou tentando fazer. Já usei o mesmos comandos abaixo para realizar interpolações de variáveis de um lago por ordinary kriging com 240 pontos, sem problemas.
Agora estou tentando interpolar a batimetria do lago, para tanto criei pontos no limite do lago com valor zero, cerca de 1200 pontos. Usando a mesma metodologia, a predição “predict( )" para o grid passa de 3 horas e não tive paciência para esperar mais, visto que as outras vezes que usei a predição para o grid demorou cerca de 2 minutos para 240 pontos interpolados. Imagino que a quantidade de pontos seja o problema.
Alguém teria alguma sugestão de caminho que tivesse menos custo computacional, pois acho que a quantidade de pontos seja o problema.
Agradeço desde já a atenção.
Segue abaixo os comandos utilizados:
library(geoR) library(raster) library(gstat) library(rgdal) interpol_pcsand<-data_interpol[,c(1,2,6)]
coordinates(interpol_pcsand)<- ~long+lat
#setting the CRS for the geospatial data, same as the mask and grid
projection(interpol_pcsand)<-'+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs '
#search for duplicated coords
dup.coords(interpol_pcsand$coords)
#Creating a gstat object
g_pcsand <- gstat(id="pcsand", formula=pcsand ~ 1, data=interpol_pcsand)
# variogram with eye ball fitting
v.eye_sand <- eyefit(variog(as.geodata(interpol_pcsand["pcsand"]), max.dist = 20000))
#saving a readable file for gstat
ve.fit <- as.vgm.variomodel(v.eye[[1]])
#updating the gstat object
g_pcsand <- gstat(g_pcsand, id="pcsand", model=ve.fit_pcsand )
#predicting to new data # grid4 diemensions - #features : 1630464 #extent : 468333.3, 499983.3, 6634933, 6681223 (xmin, xmax, ymin, ymax) #coord. ref. : +proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs
p_pcsand <- predict(g_pcsand, model=ve.fit_pcsand , newdata=grid4)
#write a raster with the interpolated matrix rasterpcsand_pred <- rasterFromXYZ(as.data.frame(p_pcsand)[, c("x", "y", "data.pred")])
#Create raster tif file
writeRaster(raster, filename="raster.tif", overwrite=TRUE)
_________________________________ Thiago Cesar Lima Silveira PhD candidate - Zoology PUCRS - Brazil Visitor PhD. student - School of Ocean Sciences - Bangor University Building Westbury Mouth Room Nautilus 327 Menai Bridge - Isle of Anglesey - UK Telephone.: 07843 115244 e-mail: thiago.cesar@acad.pucrs.br CV: http://lattes.cnpq.br/5960267776845701
será que você está fazendo global kriging? não vi no seus comandos um limite para número máximo de observações na vizinhança. Em 17 de junho de 2014 09:22, Thiago Cesar Lima Silveira < thiagoclsilveira@yahoo.com.br> escreveu:
Bom dia pessoal, Espero que possam me ajudar.
Entro em contato para pedir uma ajuda em uma interpolação que estou tentando fazer. Já usei o mesmos comandos abaixo para realizar interpolações de variáveis de um lago por ordinary kriging com 240 pontos, sem problemas.
Agora estou tentando interpolar a batimetria do lago, para tanto criei pontos no limite do lago com valor zero, cerca de 1200 pontos. Usando a mesma metodologia, a predição “predict( )" para o grid passa de 3 horas e não tive paciência para esperar mais, visto que as outras vezes que usei a predição para o grid demorou cerca de 2 minutos para 240 pontos interpolados. Imagino que a quantidade de pontos seja o problema.
Alguém teria alguma sugestão de caminho que tivesse menos custo computacional, pois acho que a quantidade de pontos seja o problema.
Agradeço desde já a atenção.
Segue abaixo os comandos utilizados:
library(geoR) library(raster) library(gstat) library(rgdal) interpol_pcsand<-data_interpol[,c(1,2,6)]
coordinates(interpol_pcsand)<- ~long+lat
#setting the CRS for the geospatial data, same as the mask and grid
projection(interpol_pcsand)<-'+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs '
#search for duplicated coords
dup.coords(interpol_pcsand$coords)
#Creating a gstat object
g_pcsand <- gstat(id="pcsand", formula=pcsand ~ 1, data=interpol_pcsand)
# variogram with eye ball fitting
v.eye_sand <- eyefit(variog(as.geodata(interpol_pcsand["pcsand"]), max.dist = 20000))
#saving a readable file for gstat
ve.fit <- as.vgm.variomodel(v.eye[[1]])
#updating the gstat object
g_pcsand <- gstat(g_pcsand, id="pcsand", model=ve.fit_pcsand )
*#predicting to new data* *# grid4 diemensions - * *#features : 1630464 * *#extent : 468333.3, 499983.3, 6634933, 6681223 (xmin, xmax, ymin, ymax)* *#coord. ref. : +proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs * *>p_pcsand <- predict(g_pcsand, model=ve.fit_pcsand , newdata=grid4) *
#write a raster with the interpolated matrix rasterpcsand_pred <- rasterFromXYZ(as.data.frame(p_pcsand)[, c("x", "y", "data.pred")])
#Create raster tif file
writeRaster(raster, filename="raster.tif", overwrite=TRUE)
_________________________________ Thiago Cesar Lima Silveira PhD candidate - Zoology PUCRS - Brazil Visitor PhD. student - School of Ocean Sciences - Bangor University Building Westbury Mouth Room Nautilus 327 Menai Bridge - Isle of Anglesey - UK Telephone.: 07843 115244 e-mail: thiago.cesar@acad.pucrs.br CV: http://lattes.cnpq.br/5960267776845701
_______________________________________________ R-br mailing list R-br@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.
Thiago, boa tarde! Me disponho a dar uma olhada, mas vou precisar dos objetos interpol_pcsand (ou data_interpol) e grid4. Me chamou a atenção o fato de aparecer a menção de 1.630.464 features no seu objeto grid4. Se for isso mesmo são mais de 1200x1200pontos ao invés de 1200 pontos como mencionado. Além disso, talvez você possa utilizar IDW para batimetria. A não ser que você esteja avaliando os parâmetros do modelo ajustado. Atte., Éder Comunello <c <comunello.eder@gmail.com>omunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
Olá Éder, Obrigado pela disponibilidade! Quando mencionei 1200 pontos, me referi aos pontos reais de batimetria amostrados. O grid4 que me refiro é isso mesmo: Features : 1630464 extent : 468333.3, 499983.3, 6634933, 6681223 (xmin, xmax, ymin, ymax) Resolution: 30 x 30 Não cheguei a usar o IDW para batimetria. É a primeira vez que tento construir mapas interpolados no R e a minha intenção é desenvolver uma linha de trabalho para gerar os mapas e avaliar a qualidade da interpolação. Segue abaixo como construí o gri4, os dados da batimetria te mando para o teu e-mail pessoal. ############# # Grid4
grid <- raster( ) grid <- raster(ncol=1544, nrow=1056, xmn=468318.3, xmx=499998.3, ymn=6634918, ymx=6681238) projection (grid) <- '+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs ‘ res(grid) <- 30 grid4 <- as(grid, ‘SpatialPoints’)
Grande abraço, Thiago _________________________________ Thiago Cesar Lima Silveira PhD candidate - Zoology PUCRS - Brazil Visitor PhD. student - School of Ocean Sciences - Bangor University Building Westbury Mouth Room Nautilus 327 Menai Bridge - Isle of Anglesey - UK Telephone.: 07843 115244 e-mail: thiago.cesar@acad.pucrs.br CV: http://lattes.cnpq.br/5960267776845701 On Jun 18, 2014, at 7:39 PM, Éder Comunello <comunello.eder@gmail.com> wrote:
Thiago, boa tarde!
Me disponho a dar uma olhada, mas vou precisar dos objetos interpol_pcsand (ou data_interpol) e grid4. Me chamou a atenção o fato de aparecer a menção de 1.630.464 features no seu objeto grid4. Se for isso mesmo são mais de 1200x1200pontos ao invés de 1200 pontos como mencionado.
Além disso, talvez você possa utilizar IDW para batimetria. A não ser que você esteja avaliando os parâmetros do modelo ajustado.
Atte.,
Éder Comunello <comunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
_______________________________________________ R-br mailing list R-br@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.
Thiago, bom dia! Olhei seus dados e pude notar que na verdade você tem 236 pontos cotados. Os 1263 restantes são na verdade da linha de borda, transformados em pontos e que receberam o valor zero. Na figura anexa, eu isolei os pontos cotados e usei os demais pra restaurar a linha de borda. O uso de linhas de borda (um tipo de breaklines) é uma prática comum quando você utiliza grades triangulares (TIN) pra modelar tridimensionalmente uma área. A partir da grade triangular você interpola por outro método pra obter uma grade regular. Talvez seja o procedimento mais adequado pro seu caso. Ao optar por outro método que não o uso de TIN, eu particularmente, nunca tinha visto adicionar os pontos da linha de borda aos dados para interpolar e, num primeiro momento, acredito que seria melhor interpolar os 236 pontos e utilizar a linha de borda apenas pra recortar a superfície gerada. Mas precisaria estudar melhor sobre a questão. Ainda pretendo seguir a análise do seu código, mas queria adiantar esses pontos. Atte., Éder Comunello <c <comunello.eder@gmail.com>omunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W] Em 19 de junho de 2014 04:39, Thiago Cesar Lima Silveira < thiagoclsilveira@yahoo.com.br> escreveu:
Olá Éder,
Obrigado pela disponibilidade!
Quando mencionei 1200 pontos, me referi aos pontos reais de batimetria amostrados.
O grid4 que me refiro é isso mesmo: Features : 1630464
extent : 468333.3, 499983.3, 6634933, 6681223 (xmin, xmax, ymin, ymax) Resolution: 30 x 30
Não cheguei a usar o IDW para batimetria.
É a primeira vez que tento construir mapas interpolados no R e a minha intenção é desenvolver uma linha de trabalho para gerar os mapas e avaliar a qualidade da interpolação.
Segue abaixo como construí o gri4, os dados da batimetria te mando para o teu e-mail pessoal.
############# # Grid4
grid <- raster( ) grid <- raster(ncol=1544, nrow=1056, xmn=468318.3, xmx=499998.3, ymn=6634918, ymx=6681238) projection (grid) <- '+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs ' res(grid) <- 30 grid4 <- as(grid, 'SpatialPoints')
Grande abraço,
Thiago
_________________________________ Thiago Cesar Lima Silveira PhD candidate - Zoology PUCRS - Brazil Visitor PhD. student - School of Ocean Sciences - Bangor University Building Westbury Mouth Room Nautilus 327 Menai Bridge - Isle of Anglesey - UK Telephone.: 07843 115244 e-mail: thiago.cesar@acad.pucrs.br CV: http://lattes.cnpq.br/5960267776845701
On Jun 18, 2014, at 7:39 PM, Éder Comunello <comunello.eder@gmail.com> wrote:
Thiago, boa tarde!
Me disponho a dar uma olhada, mas vou precisar dos objetos interpol_pcsand (ou data_interpol) e grid4. Me chamou a atenção o fato de aparecer a menção de 1.630.464 features no seu objeto grid4. Se for isso mesmo são mais de 1200x1200pontos ao invés de 1200 pontos como mencionado.
Além disso, talvez você possa utilizar IDW para batimetria. A não ser que você esteja avaliando os parâmetros do modelo ajustado.
Atte.,
Éder Comunello <c <comunello.eder@gmail.com>omunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
_______________________________________________ R-br mailing list R-br@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.
_______________________________________________ R-br mailing list R-br@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.
Bom dia Éder, Obrigado pelo retorno. Estou estudando as metodologias de interpolação não havia pensado na alternativa que mencionaste a respeito do TIN. Somente me pareceu lógico criar uma linha de borda com uma informação que sei que é zero para a batimetria. A partir disto interpolar por kriging os valores. No entanto, confesso que preciso de um caminho para construir um mapa e que também consiga avaliar a qualidade da interpolação. Estou estudando pelos exemplos apresentados pelo Bivand (2008) - Applied Spatial Data Analysis with R. Grande abraço e obrigado mais uma vez. Thiago On Jun 20, 2014, at 12:39 PM, Éder Comunello <comunello.eder@gmail.com> wrote:
Thiago, bom dia!
Olhei seus dados e pude notar que na verdade você tem 236 pontos cotados. Os 1263 restantes são na verdade da linha de borda, transformados em pontos e que receberam o valor zero. Na figura anexa, eu isolei os pontos cotados e usei os demais pra restaurar a linha de borda.
O uso de linhas de borda (um tipo de breaklines) é uma prática comum quando você utiliza grades triangulares (TIN) pra modelar tridimensionalmente uma área. A partir da grade triangular você interpola por outro método pra obter uma grade regular. Talvez seja o procedimento mais adequado pro seu caso.
Ao optar por outro método que não o uso de TIN, eu particularmente, nunca tinha visto adicionar os pontos da linha de borda aos dados para interpolar e, num primeiro momento, acredito que seria melhor interpolar os 236 pontos e utilizar a linha de borda apenas pra recortar a superfície gerada. Mas precisaria estudar melhor sobre a questão.
Ainda pretendo seguir a análise do seu código, mas queria adiantar esses pontos.
Atte.,
Éder Comunello <comunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
Em 19 de junho de 2014 04:39, Thiago Cesar Lima Silveira <thiagoclsilveira@yahoo.com.br> escreveu: Olá Éder,
Obrigado pela disponibilidade!
Quando mencionei 1200 pontos, me referi aos pontos reais de batimetria amostrados.
O grid4 que me refiro é isso mesmo: Features : 1630464
extent : 468333.3, 499983.3, 6634933, 6681223 (xmin, xmax, ymin, ymax) Resolution: 30 x 30
Não cheguei a usar o IDW para batimetria.
É a primeira vez que tento construir mapas interpolados no R e a minha intenção é desenvolver uma linha de trabalho para gerar os mapas e avaliar a qualidade da interpolação.
Segue abaixo como construí o gri4, os dados da batimetria te mando para o teu e-mail pessoal.
############# # Grid4
grid <- raster( ) grid <- raster(ncol=1544, nrow=1056, xmn=468318.3, xmx=499998.3, ymn=6634918, ymx=6681238) projection (grid) <- '+proj=utm +zone=22 +south +ellps=aust_SA +units=m +no_defs ‘ res(grid) <- 30 grid4 <- as(grid, ‘SpatialPoints’)
Grande abraço,
Thiago
_________________________________ Thiago Cesar Lima Silveira PhD candidate - Zoology PUCRS - Brazil Visitor PhD. student - School of Ocean Sciences - Bangor University Building Westbury Mouth Room Nautilus 327 Menai Bridge - Isle of Anglesey - UK Telephone.: 07843 115244 e-mail: thiago.cesar@acad.pucrs.br CV: http://lattes.cnpq.br/5960267776845701
On Jun 18, 2014, at 7:39 PM, Éder Comunello <comunello.eder@gmail.com> wrote:
Thiago, boa tarde!
Me disponho a dar uma olhada, mas vou precisar dos objetos interpol_pcsand (ou data_interpol) e grid4. Me chamou a atenção o fato de aparecer a menção de 1.630.464 features no seu objeto grid4. Se for isso mesmo são mais de 1200x1200pontos ao invés de 1200 pontos como mencionado.
Além disso, talvez você possa utilizar IDW para batimetria. A não ser que você esteja avaliando os parâmetros do modelo ajustado.
Atte.,
Éder Comunello <comunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
_______________________________________________ R-br mailing list R-br@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.
_______________________________________________ R-br mailing list R-br@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.
_______________________________________________ R-br mailing list R-br@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.
participantes (4)
-
Luis Paulo Braga -
Paulo Justiniano -
Thiago Cesar Lima Silveira -
Éder Comunello