Downloading raster data into python from postgis using psycopg2

Question

I have raster data in a postgres table that I want to get into python as a numpy array. I am using psycopg2 to connect to the db. I am able to download the data but it comes back as a string (probably a serialized binary).

Does any one know how to take this string and convert to a numpy array?

I explored other options to download the raster such as use st_astiff and encode to download the hex file and use xxd but that did not work. I keep on getting the error 'rt_raster_to_gdal: Could not load the output GDAL driver' and I don't have permissions to set the environment variables to be able to turn on the drivers.

TL, DR: want to import raster data into a numpy array (using python).

Answer 1

rt_raster_to_gdal: Could not load the output GDAL driver

As for the first error with ST_AsTIFF, you need to enable your GDAL drivers, which by default are not enabled for PostGIS 2.1. See the manual on ways to do this. For instance, I have an environment variable set up on a Windows computer with:

POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED USGSDEM AAIGrid

which can be confirmed with PostGIS with:

SELECT short_name, long_name
FROM ST_GDALDrivers();

---

PostGIS to Numpy

You can export the output to a virtual memory GeoTIFF file for GDAL to read into a Numpy array. For hints on virtual files used in GDAL, see this blog post.

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

Shows a rasterised buffered point.

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

Note that I used a 'GTiff' format in the example, but other formats might be better suited. For instance, if you have a large raster that needs to be transferred across a slow internet connection, try using 'PNG' to compress it.

Answer 2

I think the question was whether you can read from postgis raster tables WITHOUT gdal drivers enabled. As all things Python, you can!

Make sure you select your raster result as WKBinary:

select St_AsBinary(rast)...

Use the script below to decypher WKBinary into a python image format. I prefer opencv, because it handles arbitrary number of image bands, but one can use PIL/low if 1 or 3 bands are more usual.

I only handle byte imagery for now, but it is relatively trivial to expand to other datatypes.

Hope this is useful.

import struct
import numpy as np
import cv2

# Function to decypher the WKB header
def wkbHeader(raw):
    # See 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]

    return header

# Function to decypher the WKB raster data 
def wkbImage(raw):
    h = wkbHeader(raw)
    img = [] # array to store image bands
    offset = 61 # header raw length in bytes
    for i in range(h['nbands']):
        # Determine pixtype for this band
        pixtype = struct.unpack('B', raw[offset])[0]>>4
        # For now, we only handle unsigned byte
        if 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'])))))
            offset = offset + 2 + h['width']*h['height']
        # to do: handle other data types 

    return cv2.merge(tuple(img))