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.

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)
share|improve this question
    
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? –  Jay Guarneri Mar 14 '13 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 '13 at 14:32
add comment

1 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

share|improve this answer
    
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 '13 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 '13 at 16:40
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.