Exemplo: chuva no Paraná e Krigagem
Exemplo: Banco de Dados de Chuva no PR
Banco existente no SGBD
(banco geográfico no formato TerraLib)
1. Conectando banco, visualizando no R e criando vistasno TV
Fazer no TV:
- abrir o terraView e conectar ao banco "chuvaPR"
- se existir alguma vista já criada, visualize!
require(aRT)
Estabelecendo conexão com o DBMS (SGBD)
con <- openConn() con
Listando os bancos existentes no DBMS
showDbs(con)
Conectando a um banco da dados
pr <- openDb(con, "chuvaPR", update=T)
Visualizando informações deste banco
pr
Abrindo o layer com contorno do estado
pr.l_contorno <- openLayer(pr, "contorno") pr.l_contorno visualizando dado do layer plot(pr.l_contorno)
Abrindo outro layer, das estações meteorológicas e dados de precipitação
pr.l_dados <- openLayer(pr, "dados") pr.l_dados
Visualizando, sobrepondo ao gráfico anterior
plot(pr.l_dados, add=T)
Apagando uma "vista" pré-existente do banco
deleteView(pr, "Parana") type Y at the R prompt
Fazer no TV:
- reconectar o banco o TV e ver que a vista "sumiu"!!!
Criando uma vista/tema com o contorno …
tema_pol <- createTheme(pr.l_contorno, "contorno", view="Parana") setVisual(tema_pol, visualPolygons())
… e adicionando os pontos à vista já criada
tema_dados <- createTheme(pr.l_dados, "estacoes", view = "Parana") setVisual(tema_dados, visualPoints(pch=21, color="black", size=5))
Fazer no TV:
- reconectar o banco o TV e ver que "apareceu" a vista agora chamada "parana" !!!
2. Operações da TL via aRT
Ilustrando como (algumas) operações da terraLib podem ser usadas
getProj(pr.l_contorno) getMetric(pr.l_contorno, "area") getMetric(pr.l_contorno, "length") getDistance(pr.l_contorno, id=as.character(c(1,2)), pr.l_dados) getDistance(pr.l_contorno, id=as.character(c(1,3)), pr.l_dados) getDistance(pr.l_contorno, id=as.character(c(2,3)), pr.l_dados)
3. Trazendo e manipulando dados e fazendo análises no R
Abrindo o tema e trazendo dados para um objeto do R
Duas possibilidades:
- via "theme"
pr th_dados <- openTheme(pr, "estacoes") th_dados dados <- getData(th_dados) head(dados)
- via tabela do banco
pr.l_dados tb_dados <- openTable(pr.l_dados, "t_dados") tb_dados dat <- getData(tb_dados) dat[1:5,] pts <- getPoints(pr.l_dados) dat[1:5,]
Trazendo o polígono com a borda do espaço para um objeto do R
pol <- getPolygons(pr.l_contorno) pol plot(pol)
Fazendo agora uma interpolação dos dados de chuva via geoestatística/krigagem
- Carregando o pacote geoR
require(geoR)
- Convertendo dados os dados para o formato "geodata" da geoR (conveniente)
geo <- as.geodata(dados, data.col=2) geo$borders <- pol@polygons[[1]]@Polygons[[1]]@coords plot(geo)
- Estimando parâmetros
ml <- likfit(geo, trend="1st", ini=c(1000, 100)) um pouco demorado definindo e visualizando grid de predição loc0 <- pred_grid(geo$borders, by=10) points(geo, bor=borders) points(loc0, pch=".", col=2) veja gráfico!
- e fazendo krigagem
kc <- krige.conv(geo, loc=loc0, krige=krige.control(obj=ml), bor=geo$borders)
- visualizando no R
image(kc, col=terrain.colors(15), coords=parana$coords) ver gráfico
Preparing the predictions for aRT to transfer as a raster to the DBMS
georpred <- .prepare.graph.kriging(locations=loc0, borders=parana$borders, values=kc$pred) names(georpred)[3] <- "z"
Criando um novo layer no banco para armazenar o grid de predição
pr.l_pred <- createLayer(pr, "predicao") addRaster(pr.l_pred, georpred)
Checking the current status of the DB
pr
Fazer no TV:
- reconectar o banco o TV e ver que um novo layer foi criado
Optional : setting to a TV view
th <- createTheme(pr.l_pred, "raster", view="parana") setVisual(th, visualRaster(color = terrain.colors(15)), mode="r") pr
Fazer no TV:
- reconectar o banco o TV e ver que um novo tema foi adicionado a vista
- visualize!!!!
4. Adicionando uma nova coluna (com predições) à tabela de dados
Para exemplificar isto vamos fazer a krigagem nos pontos onde há dados
kc0 <- krige.conv(geo, loc=geo$coords, krige=krige.control(obj=ml), borders=NULL) prs <- data.frame(id=rownames(geo$coords), pred=kc0$pred) createColumn(tb_dados, "preds", type="numeric") updateColumns(tb_dados, prs)
veja o novo status da tabela
tb_dados
Fazer no TV:
- reconectar o banco e ver que há agora uma nova tabela no layer de dados
5. Some plots
Plotting directly from the data-base (not using R objects)
plot(pr.l_pred, col=terrain.colors(15)) plot(pr.l_dados, add=T) plot(pr.l_contorno, add=T) veja a figura!
Deletando um layer
deleteLayer(con, "pr.l_pred")