Download de dados rasterizados em python a partir do postgis usando o psycopg2

13

Eu tenho dados rasterizados em uma tabela do postgres que eu quero entrar no python como uma matriz numpy. Estou usando o psycopg2 para conectar ao banco de dados. Consigo baixar os dados, mas eles retornam como uma string (provavelmente um binário serializado).

Alguém sabe como pegar essa string e converter em uma matriz numpy?

Eu explorei outras opções para baixar o raster, como usar st_astiff e codificar para baixar o arquivo hexadecimal e usar o xxd, mas isso não funcionou. Continuo recebendo o erro 'rt_raster_to_gdal: não foi possível carregar o driver GDAL de saída' e não tenho permissões para definir as variáveis ​​de ambiente para poder ativar os drivers.

TL, DR: deseja importar dados de varredura para uma matriz numpy (usando python).

Mayank Agarwal
fonte

Respostas:

14

rt_raster_to_gdal: Não foi possível carregar o driver GDAL de saída

Quanto ao primeiro erro com ST_AsTIFF , você precisa habilitar os drivers GDAL, que por padrão não estão habilitados para o PostGIS 2.1. Veja o manual sobre maneiras de fazer isso. Por exemplo, eu tenho uma variável de ambiente configurada em um computador Windows com:

POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED USGSDEM AAIGrid

o que pode ser confirmado com o PostGIS com:

SELECT short_name, long_name
FROM ST_GDALDrivers();

PostGIS para Numpy

Você pode exportar a saída para um arquivo GeoTIFF de memória virtual para o GDAL ler em uma matriz Numpy. Para obter dicas sobre arquivos virtuais usados ​​no GDAL, consulte esta postagem do blog .

import os
import psycopg2
from osgeo import gdal

# Adjust this to connect to a PostGIS database
conn = psycopg2.connect(...)
curs = conn.cursor()

# Make a dummy table with raster data
curs.execute("""\
    SELECT ST_AsRaster(ST_Buffer(ST_Point(1, 5), 10), 10, 10, '8BUI', 1) AS rast
    INTO TEMP mytable;
""")

# Use a virtual memory file, which is named like this
vsipath = '/vsimem/from_postgis'

# Download raster data into Python as GeoTIFF, and make a virtual file for GDAL
curs.execute("SELECT ST_AsGDALRaster(rast, 'GTiff') FROM mytable;")
gdal.FileFromMemBuffer(vsipath, bytes(curs.fetchone()[0]))

# Read first band of raster with GDAL
ds = gdal.Open(vsipath)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()

# Close and clean up virtual memory file
ds = band = None
gdal.Unlink(vsipath)

print(arr)  # this is a 2D numpy array

Mostra um ponto em buffer rasterizado.

[[0 0 0 1 1 1 1 0 0 0]
 [0 1 1 1 1 1 1 1 1 0]
 [0 1 1 1 1 1 1 1 1 0]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [0 1 1 1 1 1 1 1 1 0]
 [0 1 1 1 1 1 1 1 1 0]
 [0 0 0 1 1 1 1 0 0 0]]

Observe que usei um formato 'GTiff' no exemplo, mas outros formatos podem ser mais adequados. Por exemplo, se você tem uma grande varredura que precisa ser transferida através de uma conexão lenta à Internet, tente usar 'PNG' para compactá-la.

Mike T
fonte
Isso é muito útil.
John Powell
Muito útil. obrigado! Ainda estou correndo para esse problema: ERRO: rt_raster_to_gdal: Não foi possível carregar o driver GDAL de saída, mas acho que tenho uma solução alternativa para isso. obrigado novamente!
Mayank Agarwal
@MayankAgarwal resposta atualizada para o erro rt_raster_to_gdal.
Mike T
6

Acho que a pergunta era se você pode ler as tabelas raster postgis SEM drivers gdal habilitados. Como todas as coisas Python, você pode!

Certifique-se de selecionar o resultado da varredura como WKBinary:

selecione St_AsBinary (rast) ...

Use o script abaixo para decifrar o WKBinary em um formato de imagem python. Eu prefiro o opencv, porque ele lida com um número arbitrário de bandas de imagem, mas pode-se usar PIL / low se 1 ou 3 bandas forem mais comuns.

Por enquanto, só manejo imagens de bytes, mas é relativamente trivial expandir para outros tipos de dados.

Espero que isso seja útil.

estrutura de importação
importar numpy como np
importação cv2

# Função para decifrar o cabeçalho WKB
def wkbHeader (bruto):
    # Consulte http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat

    header = {}

    cabeçalho ['endianess'] = struct.unpack ('B', bruto [0]) [0]
    cabeçalho ['versão'] = struct.unpack ('H', bruto [1: 3]) [0]
    cabeçalho ['nbands'] = struct.unpack ('H', bruto [3: 5]) [0]
    cabeçalho ['scaleX'] = struct.unpack ('d', bruto [5:13]) [0]
    cabeçalho ['scaleY'] = struct.unpack ('d', bruto [13:21]) [0]
    cabeçalho ['ipX'] = struct.unpack ('d', bruto [21:29]) [0]
    cabeçalho ['ipY'] = struct.unpack ('d', bruto [29:37]) [0]
    cabeçalho ['skewX'] = struct.unpack ('d', bruto [37:45]) [0]
    cabeçalho ['skewY'] = struct.unpack ('d', bruto [45:53]) [0]
    cabeçalho ['srid'] = struct.unpack ('i', bruto [53:57]) [0]
    cabeçalho ['largura'] = struct.unpack ('H', bruto [57:59]) [0]
    cabeçalho ['altura'] = struct.unpack ('H', bruto [59:61]) [0]

    retornar cabeçalho

# Função para decifrar os dados rasterizados WKB 
def wkbImage (bruto):
    h = wkbHeader (bruto)
    img = [] # array para armazenar faixas de imagens
    offset = 61 # comprimento bruto do cabeçalho em bytes
    para i no intervalo (h ['nbands']):
        # Determinar pixtype para esta banda
        pixtype = struct.unpack ('B', bruto [deslocamento]) [0] >> 4
        # Por enquanto, lidamos apenas com bytes não assinados
        se pixtype == 4:
            banda = np.frombuffer (bruto, dtype = 'uint8', contagem = h ['largura'] * h ['altura'], deslocamento = deslocamento + 1)
            img.append ((np.reshape (banda, ((h ['altura'], h ['largura']))))))
            deslocamento = deslocamento + 2 + h ['width'] * h ['height']
        # to do: manipular outros tipos de dados 

    retornar cv2.merge (tupla (img))

GGL
fonte
Isso é muito útil. Eu tenho tido muitos problemas com o gdal em um ambiente conda, mas essa abordagem funcionou pela primeira vez e é bom poder mergulhar um pouco na estrutura também.
John Powell