Dissolving shapefiles using shapely and psycopg2

Geospatial data is commonly stored in PostgreSQL which, through PostGIS, provides a great set of tools for geospatial analysis. However, if this analysis is part of a larger workflow chain it can sometimes be easier to work with the excellent shapley package that provides a easy set of calls for advanced geospatial analysis. The script below looks at how to extract data from PostgreSQL, insert into shapely and run a common geospatial query: dissolve (cascaded union).


#!/usr/bin/python
import psycopg2
import shapely.wkb as wkb
import shapely.wkt as wkt 
from shapely.ops import cascaded_union
import json
 
def main():
    
    #Define our connection string
    conn_string = "host='12.34.56.78' dbname='nathan_local' user='postgres' password='XXXXXXXXX'"
 
    # print the connection string
    print("Connecting to database %s" % (conn_string))
 
    # get a connection
    conn = psycopg2.connect(conn_string)
 
    # conn.cursor will return a cursor object, you can use this cursor to perform queries.
    cursor = conn.cursor()
    print("Connected!\n")

    # use the cursor object to extract data from postgreSQL
    cursor.execute("SELECT tile_name, st_asewkt(wkb_geometry) FROM test.tablename")

    tiles = list(cursor.fetchall()) #fetchall returns all matching records

    coverage = {}
    for name, geometry in tiles:
        coverage[name] = wkt.loads(geometry.split(';')[1])

    geometry = cascaded_union(coverage.values())
	
    with open('/file/path/test.json', 'w') as outfile:
        json.dump(mapping(geometry), outfile)
	 
if __name__ == "__main__":
	main()

Advertisements