diff --git a/alg/contour.cpp b/alg/contour.cpp index 7a96c1397932..23b5b4dd542e 100644 --- a/alg/contour.cpp +++ b/alg/contour.cpp @@ -196,7 +196,7 @@ struct PolygonContourWriter std::unique_ptr currentGeometry_ = {}; OGRPolygon *currentPart_ = nullptr; OGRContourWriterInfo *poInfo_ = nullptr; - double currentLevel_; + double currentLevel_ = 0; double previousLevel_ = 0; }; @@ -749,6 +749,33 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, } } } + + // Largest requested level (levels are sorted) + const double maxLevel{fixedLevels.back()}; + + // If the maximum raster value is smaller than the last requested + // level, select the requested level that is just above the + // maximum raster value + if (maxLevel > dfMaximum) + { + for (size_t i = fixedLevels.size() - 1; i > 0; --i) + { + if (fixedLevels[i] <= dfMaximum) + { + dfMaximum = fixedLevels[i + 1]; + break; + } + } + } + + // If the maximum raster value is equal to the last requested + // level, add a small value to the maximum to avoid skipping the + // last level polygons + if (maxLevel == dfMaximum) + { + dfMaximum = std::nextafter( + dfMaximum, std::numeric_limits::infinity()); + } } PolygonContourWriter w(&oCWI, dfMinimum); diff --git a/alg/marching_squares/segment_merger.h b/alg/marching_squares/segment_merger.h index a7d095f1c023..fb16123711ee 100644 --- a/alg/marching_squares/segment_merger.h +++ b/alg/marching_squares/segment_merger.h @@ -81,12 +81,12 @@ template struct SegmentMerger } } - void addBorderSegment(int levelIdx, const Point &start, const Point &end) + void addSegment(int levelIdx, const Point &start, const Point &end) { addSegment_(levelIdx, start, end); } - void addSegment(int levelIdx, const Point &start, const Point &end) + void addBorderSegment(int levelIdx, const Point &start, const Point &end) { addSegment_(levelIdx, start, end); } diff --git a/autotest/alg/contour.py b/autotest/alg/contour.py index de1f67e2566a..29d3f4c6c4bf 100755 --- a/autotest/alg/contour.py +++ b/autotest/alg/contour.py @@ -239,7 +239,7 @@ def test_contour_real_world_case(): ("1,10,20,25,30", [1, 10, 20, 25], [10, 20, 25, 30]), ("10,20,25,30", [1, 10, 20, 25], [10, 20, 25, 30]), ("10,20,24", [1, 10, 20, 24], [10, 20, 24, 25]), - ("10,20,25", [1, 10, 20], [10, 20, 25]), + ("10,20,25", [1, 10, 20, 25], [10, 20, 25, 25]), ("0,10,20", [0, 10, 20], [10, 20, 25]), ], ) @@ -290,8 +290,12 @@ def test_contour_3(input_tif, tmp_path, fixed_levels, expected_min, expected_max i = 0 for feat in lyr: - assert feat.GetField("elevMin") == expected_min[i], i - assert feat.GetField("elevMax") == expected_max[i], i + assert feat.GetField("elevMin") == pytest.approx( + expected_min[i], abs=precision / 2 * 1.001 + ), i + assert feat.GetField("elevMax") == pytest.approx( + expected_max[i], abs=precision / 2 * 1.001 + ), i envelope = feat.GetGeometryRef().GetEnvelope() for j in range(4):