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.
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()
Answer
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
No comments:
Post a Comment