Converting PostGIS table to Shapefile in Python?

Question

I want to convert a PostGIS table into a shapefile (without using pgsql2shp).

In order to create a geometry in the shapefile I have to give the Xmin,Ymin and Xmax,Ymax, and the geometry which I have in my PostGIS table is an irregularly shaped one (I can get the exterior using bounding box but that will include some extra area more than my area of ineterest). Is there any method by which I can get the task done?

I want to do the thing programmatically and using Python.

Answer 1

If you want to do it programmatically with Python then probably GDAL/OGR is the best way to do it. Look at the sample example code that copies on table from PostgreSQL into SHP file. Example is not perfect but you can easily modify it to fit your needs.

import os
os.environ['PATH'] = "c:\\Program Files\\GDAL\\bin" + ';' + os.environ['PATH']
os.environ['GDAL_DRIVER_PATH'] = "c:\\Program Files\\GDAL\\bin\\gdal\\plugins-optional"
os.environ['GDAL_DATA'] = "c:\\Program Files\\GDAL\\bin\\gdal-data"
import ogr

conn=ogr.Open("PG: host=192.168.5.3 dbname=some_database user=postgres password=xxxx")
if conn is None:
    print 'Could not open a database or GDAL is not correctly installed!'
    sys.exit(1)

output = "d:\\points.shp"

# Schema definition of SHP file
out_driver = ogr.GetDriverByName( 'ESRI Shapefile' )
out_ds = out_driver.CreateDataSource(output)
out_srs = None
out_layer = out_ds.CreateLayer("point", out_srs, ogr.wkbPoint)
fd = ogr.FieldDefn('name',ogr.OFTString)
out_layer.CreateField(fd)


layer = conn.GetLayerByName("point_data")
#layer = conn.ExecuteSQL(sql)

feat = layer.GetNextFeature()
while feat is not None:
    featDef = ogr.Feature(out_layer.GetLayerDefn())
    featDef.SetGeometry(feat.GetGeometryRef())
    featDef.SetField('name',feat.TITLE)
    out_layer.CreateFeature(featDef)
    feat.Destroy()
    feat = layer.GetNextFeature()

conn.Destroy()
out_ds.Destroy()

If your programing environment is Windows, then you can download GDAL/OGR from here. Some good starting materials you can find here. Hope it helps.

Answer 2

I can advise you to take a look at GDAL/OGR library: Link

Answer 3

When using the gdal python bindings you can use the ogr2ogr port.

from osgeo_utils.samples.ogr2ogr import main
 main_argument_list = [
        'nothing at all', # This value should contain the name of the calling script but is not used so you can simply ignore it
        destination_source,
        source,
        layer,
        '-f',
        destination_format,
    ]
main(main_argument_list)

This is not (as they say) clean or pythonic. The port is made to be used as CLI tool not as a programmable interface. However this works fine. You could create your own port that does respect the pythonic rules of the wild west. But there is also Fiona that does implement a pythonic programmable interface.

However, I prefer to keep close to the original which is a c++ program. And use the ogr2ogr port directly like I would use the cli. This way I can translate CLI usage directly to python scripts. However Fiona would be a more clean approach.