16

Excluding ArcPy, is there any python library which can do geoprocessing, such as buffer/intersect, with shapefiles?

closed as too broad by PolyGeo May 2 at 20:40

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

17

The Python GDAL/OGR Cookbook has some sample code to Buffer a Geometry.

from osgeo import ogr

wkt = "POINT (1198054.34 648493.09)"
pt = ogr.CreateGeometryFromWkt(wkt)
bufferDistance = 500
poly = pt.Buffer(bufferDistance)
print "%s buffered by %d is %s" % (pt.ExportToWkt(), bufferDistance, poly.ExportToWkt())

and to Calculate intersection between two Geometries

from osgeo import ogr

wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

intersection = poly1.Intersection(poly2)

print intersection.ExportToWkt()

The geometries can be read and written to shapefiles and a variety of other formats.

14

To simplify, Shapely: manual allows all geometry processing of PostGIS in Python.

The first premise of Shapely is that Python programmers should be able to perform PostGIS type geometry operations outside of an RDBMS...

The first example of PolyGeo

from shapely.geometry import Point, LineString, Polygon, mapping
from shapely.wkt import loads  
pt = Point(1198054.34,648493.09)
# or
pt = loads("POINT (1198054.34 648493.09)")
bufferDistance = 500
poly = pt.buffer(bufferDistance)
print poly.wkt
'POLYGON ((1198554.3400000001000000 648493.0899999999700000, 1198551.9323633362000000 
# GeoJSON
print mapping(poly)
{'type': 'Polygon', 'coordinates': (((1198554.34, 648493.09), (1198551.9323633362, 648444.0814298352), (1198544.7326402017, 648395.544838992), ....}

The polygon's example from PolyGeo:

poly1 = Polygon([(1208064.271243039,624154.6783778917), (1208064.271243039,601260.9785661874), (1231345.9998651114,601260.9785661874),(1231345.9998651114,624154.6783778917),(1208064.271243039,624154.6783778917)])    
poly2 = loads("POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

intersection = poly1.intersection(poly2)
print intersection.wkt
print mapping(intersection) -> GeoJSON

The second premise is that the persistence, serialization, and map projection of features are significant, but orthogonal problems. You may not need a hundred GIS format readers and writers or the multitude of State Plane projections, and Shapely doesn’t burden you with them.

So you combine it with other Python modules to read or write shapefiles and manipulate projections as osgeo.ogr, Fiona or PyShp.
Searching in Gis StackExchange, you can find many examples but I give you another one to illustrate the combination of shapely and Fiona and the use of the shapely functions intersection() and buffer() (This could have been done with PyShp).

Given two polyline shapefiles:

enter image description here

Compute the intersection (function intersection() of shapely)

from shapely.geometry import Point, Polygon, MultiPolygon, MumtiPoint, MultiLineString,shape, mapping
import fiona
# read the shapefiles and transform to MultilineString shapely geometry (shape())
layer1 = MultiLineString([shape(line['geometry']) for line in fiona.open('polyline1.shp')])  
layer2 = MultiLineString([shape(line['geometry']) for line in fiona.open('polyline2.shp')])
points_intersect = layer1.intersection(layer2)

Save the result as a new shapefile

# schema of the new shapefile
schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
# write the new shapefile (function mapping() of shapely)
with fiona.open('intersect.shp','w','ESRI Shapefile', schema) as e:
  e.write({'geometry':mapping(points_intersect), 'properties':{'test':1}})

Result:

enter image description here

Buffer individual points (function buffer() of shapely)

 # new schema
 schema = {'geometry': 'Polygon','properties': {'test': 'int'}}
 with fiona.open('buffer.shp','w','ESRI Shapefile', schema) as e:
     for point in points:
          e.write({'geometry':mapping(point.buffer(300)), 'properties':{'test':1}})

Result

enter image description here

Buffer the MultiPoint geometry

schema = {'geometry': 'MultiPolygon','properties': {'test': 'int'}}
points.buffer(300)
with fiona.open('buffer2.shp','w','ESRI Shapefile', schema) as e:
     e.write({'geometry':mapping(points.buffer(300)), 'properties':{'test':1}})

enter image description here

9

Shapely gives python access to GEOS which can do buffers/intersects/etc. GEOS is the library most OSGeo programs use to perform those operations.

9

Here is my list of Python geoprocessing software.

  • Shapely, python
  • OGR, python
  • QGIS, pyqgis, python
  • SagaGIS, python
  • Grass, python
  • spatialite, pyspatialite, python
  • PostreSQL/PostGIS, Psycopg, python
  • R Project, rpy2, python
  • Whitebox GAT, python -GeoScript, jython
1

My 'go to' geo-processing library is the 'Remote Sensing and GIS Library' (RSGISLib). Its easy to install and use and the documentation is really good. It has functionality for vector and raster processing - I very rarely have to go near a gui. It can be found here: http://rsgislib.org.

An example in this instance is:

rsgislib.vectorutils.buffervector(inputvector, outputvector, bufferDist, force)

A command to buffer a vector by a specified distance.

Where:

  • inputvector is a string containing the name of the input vector
  • outputvector is a string containing the name of the output vector
  • bufferDist is a float specifying the distance of the buffer, in map units
  • force is a bool, specifying whether to force removal of the output vector if it exists

Example:

from rsgislib import vectorutils
inputVector = './Vectors/injune_p142_stem_locations.shp'
outputVector = './TestOutputs/injune_p142_stem_locations_1mbuffer.shp'
bufferDist = 1
vectorutils.buffervector(inputVector, outputVector, bufferDist, True)

Not the answer you're looking for? Browse other questions tagged or ask your own question.