Search code examples
pythonpostgresqlgispostgisshapely

Importing a PostGIS geometry type into Python as a geometry type from Shapely?


So I have a situation where I have a ton of Linestrings of a broken up route, where I need to Union them together using Shapely's LineMerge or Union OR PostGIS ST_Union.

My idea right now is to use Shapely to import the Linestrings as Geometry types. Union them or merge them using Shapely, and then export back to a results table in the database.

However, the geometry type in the PostGIS database is just a bunch of gibberish. Like...

01020000020e61000....

How can I translate this from the database to a Python geometry type using Shapely, do some manipulations, and then export it back to a database?

Currently this is my code, its just importing that geom object string from the database right now and throwing errors because it isn't a Geometry type.

def create_shortest_route_geom(shortest_routes):
    conn = connect_to_database()
    cur = conn.cursor()
    shortest_route_geoms = []
    for route in shortest_routes:
        source = str(int(route[1]))
        target = str(int(route[2]))
        query = 'SELECT the_geom FROM public.ways WHERE target_osm = ' + target + ' AND source_osm = ' + source + ' OR target_osm = ' + source + ' AND source_osm = ' + target + ';'
        cur.execute(query)
        total_geom = cur.fetchone()
        for index, node in enumerate(route):
            try:
                source = str(int(node))
                target = str(int(route[index + 1]))
                query = 'SELECT the_geom FROM public.ways WHERE target_osm = ' + target + ' AND source_osm = ' + source + ' OR target_osm = ' + source + ' AND source_osm = ' + target + ';'
                cur.execute(query)
                geom = cur.fetchone()
                query = "SELECT ST_Union("+str(geom[0])+","+str(total_geom[0])+")"
                cur.execute(query)
                total_geom = cur.fetchone()
            except IndexError:
                print "Last element"
        shortest_route_geoms.insert(total_geom)
    return shortest_route_geoms

EDIT: I may have found my answer here, looking more into it and will update my question with an answer if I figure this out.


Solution

  • Shapely already has libraries for this specific problem.

    PostGIS stores the geometries as HEX values. Use Shapely's loads function to load this with the parameter of hex=True.

    Specifically...

    geom = shapely.wkb.loads(hex_geom[0], hex=True)
    

    Don't use PostGIS ST_Union, cause you have to dump and load over and over for this to work. Shapely has that configured as well with linemerge