30. Rasters¶
PostGIS soporta otro tipo de dato espacial llamado raster. Los datos raster, al igual que los datos de geometría, usan coordenadas cartesianas y un sistema de referencia espacial. Sin embargo, en lugar de datos vectoriales, los datos raster se representan como una matriz n-dimensional compuesta de píxeles y bandas. Las bandas definen el número de matrices que tienes. Cada píxel almacena un valor correspondiente a cada banda. Así, un raster de 3 bandas como una imagen RGB, tendría 3 valores para cada píxel correspondientes a las bandas Rojo-Verde-Azul.
Aunque imágenes bonitas como las que ves en la pantalla de tu televisor son raster, los raster pueden no ser tan emocionantes de ver. En pocas palabras, un raster es una matriz, fijada en un sistema de coordenadas, que tiene valores que pueden representar lo que quieras que representen.
Dado que los raster viven en el espacio cartesiano, los raster pueden interactuar con geometrías. PostGIS ofrece muchas funciones que toman como entrada tanto raster como geometrías. Muchas operaciones aplicadas a raster resultarán en geometrías. Algunas comunes son ST_Polygon, ST_Envelope, ST_ConvexHull, y ST_MinConvexHull como se muestra a continuación. Cuando conviertes un raster en una geometría, lo que se genera es el ST_ConvexHull del raster.
El formato raster se usa comúnmente para almacenar datos de elevación, datos de temperatura, datos satelitales y datos temáticos que representan cosas como contaminación ambiental, densidad poblacional y ocurrencias de riesgos ambientales. Puedes usar raster para almacenar cualquier dato numérico que tenga una ubicación de coordenadas significativa. La única restricción es que para todos los datos en una banda específica los tipos de datos numéricos deben ser los mismos.
Aunque los datos raster pueden crearse desde cero en PostGIS, un enfoque más común es cargar datos raster desde varios formatos usando la herramienta de línea de comandos raster2pgsql incluida con PostGIS. Antes de todo eso, debes habilitar el soporte raster en tu base de datos ejecutando el comando:
CREATE EXTENSION postgis_raster;
30.1. Creando raster desde geometrías¶
Comenzaremos creando datos raster a partir de datos vectoriales, y luego pasaremos al enfoque más interesante de cargar datos desde una fuente raster. Verás que los datos raster están disponibles en abundancia y, a menudo, de forma gratuita en varios sitios gubernamentales.
Comenzaremos convirtiendo algunas geometrías en raster usando la función ST_AsRaster de la siguiente manera.
CREATE TABLE rasters (name varchar, rast raster);
INSERT INTO rasters(name, rast)
SELECT f.word, ST_AsRaster(geom, width=>150, height=>150)
FROM (VALUES ('Hello'), ('Raster') ) AS f(word)
, ST_Letters(word) AS geom;
CREATE INDEX ix_rasters_rast
ON rasters USING gist(ST_ConvexHull(rast));
El ejemplo anterior crea una tabla (rasters) a partir de geometrías formadas de letras usando la función ST_Letters de PostGIS 3.2+. Los raster, al igual que las geometrías, pueden aprovechar los índices espaciales. El índice espacial usado para raster es un índice funcional que indexa el convexhull geométrico del raster.
Puedes ver algunos metadatos útiles de tus raster con la siguiente consulta que utiliza la función raster de postgis ST_Count para contar el número de píxeles que tienen datos y la función ST_MetaData para proporcionar todo tipo de información de fondo útil para nuestros raster.
SELECT name, ST_Count(rast) As num_pixels, md.*
FROM rasters, ST_MetaData(rast) AS md;
name | num_pixels | upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
--------+------------+------------+-------------------+-------+--------+--------------------+---------------------+-------+-------+------+----------
Hello | 13926 | 0 | 77.10000000000001 | 150 | 150 | 1.226888888888889 | -0.5173333333333334 | 0 | 0 | 0 | 1
Raster | 11967 | 0 | 75.4 | 150 | 150 | 1.7226319023207244 | -0.5086666666666667 | 0 | 0 | 0 | 1
(2 rows)
Nota
Hay dos niveles de funciones raster. Hay funciones como ST_MetaData que trabajan a nivel de raster y hay funciones como la función ST_Count y la función ST_BandMetaData que trabajan a nivel de banda. La mayoría de las funciones en postgis raster que trabajan a nivel de banda, funcionan con solo una banda a la vez, y asumen que la banda que quieres es la 1.
Si tienes un raster multibanda, y necesitas contar los valores de píxeles no nulos en una banda distinta de la 1, deberías especificar explícitamente el número de banda como sigue ST_Count(rast,2).
Observa cómo todos los raster tienen una dimensión de 150x150. Esto no es ideal. Esto significa que para forzarlo, nuestros raster están deformados de diversas formas. Si tan solo pudiéramos ver lo feo de los raster frente a nosotros.
30.2. Cargando raster usando raster2pgsql¶
raster2pgsql es una herramienta de línea de comandos que a menudo viene incluida con PostGIS. Si estás en Windows y usaste el instalador PostGIS Bundle con stackbuilder, encontrarás raster2pgsql.exe en la carpeta C:\Program Files\PostgreSQL\15\bin donde el 15 debe ser reemplazado con la versión de PostgreSQL que estés usando.
Si estás usando Postgres.App, encontrarás raster2pgsql entre las otras Postgres.app CLI Tools.
En Ubuntu y Debian, necesitarás
apt install postgis
tener instaladas las herramientas de línea de comandos de PostGIS. Esto también puede instalar una versión adicional de PostgreSQL. Puedes ver una lista de clusters en Debian/Ubuntu usando el comando pg_lsclusters y eliminarlos usando el comando pg_dropcluster.
Para este y los siguientes ejercicios, usaremos nyc_dem.tif que se encuentra en el archivo PG Raster Workshop Dataset https://postgis.net/stuff/workshop-data/postgis_raster_workshop.zip. Para algunos ejemplos de geometría/raster, también usaremos datos de NYC cargados de capítulos anteriores. En lugar de cargar el tif, puedes restaurar el nyc_dem.backup incluido en el archivo zip en tu base de datos usando la herramienta de línea de comandos pg_restore o el menú Restore de pgAdmin.
Nota
Estos datos raster provienen de NYC DEM 1-foot Integer que es un tif DEM de 3GB que representa elevación relativa al nivel del mar con edificios y superficies de agua eliminados. Luego creamos una versión de menor resolución de este.
La herramienta raster2pgsql es similar a shp2gpsql excepto que en lugar de cargar archivos shapefiles de ESRI en tablas de geometría/geografía de PostGIS, carga cualquier formato raster soportado por GDAL en tablas raster. Al igual que shp2pgsql puedes pasarle un identificador de referencia espacial (SRID) de la fuente. A diferencia de shp2pgsql puede inferir el sistema de referencia espacial de los datos de origen si tus datos contienen metadatos adecuados.
Para ver todas las posibles opciones ofrecidas consulta raster2pgsql options.
Algunas otras opciones notables que ofrece raster2pgsql que no cubriremos son:
Capacidad para indicar el SRID de la fuente. En su lugar, confiaremos en la capacidad de raster2pgsql para adivinar.
Capacidad para establecer el valor nodata, cuando no se especifica, raster2pgsql intenta inferirlo del archivo.
Capacidad para cargar rasters fuera de la base de datos.
Para cargar todos los archivos tif en nuestra carpeta y también crear vistas generales (overviews), ejecutaríamos lo siguiente.
raster2pgsql -d -e -l 2,3 -I -C -M -F -Y -t 256x256 *.tif nyc_dem | psql -d nyc
-d para eliminar las tablas si ya existen
El comando anterior usa -e para cargar inmediatamente en lugar de confirmar en una transacción
-C establece restricciones raster, esto es útil para que raster_columns muestre información. Puedes querer combinarlo con -x para excluir la restricción de extensión, que es una restricción lenta de verificar y también dificulta futuras cargas en la tabla.
-M para hacer vacuum y analyze después de cargar, para mejorar las estadísticas del planificador de consultas
-Y para usar copy en lotes de 50. Si estás ejecutando PostGIS 3.3 o superior, puedes usar -Y 1000 para que copy se haga en lotes de 1000, o incluso un número mayor. Esto será más rápido, pero usará más memoria.
-l 2,3 para crear tablas de vista general: o_2_ncy_dem y o_3_nyc_dem. Esto es útil para visualizar datos.
-I para crear un índice espacial
-F para añadir el nombre de archivo, si solo tienes un archivo tif, esto no tiene mucho sentido. Si tuvieras varios, sería útil para saber de qué archivo proviene cada fila.
-t para establecer el tamaño de bloque. Nota: si no estás seguro del mejor tamaño a usar, usa -t auto y raster2pgql usará el mismo tamaño de teselado que el tif original. La salida te dirá cuál es el tamaño de bloque elegido. Cancela si parece enorme o extraño. El archivo original tenía un tamaño de 300x7, lo cual no es ideal.
Usa psql para ejecutar el sql generado contra la base de datos. Si quieres volcarlo a un archivo en su lugar, usa > nyc_dem.sql
Para este ejemplo, solo tenemos un archivo tif, así que podríamos especificar el nombre completo del archivo, en lugar de *.tif. Si los archivos no están en tu directorio actual, también puedes especificar una ruta de carpeta con *.tif.
Nota
Si estás en Windows y necesitas referenciar la carpeta, asegúrate de incluir la letra de la unidad, como C:/workshop/*.tif
A menudo escucharás en la jerga de PostGIS los términos raster tile y raster usados de manera algo intercambiable. Un raster tile realmente corresponde a un raster particular en una columna raster que es un subconjunto de un raster más grande, como estos datos DEM de NYC que acabamos de cargar. Esto se debe a que cuando los raster se cargan en PostGIS desde archivos raster grandes, se cortan en muchas filas para hacerlos manejables. Cada raster en cada fila es entonces parte de un raster más grande. Cada tile cubre un área del mismo tamaño indicada por el blocksize que especificaste. Los raster están tristemente limitados por el límite de 1GB de PostgreSQL TOAST y también por el proceso lento de «detoasting», por lo que necesitamos cortarlos para lograr un rendimiento decente o incluso para poder almacenarlos.
30.3. Visualizando raster en el navegador¶
Aunque pgAdmin y psql no tienen aún un mecanismo para ver raster de PostGIS, tenemos un par de opciones. Para raster pequeños la forma más fácil es exportarlos a un formato raster amigable para la web como PNG usando funciones incluidas de PostGIS raster como ST_AsPNG o ST_AsGDALRaster listadas en PostGIS Raster output functions. A medida que tus raster sean más grandes, querrás pasar a una herramienta como QGIS para verlos en todo su esplendor o la familia de herramientas de línea de comandos GDAL como gdal_translate para exportarlos a otros formatos raster. Recuerda, sin embargo, que los raster de PostGIS están construidos para análisis, no para generar imágenes bonitas para que las mires.
Un detalle: de manera predeterminada todos los tipos de salida raster están deshabilitados. Para utilizarlos, necesitarás habilitar drivers, todos o un subconjunto, como se detalla en Enable GDAL Raster drivers
SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
Si no quieres tener que hacer esto para cada conexión, puedes establecerlo a nivel de base de datos usando:
ALTER DATABASE nyc SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
Cada nueva conexión a la base de datos usará esa configuración.
Ejecuta la siguiente consulta y copia y pega la salida en la barra de direcciones de tu navegador web.
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Hello';
Para los raster creados hasta ahora, no especificamos el número de bandas ni siquiera definimos su relación con la Tierra. Como tal, nuestros raster tienen un sistema de referencia espacial desconocido (0).
Puedes pensar en el exoesqueleto de un raster como una geometría. Una matriz encerrada en un envolvente geométrico. Para hacer análisis útil, necesitamos georreferenciar nuestros raster, lo que significa que queremos que cada píxel (rectángulo) represente alguna parcela significativa de espacio.
El ST_AsRaster tiene muchas representaciones sobrecargadas. El ejemplo anterior usó la implementación más simple y aceptó los argumentos predeterminados que son 8BUI y 1 banda, con nodata siendo 0. Si necesitas usar las otras variantes, deberías usar la sintaxis de llamada con argumentos nombrados para no caer accidentalmente en la variante equivocada de la función o recibir errores de function is not unique.
Si comienzas con una geometría que tiene un sistema de referencia espacial, terminarás con un raster con el mismo sistema de referencia espacial. En este siguiente ejemplo, colocaremos nuestras palabras en Nueva York en colores brillantes y alegres. También usaremos escala de píxel en lugar de ancho y alto para que el tamaño de nuestros píxeles raster represente 1 metro x 1 metro de espacio.
INSERT INTO rasters(name, rast)
SELECT f.word || ' in New York' ,
ST_AsRaster(geom,
scalex => 1.0, scaley => -1.0,
pixeltype => ARRAY['8BUI', '8BUI', '8BUI'],
value => CASE WHEN word = 'Hello' THEN
ARRAY[10,10,100] ELSE ARRAY[10,100,10] END,
nodataval => ARRAY[0,0,0], gridx => NULL, gridy => NULL
) AS rast
FROM (
VALUES ('Hello'), ('Raster') ) AS f(word)
, ST_SetSRID(
ST_Translate(ST_Letters(word),586467,4504725), 26918
) AS geom;
Si luego miramos esto, veremos una geometría coloreada no deformada.
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Hello in New York';
Repetir para Raster:
SELECT 'data:image/png;base64,' ||
encode(ST_AsPNG(rast),'base64')
FROM rasters
WHERE name = 'Raster in New York';
Lo más revelador es que si volvemos a ejecutar el
SELECT name, ST_Count(rast) As num_pixels, md.*
FROM rasters, ST_MetaData(rast) AS md;
Observa los metadatos de las entradas de Nueva York. Tienen el sistema de referencia espacial state plane meter de Nueva York. También tienen la misma escala. Como cada unidad es de 1x1 metro, el ancho de la palabra Raster es ahora mayor que Hello.
name | num_pixels | upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
-------------------+------------+------------+-------------------+-------+--------+--------------------+---------------------+-------+-------+-------+----------
Hello | 13926 | 0 | 77.10000000000001 | 150 | 150 | 1.226888888888889 | -0.5173333333333334 | 0 | 0 | 0 | 1
Raster | 11967 | 0 | 75.4 | 150 | 150 | 1.7226319023207244 | -0.5086666666666667 | 0 | 0 | 0 | 1
Hello in New York | 8786 | 586467 | 4504802.1 | 184 | 78 | 1 | -1 | 0 | 0 | 26918 | 3
Raster in New York | 10544 | 586467 | 4504800.4 | 258 | 76 | 1 | -1 | 0 | 0 | 26918 | 3
(4 rows)
30.4. Tablas de catálogo espacial de raster¶
Al igual que con los tipos geometry y geography, raster tiene un conjunto de catálogos que muestran todas las columnas raster en tu base de datos. Estos son raster_columns and raster_overviews.
30.4.1. raster_columns¶
La vista raster_columns es la hermana de geometry_columns y geography_columns, proporcionando básicamente los mismos datos y más, pero para columnas raster.
SELECT *
FROM raster_columns;
Explora la tabla y encontrarás esto:
r_table_catalog | r_table_schema | r_table_name | r_raster_column | srid | scale_x | scale_y | blocksize_x | blocksize_y | same_alignment | regular_blocking | num_bands | pixel_types | nodata_values | out_db | extent | spatial_index
----------------+----------------+--------------+-----------------+------+---------+---------+-------------+-------------+----------------+------------------+-----------+-------------+---------------+--------+--------+---------------
nyc | public | rasters | rast | 0 | | | | | f | f | | | | | | t
nyc | public | nyc_dem | rast | 2263 | 10 | -10 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
nyc | public | o_2_nyc_dem | rast | 2263 | 20 | -20 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
nyc | public | o_3_nyc_dem | rast | 2263 | 30 | -30 | 256 | 256 | t | f | 1 | {16BUI} | {NULL} | {f} | | t
(4 rows)
una fila decepcionante de datos, en gran parte sin llenar, para la tabla rasters.
A diferencia de geometry y geography, raster no soporta modificadores de tipo, porque el espacio de modificador de tipo es demasiado limitado y hay propiedades más críticas que las que caben en un modificador de tipo.
Raster en su lugar se basa en restricciones, y lee estas restricciones como parte de la vista.
Mira las otras filas de las tablas que cargamos usando raster2pgsql. Debido a que usamos la opción -C, raster2pgsql añadió restricciones para el srid y otra información que pudo leer del tif o que pasamos nosotros. Las tablas de overview generadas con la opción -l o_2_nyc_dem y o_3_nyc_dem también aparecen.
Intentemos añadir algunas restricciones a nuestra tabla.
SELECT AddRasterConstraints('public'::name, 'rasters'::name, 'rast'::name);
Y serás bombardeado con un montón de avisos sobre cómo tus datos raster son un desastre y nada puede ser restringido. Si vuelves a mirar raster_columns, todavía la misma decepcionante historia de muchas filas en blanco para rasters.
Para que las restricciones puedan aplicarse, todos los raster en tu tabla deben ser restringibles por al menos una regla.
Podemos quizás hacer esto, digamos simplemente que todos nuestros datos están en el sistema de coordenadas state plane de Nueva York.
UPDATE rasters SET rast = ST_SetSRID(rast,26918)
WHERE ST_SRID(rast) <> 26918;
SELECT AddRasterConstraints('public'::name, 'rasters'::name, 'rast'::name);
SELECT r_table_name AS t, r_raster_column AS c, srid,
blocksize_x AS bx, blocksize_y AS by, scale_x AS sx, scale_y AS sy,
ST_AsText(extent) AS e
FROM raster_columns
WHERE r_table_name = 'rasters';
Ah, progreso:
t | c | srid | bx | by | sx | sy | e
----------+------+-------+-----+-----+----+----+------------------------------------------
rasters | rast | 26918 | 150 | 150 | | | POLYGON((0 -0.90000000000..
(1 row)
Cuanto más puedas restringir todos tus raster, más columnas verás llenas y también más operaciones podrás realizar en todos los tiles de tu raster. Ten en cuenta que en algunos casos puede que no quieras aplicar todas las restricciones.
Por ejemplo, si planeas cargar más datos en tu tabla raster, querrás omitir la restricción de extensión ya que eso requeriría que todos los raster estén dentro de la extensión de la restricción.
30.4.2. raster_overviews¶
Las columnas de raster overview aparecen tanto en el meta-catálogo raster_columns como en otro meta-catálogo llamado raster_overviews. Los overviews se usan principalmente para acelerar la visualización en niveles de zoom altos. También pueden usarse para análisis rápidos y aproximados, proporcionando estadísticas menos precisas, pero a una velocidad mucho mayor que al aplicar sobre la tabla raster original.
Para inspeccionar los overviews, ejecuta:
SELECT *
FROM raster_overviews;
y verás la salida:
o_table_catalog | o_table_schema | o_table_name | o_raster_column | r_table_catalog | r_table_schema | r_table_name | r_raster_column | overview_factor
----------------+----------------+--------------+-----------------+-----------------+----------------+--------------+-----------------+-----------------
nyc | public | o_2_nyc_dem | rast | nyc | public | nyc_dem | rast | 2
nyc | public | o_3_nyc_dem | rast | nyc | public | nyc_dem | rast | 3
(2 rows)
La tabla raster_overviews solo te proporciona el overview_factor y el nombre de la tabla padre. Toda esta información es algo que podrías haber deducido tú mismo por la convención de nombres de raster2pgsql para overviews.
El overview_factor te dice en qué resolución está la fila con respecto a su padre. Un overview_factor de 2 significa que 2x2 = 4 tiles pueden caber en un tile overview_2. De manera similar, un overview_factor de 3 significa que 2x2x2 = 8 tiles del original pueden agruparse en un tile overview_3.
30.5. Funciones raster comunes¶
La extensión postgis_raster tiene más de 100 funciones para elegir. La funcionalidad raster de PostGIS fue diseñada siguiendo el soporte de geometría de PostGIS. Encontrarás una superposición de funciones entre raster y geometría donde tiene sentido. Algunas comunes que usarás y que tienen equivalentes en el mundo de geometría son: ST_Intersects, ST_SetSRID, ST_SRID, ST_Union, ST_Intersection, y ST_Transform.
Además de esas funciones superpuestas, soporta el operador && de solapamiento entre rasters y entre un raster y una geometría. También ofrece muchas funciones que trabajan en conjunto con geometría o que son muy específicas de raster.
Necesitas una función como ST_Union para reconstituir una región. Debido a que el rendimiento se vuelve lento mientras más píxeles necesita analizar una función, necesitas una función de acción rápida ST_Clip para recortar los raster solo a las porciones de interés para tu análisis.
Finalmente necesitas ST_Intersects o && para acercarte a los tiles raster que contienen tus áreas de interés. El operador && es un proceso más rápido que ST_Intersects. Ambos pueden aprovechar los índices espaciales raster. Cubriremos estas funciones esenciales primero antes de pasar a otras secciones donde las usaremos en conjunto con otras funciones raster y de geometría.
30.5.1. Uniendo rasters con ST_Union¶
La función ST_Union para raster, al igual que su equivalente en geometría ST_Union, agrega un conjunto de raster en un solo raster. Sin embargo, al igual que con geometría, no todos los raster pueden combinarse, pero las reglas para unir raster son más complicadas que las reglas de geometría. En el caso de geometrías, todo lo que necesitas es tener el mismo sistema de referencia espacial, pero para raster eso no es suficiente.
Si intentaras lo siguiente:
SELECT ST_Union(rast)
FROM rasters;
Serías castigado inmediatamente con un error:
ERROR: rt_raster_from_two_rasters: The two rasters provided do not have the same alignment SQL state: XX000
¿Qué es esto del mismo alineamiento, que impide unir tus preciados raster?
Para que los raster se combinen, deben estar en la misma grilla por así decirlo. Esto significa que deben tener los mismos tamaños de píxel, la misma orientación (la inclinación), el mismo sistema de referencia espacial, y sus píxeles no deben cortarse entre sí, es decir, deben compartir la misma grilla de píxeles en el mundo.
Si intentas la misma consulta, pero solo con palabras que colocamos cuidadosamente en Nueva York.
De nuevo, el mismo error. Estos son el mismo sistema de referencia espacial, los mismos tamaños de píxel, y aun así no es suficiente. Porque sus grillas no coinciden.
Podemos arreglar esto ajustando ligeramente las coordenadas superiores izquierdas en y y luego probando de nuevo. Si nuestras grillas comienzan a nivel entero, ya que nuestros tamaños de píxel son enteros, entonces los píxeles no se cortarán entre sí.
UPDATE rasters SET rast = ST_SetUpperLeft(rast,
ST_UpperLeftX(rast)::integer,
ST_UpperLeftY(rast)::integer)
WHERE name LIKE '%New York';
SELECT ST_Union(rast ORDER BY name)
FROM rasters
WHERE name LIKE '%New York%';
Voilà, funcionó, y si lo viéramos, veríamos algo como esto:
Nota
Si alguna vez no tienes claro por qué tus raster no tienen el mismo alineamiento, puedes usar la función ST_SameAlignment, que comparará 2 raster o un conjunto de raster y te dirá si tienen el mismo alineamiento. Si tienes los avisos habilitados, el NOTICE te dirá qué está mal con los raster en cuestión. La función ST_NotSameAlignmentReason, en lugar de solo un aviso, devolverá la razón. Sin embargo, solo funciona con dos raster a la vez.
Una forma importante en que la función ST_Union(raster) difiere de la función ST_Union(geometry) es que permite un argumento llamado uniontype. Este argumento por defecto se establece en LAST si no lo especificas, lo que significa tomar los valores de píxeles del LAST raster en las ocasiones donde los valores de píxeles raster se superponen. Como regla general, los píxeles en una banda que están marcados como no-data se ignoran.
Al igual que con la mayoría de los agregados en PostgreSQL, puedes poner una cláusula ORDER BY como parte de la llamada a la función, como se hace en el ejemplo anterior. Especificar el orden te permite controlar qué raster tiene prioridad. Así que en nuestro ejemplo anterior, Raster venció a Hello porque Raster es alfabéticamente el último.
Observa, si cambias el orden:
SELECT ST_Union(rast ORDER BY name DESC)
FROM rasters
WHERE name LIKE '%New York%';
Entonces Hello vence a Raster porque Hello ahora es el último superpuesto.
El tipo de unión FIRST es lo contrario de LAST.
Pero en ocasiones, LAST puede no ser la operación correcta. Supongamos que nuestros raster representan dos conjuntos diferentes de observaciones de dos dispositivos distintos. Estos dispositivos miden lo mismo, y no estamos seguros de cuál es correcto cuando se cruzan, por lo que en su lugar nos gustaría tomar el MEAN de los resultados. Haríamos esto:
SELECT ST_Union(rast, 'MEAN')
FROM rasters
WHERE name LIKE '%New York%';
Voilà, funcionó, y si lo viéramos, veríamos algo como esto:
Así que en lugar de imponerse, tenemos una mezcla de las dos fuerzas. En el caso del tipo de unión MEAN, no tiene sentido especificar el orden, porque el resultado sería el promedio de los valores de píxeles superpuestos.
Ten en cuenta que para geometrías, dado que las geometrías son vectores y, por lo tanto, no tienen valores aparte de estar o no estar, realmente no hay ambigüedad sobre cómo combinar dos vectores cuando se intersectan.
Otra característica del raster ST_Union que pasamos por alto, es esta idea de si deberías devolver todas las bandas o solo algunas bandas. Cuando no especificas qué bandas unir, ST_Union combinará los números de las mismas bandas y usará la estrategia de unión LAST. Si tienes múltiples bandas, esto puede no ser lo que quieras hacer. Quizás solo quieras unir la segunda banda. En este caso, la banda verde, y quieres el conteo de valores de píxel.
SELECT ST_BandPixelType(ST_Union(rast, 2, 'COUNT'))
FROM rasters
WHERE name LIKE '%New York%';
st_bandpixeltype
------------------
32BUI
(1 row)
Nota que en el caso del tipo de unión COUNT, que cuenta el número de píxeles rellenados y devuelve ese valor, el resultado es siempre un 32BUI similar a cómo cuando haces un COUNT en SQL, el resultado siempre es un bigint, para acomodar conteos grandes.
En otros casos, el tipo de píxel de la banda no cambia y se establece en el valor máximo o se redondea si las cantidades exceden los límites del tipo. ¿Por qué alguien querría contar píxeles que se intersectan en una ubicación? Bueno, supongamos que cada uno de tus raster representa escuadrones de policía e incidentes de arrestos en las áreas. Cada valor podría representar una razón diferente de arresto. Estás haciendo estadísticas sobre cuántos arrestos en cada región, por lo tanto solo te importa el conteo de arrestos.
O quizás quieras usar todas las bandas, pero quieres diferentes estrategias.
SELECT ST_Union(rast, ARRAY[(1, 'MAX'),
(2, 'MEAN'),
(3, 'RANGE')]::unionarg[])
FROM rasters
WHERE name LIKE '%New York%';
Usar la variante unionarg[] de la función ST_Union, también te permite reorganizar el orden de las bandas.
30.5.2. Recortando raster con ayuda de ST_Intersects¶
La función ST_Clip es una de las funciones más utilizadas para raster en PostGIS. La razón principal es que mientras más píxeles necesites inspeccionar o realizar operaciones, más lento será tu procesamiento. ST_Clip recorta tu raster solo al área de interés, para que puedas aislar tus operaciones únicamente a esa área.
Esta función es también especial en que utiliza el poder de la geometría para ayudar al análisis raster. Para reducir el número de píxeles que ST_Union tiene que manejar, cada raster se recorta primero al área de interés.
SELECT ST_Union( ST_Clip(r.rast, g.geom) )
FROM rasters AS r
INNER JOIN
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ) AS g(geom)
ON ST_Intersects(r.rast, g.geom)
WHERE r.name LIKE '%New York%';
Este ejemplo muestra varias funciones trabajando en conjunto. La función ST_Intersects empleada es la que viene con postgis_raster y puede intersectar 2 raster o un raster y una geometría. Similar a la función de geometría ST_Intersects, el raster ST_Intersects puede aprovechar los índices espaciales en las tablas raster o de geometría.
Nota
De manera predeterminada ST_Clip dejará fuera píxeles donde el centroide del píxel no intersecta la geometría. Esto puede ser molesto para píxeles grandes y puede que prefieras en su lugar incluir un píxel si cualquier parte del píxel toca la geometría. Introducido en PostGIS 3.5, está el argumento touched. Reemplaza tu ST_Clip(r.rast, g.geom) con ST_Clip(r.rast, g.geom, touched => true) y voilà, cualquier píxel que intersecte tu geometría de cualquier forma será incluido.
30.6. Convirtiendo raster en geometrías¶
Los raster pueden convertirse fácilmente en geometrías.
30.6.1. El polígono de un raster con ST_Polygon¶
Comencemos con nuestro ejemplo anterior, pero convirtámoslo en un polígono usando la función ST_Polygon.
SELECT ST_Polygon(ST_Union( ST_Clip(r.rast, g.geom) ))
FROM rasters AS r
INNER JOIN
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ) AS g(geom)
ON ST_Intersects(r.rast, g.geom)
WHERE r.name LIKE '%New York%';
Si haces clic en el visor de geometría en pgAdmin, puedes ver esto en todo su esplendor sin ningún truco.
ST_Polygon considera todos los píxeles que tienen valores (no no-data) en una banda particular, y los convierte en geometría. Al igual que muchas otras funciones en raster, ST_Polygon solo considera una banda. Si no se especifica ninguna banda, considerará solo la primera banda.
30.6.2. Los rectángulos de píxel de un raster con ST_PixelAsPolygons¶
Otra función popularmente utilizada es la función ST_PixelAsPolygons. Rara vez deberías usar ST_PixelAsPolygons en un raster grande sin antes recortarlo porque terminarás con millones de filas, una por cada píxel.
ST_PixelAsPolygons devuelve una tabla que consiste en geom, val, x y y. Donde x es el número de columna, e y es el número de fila en el raster.
ST_PixelAsPolygons, similar a otras funciones raster, trabaja en una banda a la vez y trabaja sobre la banda 1 si no se especifica ninguna banda. También, por defecto, devuelve solo los píxeles que tienen valores.
SELECT gv.*
FROM rasters AS r
CROSS JOIN LATERAL ST_PixelAsPolygons(rast) AS gv
WHERE r.name LIKE '%New York%'
LIMIT 10;
La salida es:
y si inspeccionamos usando el visor de geometría, veríamos:
Si queremos todos los píxeles de todas nuestras bandas, necesitaríamos hacer algo como lo siguiente. Nota las diferencias en este ejemplo con respecto al anterior.
1. Setting exclude_nodata_value to make sure all pixels are returned so that our sets of calls return the same number of rows. The rows out of the function will be naturally in the same order.
2. Using the PostgreSQL ROWS FROM constructor , and aliasing each set of columns from our function output with names. So for example the band 1 columns (geom, val, x, y) are renamed to g1, v1, x1, x2
SELECT pp.g1, pp.v1, pp.v2, pp.v3
FROM rasters AS r
CROSS JOIN LATERAL
ROWS FROM (
ST_PixelAsPolygons(rast, 1, exclude_nodata_value => false ),
ST_PixelAsPolygons(rast, 2, exclude_nodata_value => false),
ST_PixelAsPolygons(rast, 3, exclude_nodata_value => false )
) AS pp(g1, v1, x1, y1,
g2, v2, x2, y2,
g3, v3, x3, y3 )
WHERE r.name LIKE '%New York%'
AND ( pp.v1 = 0 OR pp.v2 > 0 OR pp.v3 > 0) ;
Nota
Usamos CROSS JOIN LATERAL en estos ejemplos porque queríamos ser explícitos en lo que estamos haciendo. Como estas son todas funciones que devuelven conjuntos, puedes reemplazar CROSS JOIN LATERAL con , como abreviatura. Usaremos , en el siguiente conjunto de ejemplos
30.6.3. Guardando polígonos con ST_DumpAsPolygons¶
Raster también introduce un tipo compuesto adicional llamado geomval. Considera un geomval como la descendencia de una geometría y un raster. Contiene una geometría y contiene un valor de píxel.
Encontrarás varias funciones raster que devuelven geomvals.
Una función comúnmente usada que produce geomvals es ST_DumpAsPolygons, que devuelve un conjunto de píxeles contiguos con el mismo valor como un polígono. Nuevamente, por defecto, solo revisará la banda 1 y excluirá los valores no-data a menos que los sobrescribas. Este ejemplo selecciona solo polígonos de la banda 2. También puedes aplicar filtros a los valores. Para la mayoría de los casos de uso, ST_DumpAsPolygons es una mejor opción que ST_PixelAsPolygons, ya que devolverá muchas menos filas.
Esto producirá 6 filas, y devolverá polígonos correspondientes a las letras en «Raster».
SELECT gv.geom , gv.val
FROM rasters AS r,
ST_DumpAsPolygons(rast, 2) AS gv
WHERE r.name LIKE '%New York%'
AND gv.val = 100;
Observa que no devuelve una geometría única, porque encuentra conjuntos contiguos de píxeles con el mismo valor que forman un polígono. Aunque todos estos valores son iguales, no son contiguos.
Un enfoque común para producir geometrías más complejas es agrupar por los valores y unir.
SELECT ST_Union(gv.geom) AS geom , gv.val
FROM rasters AS r,
ST_DumpAsPolygons(rast, 2) AS gv
WHERE r.name LIKE '%New York%'
GROUP BY gv.val;
Esto te devolverá 2 filas correspondientes a las palabras «Raster» y «Hello».
30.7. Estadísticas¶
Lo más importante que debes entender sobre los raster es que son herramientas estadísticas para almacenar datos en matrices, que puede que logres hacer ver bonitos en una pantalla.
Puedes encontrar un menú de estas funciones estadísticas en Raster Band Statistics.
30.7.1. ST_SummaryStatsAgg y ST_SummaryStats¶
¿Quieres todas las estadísticas para un conjunto de raster? Usa la función ST_SummaryStatsAgg.
Esta consulta toma unos 10 segundos y te da un resumen de toda la tabla:
SELECT (ST_SummaryStatsAgg(rast, 1, true, 1)).* AS sa
FROM o_3_nyc_dem;
Salidas:
count | sum | mean | stddev | min | max
-----------+------------+------------------+------------------+-----+-----
246794100 | 4555256024 | 18.4577184948911 | 39.4416860598687 | 0 | 411
(1 row)
Lo cual nos dice que tenemos muchos píxeles y nuestra elevación máxima es de 411 pies.
Si has construido overviews y solo necesitas una estimación aproximada de tus mínimos, máximos y medias usa uno de tus overviews. Esta siguiente consulta devuelve aproximadamente los mismos valores de mínimos, máximos y medias que la anterior, pero en alrededor de 1 segundo en lugar de 10.
SELECT (ST_SummaryStatsAgg(rast, 1, true, 1)).* AS sa
FROM o_3_nyc_dem ;
Ahora, con esta información, podemos hacer más preguntas.
30.7.2. ST_Histogram¶
Generalmente no querrás estadísticas para toda tu tabla, sino solo estadísticas para un área en particular, en ese caso, también querrás emplear a nuestros viejos amigos ST_Intersects y ST_Clip. Si también necesitas una función estadística raster que no tenga una versión agregada, deberás llevar ST_Union contigo.
Para el siguiente ejemplo usaremos una función de estadísticas diferente ST_Histogram que no tiene equivalente agregado, y para esta variante en particular, es una función que devuelve un conjunto. Estamos usando la misma área de interés que en algunos ejemplos anteriores, pero también necesitamos emplear la geometría ST_Transform para transformar nuestra geometría en metros state plane de NY a nuestros raster en pies state plane de NYC. Casi siempre es más eficiente transformar la geometría en lugar del raster, y definitivamente si tu geometría es solo una.
SELECT (ST_Quantile( ST_Union( ST_Clip(r.rast, g.geom) ), ARRAY[0.25,0.50,0.75, 1.0] )).*
FROM nyc_dem AS r
INNER JOIN
ST_Transform(
ST_Buffer(ST_Point(586598, 4504816, 26918), 100 ),
2263) AS g(geom)
ON ST_Intersects(r.rast, g.geom);
la consulta anterior se completa en menos de 60ms y produce:
quantile | value
----------+-------
0.25 | 52
0.5 | 57
0.75 | 68
1 | 78
(4 rows)
30.8. Creando raster derivados¶
PostGIS raster viene con varias funciones para editar raster. Estas funciones se usan tanto para editar como para crear conjuntos de datos raster derivados. Encontrarás estas listadas en Raster Editors y Raster Management.
30.8.1. Transformando raster con ST_Transform¶
La mayoría de nuestros datos están en NY State Plane metros (SRID: 26918), sin embargo nuestro conjunto de datos DEM raster está en NY State Plane pies (SRID: 2263). Para un flujo de trabajo menos engorroso, necesitamos que nuestros conjuntos de datos principales estén en el mismo sistema de referencia espacial.
La función raster ST_Transform es la más adecuada para este trabajo.
Para crear un nuevo conjunto de datos nyc dem en NY State Plane metros, haremos lo siguiente:
CREATE TABLE nyc_dem_26918 AS
WITH ref AS (SELECT ST_Transform(rast,26918) AS rast
FROM nyc_dem LIMIT 1)
SELECT r.rid, ST_Transform(r.rast, ref.rast) AS rast, r.filename
FROM nyc_dem AS r, ref;
Lo anterior en mi sistema tardó alrededor de 1.5 minutos. Para un conjunto de datos más grande tardaría mucho más.
Los ejemplos mencionados usaron dos variantes de la función raster ST_Transform. La primera fue para obtener un raster de referencia que se usará para transformar los otros tiles raster para garantizar que todos los tiles tengan el mismo alineamiento. Observa que la segunda variante de ST_Transform utilizada ni siquiera toma un SRID de entrada. Esto se debe a que el SRID y toda la escala de píxeles y tamaños de bloque se leen del raster de referencia. Si usaras la forma ST_Transform(rast, srid), entonces todos tus raster podrían salir con diferentes alineamientos, haciendo imposible aplicar una operación como ST_Union sobre ellos.
El único problema con el enfoque ST_Transform mencionado es que cuando transformas, lo transformado a menudo existe en otros tiles. Si miraras la salida anterior con suficiente detalle sacando el convex hull de los raster, en el siguiente ejemplo verías solapamientos molestos alrededor de los bordes.
SELECT rast::geometry
FROM nyc_dem_26918
ORDER BY rid
LIMIT 100;
visto en pgAdmin se vería algo como:
30.8.2. Usando ST_MakeEmptyCoverage para crear raster con tiles uniformes¶
Un mejor enfoque, aunque un poco más lento, es definir tu propia estructura de tiles de cobertura desde cero usando ST_MakeEmptyCoverage y luego encontrar los tiles que se intersectan para cada nuevo tile, y hacer ST_Union de estos y luego usar ST_Transform(ref, ST_Union…) para crear cada tile.
Para esto estaremos usando varias funciones que aprendimos anteriormente.
DROP TABLE IF EXISTS nyc_dem_26918;
CREATE TABLE nyc_dem_26918 AS
SELECT ROW_NUMBER() OVER(ORDER BY t.rast::geometry) AS rid,
ST_Union(ST_Clip( ST_Transform( r.rast, t.rast, 'Bilinear' ), t.rast::geometry ), 'MAX') AS rast
FROM (SELECT ST_Transform(
ST_SetSRID(ST_Extent(rast::geometry),2263)
, 26918) AS geom
FROM nyc_dem
) AS g, ST_MakeEmptyCoverage(tilewidth => 256, tileheight => 256,
width => (ST_XMax(g.geom) - ST_XMin(g.geom))::integer,
height => (ST_YMax(g.geom) - ST_YMin(g.geom))::integer,
upperleftx => ST_XMin(g.geom),
upperlefty => ST_YMax(g.geom),
scalex => 3.048,
scaley => -3.048,
skewx => 0., skewy => 0., srid => 26918) AS t(rast)
INNER JOIN nyc_dem AS r
ON ST_Transform(t.rast::geometry, 2263) && r.rast
GROUP BY t.rast;
Repitiendo el ejercicio:
SELECT rast::geometry
FROM nyc_dem_26918
ORDER BY rid
LIMIT 100;
visto en pgAdmin ya no tenemos solapamientos:
En mi sistema esto tomó ~10 minutos y devolvió 3879 filas. Después de la creación de la tabla, querrás hacer lo habitual de añadir un índice espacial, clave primaria y restricciones como sigue:
ALTER TABLE nyc_dem_26918
ADD CONSTRAINT pk_nyc_dem_26918 PRIMARY KEY(rid);
CREATE INDEX ix_nyc_dem_26918_st_convexhull_gist
ON nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
SELECT AddRasterConstraints('nyc_dem_26918'::name, 'rast'::name);
ANALYZE nyc_dem_26918;
Lo cual debería tomar menos de 2 minutos para este conjunto de datos.
30.8.3. Creando tablas de overview con ST_CreateOverview¶
Al igual que con nuestro conjunto de datos original, sería útil tener tablas de overview para acelerar el rendimiento de algunas operaciones. ST_CreateOverview es una función adecuada para ese propósito. Puedes usar ST_CreateOverview para crear overviews también si olvidaste crear una durante la carga con raster2pgsql o decidiste que necesitas más overviews.
Crearemos overviews de nivel 2 y 3 como habíamos hecho con nuestro original usando este código.
SELECT ST_CreateOverview('nyc_dem_26918'::regclass, 'rast', 2);
SELECT ST_CreateOverview('nyc_dem_26918'::regclass, 'rast', 3);
Este proceso desafortunadamente toma tiempo, y más tiempo cuanto más filas tengas, así que sé paciente. Para este conjunto de datos tomó alrededor de 3-5 minutos para el overview factor 2 y 1 minuto para el overview factor 3.
La función ST_CreateOverview también añade las restricciones necesarias para que las columnas aparezcan con todo detalle en los catálogos raster_columns y raster_overviews. Sin embargo, no añade índices ni tampoco una columna rid. La columna rid probablemente no es necesaria a menos que necesites una clave primaria para editar. Probablemente quieras un índice que puedes crear con lo siguiente:
CREATE INDEX ix_o_2_nyc_dem_26918_st_convexhull_gist
ON o_2_nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
CREATE INDEX ix_o_3_nyc_dem_26918_st_convexhull_gist
ON o_3_nyc_dem_26918 USING gist( ST_ConvexHull(rast) );
Nota
ST_CreateOverview tiene un argumento opcional para indicar el método de muestreo. Si no se especifica, usa el valor por defecto NearestNeighbor que generalmente es el más rápido de calcular pero puede no ser el ideal. Los métodos de remuestreo están fuera del alcance de este taller.
30.9. La intersección de rasters y geometrías¶
Hay un par de funciones comúnmente usadas para calcular intersecciones de rasters y geometrías. Ya hemos visto ST_Clip en acción que devuelve la intersección de un raster y una geometría como un raster, pero hay otras. Para datos puntuales, la más comúnmente usada es ST_Value y luego está ST_Intersection que tiene varias sobrecargas, algunas devolviendo rasters y otras devolviendo un conjunto de geomval.
30.9.1. Valores de píxel en un Punto geométrico¶
Si necesitas devolver valores de raster basados en la intersección de varios puntos geométricos ad-hoc, usarás ST_Value o su pariente cercano ST_NearestValue.
SELECT g.geom, ST_Value(r.rast, g.geom) AS elev
FROM nyc_dem_26918 AS r
INNER JOIN
(SELECT id, geom
FROM nyc_homicides
WHERE weapon = 'gun') AS g
ON r.rast && g.geom;
Este ejemplo tarda alrededor de 1 segundo en devolver 2444 filas. Si usaras ST_Intersects en lugar de &&, el proceso tomaría alrededor de 3 segundos. La razón por la que ST_Intersects es más lento es porque realiza una verificación adicional en algunos casos, una verificación píxel por píxel. Si esperas que todos tus puntos estén representados con datos en tu conjunto raster y tus raster representan una cobertura (un conjunto contiguo de tiles raster no superpuestos), entonces && es generalmente una opción más rápida.
Si tus datos raster no están densamente poblados o tienes rasters superpuestos (por ejemplo, representan diferentes observaciones en el tiempo), o están inclinados (no alineados a los ejes), entonces hay una ventaja en que ST_Intersects elimine los falsos positivos.
30.9.2. Estilo raster de ST_Intersection¶
Así como puedes calcular la intersección de dos geometrías usando ST_Intersection, puedes calcular la intersección de dos raster o un raster y una geometría usando raster ST_Intersection.
Lo que obtienes de esta función son dos tipos diferentes de cosas:
Intersecta una geometría con un raster, y obtienes un conjunto de descendientes geomval. Quizás uno, pero la mayoría de las veces muchos.
Intersecta 2 rasters y obtienes un único raster de vuelta.
La regla de oro tanto para la intersección raster como para la intersección de geometrías es que ambas partes involucradas deben tener el mismo sistema de referencia espacial. Para raster/raster, también deben tener el mismo alineamiento.
Aquí hay un ejemplo que responde a una pregunta que podrías haber tenido curiosidad. Si agrupamos nuestras elevaciones en 5 intervalos de valores de elevación, ¿qué rango de elevación resulta en la mayor cantidad de muertes por arma de fuego? Sabemos, basándonos en nuestras estadísticas resumidas anteriores, que 0 es el valor más bajo y 411 es el valor más alto de elevación en nuestro conjunto de datos nyc dem, así que usamos eso como valor mínimo y máximo para nuestra llamada a width_bucket.
SELECT ST_Transform(ST_Union(gv.geom),4326) AS geom ,
MIN(gv.val) AS min_elev, MAX(gv.val) AS max_elev,
count(g.id) AS count_guns
FROM nyc_dem_26918 AS r
INNER JOIN nyc_homicides AS g
ON ST_Intersects(r.rast, g.geom)
CROSS JOIN
ST_Intersection( g.geom,
ST_Clip(r.rast,ST_Expand(g.geom, 4) )
) AS gv
WHERE g.weapon = 'gun'
GROUP BY width_bucket(gv.val, 0, 411, 5)
ORDER BY width_bucket(gv.val, 0, 411, 5);
¿Existe una correlación importante entre homicidios con armas de fuego y elevación? Probablemente no.
Veamos la intersección raster/raster:
SELECT ST_Intersection(r1.rast, 1, r2.rast, 1, 'BAND1')
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
Lo que obtenemos son dos filas con NULLs, y si tienes PostgreSQL configurado para mostrar avisos, verás:
NOTICE: The two rasters provided do not have the same alignment. Returning NULL
Para solucionar esto, podemos alinear uno con el otro en el momento de la salida usando ST_Resample.
SELECT ST_Intersection(r1.rast, 1,
ST_Resample( r2.rast, r1.rast ), 1,
'BAND1')
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
Vamos también a resumirlo en un único registro estadístico
SELECT (
ST_SummaryStatsAgg(
ST_Intersection(r1.rast, 1,
ST_Resample( r2.rast, r1.rast ), 1, 'BAND1'),
1, true)
).*
FROM nyc_dem_26918 AS r1
INNER JOIN
rasters AS r2 ON ST_Intersects(r1.rast,1, r2.rast, 1);
lo cual produce:
count | sum | mean | stddev | min | max
-------+-------+-----------------+------------------+-----+-----
2075 | 99536 | 47.969156626506 | 9.57974836865737 | 33 | 62
(1 row)
30.10. Funciones de álgebra de mapas¶
El álgebra de mapas es la idea de que puedes hacer matemáticas sobre tus valores de píxeles. Las funciones ST_Union y ST_Intersection cubiertas anteriormente son un caso especial rápido de álgebra de mapas. Luego está la familia de funciones ST_MapAlgebra que te permiten definir tus propias matemáticas locas, pero a costa del rendimiento.
La gente tiene la costumbre de saltar a ST_MapAlgebra, probablemente porque el nombre suena tan genial y sofisticado. ¿Quién no querría decirle a sus amigos? Estoy usando una función llamada ST_MapAlgebra. Mi consejo: explora otras funciones antes de sacar esa escopeta. Tu vida será más simple, tu rendimiento será 100 veces mejor y tu código será más corto.
Antes de mostrar ST_MapAlgebra, exploraremos otras funciones que encajan en la familia de funciones de Map Algebra y que generalmente tienen mejor rendimiento que ST_MapAlgebra.
30.10.1. Reclasificar tu raster usando ST_Reclass¶
Una función de álgebra de mapas a menudo pasada por alto es la función ST_Reclass, que espera en segundo plano a que alguien descubra el poder y la velocidad que puede ofrecer.
¿Qué hace ST_Reclass? Como su nombre lo indica, reclasifica tus valores de píxeles basándose en un álgebra de rangos minimalista.
Revisitemos nuestros NYC Dems. Tal vez solo nos importe clasificar nuestras elevaciones como 1) baja, 2) media, 3) alta, y 4) muy alta. No necesitamos 411 valores, solo necesitamos 4. Dicho esto, hagamos algo de reclasificación.
El esquema de clasificación está gobernado por la reclass expression.
WITH r AS ( SELECT ST_Union(newrast) As rast
FROM nyc_dem_26918 AS r
INNER JOIN ST_Buffer(ST_Point(586598, 4504816, 26918), 1000 ) AS g(geom)
ON ST_Intersects( r.rast, g.geom )
CROSS JOIN ST_Reclass( ST_Clip(r.rast,g.geom), 1,
'[0-10):1, [10-50):2, [50-100):3,[100-:4','4BUI',0) AS newrast
)
SELECT SUM(ST_Area(gv.geom)::numeric(10,2)) AS area, gv.val
FROM r, ST_DumpAsPolygons(rast) AS gv
GROUP BY gv.val
ORDER BY gv.val;
Lo cual produciría:
area | val
------------+-----
6754.04 | 1
1753117.51 | 2
1355232.37 | 3
1848.75 | 4
(4 rows)
Si este fuera un esquema de clasificación que prefiriéramos, podríamos crear una nueva tabla usando ST_Reclass para recalcular cada tile.
30.10.2. Coloreando tus rasters con ST_ColorMap¶
La función ST_ColorMap es otra función de tipo mapalgebra que reclasifica tus valores de píxeles. Excepto que es creadora de bandas. Convierte un raster de una sola banda, como nuestros Dems, en un raster de 3 o 4 bandas visualmente presentable.
Podrías usar uno de los colormaps incorporados como se muestra a continuación si no quieres complicarte creando uno.
SELECT ST_ColorMap( ST_Union(newrast), 'bluered') As rast
FROM nyc_dem_26918 AS r
INNER JOIN
ST_Buffer(
ST_Point(586598, 4504816, 26918), 1000
) AS g(geom)
ON ST_Intersects( r.rast, g.geom)
CROSS JOIN ST_Clip(rast, g.geom) AS newrast;
Lo cual se vería así:
Cuanto más azul es el color, menor es la elevación, y cuanto más rojo es el color, mayor es la elevación.
