1

TL;DR

I want to decode a WKB string which comes from a PostGIS database using the Python standard library. The WKB string seems to actually embed the SRID information, but I'm not 100% sure and shapely in agnostic to that.

Details

I have a PostGIS point:

SELECT ST_AsText(ST_SetSrid(ST_MakePoint(2670599, 1207749), 2056))

which one can load in Python this way using a psycopg2 connector:

sql = """SELECT ST_SetSrid(ST_MakePoint(2670599, 1207749), 2056)"""

with psycopg2.connect(**DATABASE_PARAMETERS) as conn:
    with conn.cursor() as curs:
        curs.execute(sql)
        pt = curs.fetchone()

point = pt[0]

Now point looks like '010100002008080000000000800360444100000000C56D3241' which is a WKB string.

Great. But shapely is agnostic to the SRS information, so this:

from shapely import wkb

shape = wkb.loads(point, hex=True)
wkb.dumps(shape, hex=True)

prints '0101000000000000800360444100000000C56D3241', which is not equal to the first string.
And if one need to dump this, one has to manually specify the SRID this way:

wkb.dumps(shape, hex=True, srid=2056)

which prints '010100002008080000000000800360444100000000C56D3241' and which is now equal to the first string:

>>> point == wkb.dumps(shape, hex=True, srid=2056)
True

So, is the difference an actual encoding of the SRID information or not?

If yes, how to decode it properly using only the Python standard library?

I know it can be easily done using, e.g. GeoPandas:

import geopandas as gpd
sql = """SELECT ST_SetSrid(ST_MakePoint(2670599, 1207749), 2056) as geom"""  # <-- give it a name so that GeoPandas doesn't compain.
with psycopg2.connect(**DATABASE_PARAMETERS) as conn:
    with conn.cursor() as curs:
        curs.execute(sql)
        pt = curs.fetchone()

G = gpd. gpd.read_postgis(sql, connector())

Then:

>>> GG.crs

<Projected CRS: EPSG:2056>
Name: CH1903+ / LV95
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Liechtenstein; Switzerland.
- bounds: (5.96, 45.82, 10.49, 47.81)
Coordinate Operation:
- name: Swiss Oblique Mercator 1995
- method: Hotine Oblique Mercator (variant B)
Datum: CH1903+
- Ellipsoid: Bessel 1841
- Prime Meridian: Greenwich

I'm simply curious about this, because for the moment, my workaround to avoid loading an other library, is to query it separately, which works nicely:

sql = """SELECT ST_SetSrid(ST_MakePoint(2670599, 1207749), 2056) as geom,
         ST_SRID(ST_SetSrid(ST_MakePoint(2670599, 1207749), 2056)) as srid"""

with psycopg2.connect(**DATABASE_PARAMETERS) as conn:
    with conn.cursor() as curs:
        curs.execute(sql)
        pt = curs.fetchone()

point = pt[0]
srid = pt[1] # <- this prints the integer: 2056
3

From the QuickWKT QGIS plugin, we can see the code (lightly modified) :

import binascii


def decodeBinary(wkb):
    """Decode the binary wkb and return as a hex string"""
    value = binascii.a2b_hex(wkb)
    value = value[::-1]
    value = binascii.b2a_hex(value)
    return value.decode("UTF-8")

def encodeBinary(value):
    wkb = binascii.a2b_hex("%08x" % value)
    wkb = wkb[::-1]
    wkb = binascii.b2a_hex(wkb)
    return wkb.decode("UTF-8")


wkb = "010100002008080000000000800360444100000000C56D3241"
SRID_FLAG = 0x20000000
srid = ""
geomType = int("0x" + decodeBinary(wkb[2:10]), 0)
if geomType & SRID_FLAG:
    srid = int("0x" + decodeBinary(wkb[10:18]), 0)
    # String the srid from the wkb string
    wkb = wkb[:2] + encodeBinary(geomType ^ SRID_FLAG) + wkb[18:]

1
  • Hmm, great! I guess the index position such as 2:10 or 10:18 are defined somewhere? Probably in an OGC document? Probably this one: ogc.org/standards/sfa – s.k Mar 24 at 20:06

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.