[gdal_contour] Fix regression when fixed level == raster max (#10885)

Fix #10854
OSGeo · Sep 28, 2024 · d04e73d · d04e73d
1 parent 68c4838
commit d04e73d
Show file tree

Hide file tree

Showing 3 changed files with 37 additions and 6 deletions.
diff --git a/alg/contour.cpp b/alg/contour.cpp
@@ -196,7 +196,7 @@ struct PolygonContourWriter
     std::unique_ptr<OGRMultiPolygon> 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<double>::infinity());
+                }
             }
 
             PolygonContourWriter w(&oCWI, dfMinimum);

diff --git a/alg/marching_squares/segment_merger.h b/alg/marching_squares/segment_merger.h
@@ -81,12 +81,12 @@ template <typename LineWriter, typename LevelGenerator> 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
@@ -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):