3

I am trying to combine polygons, using Qgis 1.8, and every time I run this script in the python console, I get a segmentation fault. Here is the script:

from PyQt4.QtCore import *

iface = qgis.utils.iface;
legend = iface.legendInterface()
poly_layer = iface.activeLayer()

MapLayer = QgsMapLayerRegistry.instance()

#######NEW POLY LAYER TO CREATE#########
poly_provider = poly_layer.dataProvider()

fields = { 0 : QgsField("STAT_NO", QVariant.String) }

writer = QgsVectorFileWriter("/localstore/data/shps/contrib_area.shp",     "CP1250", fields, poly_provider.geometryType(), poly_layer.crs(), "ESRI Shapefile")

if writer.hasError() != QgsVectorFileWriter.NoError:
  print "Error when creating shapefile: ", writer.hasError()



#######POINT PROCESSING########
#get the point layer that contains the points contained by the polygons
point_layer = None

for layer_id, layer in MapLayer.mapLayers().iteritems():
  if(layer.name() == "slope_and_catchment"):
    point_layer = MapLayer.mapLayer(layer_id)
    break

point_provider = point_layer.dataProvider()
point_id_index = point_provider.fieldNameIndex("SITE_NO")
point_link_index = point_provider.fieldNameIndex("pw_comb")
point_provider.select([point_id_index,point_link_index])
point_feat = QgsFeature()

intersects = []

point_feat_dict = {}
while point_provider.nextFeature(point_feat):
  point_attrmap = point_feat.attributeMap()
  point = point_attrmap[point_link_index].toString()
  intersects.append(point)
  point_feat_dict[point] = point_attrmap[point_id_index].toString()

for point in intersects:
  print point


########LINE PROCESSING######
#get the line layer that contains the stream network
line_layer = None

for layer_id, layer in MapLayer.mapLayers().iteritems():
  if(layer.name() == "projectednet"):
    line_layer = MapLayer.mapLayer(layer_id)
    break

line_provider = line_layer.dataProvider()
line_index = line_provider.fieldNameIndex("LINKNO")
line_US1_index = line_provider.fieldNameIndex("USLINKNO1")
line_US2_index = line_provider.fieldNameIndex("USLINKNO2")

line_provider.select([line_index,line_US1_index,line_US2_index])
line_feat = QgsFeature()


line_feat_dict = {}
#loop over all features and store relevant information
while line_provider.nextFeature(line_feat):
  line_feat_id = line_feat.id()
  line_attrmap = line_feat.attributeMap()
  link = line_attrmap[line_index].toString()
  us1 = line_attrmap[line_US1_index].toString()
  us2 = line_attrmap[line_US2_index].toString()
  line_feat_dict[link] = [line_feat_id, link, us1, us2] 

#all_lines = line_feat_dict.keys()

#line_feat_dict is now a map from link names to feature ids, and the ustream links


########POLYGON PROCESSING#######
#set up attribute information
poly_id_index = poly_provider.fieldNameIndex("DN")
#get all features allong with the attributes
poly_provider.select([poly_id_index])
#make a feature variable to store features in
poly_feat = QgsFeature()

#######NOW WE NEED TO GET ALL CATCHMENTS THAT SHARE LINKS OUTLINED IN STREAM NETWORK########
poly_feat_dict = {}
#get a list of all related polygons:
while poly_provider.nextFeature(poly_feat):
  poly_id = poly_feat.id()
  poly_attrmap = poly_feat.attributeMap()
  poly = poly_attrmap[poly_id_index].toString()
  poly_feat_dict[poly] = poly_id

#poly_feat_dict now maps the link number to the feature id

#start traversing the network and getting geometries
geoms = []

def traverse(link1, link2):
  #work down one side of the "tree"
  if link1 != "-1":
    feat = QgsFeature()
    poly_provider.featureAtId( poly_feat_dict[link1], feat )
    geoms.append( feat.geometry() )
    #traverse this links uplinks
    traverse( line_feat_dict[link1][2], line_feat_dict[link1][3] )
  else:
    #have hit headwater
    return
  if link2 != "-1":
    feat = QgsFeature()
    poly_provider.featureAtId( poly_feat_dict[link2], feat )
    geoms.append( feat.geometry() )
    #traverse this links uplinks
    traverse( line_feat_dict[link2][2], line_feat_dict[link2][3] )
  else:
    #have hit headwater
    return


#process the catchment up to station  
for point in intersects:
  #make sure the geoms list is clear for each point
  geoms = []
  #point is the starting link number
  #get this initial geometry from poly_feat_dict
  feat = QgsFeature()
  poly_provider.featureAtId( poly_feat_dict[point], feat )
  geoms.append( feat.geometry() )

  #get the corresponding upstream links
  up1 = line_feat_dict[point][2]
  up2 = line_feat_dict[point][3]
  #traverse the stream network links till we can get all geometries
  traverse(up2, up2)
  #have all geometries, join them
  union_geoms = geoms[0]
  for geometry in geoms:
    union_geoms = union_geoms.combine(geometry)

  fet = QgsFeature()
  fet.setGeometry(union_geoms)
  fet.addAttribute(0, QVariant( point_feat_dict[point] ))
  writer.addFeature(fet)


