2

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)
2
  • Welcome to the site! It sounds like an interesting problem. Could you explain what isn't working about the script? What error messages are you getting (if any)? How does your output differ from what you expect? Mar 14, 2013 at 14:24
  • Hi! @JayGuarneri:) Thanks! The script is basically unfinished. I still have no lines for the loop searching for the intersections and dangles. Also I do not know how to write intersection and dangle results to a new geometry file, especially since org intersect, within and other geometry calls return binary outputs.
    – Spencer
    Mar 14, 2013 at 14:32

1 Answer 1

1

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

2
  • thanks @Pablo! I actually considered shapely but I was not that familiar. But I will try your suggestions. fingers crossed! hope it works :)
    – Spencer
    Mar 14, 2013 at 19:16
  • @Spencer, please post some feedback if you make any progress. And remember, StackExchange do encourage users to answer their own questions. (blog.stackoverflow.com/2011/07/…)
    – Pablo
    Mar 15, 2013 at 16:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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