0

I'm facing a problem with a process which should be straight forward, namely transforming a bunch of lines sharing their end nodes to a polygon.

The thing is, it just does not work. I tried a dozen of methods, went through validation and cleaning steps to make sure the data is OK, read dozens of online contributions, tried them all, and so on and so on, but there seems to be something intrinsically wrong with the data that makes it impossible to fulfill this simple task.

I could describe all the variations I went through, but I decided to share the shapefile for whoever is available to check what the problem with the data is.

You'll find it here

I would be very pleased to learn from you, what the clou is. Even some colleagues with many years of experience with QGIS and GIS in general were not able to find the problem.

ADDED 22.03.17

So, since I have been asked to be more specific with my problems, here it goes.

History I digitalized several features of satellite images using a 'line' layer. Since these features represent the border of management units and I want to transform the layer to a polygon layer, I made sure, that the last node of all lines coincide with the first node of the next line.

So, I end up with a layer with 213 individual lines (i.e.213 attributes), which looks like this (by the way, I realy made sure, that all end nodes coincide. original layer

I will call it START_LAYER from now on.

My goal: Transform these lines into a polygon.

So, that simple tasks turns out to be impossible, and I do not understand why.

So, what did I do?

Variant 1:

Vector > Geometry tools > Lines to polygons

The result looks like this bad result 1

I thought, Ok, sure, I have all individual lines, sure it will not work, so I went over to...

Variant 2:

... where I marked all lines (attributes) and merged them using the Merge selected features tool.

I end up with one attribute (I'll refer to it as LAYER_2) Go back to Var. 1 and create the polygon. The result looked like the first one, but the result was different (two different colors in the next picture). Any how, it was not what I was looking for. enter image description here

Ok, it has to work in a different way.

Variant 3:

I went to to do the Geometry Check (tried different parameters, the result might be different, but the following problematic seems to be the same) and tested it on LAYER_2.

Ok, 1) duplicate nodes, that makes sense, since each original line contributes with a node in the intersection point (just strange, having 213 initial lines, why only 51 double nodes?) and 2) Error = Geometry type (MultiLineString) - what ever that means (I still did not understand it).

OK, first get rid of the double nodes. Using the automated fix works for some, but not all.

enter image description here

Mmm, ok, so I'll just do it manually, which however is not really possible, since I just have one attribute.

But well, if I address then the second error "Geometry type (MultiLineString)" by running the proposed solution, several attributes are recreated, funny enough, more than I had initially. (apparently new double nodes appear which did not exist before)

Now I go and remove the double nodes by hand. (LAYER_3)

Go through Var. 1. and Var. 2. > nothing new.

Variant 4: Take all the variants attempted until now an play around with the function Merge Lines available in the menu point Vector. THe result, nothing new.

Note: By the way, I also cleaned my layer also via Check validation, but in the meantime it is clear, that the result has nothing to do with double nodes or similar stuff.

Variant 5: Go through the process using Polygonize. It doesn't matter what I do here (i.e. which way I go until I go for the last step to generate the polygon), the tool just shows enter image description here

I did try to recreate the process which lead me to the START_LAYER, i.e. creating lines with common end nodes and then I apply Polygonize and it works perfectly. Not with the START_LAYER.

Variant 6: Something drastic:

Extract Points

and then use the plugin

Points2One

to generate the polygon form the points.

The result looks like this. enter image description here

I try it only for one segment

enter image description here

I tried the plug-in for a data freshly created first lines, then node extraction, then polygon creation with the Points2One plug-in and it works fine.

It is as if some intrinsic value of the nodes is messing up every step I try. Some time-stamp, etc. I do not know, I really do not understand the system in depth.

Conclusion

I did try other stuff, which I will not explain now , I think the point is more than clear, probably I was redundant with some explanations.

What I really want to use is the Polygonize which normally works sooo well with fresh data.

So that is why I shared the data so you can try to understand what is wrong (or what I AM doing wrong) with it.


The full data set can be downloaded here:. This line shapefile should be free of double nodes or some other type of error according to the validation functions in QGIS.

  • 2
    What process are you following? What results are you getting? You have a gap of 21 metres at 674970 210391 which makes a leaky polygon. – Michael Stimson Mar 21 '17 at 0:52
  • 1
    Please describe what the problem is. How is the operation failing? What is your question? – Aaron Mar 21 '17 at 1:50
  • I voted +1 and for reopening this question because there are more than two detected problems (see answer an comment below). However, its geometry validation is OK. Interesting. – xunilk Mar 21 '17 at 12:14
  • 1
    A good question should show some degree of research and attempt. You have said you've tried different things, but as we don't know what you've tried we're likely to go through many of those same things as well. To help us help you, tell us what you've tried and what happens when you try them. This can help narrow the focus onto things you haven't tried, or perhaps that you've tried the wrong way. Also a question should have a question - a specific, tested problem that you need help with. Don't just give some data and tell people to have a look and determine the problem for themselves. – Midavalo Mar 21 '17 at 14:34
  • I tried to describe in the main description what I tried. It is still not everything, but this might be a hint. if you have questions, please let me know. – Patribus Mar 22 '17 at 22:26
