1

I'm currently working to find the intersections between linears in a shapefile. I'm a newbie to the GDAL/OGR library. I've done some googling and came up with what seems to be a useful example(see below). There is a problem though. I don't see the method ExportToWkt(). My intellisense(or whatever it's called in IDLE) for the variable "intersection"doesn't have that option. Can anybody tell me what I'm doing wrong. Or just give me some general guidance to comparing two features in a shapefile and finding the intersection(if they intersect at all). Thanks.

enter image description here

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()
2

use the dir() function

dir(intersection)
['AddGeometry', 'AddGeometryDirectly', 'AddPoint', 'AddPoint_2D', 'Area', 'AssignSpatialReference', 'Boundary', 'Buffer', 'Centroid', 'Clone', 'CloseRings', 'Contains', 'ConvexHull', 'Crosses', 'Destroy', 'Difference', 'Disjoint', 'Distance', 'Empty', 'Equal', 'Equals', 'ExportToGML', 'ExportToJson', 'ExportToKML', 'ExportToWkb', 'ExportToWkt', 'FlattenTo2D', 'GetArea', 'GetBoundary', 'GetCoordinateDimension', 'GetDimension', 'GetEnvelope', 'GetEnvelope3D', 'GetGeometryCount', 'GetGeometryName', 'GetGeometryRef', 'GetGeometryType', 'GetPoint', 'GetPointCount', 'GetPoint_2D', 'GetPoints', 'GetSpatialReference', 'GetX', 'GetY', 'GetZ', 'Intersect', 'Intersection', 'Intersects', 'IsEmpty', 'IsRing', 'IsSimple', 'IsValid', 'Length', 'Overlaps', 'PointOnSurface', 'Segmentize', 'SetCoordinateDimension', 'SetPoint', 'SetPoint_2D', 'Simplify', 'SimplifyPreserveTopology', 'SymDifference', 'SymmetricDifference', 'Touches', 'Transform', 'TransformTo', 'Union', 'UnionCascaded', 'Within', 'WkbSize', '__class__', '__del__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattr__', '__getattribute__', '__hash__', '__init__', '__iter__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '__swig_destroy__', '__swig_getmethods__', '__swig_setmethods__', '__weakref__', 'next', 'this']

And you have all the class/function of the module. You may find what you are looking for with a little function adapted from a script Script de Python para filtrar por patrón de texto los métodos de Clases en PyQGIS de José Guerrero (already used in Iterating over selected features in QGIS Processing)

import re # regular expression
def get_patt(keyword, L):
    return [item for item in dir(L) if re.search(keyword,item)]

print get_patt(("Export",dir(intersection))
['ExportToGML', 'ExportToJson', 'ExportToKML', 'ExportToWkb', 'ExportToWkt']

so

print intersection.ExportToWkt()
POLYGON ((1208064.271243039052933 614453.958118694950826,1208064.271243039052933 624154.678377891657874,1219317.106743707787246 624154.678377891657874,1219317.106743707787246 614453.958118694950826,1208064.271243039052933 614453.958118694950826))
# or
print intersection.ExportToJson() # GeoJSON format
{ "type": "Polygon", "coordinates": [ [ [ 1208064.271243039052933, 614453.958118694950826 ], [ 1208064.271243039052933, 624154.678377891657874 ], [ 1219317.106743707787246, 624154.678377891657874 ], [ 1219317.106743707787246, 614453.958118694950826 ], [ 1208064.271243039052933, 614453.958118694950826 ] ] ] }

are correct

Rather than using ogr, use the Shapely module, it is easier.

from shapely.wkt import loads
from shapely.geometry import polygon, mapping
poly1 = loads(wkt1)
poly2 = loads(wkt2)
poly1.intersects(poly2)
True
print poly1.intersection(poly2).wkt
POLYGON ((1208064.271243039 614453.958118695, 1208064.271243039 624154.6783778917, 1219317.106743708 624154.6783778917, 1219317.106743708 614453.958118695, 1208064.271243039 614453.958118695))
# GeoJSON
print mapping(poly1.intersection(poly2))
{'type': 'Polygon', 'coordinates': (((1208064.2712430391, 614453.95811869495), (1208064.2712430391, 624154.67837789166), (1219317.1067437078, 624154.67837789166), (1219317.1067437078, 614453.95811869495), (1208064.2712430391, 614453.95811869495)),)}

And if you need to work with shapefiles, use the Fiona module (easier and also based on GDAL/OGR, with many examples in GIS stackexchange)

If you want to continue with osgeo/ogr, look at the Python GDAL/OGR Cookbook

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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