A funcionalidade Raster no PostGIS faz parte da extensão principal desde que foi introduzida. Com o lançamento do PostGIS 3, se você quiser a funcionalidade Raster, você precisará instalar a extensão principal (postgis) e também a extensão postgis_raster. Veja:

CREATE EXTENSION postgis;
CREATE EXTENSION postgis_raster;

A separação da funcionalidade Raster permite que se crie os pacotes do PostGIS de uma forma simplificada, sem a necessidade das dependências da extensão Raster, que incluem a biblioteca GDAL, um tanto quanto pesada. A funcionalidade Raster permanece intacta, no entanto, e você ainda pode fazer coisas bacanas com ela.

Por exemplo, faça o download e importe um arquivo de “modelo de elevação digital”. No meu caso, vou utilizar um arquivo da cidade de Vancouver.

Para baixar essa imagem eu fiz o seguinte:

createdb yvr_raster
psql -c 'CREATE EXTENSION postgis' yvr_raster
psql -c 'CREATE EXTENSION postgis_raster' yvr_raster

1. Criei uma base de dados e instalei nela as extensões postgis e postgis_raster

wget https://pub.data.gov.bc.ca/datasets/175624/92g/092g06_e.dem.zip
unzip 092g06_e.dem.zip

2. Baixei os dados de DEM no site https://pub.data.gov.bc.ca/datasets/175624/

raster2pgsql -I -F -s 4269 -t 56x56 092g06_e.dem dem092g06e | psql yvr_raster

3. Carreguei os dados no banco.

Após o carregamento dos dados, temos uma tabela com o nome dem092g06e de elevação de 56 por 56 pixels. Se você mapear a extensão, eles se parecerão com isso:



Imagine uma elevação de 30 metros no nível do solo (em um caso extremo, a Groenlândia e a Antártica seriam 68 metros). Quanto de Vancouver ficaria debaixo d’água? É principalmente um lugar montanhoso. Vamos descobrir:

CREATE TABLE poly_30 AS 
  SELECT (
   ST_DumpAsPolygons(
    ST_SetBandNoDataValue(
     ST_Reclass(
      ST_Union(rast), 
      '0-30:1-1, 31-5000:0-0', '2BUI'),
     0))).*
FROM dem092g06e d

Existem muitas funções aninhadas, portanto, lendo de dentro para fora, nós temos:

1. União de todas os chips em uma grande imagem Raster
2. Reclassificação de todos os valores de 0 a 30 para 1 e todos os valores mais altos para 0
3. Definição do valor “nodata” como 0, não nos importamos com coisas que estão acima do nosso limite
4. Criação de um polígono vetorial para cada valor no raster (existe apenas um valor: “1”)

O resultado fica assim:



Podemos pegar os dados de construção em Vancouver e verificar quantos edifícios ficarão embaixo da água:

wget ftp://webftp.vancouver.ca/OpenData/shape/building_footprints_2009_shp.zip
unzip building_footprints_2009_shp.zip
shp2pgsql -I -s 26910 -D building_footprints_2009 buildings | psql yvr_raster

Antes comparar os edifícios com a zona de inundação, precisamos colocá-los na mesma projeção que a zona de inundação (SRID 4269).

ALTER TABLE buildings
ALTER COLUMN geom 
TYPE geometry(MultiPolygonZ, 4269)
USING ST_Transform(geom, 4269)

Agora podemos encontrar edifícios inundados:

CREATE TABLE buildings_30_poly AS
  SELECT b.* 
    FROM buildings b
    JOIN poly_30 p
      ON ST_Intersects(p.geom, b.geom)
    WHERE ST_Intersects(p.geom, ST_Centroid(b.geom))

Existe outra maneira de encontrar edifícios abaixo de 30 metros, sem ter que construir um polígono, assim:

CREATE TABLE buildings_30_rast AS
  SELECT b.*
    FROM buildings b
    JOIN dem092g06e d
      ON ST_Intersects(b.geom, d.rast)
    WHERE ST_Value(d.rast, ST_Centroid(b.geom)) < 30;

Como a construção de polígonos a partir da Raster pode ser cara, unir as camadas raster e a vetorial é geralmente a maneira que queremos realizar essa análise.



Para trazer dados contínuos (elevações, temperaturas, resultados do modelo) para a análise GIS, a extensão postgis_raster pode ser bastante útil.

Fonte: Clever Elephant Blog