3

Since QGIS is not able to solve polygonizing in an easy way, I use the Openstreetmap editor JOSM for such things. In QGIS, convert your data into a GPX file setting the FORCE_GPX_TRACKoption to YES.

In JOSM, load the GPX file, and convert it to a data layer. Then select the lines that should build a polygon, and create a new relation of type multipolygon from it. You might add landuse=forest to get a reasonable tagging.

In the relation editor, you can order the elements. In most cases, they will not yet chain. You have to merge the end points manually, then re-sort.

In some cases, you have to split lines where three polygons share one point, but one line continues. At least two points were too close to each other, so I had to merge them too. And there were still some gaps, and dangling line ends.

If the lines close, you get a loop in the relation editor. You might run the check tool of JOSM as well, before saving to a local file. DO NOT UPLOAD ON EXITING JOSM.

Then you can load the .osm file back into QGIS with Add Vector Layer or the QuickOSM plugin.

  • Thanks, I will have a look and report about the results. – Patribus Mar 24 '17 at 14:44
  • Dear AndreJ, I tested your approach and it worked, thanks! It is an effective way of getting the results I needed, but it is far from being efficient when you have that many lines and nodes to edit. But well, I got my result. As soon as I find the time, I will also check the other methods. – Patribus Mar 28 '17 at 12:55
  • You are right about the efficiency. I tend to digitize polygons as polygons, not lines. You have to catch every line twice, but for the second run I turn off the background map, zoom a bit out, and just snap to the nodes of the previous polygon. – AndreJ Mar 28 '17 at 15:59
5

The Polyline with FID=4 has 3 parts. There is a 22 meter gap between part 1 and part 2. This is the only gap in the data.

enter image description here

  • 1
    I found another miniscule gap (less than a millimetre) at 680651.71547115,220922.242444406). Depending on the tolerance of the tool, it may fail on this one too. – Fezter Mar 21 '17 at 1:29
  • dear @klewis, with which software are you looking at the data? – Patribus Mar 22 '17 at 20:32
  • @Patribus, I'm using ArcGIS ArcMap 10.5. The window with coordinate values is the Edit Sketch Properties Window. – klewis Mar 22 '17 at 22:48
  • Ok, thanks, i just installed ArcGis two days ago to test the transformation from lines to polygon.Even ArcGis has problems and the result shows some problems I could share with you if you want. I will try to use ArcGis to search for further problems with the data – Patribus Mar 23 '17 at 9:36
  • Dear @Fezter, I did upload new data which should be clean from errors. The transformation to polygon still does not work for me. Maybe you could a have a look into it now. (The file is at the end of the post) – Patribus Mar 23 '17 at 13:59
2

I downloaded test shapefile and it originally looks like:

enter image description here

and after using "Lines to polygon" method of Processing tool box of QGIS looks as:

enter image description here

Obviously, there are problems with it.

On the other hand, when two consecutive features are selected at attributes table they are no showed in this way at the Map Canvas.

enter image description here

As @klewis pointed out that "There is a 22 meter gap between part 1 and part 2", I used next PyQGIS code to evaluate how many gaps are in the geometry.

import itertools

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

n = len(feats)

for i, j in itertools.combinations(range(n), 2):
    distance = feats[i].geometry().distance(feats[j].geometry())
    if distance < 50:
        print i, j, distance

After running it at the Python Console of QGIS, I got:

0 3 0.0
0 5 6.32561007326e-08
1 2 0.0
1 4 2.77618203084e-08
2 3 0.0
4 6 3.69805894613e-10
5 6 5.00951291099e-08

suggesting only minuscules gaps and not 22 meters gap. So, there are features that are intrinsically multiparts.

By using "Multipart to singlepart" method of Processing tool box of QGIS this was corroborated and above code was run again to produce next result:

0 3 0.0
0 7 6.32561007326e-08
1 2 0.0
1 6 2.77618203084e-08
2 3 0.0
4 5 6.09094222774e-08
4 8 3.69805894613e-10
5 6 22.0495367734
7 9 5.00951291099e-08
8 9 9.68909528121e-08

Now, it can be observed 22 meters gap.

By using "v.clean" method of Processing tool box of QGIS with snap option and threshold 30 m, I got a layer without gaps; as it can be observed at a new execution of above code:

0 3 0.0
0 7 0.0
1 2 0.0
1 6 0.0
2 3 0.0
4 5 0.0
4 8 0.0
5 6 0.0
7 9 0.0
8 9 0.0

However, gaps are not the cause of problems. So, I used next code:

import itertools

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

selected_feats = []

