Wednesday, 2 March 2016

pyqgis - Get intersection in polygons with holes or multi parts


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):


enter image description here


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):


enter image description here


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...