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.