22 Oct 2012 Mapas de densidad con R y PostgreSQL/PostGIS

En este artículo vamos a ver cómo crear con SQL, mapas de densidad a partir de datos almacenados en una base de datos PostgreSQL/PostGIS . Para ello vamos a utilizar las funciones de R (http://www.r-project.org/ ) mediante el lenguaje PL/R. Este lenguaje soportado por PostgreSQL permite la creación de stored procedures dentro de nuestras bases de datos utilizando las funciones estadísticas de R.

Software

Para la elaboración de este mapa de densidad hemos utilizado un PC con windows 7 32 bit. Las versiones de software utilizadas son:

Postgresql 9.1.6 + PostGIS 2.0. + R 2.15.1

Una vez instalado R y siguiendo las indicaciones que encontraréis aquí (http://www.bostongis.com/PrinterFriendly.aspx?content_name=postgresql_plr_tut01 ) definimos una nueva variable de entorno (R_HOME) en nuestro sistema operativo:

R_HOME debe apuntar a la carpeta de la instalación de R . En nuestro caso C:Program FilesRR-2.15.1bin

A continuación modificamos la variable de entorno PATH para añadir la ruta a la carpeta donde se encuentran los archivos binarios de R. En nuestro caso C:Program FilesRR-2.15.1bini386 (versiones de R anteriores a la 2.12 deben apuntar directamente a la carpeta bin)

Una vez instalado R debemos instalar la librería spatstat junto con sus librerías dependientes. Esta librería contiene las funciones necesarias para la creación de mapas de densidad.

Desde una consola de R tecleamos: install.packages(“spatstat”, dependencies = TRUE)

Instalación de PL/R

Para la instalación de PL/R hemos descargado el fichero que encontraréis en http://www.joeconway.com/plr/ correspondiente a la versión de PostgreSQL 9.1 (a fecha de hoy PL/R está disponible hasta la versión 9.1 de PostgreSQL). Este fichero contiene dos ficheros importantes: plr.dll y plr.sql

Copiamos el fichero plr.dll dentro de la carpeta lib que aparece en la carpeta de instalación de PostgreSQL. En nuestro ejemplo C:Archivos de ProgramaPostgreSQLlib

Finalmente copiamos el fichero plr.sql a la carpeta que resulte de nuestro agrado. Por ejemplo, C:Archivos de ProgramaPostgreSQLsharecontrib

Datos.

Para la elaboración del mapa de densidad hemos utilizado la siguiente tabla.

Create table flickr
“user” character varying(15),
url character varying(255),
the_geom geometry
);

Esta tabla contiene información de las imágenes de la plataforma flickr geolocalizadas en Barcelona. La siguiente imagen muestra su disposición en el espacio.

Función PL/R.

Utilizando el lenguaje PL/R, definimos la función f_density como:

CREATE OR REPLACE FUNCTION f_density(text, text, text, float)
RETURNS text AS
$BODY$
–La función f_density recibe 4 parámetros. La columna que contiene la geometría, la tabla que contiene la geometría, la ruta a la carpeta donde se creará la imagen y el grado de densidad que queremos aplicar.

–Cargamos las funciones necesarias
library(spatstat)

–Comprobamos que las geometrías son de tipo POINT y que aparecen en la tabla geometry_columns

geometry <- pg.spi.exec(sprintf(«select type from geometry_columns where f_table_name=’%2$s’ and f_geometry_column=’%1$s'» ,arg1, arg2))
if (geometry$type != ‘POINT’)
{
return(‘Invalid geometry type’)
}

–Descomponemos las coordenadas de los puntos en las variables latitude y longitude
dataset <- pg.spi.exec(sprintf(‘select st_x(%1$s) as x, st_y(%1$s) as y from %2$s’,arg1, arg2))
latitude <- dataset$x
longitude <- dataset$y

–Calculamos los rangos de cada variable para definir un ppp (planar point pattern) sobre el que poder calcular la función density
xrange <-c(min(longitude), max(longitude))
yrange <-c(min(latitude), max(latitude))
d <- ppp(latitude, longitude, yrange, xrange)

–Creamos nuestra paleta. Valores altos del parámetro bias otorgan rangos mayores para los valores más altos
Lab.palette <- colorRampPalette(c(«white», «yellow», «orange», «#ff4100», «redo»), bias=2)

–Preparamos la salida en formato png
png(arg3, width=500,height=500)

–Finalmente creamos nuestro mapa de densidad sobre el ppp utilizando la paleta anterior.
plot(density(d,arg4), col=Lab.palette(256))

–Cerramos el fichero
dev.off()

–Devolvemos un mensaje de salida
print(«done»)
$BODY$
LANGUAGE plr VOLATILE COST 100;

En este momento ya podemos crear nuestros mapas de densidad ejecuntando la siguiente consulta SQL:

select f_density(‘the_geom’,’flickr’,’c:/temp/mapadensidad.png’, 0.002);

Y obtendremos la siguiente imagen: