Mengunduh data raster ke python dari postgis menggunakan psycopg2

13

Saya punya data raster di tabel postgres yang ingin saya masuki ke python sebagai array yang numpy. Saya menggunakan psycopg2 untuk terhubung ke db. Saya dapat mengunduh data tetapi kembali sebagai string (mungkin biner serial).

Apakah ada yang tahu cara mengambil string ini dan mengonversi ke array numpy?

Saya menjelajahi opsi lain untuk mengunduh raster seperti menggunakan st_astiff dan menyandikan untuk mengunduh file hex dan menggunakan xxd tetapi itu tidak berhasil. Saya terus mendapatkan kesalahan 'rt_raster_to_gdal: Tidak dapat memuat driver GDAL keluaran' dan saya tidak memiliki izin untuk mengatur variabel lingkungan untuk dapat mengaktifkan driver.

TL, DR: ingin mengimpor data raster ke array numpy (menggunakan python).

Mayank Agarwal
sumber

Jawaban:

14

rt_raster_to_gdal: Tidak dapat memuat driver GDAL keluaran

Adapun kesalahan pertama dengan ST_AsTIFF , Anda harus mengaktifkan driver GDAL Anda, yang secara default tidak diaktifkan untuk PostGIS 2.1. Lihat manual tentang cara melakukan ini. Sebagai contoh, saya memiliki variabel lingkungan yang diatur di komputer Windows dengan:

POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED USGSDEM AAIGrid

yang dapat dikonfirmasi dengan PostGIS dengan:

SELECT short_name, long_name
FROM ST_GDALDrivers();

PostGIS ke Numpy

Anda dapat mengekspor output ke file GeoTIFF memori virtual untuk GDAL untuk dibaca ke array Numpy. Untuk petunjuk tentang file virtual yang digunakan dalam GDAL, lihat posting blog ini .

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

Menunjukkan titik buffer yang dirasterisasi.

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

Perhatikan bahwa saya menggunakan format 'GTiff' dalam contoh, tetapi format lain mungkin lebih cocok. Misalnya, jika Anda memiliki raster besar yang perlu ditransfer melalui koneksi internet yang lambat, coba gunakan 'PNG' untuk mengompresnya.

Mike T
sumber
Itu sangat membantu.
John Powell
Sangat membantu. Terima kasih! Saya masih mengalami masalah ini bahwa: ERROR: rt_raster_to_gdal: Tidak dapat memuat driver GDAL output tapi saya pikir saya punya solusi untuk itu. Terima kasih lagi!
Mayank Agarwal
@MayankAgarwal memperbarui jawaban untuk kesalahan rt_raster_to_gdal.
Mike T
6

Saya pikir pertanyaannya adalah apakah Anda dapat membaca dari tabel raster postgis TANPA driver gdal diaktifkan. Seperti semua hal Python, Anda bisa!

Pastikan Anda memilih hasil raster Anda sebagai WKBinary:

pilih St_AsBinary (rast) ...

Gunakan skrip di bawah ini untuk mendekripsi WKBinary menjadi format gambar python. Saya lebih suka opencv, karena menangani jumlah band gambar yang sewenang-wenang, tetapi orang dapat menggunakan PIL / rendah jika 1 atau 3 band lebih biasa.

Saya hanya menangani citra byte untuk saat ini, tetapi relatif sepele untuk memperluas ke tipe data lain.

Semoga ini bermanfaat.

impor struct
impor numpy sebagai np
impor cv2

# Berfungsi untuk mendekripsi header WKB
def wkbHeader (mentah):
    # Lihat http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat

    header = {}

    header ['endianess'] = struct.unpack ('B', raw [0]) [0]
    tajuk ['versi'] = struct.unpack ('H', mentah [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]

    kembali header

# Berfungsi untuk mendekripsi data raster WKB 
def wkbImage (mentah):
    h = wkbHeader (mentah)
    img = [] # array untuk menyimpan pita gambar
    offset = 61 # panjang raw header dalam byte
    untuk saya dalam jangkauan (h ['nbands']):
        # Tentukan pixtype untuk band ini
        pixtype = struct.unpack ('B', raw [offset]) [0] >> 4
        # Untuk saat ini, kami hanya menangani byte yang tidak ditandai
        jika pixtype == 4:
            band = np.frombuffer (raw, dtype = 'uint8', hitung = h ['lebar'] * h ['tinggi'], offset = offset + 1)
            img.append ((np.reshape (band, ((h ['tinggi'], h ['lebar'])))))))
            offset = offset + 2 + h ['lebar'] * h ['tinggi']
        # yang harus dilakukan: menangani tipe data lainnya 

    return cv2.merge (tuple (img))

GGL
sumber
Itu sangat membantu. Saya telah mengalami banyak masalah dengan gdal di lingkungan conda, tetapi pendekatan ini bekerja pertama kali, dan itu bagus untuk dapat mempelajari struktur sedikit juga.
John Powell