Cómo extraer el valor de un píxel a una geometría con PostGIS

Convertir entre el formato ráster y vectorial permite emplear ambos tipos de datos cuando nos encontramos inmersos en un proyecto SIG. Esto además permite emplear las diferentes herramientas de análisis específicas para cada uno de estos tipos de datos geográficos.

Pero, en ocasiones, puede que lo que nos interese es obtener los valores de los píxeles de un ráster en una geometría de tipo vectorial, bien sea un polígono, una línea o un punto, puesto que el resto no nos interesa y queremos desecharlos. Con PostGIS es posible extraer los valores de un ráster a una tabla con geometrías (archivo vectorial) almacenada en la base de datos.

NDVI transectos

En este tutorial, vamos a comentarte cómo conseguir esos valores explotando la información contenida, por ejemplo, en un ráster con el índice de vegetación NDVI calculado. Para ello, crearemos dos transectos para obtener los valores del NDVI en dos zonas diferenciadas, sobre una masa de agua y sobre un espacio de vegetación vigorosa.

Creación de una tabla

En primer lugar, y tras importar a la base de datos el dato ráster, abriríamos pgAdmin y nos conectaríamos a la base de datos que lo contiene para después abrir el constructor de consultas.

Constructor consultas

Creamos una tabla denominada, por ejemplo, transectos que tenga un campo que se llame cobertura con una longitud de 20 caracteres.

CREATE TABLE "transectos"

(gid serial PRIMARY KEY, 

"cobertura" varchar(20));

Create table

El siguiente paso es añadir la columna de geometría, SRID, tipo y dimensión:

SELECT AddGeometryColumn ('transectos', 'geom', '32629', 'LINESTRING', 2, false);

Add geometry

Introducción de geometrías

Ahora es el momento de introducir las geometrías que deseamos para ambos transectos:

INSERT INTO transectos (gid, geom, cobertura) VALUES (1, ST_GeomFromText('LINESTRING (499601.5257 4759846.214, 500341.0252 4760513.917)',32629), 'agua');

INSERT INTO transectos (gid, geom, cobertura) VALUES (2, ST_GeomFromText('LINESTRING (499578.4553 4762324.24, 499803.4774 4761983.772)',32629), 'cultivos');

Insert into

Como hemos comentado, vamos a extraer el dato de dos zonas diferenciadas, masas de agua y zona de cultivos. Una recomendación es que crees una conexión a PostGIS desde QGIS para visualizar los datos que almacenas en tus bases de datos.

Si no sabes como conectar con PostGIS desde QGIS te enseñamos en nuestro tutorial Cómo conectar con PostGIS desde QGIS donde, además, encontrarás un videotutorial.

transectos

Extracción del valor del píxel en el centroide

Una vez tenemos almacenadas las geometrías sobre las que queremos extraer la información del NDVI vamos a extraer el valor del píxel en el centroide de los mismos.

Es evidente que, para cruzar la información de ambas tablas en PostGIS requeriremos emplear la función ST_Intersect como herramienta de geoprocesamiento.

SELECT transectos.gid, transectos.cobertura,

ST_Value (ndvi_lugo_tm.rast, ST_Centroid (transectos.geom), false)

AS ndvi_value

FROM transectos, ndvi_lugo_tm

WHERE ST_Intersects (ndvi_lugo_tm.rast, transectos.geom);

El resultado sería el siguiente:

centroide

Extracción del valor de los píxeles a lo largo de la geometría

Aunque en nuestro ejemplo hemos calculado los valores ráster en el centroide del transecto es posible calcular dichos valores a lo largo del mismo, la consulta sería, por ejemplo:

SELECT tr.gid, tr.cobertura,

ST_Value (ndvi.rast, ST_Centroid ((ST_Intersection (ndvi.rast, tr.geom)).geom), false)

AS ndvi_value

FROM transectos AS tr, ndvi_lugo_tm AS ndvi

WHERE ST_Intersects (ndvi.rast, tr.geom) AND tr.cobertura = 'cultivos'

El resultado sería el siguiente:

linestring

Podríamos pensar que por qué es necesario obtener el centroide de la geometría que resulta de intersectar el transecto vectorial con la capa ráster. La razón sería que dicha intersección no es una serie de puntos, sino segmentos de línea. Para explicar esto, podríamos preguntar a PostGIS cuál es la geometría que resulta de intersectar el transecto vectorial y la capa ráster, para ello podemos emplear la siguiente consulta:

SELECT tr.gid, tr.cobertura,

ST_AsText ((ST_Intersection (ndvi.rast, tr.geom)).geom)

AS geometría,

ST_Value(ndvi.rast, ST_Centroid(

(ST_Intersection (ndvi.rast, tr.geom)).geom), false) AS valor_ndvi

FROM transectos AS tr, ndvi_lugo_tm AS ndvi

WHERE ST_Intersects (ndvi.rast, tr.geom) AND tr.cobertura = 'agua'

El resultado sería el siguiente:

consulta geometría

La explicación es que antes de realizar la intersección, PostGIS convierte las celdas del ráster en polígonos, por lo que al intersectar estos polígonos con el transecto se obtiene como resultado una serie de segmentos de línea.

Por otro lado, la función ST_Value retorna el valor de una banda específica en un punto geométrico, por lo que debe de pasarse como argumento una geometría de punto y no una de línea. Por tanto, en la consulta correspondiente es necesario obtener el centroide de cada segmento de línea que conforma la intersección.

Extracción del valor de los píxeles por bandas

Por otro lado, también podemos obtener los valores ráster en una serie de bandas de nuestro interés. Por ejemplo, vamos a conocer los valores en las bandas 4 y 5 asociados al centroide de nuestros transectos:

SELECT transectos.gid, transectos.cobertura,

ST_Value (lugo_tm.rast, 4, ST_Centroid (transectos.geom), true) AS banda4,

ST_Value (lugo_tm.rast, 5, ST_Centroid (transectos.geom), true) AS banda5

FROM transectos, lugo_tm

WHERE ST_Intersects (lugo_tm.rast, transectos.geom);

El resultado sería el siguiente:

valor por bandas

Si quieres aprender a trabajar con Bases de Datos Espaciales de la mano de un tutor cualificado, inscribete ya en nuestro curso online en Bases de Datos Espaciales: PostGIS.

3 comentarios en «Cómo extraer el valor de un píxel a una geometría con PostGIS»

  1. Hola ALC,

    Supongo que te refieres a convertir de formato JSON a shapefile (.shp). Esto puedes hacerlo fácilmente con QGIS, cargando el JSON en su interfaz y guardándolo después desde sus propiedades (botón derecho sobre el JSON/opción Guardar como….) con el tipo .shp

    Un saludo!

    • Gracias por tu respuesta. En realidad es un servicio Json, desde una web la cual contiene la latitud y longitud de vehículos, y cada cierto tiempo se esta actualizando ese Json. La duda es, cual será el mejor metodo para pasar esa latitud y longitud del servicio Json a un mapa, como OpenStreet. Estuve leyendo un poco y existen los marcadores de OpenLayers, pero no me queda claro, no hay un ejemplo básico para entender.

  2. Hola, quisiera preguntar si existe alguna forma o método explicado de pasar de un json que muestra medios de transporte en tiempo real, a un mapa, el json contiene las latitudes, longitudes.

    Muchas Gracias

Los comentarios están cerrados.