1

I have a problem with my code:

driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open('merge.shp', 0)
if dataSource is None:
    print 'Could not open file'
    sys.exit(1)
layer = dataSource.GetLayer() #layer = Shapefile

numFeatures = layer.GetFeatureCount()
i = 0
intersec_count=0
while i < (numFeatures-2):
    feature1 = layer.GetFeature(i)
    geom1 = feature1.GetGeometryRef()
    j = i+1
    while j < numFeatures:
        feature2 = layer.GetFeature(j)
        geom2 = feature2.GetGeometryRef()   
           if geom1.Distance(geom2)<0.0001 and geom1.Distance(geom2)!=0:
              buff = geom1.Buffer(0.0001)
              linkline = geom2.Intersection(buff)
              linkpoint = linkline.GetPoint(0)
              geom1.AddPoint(linkpoint[0],linkpoint[1])
              feature1.SetGeometry(geom1)
        j = j+1
    i = i+1
dataSource.SyncToDisk()
dataSource = None

When running this script the folloging Message appears: "Segmentation fault (core dumped)"

All I'm trying to do, is to add a Point (linkpoint) to a LineString (geom1) and save these changes.

I hope someone can help me :)

3
  • 1
    I do not understand what you're trying to do
    – gene
    Commented Nov 21, 2014 at 18:56
  • Segmentation fault is a memory error. How many records are in your shapefile?
    – andytilia
    Commented Nov 21, 2014 at 20:41
  • Any chance your geometries contain 4d(XYZM) coordinates? I've seen a lot of segfaults due to that. Commented Nov 21, 2014 at 21:19

1 Answer 1

5

For this, you need to create new ogr features and save them to the layer and not simply create geometries ( geom1.AddPoint(linkpoint[0],linkpoint[1]), feature1.SetGeometry(geom1)), with dataSource.SyncToDisk() at the end (SyncToDisk() might be helpful to ensure that a particular feature is flushed to disk but it is not necessary here, look at the script).

# open the shapefile specifying that the data can be updated
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
datasource = driver.Open('merge.shp', 1) # 1 -> the layer can be updated
layer = datasource.GetLayer()
featureDefn = layer.GetLayerDefn()
numFeatures = layer.GetFeatureCount()
i = 0
while i < (numFeatures-2):
    feature1 = layer.GetFeature(i)
    geom1 = feature1.GetGeometryRef()
    j = i+1
    while j < numFeatures:
         feature2 = layer.GetFeature(j)
         geom2 = feature2.GetGeometryRef()  
         if geom1.Distance(geom2)<0.0001 and geom1.Distance(geom2)!=0:
             buff = geom1.Buffer(0.0001)
             linkline = geom2.Intersection(buff)
             linkpoint = linkline.GetPoint(0)
             # create a new ogr geometry with wkt here
             wkt = 'POINT(%f %f)' % (linkpoint[0],linkpoint[1])
             point = ogr.CreateGeometryFromWkt(wkt)  
             # write the new feature  
             newFeature = ogr.Feature(featureDefn)      
             newFeature.SetGeometryDirectly(point)
             layer.CreateFeature(newFeature)
             # here you can place layer.SyncToDisk() if you want
             newFeature.Destroy()
         j = j+1
    i = i+1
datasource.Destroy() # or datasource = None

But instead of using nested while loops use a for loop or better, the itertools standard module.

import itertools
...
datasource = driver.Open('merge.shp', 1)
layer = datasource.GetLayer()
featureDefn = layer.GetLayerDefn()
for feature1,feature2 in  itertools.combinations(layer, 2):
     geom1 = feature1.GetGeometryRef()
     geom2 = feature2.GetGeometryRef()
     if geom1.Distance(geom2)<0.0001 and geom1.Distance(geom2)!=0:
         outFeature = ogr.Feature(featureDefn)
         buff = geom1.Buffer(0.0001)
         linkline = geom2.Intersection(buff)
         linkpoint = linkline.GetPoint(0)
         wkt = 'POINT(%f %f)' % (linkpoint[0],linkpoint[1])
         point = ogr.CreateGeometryFromWkt(wkt)         
         outFeature.SetGeometryDirectly(point)
         layer.CreateFeature(outFeature)
         outFeature.Destroy()

datasource.Destroy() 

Your Answer

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

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