4

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?

4
  • 1
    Why are you creating polygons from the contourlines? Contourlines are, as the name says, lines - not polygons. Polygons would indicate steps in the landscape (since the whole polygon is supposed to have the same height), while contourlines just indicate "here the elevation is like X". – Erik Feb 7 '20 at 8:19
  • Also, it would help to see the size and resolution of the original elevation raster. i.e. the output of gdalinfo elev.tif` – Micha Feb 9 '20 at 6:49
  • Hi Micha. Thanks for your interest. The contours are not on elevation but on slope, and the intention is to use the polygons as a starting point in regolith-terrain mapping of the Tanami Desert in Australia. Changes in slope (breaks) are a good tool in geological mapping, and the polygons are as useful, or more useful, than contours alone. I use a bil file created in ER Mapper and/or LandSerf. The file I am attempting to polygonize is 460 x 240 km, with 90-m pixels. – Dr Richard Langford Feb 12 '20 at 4:37
  • I have used the SAGA tools with some success, but can only show the nested polygons if I dissolve the 'lines to polygons' output. I then end up with multiparts that I cannot break down, perhaps because dissolve only works with one parameter (one field in the table)...ELEV (i.e slope). Apologies for the delay...my computer packed up, possibly because of the many hours of intense processing. GDAL contour with -p is really good, but I need a bigger computer! – Dr Richard Langford Feb 12 '20 at 4:37
4

This is my option to do that. The process takes less than an hour.

I will start with the following raster:

1

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[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            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:

2


Now, I will extract the contours lines each 10 meters of elevation. But I will check on "Produce 3D vector".

3

Contours lines looks like this:

4


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:

5

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:

6

1
  • Worked perfectly on my 87000 polygons, and quick. However, contours have to start at 1, not zero. Still struggling to get zero recognized as a legitimate value. I have notionally correct polygons inside polygons, both with Z=1, but some are actually 1 and some are zero renamed as 1. My clunky process also uses 1 as the base, although I sort of fool it by renaming it as zero...not ideal. – Dr Richard Langford Feb 18 '20 at 6:11
1

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.

0

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.

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.