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?
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.tif returns:
Driver: GTiff/GeoTIFF Files: input_1.tif Size is 1801, 3601 Coordinate System is: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER, ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (10.999722222222223,52.000138888888884) Pixel Size = (0.000555555555556,-0.000277777777778) Metadata: AREA_OR_POINT=Point DTED_CompilationDate=0002 DTED_DataEdition=02 DTED_DigitizingSystem=SRTM DTED_HorizontalAccuracy=0009 DTED_HorizontalDatum=WGS84 DTED_MaintenanceDate=0000 DTED_MaintenanceDescription=0000 DTED_MatchMergeDate=0000 DTED_MatchMergeVersion=A DTED_NimaDesignator=DTED2 DTED_OriginLatitude=0510000N DTED_OriginLongitude=0110000E DTED_Producer=USCNIMA DTED_RelHorizontalAccuracy=NA DTED_RelVerticalAccuracy=0003 DTED_SecurityCode_DSI=U DTED_SecurityCode_UHL=U DTED_UniqueRef_DSI=E02 068 DTED_UniqueRef_UHL=E02 068 DTED_VerticalAccuracy_ACC=0003 DTED_VerticalAccuracy_UHL=0003 DTED_VerticalDatum=E96 Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 10.9997222, 52.0001389) ( 10d59'59.00"E, 52d 0' 0.50"N) Lower Left ( 10.9997222, 50.9998611) ( 10d59'59.00"E, 50d59'59.50"N) Upper Right ( 12.0002778, 52.0001389) ( 12d 0' 1.00"E, 52d 0' 0.50"N) Lower Right ( 12.0002778, 50.9998611) ( 12d 0' 1.00"E, 50d59'59.50"N) Center ( 11.5000000, 51.5000000) ( 11d30' 0.00"E, 51d30' 0.00"N) Band 1 Block=1801x2 Type=Int16, ColorInterp=Gray NoData Value=-32767 Unit Type: m
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.tif
The 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:
Finally, we can style the layer with that field. I have used the Assistant:
Micha's solution works well, but appears not to allow for zero (null) values to be contoured. Building on this excellent solution I now use r.contour with a specific advanced parameter.
r.contour will handle negative and positive values as well as zero, and geometry is preserved if -m is added as a v.out.ogr output layers option. -c also works, but not preferred, as both export the data as multi-features.
Not setting a minimum number of points can be a problem, so I always choose 20 or 50.
The best solution for small rasters is GDAL_contour with -amin, -amax and -p. Categorizing using min or max gives a good result. However, for large rasters it is proving to be too slow.
The most pragmatic (very fast) solution is GRASS r.reclass, with a reclass rules text file for a mixed set of values, then GDAL_polygonize, and this solved all problems with holes using other methods. You don't need the bounding box of zero values as you are going directly from the raster. Smoothing the polygons from rasters gives a perfectly acceptable result, with good geometry. Integers are needed for the rules file, but not a significant problem if you use the raster and/or field calculator to get real values after polygonizing.
A more refined and faster procedure, which does not have a problem with zero or negative is as follows:
Make sure the raster has values that can be divided into classes as integers, so multiply by 10, 100 or whatever gives you large numbers. These can be real numbers, not necessarily integers; in the next step they will be rounded.
Use GRASS r.reclass, and after go to Symbology -->Delete All, then Classify to see if it is what you want.
Finally, GRASS r.to.vect, with no name needed for attribute as the "cat" field will hold the value you need, Smooth corners selected (if you want) and Use categories selected. Put -m in v.out.ogr layer options for multipart.
Works well in Windows but of course much faster in Ubuntu.