5

I am trying to use QGIS to split a set of polygons into single parts where there is a join at a single vertex. The image shows a number of examples.

I have taken a single polygon and used 'Singlepart to multipart' to generate the example. This does not however create a single part where there are touching / shared vertices. Is there a way to further split these polygons down where only a single vertex is shared in common?

I do not want a manual process: any solution must be be able to be used as part of an automated PyQGIS script.

Example joined polygons

Test geometry WKT:

Polygon ((393796.08000000001629815 389603.25, 393784.08000000001629815 389603.25, 393784.08000000001629815 389615.25, 393780.08000000001629815 389615.25, 393780.08000000001629815 389611.25, 393776.08000000001629815 389611.25, 393776.08000000001629815 389615.25, 393780.08000000001629815 389615.25, 393780.08000000001629815 389619.25, 393784.08000000001629815 389619.25, 393784.08000000001629815 389623.25, 393792.08000000001629815 389623.25, 393792.08000000001629815 389635.25, 393796.08000000001629815 389635.25, 393796.08000000001629815 389631.25, 393800.08000000001629815 389631.25, 393800.08000000001629815 389627.25, 393804.08000000001629815 389627.25, 393804.08000000001629815 389623.25, 393808.08000000001629815 389623.25, 393808.08000000001629815 389627.25, 393812.08000000001629815 389627.25, 393812.08000000001629815 389619.25, 393820.08000000001629815 389619.25, 393820.08000000001629815 389611.25, 393824.08000000001629815 389611.25, 393824.08000000001629815 389603.25, 393800.08000000001629815 389603.25, 393800.08000000001629815 389607.25, 393796.08000000001629815 389607.25, 393796.08000000001629815 389603.25))

  • Could you include a test geometry as WKT? – user30184 Feb 7 '17 at 14:47
  • I've edited the question to include the highlighted test example as WKT. – AJ Walsh Feb 7 '17 at 15:18
  • To be sure, do you want to get two simple polygons, this small one POLYGON (( 393780.08 389615.25, 393780.08 389611.25, 393776.08 389611.25, 393776.08 389615.25, 393780.08 389615.25 )) and the rest as another polygon? – user30184 Feb 7 '17 at 15:24
  • Yes exactly that. In this example, two polygons, one small one as per your WKT and the rest as another polygon. – AJ Walsh Feb 7 '17 at 15:30
  • This is interesting. QGIS thinks that this polygon is topologically correct and therefore it is not so easy to correct. I will need to think. – user30184 Feb 7 '17 at 15:44
3

Short Answer-

For topology correction/cleaning you can use the GRASS 'v.clean' algorithm from the processing toolbox. In your case break option of the cleaning tool would be perfect for splitting polygons.

For information on how to access the processing algorithms in PyQGIS refer to - http://docs.qgis.org/2.14/en/docs/user_manual/processing/console.html#calling-algorithms-from-the-python-console


Detailed answer -

The syntax definition of v.clean -

ALGORITHM: v.clean - Toolset for cleaning topology of vector map.
input <ParameterVector>
tool <ParameterSelection>
thresh <ParameterNumber>
GRASS_REGION_PARAMETER <ParameterExtent>
GRASS_SNAP_TOLERANCE_PARAMETER <ParameterNumber>
GRASS_MIN_AREA_PARAMETER <ParameterNumber>
output <OutputVector>
error <OutputVector>
tool(Cleaning tool)
    0 - break
    1 - snap
    2 - rmdangle
    3 - chdangle
    4 - rmbridge
    5 - chbridge
    6 - rmdupl
    7 - rmdac
    8 - bpol
    9 - prune
    10 - rmarea
    11 - rmline
    12 - rmsa

For your purpose, you'll use 'break' option in the cleaning tool. The thresh, GRASS_SNAP_TOLERANCE_PARAMETER and GRASS_MIN_AREA_PARAMETER can be 0. The QGIS processing algorithms expect the extent parameter in a different format which can be obtained from QgsRectangle object -

rect = layer.extent();            #Where 'layer' is the input QgsVectorLayer object
extent = "%f , %f , %f , %f "% (rect.xMinimum(), rect.xMaximum() , rect.yMinimum(), rect.yMaximum() )

