[GIS] need help for script checking topology using ogr in python

ogrpythonshapefiletopology

My main goal is to create a new feature indicating if the line features intersect each other as well as if they are not connected to a feature.

My problem is how to identify the intersecting and not connected points and converting them to another geometry layer.

The only thing I came up with is:

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')

inDatasource = driver.Open('test.shp', 0)
inLayer = inDatasource.GetLayer()

#loop through features in the layer
feature = inLayer.GetNextFeature()
while feature:
    #condition
    feature.Destroy()
    feature = inLayer.GetNextFeature()

#check topology

outData = driver.CreateDataSource('TopologyResult.shp')
outLayer = outData.CreateLayer ('TopologyResult', geom_type = ogr.wkbPoint)

fieldDefn = ogr.FieldDefn ('id', ogr.OFTInteger)
outLayer.CreateField(fieldDefn)

featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(point)

Best Answer

Working with ogr is as hard as it is computational efficient, I don't have a full script for what you want so I will point some directions and maybe we can come up with something.

I think a good starting point would be to use osgeo.ogr.Layer.SetSpatialFilter(geom). It checks for bounding box intersection, so you would get a subset of the data that would potentially intersects a given geometry. Then you could use osgeo.ogr.Geometry.Intersects(geom) to check for the line intersection or osgeo.ogr.Geometry.Intersection(geom) to get the intersection point. It will return a new geometry, so just put it in a new feature like you thought.

An other option would be to use Shapely, shapely is more intuitive to use and claims to be more efficient than ogr python bindings. It comes with a built in topology validity check if is that what you want, check this section of the docs: http://toblerity.github.com/shapely/manual.html#spatial-analysis-methods

Related Question