Like in a title. Let's say I have normal polygon that overlaps other. And in my case, when the other polygon has holes or is splitted (like on attached image), normal methods "intersection" or "intersects" don't work.
One of my try was something like this, but it failed.
"some_feature" is that normal polygon.
# if polygon has holes, create new one without them
if len(feature.geometry().asPolygon()) > 1:
feat = QgsFeature()
geom = feature.geometry().asPolygon()
feat.setGeometry(QgsGeometry.fromPolygon([geom[0]]))
print(some_feature.geometry().intersection(feat.geometry()).area())
Below are some images, that show cases I'm talking about.
On the left is case with holes and on the right is one object as two polygons.
Does anyone have idea, how to achieve such thing?
Answer
I don't know why normal methods "intersection" or "intersects" don't work in your script for polygons with holes or multi parts but, in my script, they did. I used polygon layers of next image (where attributes table belongs to multipart layer):
My complete script is:
registry = QgsMapLayerRegistry.instance()
polygon3_test_holes = registry.mapLayersByName('polygon3_test_holes')
polygon7_test = registry.mapLayersByName('polygon7_test')
polygon_multipart = registry.mapLayersByName('polygon_multipart')
featp3_hole = polygon3_test_holes[0].getFeatures().next()
featp7_test = polygon7_test[0].getFeatures().next()
feat_multip = polygon_multipart[0].getFeatures().next()
epsg = polygon3_test_holes[0].crs().postgisSrid()
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
mem_layer = QgsVectorLayer(uri,
'polygon',
'memory')
prov = mem_layer.dataProvider()
if featp3_hole.geometry().intersects(featp7_test.geometry()):
geom = featp3_hole.geometry().intersection(featp7_test.geometry())
feat = QgsFeature()
feat.setAttributes([1])
feat.setGeometry(geom)
prov.addFeatures([feat])
if featp3_hole.geometry().intersects(feat_multip.geometry()):
geom = featp3_hole.geometry().intersection(feat_multip.geometry())
feat = QgsFeature()
feat.setAttributes([2])
feat.setGeometry(geom)
prov.addFeatures([feat])
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
After running the script at the Python Console of QGIS, I got memory layer (violet) of next image (where it was selected second feature):
Memory layer is a valid layer (it was corroborated directly at command line).
>>>layer = iface.activeLayer() #memory layer
>>>new_feat = layer.getFeatures().next()
>>>new_feat.geometry().isGeosValid()
True
Editing Note:
To validate feature geometry:
layer = iface.activeLayer()
feats = [ feat for feat in layer.getFeatures() ]
for feat in feats:
if feat.geometry().isGeosValid():
pass
else:
print "geometry is not valid for ID: ", feat.attribute("ID")
No comments:
Post a Comment