1

I have a vector layer with polygons representing administrative districts of germany. For every polygon smaller than 500 km² I would like to create bigger polygons by merging each of the small polygons with all its adjacent polygons.

You can find the dataset here under vg250_01-01.utm32s.shape.ebenen\vg250_ebenen_0101\VG250_KRS.shp

In this picture you can see the dataset with my selection of the polygons smaller than 500 km² by the processing tool "selectbyexpression" with the expression "$area < 500000000".

enter image description here

I tried the processing tool "eliminateselectedpolygons", however, it is only able to merge with either the biggest or the smallest adjacent polygon or the one with the longest shared border. Unfortunately, I wasn't able to find a tool that does the same with all adjacent polygons. I also tried to find a solution with "selectbylocation" and "dissolve", which works fine for the isolated small polygons, but I wasn'a able to make it work with the clusters of small polygons, which will end up in one big polygon.

How can I realize this preferably with QGIS? If this isn't possible with QGIS I'm also open for suggestions with other tools, e.g. R.

Improve this question
5
  • For example Schwerin is adjacent to Nordwestmecklenburg and Ludwigslust-Parchim. By your description you want to put all them together into one geometry. Is that right?
    – user30184
    Commented Feb 14, 2023 at 14:07
  • I'm not sure what you mean with "I would like to create bigger polygons by merging each of the small polygons with all its adjacent polygons." You want to blend them into eachother? Like one small polygon has two bigger neighbours on the left and right. Then both neighbours must have that small polygon inside it? Then you do realize the total sum of the area of Germany becomes much bigger than the actual size of Germany? If that is what you want I think you need a script which loops through each small polygon looking for surrounded ones. You might find your solution in a Python script.
    – ThomasHaarlem
    Commented Feb 14, 2023 at 15:33
  • @user30184 Exactly. So for your example with Schwerin I want to find all adjacent features (Nordwestmecklenburg and Ludwigslust-Parchim) and then put all of them in one geometry. The result should be one single feature with the extend of the three combined and the attributes of Schwerin. And then do the same process for every single feature smaller than 500 km².
    – zfn
    Commented Feb 14, 2023 at 16:48
  • @ThomasHaarlem Sorry. I had a lot of difficulties to explain it. In the end, the aim is to have a new shapefile with the same administrative districts, but each feature smaller 500 km² (for example A) is replaced by a union of itself (A) and its adjacent features (B+C) while keeping its original attributes (A_attributes). But only feature A is replaced. B+C remain in the shapefile. A loop sounds like a reasonable plan. Could you give me a hint how to realize this special task in a phython script?
    – zfn
    Commented Feb 14, 2023 at 17:07
  • alright, that clears it up. You just want to increase the geometry of A and keep features B and C intact (unless those are also smaller then 500 km2, then those will also unionize the geometry with A)? You do realize that A will overlap B and C? And that is no problem? I think you need Python, with Geopandas-library and maybe Shapely as well. In your Python-script you need a For-loop with if-statement construction. Unfortunatly I can not help you with that. At least not in the short-term.
    – ThomasHaarlem
    Commented Feb 16, 2023 at 10:32

1 Answer 1

Reset to default
1

I asked the question some time ago, but just recently found the time to work on that problem again.

I wrote a Python-script for QGIS that worked fine for me. I just put it here, if anybody has a similar problem and needs some inspiration for a script. The script takes the active layer in QGIS as input and adds the resulting layer to the project.

from qgis import processing

layer = iface.activeLayer()
provider = layer.dataProvider()
layer.selectByExpression("$area<500000000")
selection = layer.selectedFeatures()

vl = QgsVectorLayer("MultiPolygon?crs=epsg:25832", "temporary_polygons", "memory")
pr = vl.dataProvider()
pr.addAttributes(provider.fields())
vl.updateFields()

for feature in selection:
    layer.removeSelection()
    layer.select(feature.id())
    processing.run('native:selectbylocation', {
                'INPUT': layer,
                'PREDICATE': [4],
                'INTERSECT': QgsProcessingFeatureSourceDefinition(layer.id(), selectedFeaturesOnly=True, featureLimit=-1, geometryCheck=QgsFeatureRequest.GeometryAbortOnInvalid),
                'METHOD': 1})
    
    dis = processing.run('native:dissolve', {
                'INPUT': QgsProcessingFeatureSourceDefinition(layer.id(), selectedFeaturesOnly=True, featureLimit=-1, geometryCheck=QgsFeatureRequest.GeometryAbortOnInvalid),
                'OUTPUT': 'memory:'})
    
    fet = next(dis['OUTPUT'].getFeatures())
    fet.setAttributes(feature.attributes())
    pr.addFeatures([fet])

vl.updateExtents()
QgsProject.instance().addMapLayer(vl)
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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