Конвертировать геометрию Postgis с помощью osgeo или shapely из EWKB и обратно в EWKB

Question

Предположим, у меня есть следующая геометрия, созданная с помощью Postgis:

SELECT 'POINT(1 2 3 2)'::geometry

Я хотел бы прочитать форму EWKB и просто переписать ее обратно в EWKB (в реальном приложении я бы внес некоторые изменения в геометрию, для этого примера мы просто конвертируем геометрию)

Используя osgeo, я могу сделать следующее:

import binascii
from osgeo import ogr

# The geometry as EWKB as it is returned by postgis:
s = '01010000C0000000000000F03F000000000000004000000000000008400000000000000040'

binary_string = binascii.unhexlify(s)

g = ogr.CreateGeometryFromWkb(binary_string)
print(g.ExportToIsoWkt())
print(binascii.hexlify(g.ExportToWkb()))

Это дает:

POINT ZM (1 2 3 2)
b'0000000bb93ff0000000000000400000000000000040080000000000004000000000000000'

Как вы можете легко видеть, геометрия не идентична исходной. Как мне вернуть результат такой, что

export(g) == s

Я был бы открыт для решения, используя также shapely, однако мои тесты дали некоторые проблемы с M-координатой.

Answer 1

Я думаю, что одним из виновников является функция ST_AsEWKB сам, так как он возвращает представление WKB, включая метаданные SRID. Как указано в документации:

Спецификация WKB не включает SRID. Чтобы получить формат OGC WKB, используйте ST_AsBinary

Сейчас, ST_AsBinary выходы:

test=> SELECT ST_AsBinary(ST_GeomFromText('POINT(1 2 3 2)'));
                                 st_asbinary
------------------------------------------------------------------------------
 \x01b90b0000000000000000f03f000000000000004000000000000008400000000000000040
(1 row)

Первый байт 01 этого представления указывает на порядковый номер. Итак, здесь мы видим, что по умолчанию PostGIS возвращает данные в формате с прямым порядком байтов, тогда как по умолчанию для ExportToWkb похоже на большой порядок байтов. Тем не мение, ExportToWkb принимает параметр byte_order используя который можно переопределить это (или можно переопределить значение по умолчанию для ST_AsBinary с ST_AsBinary(ST_GeomFromText('POINT(1 2 3 2)'), 'XDR')):

import binascii
from osgeo import ogr

# The geometry as it is returned by ST_AsBinary:
s = '01b90b0000000000000000f03f000000000000004000000000000008400000000000000040'

binary_string = binascii.unhexlify(s)

g = ogr.CreateGeometryFromWkb(binary_string)

print(g.ExportToIsoWkt())
print(binascii.hexlify(g.ExportToWkb(byte_order = 1)))

это дает:

POINT ZM (1 2 3 2)
b'01b90b0000000000000000f03f000000000000004000000000000008400000000000000040'

что согласуется с выходом ST_AsBinary...