python - PostGIS: parse geometry wkb with OGR

I'm trying to pull a LineString geometry out of PostGIS and parse it with OGR (python bindinds).

from osgeo import ogr
import psycopg2

connection = psycopg2.connect("...")
cursor = connection.cursor()

query = "SELECT geom FROM points LIMIT 1"


cursor.execute(query)
row = cursor.fetchone()

wkb = row[0]
geom = ogr.CreateGeometryFromWkb(wkb)

cursor.close()
connection.close()

I've tried:

wkb = bin(int(row[0], 16))

and:

SELECT ST_AsEWKB(geom) FROM points LIMIT 1

OGR does not want to parse it. Keeps giving the following error:

ERROR 3: OGR Error: Unsupported geometry type

Answer

Internally, PostGIS stores geometries in a binary specification, but it is queried and viewed outside as a hex-encoded string. There are two popular variations of well-known binary (WKB):

• EWKB (via ST_AsEWKB) - an extended WKB specification designed by PostGIS.


• OGC WKB (via ST_AsBinary) - specified by the OGC and ISO. For a while it was 2D-only, but later extended to support Z, M and ZM geometries.

The two specifications are the same for 2D geometries, but are different for higher-order geometries with Z, M and ZM coordinates.

---

Older versions of GDAL/OGR (1.x) only understand the EWKB for 3D geometries, so for these I recommend using ST_AsEWKB. (But if you only have 2D geometries, either format are fine). For example:

import psycopg2
from osgeo import ogr

ogr.UseExceptions()    

conn = psycopg2.connect('dbname=postgis user=postgres')
curs = conn.cursor()

curs.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 0101000080000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToWkt())  # POINT (1 2 3)

curs.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")

b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 01e9030000000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
# RuntimeError: OGR Error: Unsupported geometry type

Also, note that older GDAL/OGR versions do not support M coordinates, and these will be parsed but ignored.

---

With GDAL 2.0 and more recent, ISO WKT/WKB is supported. This means that CreateGeometryFromWkb can read either WKB flavour (without specifying) and ExportToIsoWkt() shows output with a modern WKT syntax.

import psycopg2
from osgeo import ogr


ogr.UseExceptions()
conn = psycopg2.connect('dbname=postgis user=postgres')
curs = conn.cursor()

curs.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 0101000080000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToIsoWkt())  # POINT Z (1 2 3)


curs.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 01e9030000000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToIsoWkt())  # POINT Z (1 2 3)

Additionally, GDAL 2.1 or later will create/export WKT/WKB with M or ZM coordinates as expected.