Integração do PostgreSQL e R para Estimador de Densidade

Começando a aprender o versátil e completo software R, hoje consegui fazer uma pequena e simples rotina que quero compartilhar. Pessoalmente, escolhi o RStudio como interface (quase) gráfica para o R.

O objetivo é fazer de forma rápida mapas de densidade de ocorrência de um determinado evento.

É bem provável que eu não utilize um mapa como o demonstrado acima como produto final. Mas se apresentou como uma solução boa, bonita e barata para uma análise exploratória dos dados. Principalmente se eu quero fazer isto para 5 espécies e 295 saídas de campo. Não dá para fazer isto clicando no QGIS. O construtor de modelos do GRASS pode ser uma saída. Mas certamente esta solução no R foi a mais rápida.

Neste exemplo, utilizaremos o caso encalhes de tartarugas marinhas pelo litoral do Rio Grande do Sul. Nestes links, você pode ter mais informações sobre o Projeto Tartarugas no Mar e sobre o NEMA, instituição responsável pelo banco de dados de monitoramento de praia aqui utilizado.

Este procedimento será realizado em 7 passos:

1. Carregar as bibliotecas necessárias;
2. Conectar ao banco de dados;
3. Criar as variáveis que serão utilizadas;
4. Calcular o índice de densidade;
5. Montar o raster de densidade;
6. Plotar os dados e salvar a figura;
7. Salvar o raster de densidade em GeoTIFF para abrir em outro software de SIG (QGIS, GRASS, GvSIG ou qualquer outro).

Então mãos à obra. Como dados de entrada, precisamos apenas de uma série de pontos (longitude e latitude) e de um shapefile para representar a divisão política e facilitar a compreensão dos resultados.

1. Bibliotecas:

library(sqldf) ##execução de sql

library(RPostgreSQL) ##conexão com Postgres

library(KernSmooth)

library(raster) ##construção de rasters

library(maps) ##importação de shapes e muitas outras funções de mapas

2.  Conexão com o PostgreSQL:

options(sqldf.RPostgreSQL.user="seu_usuario",

        sqldf.RPostgreSQL.password="sua_senha",

        sqldf.RPostgreSQL.dbname ="nome_do_bd",

        sqldf.RPostgreSQL.host ="localhost_ou_IP_do_servidor",

        sqldf.RPostgreSQL.port =5432)

3.  Criação das variáveis:

O sqldf permite você fazer qualquer consulta em SQL dentro do BD. Neste caso, selecionei apenas as colunas lat e lon da tabela de encalhes, para os registros do ano de 2013.

encalhes<- sqldf("select lon,lat from encalhes_2002p where extract ('year' from data)= 2013")

No caso da linha de costa, iremos carregar um shapefile:

costa <-shapefile('../divisao_politica.shp', stringsAsFactors=FALSE, verbose=FALSE) 

4. Cálculo do Índice de Densidade Kernel

dens <- bkde2D(encalhes, # nome da variável com lat, lon

              bandwidth=c(0.05,0.05), #maior distância de influência de cada ponto, neste caso em graus 

              gridsize=c(2000,2000), #número de linhas e colunas (resolução) do arquivo gerado

              range.x=list(c(-53.5,-50),c(-33.8,-31.2))) #coordenadas mínimas e máximas do cálculo

5. Transformando em raster

O caso é que o comando bkde2d gera uma lista com lon, lat e índice de densidade calculado. Portanto é necessário criar um raster, que irá ocupar menos memória e será mais efetivo.

dens.raster = raster(list(x=dens$x1, y=dens$x2, z=dens$fhat))

projection(dens.raster) <- CRS("+init=epsg:4326")# Sistema de Referência

xmin(dens.raster) <- -53.5 #Coordenadas do retângulo de extensão do raster

xmax(dens.raster) <- -50

ymin(dens.raster) <- -33.8

ymax(dens.raster) <- -31.2

6. Fazendo o plot:

 plot(dens.raster)

plot(costa, bg="transparent", add=TRUE)

Note que o parâmetro "bg" define que o background dos polígonos será transparente e o "add" que este plot deverá ser adicionado ao anterior, e não simplesmente gerado em uma nova figura

7. Se quiser salvar como figura, basta adicionar 1 linha antes e uma depois:

png("../heatmap_total.png", width=1000, height=1000)

plot(dens.raster)

plot(costa, bg="transparent", add=TRUE)

dev.off()

8. Se quiser salvar o kernel como GeoTIFF para utilizar em outros programas, basta:

writeRaster(../kern.tif',"GTiff",overwrite=TRUE) 

Extra: Agora, com todos os mapinhas padronizados, fica muito fácil fazer uma animação tratando cada mapa como um frame, e montando um gif. Um comando no terminal do Ubuntu:

convert -delay 100 -loop 0 heatmap*.jpg animation.gif