Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
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 at 16:41

3 Answers 3

up vote 2 down vote accepted

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.

OGR (and GEOS) only understands the EWKB representation for higher-order geometries, so for this reason I recommend using ST_AsEWKB. For example:

from osgeo import ogr
ogr.UseExceptions()

import psycopg2
connection = psycopg2.connect('dbname=postgis user=postgres')
cursor = connection.cursor()

cursor.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
res = cursor.fetchone()
g = ogr.CreateGeometryFromWkb(bytes(res[0]))
print(g.ExportToWkt())  # POINT (1 2 3)

cursor.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")
res = cursor.fetchone()
g = ogr.CreateGeometryFromWkb(bytes(res[0]))
# RuntimeError: OGR Error: Unsupported geometry type

Also, note that OGR (and GEOS) do not support M coordinates, and these will be parsed but ignored.

share|improve this answer

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)

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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