Now you can run the algorithm in PyQGIS by importing processing and using runalg or runandload -

import processing
processing.runalg("grass7:v.clean", layer , 0 , 0 , extent , 0 , 0 , <OutputVector> , None );

I followed the above steps for the WKT provided in the question on QGIS 2.18 with GRASS 7 and it worked well.

Hope this helps!

  • Many thanks for the answer and the detailed description provided, this was really useful. The steps you suggest work correctly. – AJ Walsh Feb 9 '17 at 15:28
3

First, it feels like QGIS has a bug because it does not consider your polygon as invalid. I would consider writing a question to qgis-users mailing list and based on the feedback perhaps a bug report.

I don't know how to correct your polygons with QGIS script if program does not even think there is anything to repair. One alternative would be to use a spatial database and some SQL. First thing to do would be to run ST_MakeValid which will convert self-intersecting polygons into multipolygons. How to explode multipolygons depends on the database: PostGIS has "ST_Dump" and SpatiaLite has "ElementaryGeometries".

  1. PostGIS

SELECT ST_Dump( ST_MakeValid( ST_GeomFromText(' Polygon ((393796.08000000001629815 389603.25, 393784.08000000001629815 389603.25, 393784.08000000001629815 389615.25, 393780.08000000001629815 389615.25, 393780.08000000001629815 389611.25, 393776.08000000001629815 389611.25, 393776.08000000001629815 389615.25, 393780.08000000001629815 389615.25, 393780.08000000001629815 389619.25, 393784.08000000001629815 389619.25, 393784.08000000001629815 389623.25, 393792.08000000001629815 389623.25, 393792.08000000001629815 389635.25, 393796.08000000001629815 389635.25, 393796.08000000001629815 389631.25, 393800.08000000001629815 389631.25, 393800.08000000001629815 389627.25, 393804.08000000001629815 389627.25, 393804.08000000001629815 389623.25, 393808.08000000001629815 389623.25, 393808.08000000001629815 389627.25, 393812.08000000001629815 389627.25, 393812.08000000001629815 389619.25, 393820.08000000001629815 389619.25, 393820.08000000001629815 389611.25, 393824.08000000001629815 389611.25, 393824.08000000001629815 389603.25, 393800.08000000001629815 389603.25, 393800.08000000001629815 389607.25, 393796.08000000001629815 389607.25, 393796.08000000001629815 389603.25))') ));

  1. SpatiaLite

Step one:

SELECT ST_MakeValid( ST_GeomFromText(' Polygon ((393796.08000000001629815 389603.25, 393784.08000000001629815 389603.25, 393784.08000000001629815 389615.25, 393780.08000000001629815 389615.25, 393780.08000000001629815 389611.25, 393776.08000000001629815 389611.25, 393776.08000000001629815 389615.25, 393780.08000000001629815 389615.25, 393780.08000000001629815 389619.25, 393784.08000000001629815 389619.25, 393784.08000000001629815 389623.25, 393792.08000000001629815 389623.25, 393792.08000000001629815 389635.25, 393796.08000000001629815 389635.25, 393796.08000000001629815 389631.25, 393800.08000000001629815 389631.25, 393800.08000000001629815 389627.25, 393804.08000000001629815 389627.25, 393804.08000000001629815 389623.25, 393808.08000000001629815 389623.25, 393808.08000000001629815 389627.25, 393812.08000000001629815 389627.25, 393812.08000000001629815 389619.25, 393820.08000000001629815 389619.25, 393820.08000000001629815 389611.25, 393824.08000000001629815 389611.25, 393824.08000000001629815 389603.25, 393800.08000000001629815 389603.25, 393800.08000000001629815 389607.25, 393796.08000000001629815 389607.25, 393796.08000000001629815 389603.25))') );

Step 2:

Collect all data from the original data, with corrected polygons, into a table and continue with ElementaryGeometries function https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html.

  • Thanks for taking the time to look into this problem. The answer provided by nash sets out an approach to use within QGIS. Nevertheless, the approach you have used is is a good alternative. – AJ Walsh Feb 9 '17 at 15:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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