gdal_contour generates polygon for value not specified in fixed levels

[sirisec]

What is the bug?

When using gdal_contour and specifying fixed levels to contour, it seems to output the lowest value in the raster data regardless of the fixed levels specified in the command.

Steps to reproduce the issue

Using the supplied grib2 file (unzip it first), run the following command

gdal_contour -p -amin value -fl 20 22 22 24 26 28 30 32 34 36 38 40 44 48 52 56 60 64 68 72 76 80 ./base-ref.grib2 ./base-ref.json

The output JSON will have all the fixed levels, but will contour the lowest value in the grib file (-10) which drastically bloats the output.

base-ref.grib2.zip

Versions and provenance

Mac OS: 15.0
GDAL 3.10.0, released 2024/11/01

Additional context

No response

[rouault]

CC @elpaso

[elpaso]

@rouault I have looked to the issue and I have wrote a reproducer test and a patch to fix the issue, but the patch breaks other tests (test_contour_3) and therefore I am not convinced that this is an actual bug but perhaps a documentation issue.

Looking at the the code and the tests it seems that the polygonize option is expected to create polygons that include all pixels that are between a range of values, the fixed levels define the upper value of the ranges and the minimum value of the first range is taken from the raster's minimum value, that would be -10 in the sample attached to this issue.

If you believe that this is a bug, I can submit the patch and rewrite the test_contour_3 to behave accordingly and complete the documentation of the utility and the library function to match the new behavior.

[rouault]

Looking at the the code and the tests it seems that the polygonize option is expected to create polygons that include all pixels that are between a range of values, the fixed levels define the upper value of the ranges and the minimum value of the first range is taken from the raster's minimum value

yes, I can see it makes sense from a coding point of view, but it is also a bit counter-intuitive from a user point of view.
So I'd be tempted to change it (master only), but we should make sure that it works properly with the -i switch, that is if you have a raster with values in range [4,36] and you specify -i 10, then you would get classes for [4,10],[10,20],[20,30] and [30,36]

e.g. given test.asc with:

ncols        2
nrows        2
xllcorner    0
yllcorner    0
cellsize     1
4 15
25 36

$ gdal_contour test.asc -f GeoJSON /vsistdout/ -i 10 -p -amin min -amax max
{
"type": "FeatureCollection",
"name": "contour",
"features": [
{ "type": "Feature", "properties": { "ID": 0, "min": 4.0, "max": 10.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 0.5, 1.214285714285714 ], [ 1.045454545454545, 1.5 ], [ 1.045454545454545, 2.0 ], [ 1.0, 2.0 ], [ 0.5, 2.0 ], [ 0.0, 2.0 ], [ 0.0, 1.5 ], [ 0.0, 1.214285714285714 ], [ 0.5, 1.214285714285714 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 1, "min": 10.0, "max": 20.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 1.5, 1.261904761904762 ], [ 2.0, 1.261904761904762 ], [ 2.0, 1.5 ], [ 2.0, 2.0 ], [ 1.5, 2.0 ], [ 1.045454545454545, 2.0 ], [ 1.045454545454545, 1.5 ], [ 0.5, 1.214285714285714 ], [ 0.0, 1.214285714285714 ], [ 0.0, 1.0 ], [ 0.0, 0.738095238095238 ], [ 0.5, 0.738095238095238 ], [ 1.5, 1.261904761904762 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 2, "min": 20.0, "max": 30.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 0.954545454545455, 0.0 ], [ 0.954545454545455, 0.5 ], [ 1.5, 0.785714285714286 ], [ 2.0, 0.785714285714286 ], [ 2.0, 1.0 ], [ 2.0, 1.261904761904762 ], [ 1.5, 1.261904761904762 ], [ 0.5, 0.738095238095238 ], [ 0.0, 0.738095238095238 ], [ 0.0, 0.5 ], [ 0.0, 0.0 ], [ 0.5, 0.0 ], [ 0.954545454545455, 0.0 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 3, "min": 30.0, "max": 36.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 1.5, 0.0 ], [ 1.0, 0.0 ], [ 0.954545454545455, 0.0 ], [ 0.954545454545455, 0.5 ], [ 1.5, 0.785714285714286 ], [ 2.0, 0.785714285714286 ], [ 2.0, 0.5 ], [ 2.0, 0.0 ], [ 1.5, 0.0 ] ] ] ] } }
]
}

What is a bit tricky is when the user specify both -i and -fl. Currently this does

$ gdal_contour test.asc -f GeoJSON /vsistdout/ -i 10 -p -amin min -amax max -fl 15
{
"type": "FeatureCollection",
"name": "contour",
"features": [
{ "type": "Feature", "properties": { "ID": 0, "min": 4.0, "max": 10.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 0.5, 1.214285714285714 ], [ 1.045454545454545, 1.5 ], [ 1.045454545454545, 2.0 ], [ 1.0, 2.0 ], [ 0.5, 2.0 ], [ 0.0, 2.0 ], [ 0.0, 1.5 ], [ 0.0, 1.214285714285714 ], [ 0.5, 1.214285714285714 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 1, "min": 10.0, "max": 15.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 0.5, 0.976190476190476 ], [ 1.499999909090926, 1.5 ], [ 1.499999909090926, 2.0 ], [ 1.045454545454545, 2.0 ], [ 1.045454545454545, 1.5 ], [ 0.5, 1.214285714285714 ], [ 0.0, 1.214285714285714 ], [ 0.0, 1.0 ], [ 0.0, 0.976190476190476 ], [ 0.5, 0.976190476190476 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 2, "min": 15.0, "max": 20.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 1.5, 1.261904761904762 ], [ 2.0, 1.261904761904762 ], [ 2.0, 1.5 ], [ 2.0, 2.0 ], [ 1.5, 2.0 ], [ 1.499999909090926, 2.0 ], [ 1.499999909090926, 1.5 ], [ 0.5, 0.976190476190476 ], [ 0.0, 0.976190476190476 ], [ 0.0, 0.738095238095238 ], [ 0.5, 0.738095238095238 ], [ 1.5, 1.261904761904762 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 3, "min": 20.0, "max": 30.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 0.954545454545455, 0.0 ], [ 0.954545454545455, 0.5 ], [ 1.5, 0.785714285714286 ], [ 2.0, 0.785714285714286 ], [ 2.0, 1.0 ], [ 2.0, 1.261904761904762 ], [ 1.5, 1.261904761904762 ], [ 0.5, 0.738095238095238 ], [ 0.0, 0.738095238095238 ], [ 0.0, 0.5 ], [ 0.0, 0.0 ], [ 0.5, 0.0 ], [ 0.954545454545455, 0.0 ] ] ] ] } },
{ "type": "Feature", "properties": { "ID": 4, "min": 30.0, "max": 36.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 1.5, 0.0 ], [ 1.0, 0.0 ], [ 0.954545454545455, 0.0 ], [ 0.954545454545455, 0.5 ], [ 1.5, 0.785714285714286 ], [ 2.0, 0.785714285714286 ], [ 2.0, 0.5 ], [ 2.0, 0.0 ], [ 1.5, 0.0 ] ] ] ] } }
]
}

and I believe this should still behave the same.

So the change should only be done if "-p" and "-fl" are specified, but not when "-i" or "-e" are specified

It would be good to add numerous examples (e.g. using above example) in the documentation to explain the subtelties.