del writer

layer = QgsVectorLayer("/localstore/data/shps/contrib_area.shp", "contrib_area", "ogr")

MapLayer.addMapLayer(layer)

if iface.setActiveLayer(layer):
  print "New Layer is Active"
else:
  print "Couldn't set active layer"

The following is the GDB backtrace from the coredump:

Thread 1 (Thread 0x7ffff2f27820 (LWP 29411)):
#0  0x00007ffff3645cb0 in geos::geom::Geometry::getPrecisionModel() const () from /usr/lib64/libgeos-3.3.2.so
No symbol table info available.
#1  0x00007ffff369f876 in   geos::operation::GeometryGraphOperation::GeometryGraphOperation(geos::geom::Geometry const*, geos::geom::Geometry const*) () from /usr/lib64/libgeos-3.3.2.so
No symbol table info available.
#2  0x00007ffff36baea0 in geos::operation::overlay::OverlayOp::OverlayOp(geos::geom::Geometry const*, geos::geom::Geometry const*) () from /usr/lib64/libgeos-3.3.2.so
No symbol table info available.
#3  0x00007ffff36bb901 in geos::operation::overlay::OverlayOp::overlayOp(geos::geom::Geometry const*, geos::geom::Geometry const*, geos::operation::overlay::OverlayOp::OpCode) () from /usr/lib64/libgeos-3.3.2.so
No symbol table info available.
#4  0x00007ffff46165da in std::auto_ptr<geos::geom::Geometry> geos::geom::BinaryOp<geos::operation::overlay::overlayOp>(geos::geom::Geometry const*, geos::geom::Geometry const*, geos::operation::overlay::overlayOp) () from /usr/lib64/libgeos_c.so.1
No symbol table info available.
#5  0x00007ffff46109ae in GEOSUnion_r () from /usr/lib64/libgeos_c.so.1
No symbol table info available.
#6  0x00007ffff56b1294 in QgsGeometry::combine(QgsGeometry*) () from /opt/lib64/libqgis_core.so.1.8.0
No symbol table info available.
#7  0x00007fffe346185e in meth_QgsGeometry_combine () from /opt/share/qgis/python/qgis/core.so
No symbol table info available.
#8  0x00000034c42d59e4 in PyEval_EvalFrameEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#9  0x00000034c42d7657 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#10 0x00000034c42d7732 in PyEval_EvalCode () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#11 0x00000034c42d5c92 in PyEval_EvalFrameEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#12 0x00000034c42d7657 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#13 0x00000034c42d5aa4 in PyEval_EvalFrameEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#14 0x00000034c42d7657 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#15 0x00000034c42d5aa4 in PyEval_EvalFrameEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#16 0x00000034c42d6b8f in PyEval_EvalFrameEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#17 0x00000034c42d7657 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#18 0x00000034c426acb0 in ?? () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#19 0x00000034c4243c63 in PyObject_Call () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#20 0x00000034c42566af in ?? () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#21 0x00000034c4243c63 in PyObject_Call () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#22 0x00000034c42cfc93 in PyEval_CallObjectWithKeywords () from /usr/lib64/libpython2.6.so.1.0
No symbol table info available.
#23 0x00007fffe45bf3db in ?? () from /usr/lib64/python2.6/site-packages/sip.so
No symbol table info available.
#24 0x00007fffe39927c3 in ?? () from /usr/lib64/python2.6/site-packages/PyQt4/QtGui.so
No symbol table info available.

My question is am I doing something wrong, or is this a bug? Any thoughts would be appreciated, Thanks!

EDIT: More stack trace is available if needed. I only included the last 25 steps.

Improve this question

1 Answer 1

Reset to default
3

Seems like this is a "feature"...http://hub.qgis.org/issues/777

When the provider loops through the features, python stores the geometry, but as soon as the loop continues to the next iteration, the feature gets destroyed, and with it the geometry gets destroyed--passing invalid memory when trying to combine the geometries. So to save the geometries to combine later you need to use

geoms.append( feat.geometryAndOwnership() )

instead of

geoms.append( feat.geometry() )

when looping over the features.

Improve this answer
2
  • 1
    You can also do geoms.append( QgsFeature(feat.geometry()) ) which will make a copy of the geometry.
    – Nathan W
    May 21, 2013 at 23:04
  • Ah, good to know. I looked through the documentation on the Geometry class and was wondering how to use a copy constructor, didn't think to use the Feature constructor. Curious, do you know which is more efficient? (This has to run on over 500,000 polys) My initial guess is that getting the already created geometry would be faster than recreating the feature and capturing its geometry...
    – Nelz11
    May 22, 2013 at 2:39

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.