Mam dane rastrowe w tabeli postgres, którą chcę wprowadzić do Pythona jako tablicę numpy. Korzystam z psycopg2, aby połączyć się z bazą danych. Jestem w stanie pobrać dane, ale są one zwracane jako ciąg (prawdopodobnie serializowany plik binarny).

Czy ktoś wie, jak pobrać ten ciąg i przekonwertować na tablicę numpy?

Zbadałem inne opcje pobierania rastra, takie jak użyj st_astiff i kodowania, aby pobrać plik hex i użyć xxd, ale to nie działało. Ciągle pojawia się błąd „rt_raster_to_gdal: Nie można załadować wyjściowego sterownika GDAL” i nie mam uprawnień do ustawiania zmiennych środowiskowych, aby móc włączyć sterowniki.

TL, DR: chcesz zaimportować dane rastrowe do tablicy numpy (używając Pythona).

rt_raster_to_gdal: Nie można załadować wyjściowego sterownika GDAL

Jeśli chodzi o pierwszy błąd w ST_AsTIFF , musisz włączyć sterowniki GDAL, które domyślnie nie są włączone dla PostGIS 2.1. Zobacz, jak to zrobić. Na przykład mam zmienną środowiskową skonfigurowaną na komputerze z systemem Windows z:

POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED USGSDEM AAIGrid

co można potwierdzić za pomocą PostGIS za pomocą:

SELECT short_name, long_name
FROM ST_GDALDrivers();

PostGIS na Numpy

Możesz wyeksportować dane wyjściowe do pliku GeoTIFF pamięci wirtualnej, aby GDAL mógł je odczytać do tablicy Numpy. Wskazówki dotyczące plików wirtualnych używanych w GDAL można znaleźć w tym poście na blogu .

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

Pokazuje zrasteryzowany buforowany punkt.

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

Zauważ, że użyłem formatu „GTiff” w tym przykładzie, ale inne formaty mogą być bardziej odpowiednie. Na przykład, jeśli masz duży raster, który musi zostać przesłany przez wolne połączenie internetowe, spróbuj użyć „PNG”, aby go skompresować.

Myślę, że pytanie brzmiało, czy można czytać z tabel rastrowych Postgis BEZ włączonych sterowników gdal. Jak wszystko Python, możesz!

Upewnij się, że wybrałeś wynik rastra jako WKBinary:

wybierz St_AsBinary (rast) ...

Użyj poniższego skryptu, aby rozszyfrować WKBinary w formacie obrazu python. Wolę opencv, ponieważ obsługuje dowolną liczbę pasm obrazu, ale można użyć PIL / low, jeśli 1 lub 3 pasma są bardziej typowe.

Obecnie obsługuję tylko bajtowe obrazy, ale rozszerzenie na inne typy danych jest stosunkowo proste.

Mam nadzieję, że to się przyda.

import struct
zaimportuj numpy jako np
import cv2

# Funkcja do odszyfrowania nagłówka WKB
def wkbHeader (raw):
    # Patrz http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat

    nagłówek = {}

    nagłówek ['endianess'] = struct.unpack ('B', raw [0]) [0]
    nagłówek [„wersja”] = struct.unpack („H”, raw [1: 3]) [0]
    nagłówek ['nbands'] = struct.unpack ('H', raw [3: 5]) [0]
    nagłówek ['scaleX'] = struct.unpack ('d', raw [5:13]) [0]
    nagłówek ['scaleY'] = struct.unpack ('d', raw [13:21]) [0]
    nagłówek ['ipX'] = struct.unpack ('d', raw [21:29]) [0]
    nagłówek ['ipY'] = struct.unpack ('d', raw [29:37]) [0]
    nagłówek ['skewX'] = struct.unpack ('d', raw [37:45]) [0]
    nagłówek ['skewY'] = struct.unpack ('d', raw [45:53]) [0]
    nagłówek ['srid'] = struct.unpack ('i', raw [53:57]) [0]
    nagłówek ['szerokość'] = struct.unpack ('H', raw [57:59]) [0]
    nagłówek ['wysokość'] = struct.unpack ('H', raw [59:61]) [0]

    zwraca nagłówek

# Funkcja do odszyfrowywania danych rastrowych WKB 
def wkbImage (raw):
    h = wkbHeader (surowy)
    img = [] # tablica do przechowywania pasm obrazu
    offset = 61 # nieprzetworzona długość nagłówka w bajtach
    dla i w zakresie (h [„nbands”]):
        # Określ piksele dla tego zespołu
        pixtype = struct.unpack ('B', raw [offset]) [0] >> 4
        # Na razie obsługujemy tylko bajty niepodpisane
        jeśli pixtype == 4:
            band = np. frombuffer (raw, dtype = 'uint8', count = h ['width'] * h ['height'], offset = offset + 1)
            img.append ((np .reshape (pasmo, ((h ['wysokość'], h ['szerokość'])))))
            offset = offset + 2 + h [„szerokość”] * h [„wysokość”]
        # do zrobienia: obsługa innych typów danych 

    return cv2.merge (krotka (img))