PostGIS: analyser la géométrie wkb avec OGR

J'essaie d'extraire une LineStringgéométrie de PostGIS et de l'analyser avec OGR (bindinds python).

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

J'ai essayé:

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

et:

SELECT ST_AsEWKB(geom) FROM points LIMIT 1

OGR ne veut pas l'analyser. Continue de donner l'erreur suivante:

ERROR 3: OGR Error: Unsupported geometry type

En interne, PostGIS stocke les géométries dans une spécification binaire, mais elle est interrogée et vue à l'extérieur sous la forme d'une chaîne codée hexadécimale. Il existe deux variantes populaires de binaire bien connu (WKB) :

• EWKB (via ST_AsEWKB) - une spécification WKB étendue conçue par PostGIS .
• OGC WKB (via ST_AsBinary) - spécifié par l'OGC et l'ISO. Pendant un certain temps , il était 2D seulement, mais plus tard étendu à l' appui Z, Met ZMgéométries.

Les deux spécifications sont les mêmes pour les géométries 2D, mais sont différentes pour des géométries d'ordre supérieur avec Z, Met les ZMcoordonnées.

Les versions plus anciennes de GDAL / OGR (1.x) ne comprennent que l'EWKB pour les géométries 3D, donc je recommande d'utiliser celles-ci ST_AsEWKB. (Mais si vous n'avez que des géométries 2D, les deux formats conviennent.) Par exemple:

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

Notez également que les anciennes versions de GDAL / OGR ne prennent pas en charge les Mcoordonnées, et celles-ci seront analysées mais ignorées.

Avec GDAL 2.0 et plus récent , ISO WKT / WKB est pris en charge . Cela signifie que CreateGeometryFromWkbpeut lire soit la saveur WKB (sans spécifier) ​​et ExportToIsoWkt()affiche la sortie avec une syntaxe WKT moderne.

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)

De plus, GDAL 2.1 ou version ultérieure créera / exportera WKT / WKB avec Mou ZMcoordonne comme prévu.

Dans la base de données, les géométries sont stockées sur disque dans un format uniquement utilisé par le programme PostGIS. Pour que les programmes externes puissent insérer et récupérer des géométries utiles, ils doivent être convertis dans un format que d'autres applications peuvent comprendre. Heureusement, PostGIS prend en charge l'émission et la consommation de géométries dans un grand nombre de formats:

de Introduction à PostGIS

Avec le format WKB:

Binaire bien connu (WKB): ST_GeomFromWKB (bytea) renvoie une géométrie
ST_AsBinary (géométrie) renvoie bytea
ST_AsEWKB (géométrie) renvoie bytea

ogr reconnaît les géométries et n'est pas un résultat bytea ( 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;"

Testez avec une de mes tables:

rien:

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

et une géométrie:

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

Essayons donc:

 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

Pourquoi ?

Parce que le résultat de la requête est une chaîne:

'01010000208A7A0000DD2571661A9B10410CCD751AEBF70241'

et non un bytecode.

Vous devez décoder cette chaîne (regardez Créer une géométrie à partir de WKB dans le livre de recettes Python GDAL / OGR ).

C'est pourquoi il est beaucoup plus facile à utiliser:

1) autres formats de sortie (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) directement osgeo.ogr ( Comment convertir une table PostGIS en Shapefile en Python?, Par exemple)

Vous voudrez utiliser ST_AsBinary(geom)pour convertir votre géométrie du format interne PostGIS en WKB que vous pouvez lire avec ogr:

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

En termes Postgres, votre résultat est a bytea. La bibliothèque psycpopg2 mappera ceci à un memoryviewtype Python:

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

Il suffit de lancer votre memoryviewpour byteslire le WKB avec ogr:

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

Si vous vous souciez de la précision numérique, évitez absolument de l'utiliser ST_AsText(). Cette fonction convertit votre géométrie en WKT, tronquant vos coordonnées avec une précision qui dépend de votre version et plate-forme PostGIS.