Загрузка растровых данных в python из postgis с помощью psycopg2

13

У меня есть растровые данные в таблице postgres, которые я хочу добавить в python в виде массива. Я использую psycopg2 для подключения к БД. Я могу загрузить данные, но они возвращаются в виде строки (вероятно, сериализованного двоичного файла).

Кто-нибудь знает, как взять эту строку и преобразовать в массив Numpy?

Я исследовал другие варианты загрузки растра, такие как использование st_astiff и кодирование для загрузки шестнадцатеричного файла и использование xxd, но это не сработало. Я продолжаю получать сообщение об ошибке «rt_raster_to_gdal: не удалось загрузить выходной драйвер GDAL», и у меня нет прав для установки переменных среды, чтобы иметь возможность включать драйверы.

TL, DR: хотите импортировать растровые данные в пустой массив (используя python).

postgis gdal numpy psycopg2 postgis-raster Маянк Агарвал
источник

Ответы:

14

rt_raster_to_gdal: не удалось загрузить выходной драйвер GDAL

Что касается первой ошибки с ST_AsTIFF , вам необходимо включить драйверы GDAL, которые по умолчанию не включены для PostGIS 2.1. Смотрите руководство о том, как это сделать. Например, у меня есть переменная окружения, настроенная на компьютере Windows с:

POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED USGSDEM AAIGrid

что может быть подтверждено PostGIS с помощью:

SELECT short_name, long_name
FROM ST_GDALDrivers();

PostGIS для Numpy

Вы можете экспортировать вывод в файл GeoTIFF виртуальной памяти для чтения GDAL в массив Numpy. Для подсказок на виртуальных файлах, используемых в GDAL, см. Этот пост в блоге .

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

Показывает растеризованную буферную точку.

[[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]]

Обратите внимание, что я использовал формат 'GTiff' в примере, но другие форматы могли бы быть более подходящими. Например, если у вас большой растр, который необходимо передать по медленному интернет-соединению, попробуйте использовать «PNG» для его сжатия.

Майк Т
источник
Это очень полезно.
Джон Пауэлл
Очень полезно. Благодарность! Я все еще сталкиваюсь с этой проблемой, которая: ОШИБКА: rt_raster_to_gdal: Не удалось загрузить выходной драйвер GDAL, но я думаю, что у меня есть обходной путь для этого. Спасибо еще раз!
Маянк Агарвал
@MayankAgarwal обновил ответ для ошибки rt_raster_to_gdal.
Майк Т
6

Я думаю, что вопрос был в том, можно ли читать из таблиц постгисовых растров БЕЗ драйверов gdal. Как и все вещи Python, вы можете!

Убедитесь, что вы выбрали свой растровый результат как WKBinary:

выберите St_AsBinary (раст)

Используйте сценарий ниже, чтобы расшифровать WKBinary в формате изображения Python. Я предпочитаю opencv, потому что он обрабатывает произвольное количество полос изображения, но можно использовать PIL / low, если 1 или 3 полосы более обычны.

Пока я работаю только с байтовыми изображениями, но распространить их на другие типы данных довольно просто.

Надеюсь, это полезно. 

структура импорта
импортировать numpy как np
импорт cv2

# Функция для расшифровки заголовка WKB
def wkbHeader (raw):
    # См. Http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat

    header = {}

    header ['endianess'] = struct.unpack ('B', raw [0]) [0]
    header ['version'] = struct.unpack ('H', raw [1: 3]) [0]
    header ['nbands'] = struct.unpack ('H', raw [3: 5]) [0]
    header ['scaleX'] = struct.unpack ('d', raw [5:13]) [0]
    header ['scaleY'] = struct.unpack ('d', raw [13:21]) [0]
    header ['ipX'] = struct.unpack ('d', raw [21:29]) [0]
    header ['ipY'] = struct.unpack ('d', raw [29:37]) [0]
    header ['skewX'] = struct.unpack ('d', raw [37:45]) [0]
    header ['skewY'] = struct.unpack ('d', raw [45:53]) [0]
    header ['srid'] = struct.unpack ('i', raw [53:57]) [0]
    header ['width'] = struct.unpack ('H', raw [57:59]) [0]
    header ['height'] = struct.unpack ('H', raw [59:61]) [0]

    вернуть заголовок

# Функция для расшифровки растровых данных WKB 
def wkbImage (raw):
    h = wkbHeader (необработанный)
    img = [] # массив для хранения полос изображения
    offset = 61 # необработанная длина заголовка в байтах
    для i в диапазоне (h ['nbands']):
        # Определить pixtype для этой полосы
        pixtype = struct.unpack ('B', raw [offset]) [0] >> 4
        # На данный момент мы обрабатываем только неподписанный байт
        если pixtype == 4:
            band = np.frombuffer (raw, dtype = 'uint8', count = h ['width'] * h ['height'], offset = offset + 1)
            img.append ((np.reshape (band, ((h ['height'], h ['width'])))))
            смещение = смещение + 2 + h ['ширина'] * h ['высота']
        # делать: обрабатывать другие типы данных 

    return cv2.merge (tuple (img))
GGL
источник
Это очень полезно. У меня было много проблем с gdal в среде conda, но этот подход сработал в первый раз, и приятно иметь возможность немного углубиться в структуру.
Джон Пауэлл