0

I am trying to select lines from a shapefile using an IN statement, then dissolve these lines into one feature, and then export.

A simplified version of my script is here. It is not producing any output in the output shapefile, but I don't know why!

from osgeo import ogr
from shapely.wkb import loads
from shapely.ops import cascaded_union

# Input line shapefile
links_shapefile = workspace + r"/Input/Links.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
links_ds = driver.Open(links_shapefile)
links_layer = links_ds.GetLayer(0)

# Output shapefile
output_file = workspace + r"/Output/Output_test.shp"
output_ds = ogr.Open(output_file)
output_layer = output_ds.GetLayer(0)

# Query input shapefile to a list of link IDs
link_attributes = ['1001_10002', '1002_1003', '1003_1004']
file_name = links_shapefile[:-4].split('/')[-1]
sql = "SELECT A_B FROM " + file_name + " WHERE A_B IN " + str(tuple(link_attributes))
selection_layer = links_ds.ExecuteSQL(sql)

# Loop selected lines and append geometry to a list
selection_feature_count = selection_layer.GetFeatureCount()
lines = [] #list to store line geometry
for index in range(0, selection_feature_count):
      feature = selection_layer.GetFeature(index)
      geometry = feature.geometry()
      lines.append(loads(geometry.ExportToWkb()))

#Dissolve list of lines
u = cascaded_union(lines)

#Export dissolved lines to output shapefile
defn = output_layer.GetLayerDefn()
feat = ogr.Feature(defn)
feat.SetField('UniqID', 9999)
geom = ogr.CreateGeometryFromWkb(u.wkb)
feat.SetGeometry(geom)
output_layer.CreateFeature(feat)

Please let me know if you see any issues with this script or if there is a better way to do something! I am new to using Python.

Improve this question
4
  • 1
    Check your input file links_shapefile = "".
    – swatchai
    Commented Jul 10, 2017 at 15:32
  • Thanks @swatchai I have it correct in the code I was running but simplified it for this post. Have added a path in now so as not to confuse!
    – Andy
    Commented Jul 10, 2017 at 15:41
  • 2
    With Shapely, it is easier to use Fiona in place of ogr. For LineStrings, use Shapely unary_union ( The unary union function can operate on different geometry types, not only polygons as is the case for the older cascaded unions)
    – gene
    Commented Jul 10, 2017 at 15:49
  • Thanks, I just gave this a go (using Fiona instead of OGR) but it seems a lot slower and I have 100,000's records to loop!
    – Andy
    Commented Jul 10, 2017 at 17:57

1 Answer 1

Reset to default
1

The issue was I was using GetFeature, rather than GetNextFeature.

This looped through the features in the layer, rather than features in queried/selected layer.

Improve this answer

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.