8

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
  • 2
    Mistype on this line; make sure it is not in your code: geom = org.CreateGeometryFromWkb(wkb) (should be ogr not org). – Arthur Mar 11 '14 at 16:41
9

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.

  • 1
    For me it is the other way around. With ST_AsEWKB I get ERROR 3: OGR Error: Unsupported geometry type. But with ST_AsBinary it reads fine: MULTILINESTRING ((-4.625433 40.682732,-4.6275242 40.6820109,-4.6293233 40.681392,-4.6301239 40.681117)) – HeikkiVesanto Jan 25 '17 at 14:36
8

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)

  • +1 for the ST_AsGeoJSON(geom) solution. – Michael Jan 30 '18 at 6:03
6

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.

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.