Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
1  
Do you mean a pure python library which doesn't require any other packages such as QGIS etc? –  Dan Patterson Jan 10 at 1:48
    
Yes, I mean a pure python lib –  Mingshu Jan 13 at 15:31
add comment

4 Answers

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.

share|improve this answer
2  
"Python GDAL/OGR Cookbook" - what a perfect collection! This is getting me excited to start digging into this stuff even more. –  SaultDon Jan 10 at 6:11
add comment

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

share|improve this answer
add comment

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

share|improve this answer
    
Examples really help =) –  SaultDon Jan 10 at 1:36
    
The question did not asked for examples... –  Alexandre Neto Jan 10 at 9:46
add comment

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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