12

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).

14

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.

  • That is very helpful. – John Powell Jan 14 '15 at 6:19
  • Very helpful. thanks! I am still running into this issue that: ERROR: rt_raster_to_gdal: Could not load the output GDAL driver but I think I have a workaround for that. thanks again! – Mayank Agarwal Jan 14 '15 at 19:47
  • @MayankAgarwal updated answer for the rt_raster_to_gdal error. – Mike T Jan 14 '15 at 20:03
6

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))

  • That is very helpful. I have been having a lot of problems with gdal in a conda environment, but this approach worked first time, and it is nice to be able to delve into the structure a bit also. – John Powell May 17 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.