PostGIS: parse geometry wkb with OGR

Question

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 1

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.

Answer 2

Within the database, geometries are stored on disk in a format only used by the PostGIS program. In order for external programs to insert and retrieve useful geometries, they need to be converted into a format that other applications can understand. Fortunately, PostGIS supports emitting and consuming geometries in a large number of formats:

from Introduction to PostGIS

With the WKB format:

Well-known binary (WKB):
ST_GeomFromWKB(bytea) returns geometry
ST_AsBinary(geometry) returns bytea
ST_AsEWKB(geometry) returns bytea

ogr recognize geometries an not a bytea result (ST_AsEWKB())

# result -> bytea format:
query = "SELECT ST_AsEWKB(geom) FROM points LIMIT 1"
# result -> geometry from bytea:
query = "SELECT ST_GeomFromWKB(ST_AsEWKB(geom)) from points LIMIT 1;"

Test with one of my tables:

nothing:

query = """SELECT ST_AsText(ST_AsEWKB(geom)) from mytable;"""
cur = conn.cursor()
cur.execute(query)
row = cur.fetchone()
print row[0]
'01010000208A7A0000DD2571661A9B10410CCD751AEBF70241'

and a geometry:

query = """SELECT ST_AsText(ST_GeomFromWKB(ST_AsEWKB(geom))) from mytable;"""
# result
cur.execute(query)
row = cur.fetchone()
print row
('POINT(272070.600041 155389.38792)',)

So, let's try:

 query = """SELECT ST_AsText(ST_GeomFromWKB(ST_AsEWKB(geom))) from mytable;"""
 cur = conn.cursor()
 cur.execute(query)
 row = cur.fetchone()    
 wkb = row[0];
 geom = ogr.CreateGeometryFromWkb(wkb)
 ERROR 3: OGR Error: Unsupported geometry type

Why ?

Because the result of the query is a string:

'01010000208A7A0000DD2571661A9B10410CCD751AEBF70241'

and not a bytecode.

You need to decode this string (look at Create Geometry from WKB in the Python GDAL/OGR Cookbook ).

That is why it is much easier to use:

1) other output formats (WKT, GeoJSON, ...)

 query = """SELECT ST_AsGeoJSON(geom) from mytable;"""  
 cur.execute(query)
 row = cur.fetchone()
 point = ogr.CreateGeometryFromJson(row[0])
 print "%d,%d" % (point.GetX(), point.GetY())
 272070,155389

2) directly osgeo.ogr (How to convert PostGIS table to Shapefile in Python?, for example)

Answer 3

You'll want to use ST_AsBinary(geom) to convert your geometry from the PostGIS internal format to WKB that you can read with ogr:

cur.execute('SELECT ST_AsBinary(geom) FROM mytable LIMIT 1')
result = cur.fetchone()

In Postgres terms, your result is a bytea. The psycpopg2 library will map this to a memoryview Python type:

>>>> type(result[0])
<class 'memoryview'>

Just cast your memoryview to bytes to read the WKB with ogr:

>>>>geom = ogr.CreateGeometryFromWkb(bytes(result[0]))
<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at 0x0000000002D179F0> >

If you're concerned with numerical precision, definitely avoid using ST_AsText(). That function converts your geometry to WKT, truncating your coordinates at a precision that depends on your PostGIS version and platform.