I have Windows 7 running 3.10.2-A Coruña. The only success I have had in creating nested polygons from a raster dataset is using GDAL with -p to create polygons. However, I recently cancelled the process after it had been running for more than a day, as it seems to struggle with large datasets (around 50000 contour lines).
Ideally, I should be able to use lines to polygons (QGIS or SAGA), but the nested polygons are not shown. I have tried Dissolve on Z/ELEV, but only with partial success. v.clean in GRASS consistently fails, regardless of parameters chosen. I have tried checking geometry and topology, again with no success. What is frustrating is that contours created with GDAl cannot be converted to polygons using a separate process rather than the -p option.
Does anyone have a process that may work as well as GDAL -p, but a lot quicker/more efficiently?
[GIS] QGIS alternatives to GDAL contours to polygons
gdal-contourgrasslinepolygonqgis
Related Solutions
To illustrate a raster/image processing solution, I began with the posted image. It is of much lower quality than the original data, due to the superposition of blue dots, gray lines, colored regions, and text; and the thickening of the original red lines. As such it presents a challenge: nevertheless, we can still obtain Voronoi cells with high accuracy.
I extracted the visible parts of the red linear features by subtracting the green from the red channel and then dilating and eroding the brightest parts by three pixels. This was used as the base for a Euclidean distance calculation:
i = Import["http://i.stack.imgur.com/y8xlS.png"];
{r, g, b} = ColorSeparate[i];
string = With[{n = 3}, Erosion[Dilation[Binarize[ImageSubtract[r, g]], n], n]];
ReliefPlot[Reverse@ImageData@DistanceTransform[ColorNegate[string]]]
(All code shown here is Mathematica 8.)
Identifying the evident "ridges"--which must include all points that separate two adjacent Voronoi cells--and re-combining them with the line layer provides most of what we need to proceed:
ridges = Binarize[ColorNegate[
LaplacianGaussianFilter[DistanceTransform[ColorNegate[string]], 2] // ImageAdjust], .65];
ColorCombine[{ridges, string}]
The red band represents what I could save of the line and the cyan band shows the ridges in the distance transform. (There is still a lot of junk due to the breaks in the original line itself.) These ridges need to be cleaned and closed up through a further dilation--two pixels will do--and then we can identify the connected regions determined by the original lines and the ridges between them (some of which need explicitly to be recombined):
Dilation[MorphologicalComponents[
ColorNegate[ImageAdd[ridges, Dilation[string, 2]]]] /. {2 -> 5, 8 -> 0, 4 -> 3} // Colorize, 2]
What this has accomplished, in effect, is to identify five oriented linear features. We can see three separate linear features emanating from a point of confluence. Each has two sides. I have considered the right side of the two rightmost features as being the same, but have otherwise distinguished everything else, giving the five features. The colored areas show the Voronoi diagram from these five features.
A Euclidean Allocation command based on a layer that distinguishes the three linear features (which I did not have available for this illustration) would not distinguish the different sides of each linear feature, and so it would combine the green and orange regions flanking the leftmost line; it would split the rightmost teal feature into two; and it would combine those split pieces with the corresponding beige and magenta features on their other sides.
Evidently, this raster approach has the power to construct Voronoi tessellations of arbitrary features--points, linear pieces, and even polygons, regardless of their shapes--and it can distinguish sides of linear features.
I got a pretty good result using GRASS v.generalize from the Processing Toolbox using default values:
Some lines remain untouched.
Best Answer
This is my option to do that. The process takes less than an hour.
I will start with the following raster:
It is a SRTM tile, from latitude 51 degrees to 52 degrees, from longitude 11 degrees to 12 degrees, with a spatial resolution of 1 arc second in latitudes and 2 arc seconds in longitudes.
gdalinfo input_1.tifreturns:Elevations go from 33 to 585 meters. To extract closed contours we need some more rows and columns for the edges. We can create them with a command like:
gdal_merge.py -o to_contours.tif -ul_lr 10.9 52.1 12.1 50.9 -init 1 input_1.tifThe output looks:
Now, I will extract the contours lines each 10 meters of elevation. But I will check on "Produce 3D vector".
Contours lines looks like this:
Now, we can Polygonize those closed lines. Do not worry about the attributes, we will extract the Z from the geometry.
The polygons are all single part geometries, with holes, and look like this:
Now, it is a good time to check the validity of the geometries. I have 24123 Valid geometries there.
I will populate a new Integer type field with the Z value of the first point of the exterior ring of each polygon (all points of the exterior ring has the same Z, but we need just one), with the following expression:
z(start_point(exterior_ring($geometry)))Finally, we can style the layer with that field. I have used the Assistant: