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.

  • 1
    Check your input file links_shapefile = "". – swatchai Jul 10 '17 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 Jul 10 '17 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 Jul 10 '17 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 Jul 10 '17 at 17:57
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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