Decode SRID information from a PostGIS WKB string in Python using the standard library

Question

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

Answer 1

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:]