selected_feats.append(feats[0])

del feats[0]

k=0

while len(feats) != 0:
    for i, feat in enumerate(feats):
        d = selected_feats[k].geometry().distance(feat.geometry())
        if d == 0:
            selected_feats.append(feat)
            k += 1
            del feats[i]

lines = [ feat.geometry().asPolyline() for feat in selected_feats ]

epsg = layer.crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           '4testing_feats_ordered',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(selected_feats)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(lines[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

for ordering features in a new memory layer; as can see at next image:

enter image description here

Finally, next code produces a polygon memory layer (points clockwise) by using '4testing_feats_ordered' layer:

registry = QgsMapLayerRegistry.instance()

layer = registry.mapLayersByName("4testing_feats_ordered") 

lines = [ feat.geometry().asPolyline() for feat in layer[0].getFeatures() ]

sum = 0

for i in range(len(lines[0])-2):
    det =  (lines[0][i+1][0] - lines[0][i][0])*(lines[0][i+2][1] - lines[0][i][1])-(lines[0][i+1][1] - lines[0][i][1])*(lines[0][i+2][0] - lines[0][i][0])
    sum += det

if sum > 0:

    lines[0].reverse()

new_points = []

for point in lines[0]:
    new_points.append(point)

for i in range(1, len(lines), 1):
    if lines[i][-1] == new_points[-1]:
        lines[i].reverse()

    for point in lines[i]:
        new_points.append(point)

polygon = [[new_points]]

epsg = layer[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'malasya_polygon',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(polygon)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolygon(polygon[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

as can be observed at next image (desired polygon):

enter image description here

  • Uau, thanks! Now my next question, since I'm more or less new to QGis. Where can I run the code? An since the ordering of the features seems to be the problem in many of the reported procedures, is there another way of ordering the features that does not need coding? Some function, where you can indicate which criteria to use to reorder the features? – Patribus Mar 24 '17 at 14:46
  • You can run the code at Python Console of QGIS. However, you need to do first gap elimination with v.clean. Criteria was established on the basis of watching point ordering in each feature. They come from Calculus (when you have an inflection point in a function, point ordering changes from clockwise to counterclockwise or vice versa). Procedure is "equivalent" to exploratory data analysis in statistic. By the way, I practically know procedures with coding. – xunilk Mar 25 '17 at 8:33
1

Unfortunatly i couldn't get QGIS to do this. However ArcGIS can do it, here is the easiest way i found...

  1. Open your line file, i added a couple extra lines to your dataset to make it more interesting and more like your own scenario :-) enter image description here

  2. Dissolve all your lines to one big line.

edit: forgot to mention to close the gaps, i did this in QGIS using the topology checker and setting it to "find dangles".

  1. Create a new shapefile with a polygon that encloses all your lines, you can do this by using the draw toolbar and then using the convert graphics to feature in the drawing toolbar drop down menu:enter image description hereenter image description here
  2. Select the dissolved line
  3. Start an edit session on the polygon layer. Select the polygon layer, do not select the polygon it self...
  4. Use the advanced editing toolbar, you can find this by right-clicking on a blank spot on the toolbar area in the top and selecting it here.
  5. split the polygon with the lines :-)enter image description here
  6. DONE! you now have a polygon for each area that was spanned between the original lines. you can remove the excess outside polygon if you want to...enter image description here
  • Ok, I will try it and report to you if it worked. In ArcGis I tried the "Feature to Polygon" function, and it does work to some extend, but than looking at it in detail, the results if far from being correct. I will report on this below. – Patribus Mar 24 '17 at 14:43
0

OK, I don't have QGIS experience, so I am using other tools to look at your data. You say the end nodes now match, and that's true... but not totally true.

For example, you have one line who's start point is at:

Coordinate: 112.55618528638789, 1.996431306892094

While the line that should connect starts with:

Coordinate: 112.55618528638789, 1.9964313068920885

OK, there's only a difference of 0.0000000000000055 in latitude, but that's enough on my computer to not make a match. It depends on how picky QGIS is.

Here's what I did using FME.

enter image description here

Basically I intersected the lines (because some lines pass through a boundary) and picked out the example you showed above. All I did was round the coordinates to 4 decimal places. That removed the difference in coordinates. Then I can build that area feature.

So if you can round off the coordinates (maybe not to 4, but to a few less decimal places) I figure it will be much improved. Maybe QGIS has a precision setting that can be used???

Also, I said before that there might be an ordering problem, and there could still be depending on how QGIS tackles the problem.

For example, take the coordinates above. Both of them are the start point of a line feature. If it were digitized correctly I would expect one of them to be a start point, the other to be an end point.

So one is going in a clockwise direction, the other anti-clockwise. FME doesn't have an issue with that. QGIS may or may not, I don't know. If it does, I'd suggest trying to reverse the order so that they all travel in a clockwise direction.

I hope this helps.

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.