From 621a86a3238efc9f6016df5443509670e91a015b Mon Sep 17 00:00:00 2001 From: Alessandro Pasotti Date: Thu, 16 Jan 2025 21:06:34 +0100 Subject: [PATCH] [gdal_contour] Do not include min/max when polygonize (#11673) Fixes #11564 by not including raster min/max when polygonize and fixed levels are specified. MIN and MAX literals can be used to specify raster min/max values, e.g.: -fl MIN 10 20 30 MAX --- alg/contour.cpp | 48 +++++- alg/marching_squares/level_generator.h | 5 + alg/marching_squares/polygon_ring_appender.h | 6 +- alg/marching_squares/segment_merger.h | 38 ++++- apps/gdal_contour.cpp | 56 +++++-- autotest/alg/contour.py | 145 +++++++++++++++++- .../cpp/test_marching_squares_polygon.cpp | 94 ++++++++---- autotest/utilities/test_gdal_contour.py | 40 ++++- doc/source/programs/gdal_contour.rst | 125 ++++++++++++++- 9 files changed, 493 insertions(+), 64 deletions(-) diff --git a/alg/contour.cpp b/alg/contour.cpp index e2081cbb7121..808cf6a43270 100644 --- a/alg/contour.cpp +++ b/alg/contour.cpp @@ -141,8 +141,11 @@ struct PolygonContourWriter currentGeometry_->addGeometryDirectly(currentPart_); } - OGRPolygonContourWriter(previousLevel_, currentLevel_, - *currentGeometry_, poInfo_); + if (currentGeometry_->getNumGeometries() > 0) + { + OGRPolygonContourWriter(previousLevel_, currentLevel_, + *currentGeometry_, poInfo_); + } currentGeometry_.reset(nullptr); currentPart_ = nullptr; @@ -514,7 +517,9 @@ an averaged value from the two nearby points (in this case (12+3+5)/3). * FIXED_LEVELS=f[,f]* * * The list of fixed contour levels at which contours should be generated. - * This option has precedence on LEVEL_INTERVAL + * This option has precedence on LEVEL_INTERVAL. MIN and MAX can be used + * as special values to represent the minimum and maximum values of the + * raster. * * NODATA=f * @@ -607,7 +612,19 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, fixedLevels.resize(aosLevels.size()); for (size_t i = 0; i < fixedLevels.size(); i++) { - fixedLevels[i] = CPLAtof(aosLevels[i]); + // Handle min/max values + if (EQUAL(aosLevels[i], "MIN")) + { + fixedLevels[i] = std::numeric_limits::lowest(); + } + else if (EQUAL(aosLevels[i], "MAX")) + { + fixedLevels[i] = std::numeric_limits::max(); + } + else + { + fixedLevels[i] = CPLAtof(aosLevels[i]); + } if (i > 0 && !(fixedLevels[i] >= fixedLevels[i - 1])) { CPLError(CE_Failure, CPLE_AppDefined, @@ -712,6 +729,19 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, bool ok = true; + // Replace fixed levels min/max values with raster min/max values + if (!fixedLevels.empty()) + { + if (fixedLevels[0] == std::numeric_limits::lowest()) + { + fixedLevels[0] = dfMinimum; + } + if (fixedLevels.back() == std::numeric_limits::max()) + { + fixedLevels.back() = dfMaximum; + } + } + try { if (polygonize) @@ -719,6 +749,7 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, if (!fixedLevels.empty()) { + // If the minimum raster value is larger than the first requested // level, select the requested level that is just below the // minimum raster value @@ -781,6 +812,8 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, } // Append minimum value to fixed levels fixedLevels.push_back(dfMinimum); + // Append maximum value to fixed levels + fixedLevels.push_back(dfMaximum); } else if (contourInterval != 0) { @@ -798,6 +831,8 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, } // Append minimum value to fixed levels fixedLevels.push_back(dfMinimum); + // Append maximum value to fixed levels + fixedLevels.push_back(dfMaximum); } if (!fixedLevels.empty()) @@ -814,6 +849,11 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer, -std::numeric_limits::infinity(), dfMaximum); SegmentMerger writer( appender, levels, /* polygonize */ true); + std::vector aoiSkipLevels; + // Skip first and last levels (min/max) in polygonal case + aoiSkipLevels.push_back(0); + aoiSkipLevels.push_back(static_cast(levels.levelsCount())); + writer.setSkipLevels(aoiSkipLevels); ContourGeneratorFromRaster cg(hBand, useNoData, noDataValue, writer, levels); diff --git a/alg/marching_squares/level_generator.h b/alg/marching_squares/level_generator.h index 39310e02d5ba..a2c08ad22f46 100644 --- a/alg/marching_squares/level_generator.h +++ b/alg/marching_squares/level_generator.h @@ -118,6 +118,11 @@ class FixedLevelRangeIterator return minLevel_; } + size_t levelsCount() const + { + return count_; + } + private: const double *levels_; size_t count_; diff --git a/alg/marching_squares/polygon_ring_appender.h b/alg/marching_squares/polygon_ring_appender.h index 186140638841..d53bf13186e2 100644 --- a/alg/marching_squares/polygon_ring_appender.h +++ b/alg/marching_squares/polygon_ring_appender.h @@ -137,10 +137,14 @@ template class PolygonRingAppender void addLine(double level, LineString &ls, bool) { + auto &levelRings = rings_[level]; + if (ls.empty()) + { + return; + } // Create a new ring from the LineString Ring newRing; newRing.points.swap(ls); - auto &levelRings = rings_[level]; // This queue holds the rings to be checked std::deque queue; std::transform(levelRings.begin(), levelRings.end(), diff --git a/alg/marching_squares/segment_merger.h b/alg/marching_squares/segment_merger.h index f3557f3f4e45..9a8aa2db7df5 100644 --- a/alg/marching_squares/segment_merger.h +++ b/alg/marching_squares/segment_merger.h @@ -12,6 +12,7 @@ #ifndef MARCHING_SQUARES_SEGMENT_MERGER_H #define MARCHING_SQUARES_SEGMENT_MERGER_H +#include "cpl_error.h" #include "point.h" #include @@ -38,7 +39,7 @@ template struct SegmentMerger SegmentMerger(LineWriter &lineWriter, const LevelGenerator &levelGenerator, bool polygonize_) : polygonize(polygonize_), lineWriter_(lineWriter), lines_(), - levelGenerator_(levelGenerator) + levelGenerator_(levelGenerator), m_anSkipLevels() { } @@ -56,6 +57,13 @@ template struct SegmentMerger for (auto it = lines_.begin(); it != lines_.end(); ++it) { const int levelIdx = it->first; + + // Skip levels that should be skipped + if (std::find(m_anSkipLevels.begin(), m_anSkipLevels.end(), + levelIdx) != m_anSkipLevels.end()) + { + continue; + } while (it->second.begin() != it->second.end()) { lineWriter_.addLine(levelGenerator_.level(levelIdx), @@ -125,6 +133,23 @@ template struct SegmentMerger SegmentMerger & operator=(const SegmentMerger &) = delete; + /** + * @brief setSkipLevels sets the levels that should be skipped + * when polygonize option is set. + * @param anSkipLevels integer 0-based levels to skip. + */ + void setSkipLevels(const std::vector &anSkipLevels) + { + // Warn if polygonize is not set + if (!polygonize) + { + CPLError( + CE_Warning, CPLE_NotSupported, + "setSkipLevels is ignored when polygonize option is not set"); + } + m_anSkipLevels = anSkipLevels; + } + const bool polygonize; private: @@ -133,8 +158,12 @@ template struct SegmentMerger std::map lines_; const LevelGenerator &levelGenerator_; + // Store 0-indexed levels to skip when polygonize option is set + std::vector m_anSkipLevels; + void addSegment_(int levelIdx, const Point &start, const Point &end) { + Lines &lines = lines_[levelIdx]; if (start == end) @@ -257,11 +286,18 @@ template struct SegmentMerger typename Lines::iterator emitLine_(int levelIdx, typename Lines::iterator it, bool closed) { + Lines &lines = lines_[levelIdx]; if (lines.empty()) lines_.erase(levelIdx); // consume "it" and remove it from the list + // but clear the line if the level should be skipped + if (std::find(m_anSkipLevels.begin(), m_anSkipLevels.end(), levelIdx) != + m_anSkipLevels.end()) + { + it->ls.clear(); + } lineWriter_.addLine(levelGenerator_.level(levelIdx), it->ls, closed); return lines.erase(it); } diff --git a/apps/gdal_contour.cpp b/apps/gdal_contour.cpp index 222c810435a3..406284777525 100644 --- a/apps/gdal_contour.cpp +++ b/apps/gdal_contour.cpp @@ -42,7 +42,7 @@ struct GDALContourOptions std::string osElevAttrib; std::string osElevAttribMin; std::string osElevAttribMax; - std::vector adfFixedLevels; + std::vector aosFixedLevels; CPLStringList aosOpenOptions; CPLStringList aosCreationOptions; bool bQuiet = false; @@ -140,7 +140,7 @@ GDALContourAppOptionsGetParser(GDALContourOptions *psOptions) "integer.")); // Dealt manually as argparse::nargs_pattern::at_least_one is problematic - argParser->add_argument("-fl").scan<'g', double>().metavar("").help( + argParser->add_argument("-fl").metavar("").help( _("Name one or more \"fixed levels\" to extract.")); argParser->add_argument("-off") @@ -262,22 +262,49 @@ MAIN_START(argc, argv) CSLTokenizeString(argv[i + 1])); for (const char *pszToken : aosTokens) { - sOptions.adfFixedLevels.push_back(CPLAtof(pszToken)); + // Handle min/max special values + if (EQUAL(pszToken, "MIN")) + { + sOptions.aosFixedLevels.push_back("MIN"); + } + else if (EQUAL(pszToken, "MAX")) + { + sOptions.aosFixedLevels.push_back("MAX"); + } + else + { + sOptions.aosFixedLevels.push_back( + std::to_string(CPLAtof(pszToken))); + } } i += 1; } else { - auto isNumeric = [](const char *pszArg) -> bool + auto isNumericOrMinMax = [](const char *pszArg) -> bool { + if (EQUAL(pszArg, "MIN") || EQUAL(pszArg, "MAX")) + return true; char *pszEnd = nullptr; CPLStrtod(pszArg, &pszEnd); return pszEnd != nullptr && pszEnd[0] == '\0'; }; - while (i < argc - 1 && isNumeric(argv[i + 1])) + while (i < argc - 1 && isNumericOrMinMax(argv[i + 1])) { - sOptions.adfFixedLevels.push_back(CPLAtof(argv[i + 1])); + if (EQUAL(argv[i + 1], "MIN")) + { + sOptions.aosFixedLevels.push_back("MIN"); + } + else if (EQUAL(argv[i + 1], "MAX")) + { + sOptions.aosFixedLevels.push_back("MAX"); + } + else + { + sOptions.aosFixedLevels.push_back( + std::to_string(CPLAtof(argv[i + 1]))); + } i += 1; } } @@ -291,7 +318,7 @@ MAIN_START(argc, argv) auto argParser = GDALContourAppOptionsGetParser(&sOptions); argParser->parse_args_without_binary_name(aosArgv.List()); - if (sOptions.dfInterval == 0.0 && sOptions.adfFixedLevels.empty() && + if (sOptions.dfInterval == 0.0 && sOptions.aosFixedLevels.empty() && sOptions.dfExpBase == 0.0) { fprintf(stderr, "%s\n", argParser->usage().c_str()); @@ -465,24 +492,19 @@ MAIN_START(argc, argv) sOptions.osElevAttribMax.c_str()); char **options = nullptr; - if (!sOptions.adfFixedLevels.empty()) + if (!sOptions.aosFixedLevels.empty()) { std::string values = "FIXED_LEVELS="; - for (size_t i = 0; i < sOptions.adfFixedLevels.size(); i++) + for (size_t i = 0; i < sOptions.aosFixedLevels.size(); i++) { - const int sz = 32; - char *newValue = new char[sz + 1]; - if (i == sOptions.adfFixedLevels.size() - 1) + if (i == sOptions.aosFixedLevels.size() - 1) { - CPLsnprintf(newValue, sz + 1, "%f", sOptions.adfFixedLevels[i]); + values = values + sOptions.aosFixedLevels[i]; } else { - CPLsnprintf(newValue, sz + 1, "%f,", - sOptions.adfFixedLevels[i]); + values = values + sOptions.aosFixedLevels[i] + ","; } - values = values + std::string(newValue); - delete[] newValue; } options = CSLAddString(options, values.c_str()); } diff --git a/autotest/alg/contour.py b/autotest/alg/contour.py index e8714920fba5..bfbec65674ee 100755 --- a/autotest/alg/contour.py +++ b/autotest/alg/contour.py @@ -211,6 +211,78 @@ def test_contour_real_world_case(): # Test with -p option (polygonize) +@pytest.mark.parametrize( + "fixed_levels, expected_min, expected_max", + [ + ("10,20", [10], [20]), + ("0,20", [0], [20]), + ("20,1000", [20], [1000]), + ("20", [], []), # Nothing to do here! + ("min,20", [1], [20]), + ("min,max", [1], [25]), + ("min,25", [1], [25]), + ("min,10,max", [1, 10], [10, 25]), + ], +) +def test_contour_polygonize( + input_tif, tmp_path, fixed_levels, expected_min, expected_max +): + """This tests the min/max values of the polygonize option in simple cases + without testing the geometry itself.""" + + output_shp = str(tmp_path / "contour_polygonize.shp") + ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(output_shp) + ogr_lyr = ogr_ds.CreateLayer("contour", geom_type=ogr.wkbMultiPolygon) + field_defn = ogr.FieldDefn("ID", ogr.OFTInteger) + ogr_lyr.CreateField(field_defn) + field_defn = ogr.FieldDefn("elevMin", ogr.OFTReal) + ogr_lyr.CreateField(field_defn) + field_defn = ogr.FieldDefn("elevMax", ogr.OFTReal) + ogr_lyr.CreateField(field_defn) + + ds = gdal.Open(input_tif) + gdal.ContourGenerateEx( + ds.GetRasterBand(1), + ogr_lyr, + options=[ + "FIXED_LEVELS=" + fixed_levels, + "ID_FIELD=0", + "ELEV_FIELD_MIN=1", + "ELEV_FIELD_MAX=2", + "POLYGONIZE=TRUE", + ], + ) + ds = None + + with ogr_ds.ExecuteSQL( + "select * from contour_polygonize order by elevMin asc" + ) as lyr: + + # Get the values from the layer + values = [[], [], []] + for feat in lyr: + values[0].append(feat.GetField("elevMin")) + values[1].append(feat.GetField("elevMax")) + values[2].append(feat.GetGeometryRef().GetEnvelope()) + + assert len(values[0]) == len(expected_min), ( + values[0], + values[1], + values[2], + expected_min, + expected_max, + ) + + i = 0 + for valMin, valMax in zip(values[0], values[1]): + + assert valMin == pytest.approx( + expected_min[i], abs=precision / 2 * 1.001 + ), i + assert valMax == pytest.approx( + expected_max[i], abs=precision / 2 * 1.001 + ), i + i = i + 1 @pytest.mark.parametrize( @@ -220,13 +292,12 @@ def test_contour_real_world_case(): ("-10,0,10,20,25,30", [0, 10, 20, 25], [10, 20, 25, 30]), ("0,10,20,25,30", [0, 10, 20, 25], [10, 20, 25, 30]), ("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, 25], [10, 20, 25, 25]), - ("0,10,20", [0, 10, 20], [10, 20, 25]), + ("0,10,20,24,25", [0, 10, 20, 24], [10, 20, 24, 25]), + ("0,10,20,25", [0, 10, 20], [10, 20, 25]), ], ) def test_contour_3(input_tif, tmp_path, fixed_levels, expected_min, expected_max): + """This tests the min/max values of the polygonize option and the geometry itself.""" output_shp = str(tmp_path / "contour.shp") @@ -273,6 +344,7 @@ def test_contour_3(input_tif, tmp_path, fixed_levels, expected_min, expected_max i = 0 for feat in lyr: + assert feat.GetField("elevMin") == pytest.approx( expected_min[i], abs=precision / 2 * 1.001 ), i @@ -470,3 +542,68 @@ def test_contour_constant_raster_value(tmp_vsimem): f = lyr.GetNextFeature() assert f is None + + +############################################################################### +# Test scenario of https://github.com/OSGeo/gdal/issues/11564 + + +@pytest.mark.parametrize( + "options, polygonize, expected_elev_values", + [ + (["LEVEL_INTERVAL=10"], "TRUE", [(4, 10), (10, 20), (20, 30), (30, 36)]), + ( + ["FIXED_LEVELS=15", "LEVEL_INTERVAL=10"], + "TRUE", + [(4, 10), (10, 15), (15, 20), (20, 30), (30, 36)], + ), + ], +) +@pytest.mark.require_driver("AAIGRID") +def test_contour_lowest_fixed_value( + tmp_vsimem, options, polygonize, expected_elev_values +): + + content = """ncols 2 +nrows 2 +xllcorner 0 +yllcorner 0 +cellsize 1 +4 15 +25 36""" + + srcfilename = str(tmp_vsimem / "test.asc") + + def _create_output_ds(): + ogr_ds = ogr.GetDriverByName("Memory").CreateDataSource("") + lyr = ogr_ds.CreateLayer("contour", geom_type=ogr.wkbLineString) + lyr.CreateField(ogr.FieldDefn("ID", ogr.OFTInteger)) + lyr.CreateField(ogr.FieldDefn("ELEV_MIN", ogr.OFTReal)) + lyr.CreateField(ogr.FieldDefn("ELEV_MAX", ogr.OFTReal)) + return ogr_ds, lyr + + with gdaltest.tempfile(srcfilename, content): + + src_ds = gdal.Open(srcfilename) + ogr_ds, lyr = _create_output_ds() + + options.extend( + [ + "ID_FIELD=0", + "ELEV_FIELD_MIN=1", + "ELEV_FIELD_MAX=2", + f"POLYGONIZE={polygonize}", + ] + ) + + assert ( + gdal.ContourGenerateEx(src_ds.GetRasterBand(1), lyr, options=options) + == gdal.CE_None + ) + + # Get all elev values + elev_values = [] + for f in lyr: + elev_values.append((f["ELEV_MIN"], f["ELEV_MAX"])) + + assert elev_values == expected_elev_values, (elev_values, expected_elev_values) diff --git a/autotest/cpp/test_marching_squares_polygon.cpp b/autotest/cpp/test_marching_squares_polygon.cpp index 72d62d1580fe..6e907cbc2a2b 100644 --- a/autotest/cpp/test_marching_squares_polygon.cpp +++ b/autotest/cpp/test_marching_squares_polygon.cpp @@ -287,37 +287,75 @@ TEST_F(test_ms_polygon, four_pixels_2) // NaN NaN NaN NaN std::vector data = {155.0, 155.01, 154.99, 155.0}; - TestPolygonWriter w; - { - PolygonRingAppender appender(w); - const double levels[] = {155.0}; - FixedLevelRangeIterator levelGenerator( - levels, 1, -std::numeric_limits::infinity(), - std::numeric_limits::infinity()); - SegmentMerger, - FixedLevelRangeIterator> - writer(appender, levelGenerator, /* polygonize */ true); - ContourGenerator cg( - 2, 2, false, NaN, writer, levelGenerator); - cg.feedLine(&data[0]); - cg.feedLine(&data[2]); - } { - std::ostringstream ostr; - w.out(ostr, 155.0); - // "Polygon #0" - EXPECT_EQ(ostr.str(), - "{ { (1.4999,2) (1.4999,1.5) (0.5,0.5001) (0,0.5001) (0,1) " - "(0,1.5) (0,2) (0.5,2) (1,2) (1.4999,2) } } "); + TestPolygonWriter w; + { + PolygonRingAppender appender(w); + const double levels[] = {155.0}; + FixedLevelRangeIterator levelGenerator( + levels, 1, -std::numeric_limits::infinity(), + std::numeric_limits::infinity()); + SegmentMerger, + FixedLevelRangeIterator> + writer(appender, levelGenerator, /* polygonize */ true); + ContourGenerator cg( + 2, 2, false, NaN, writer, levelGenerator); + cg.feedLine(&data[0]); + cg.feedLine(&data[2]); + } + EXPECT_EQ(w.polygons_.size(), 2); + { + std::ostringstream ostr; + w.out(ostr, 155.0); + // "Polygon #0" + EXPECT_EQ( + ostr.str(), + "{ { (1.4999,2) (1.4999,1.5) (0.5,0.5001) (0,0.5001) (0,1) " + "(0,1.5) (0,2) (0.5,2) (1,2) (1.4999,2) } } "); + } + { + std::ostringstream ostr; + w.out(ostr, Inf); + // "Polygon #1" + EXPECT_EQ( + ostr.str(), + "{ { (1.5,2) (2,2) (2,1.5) (2,1) (2,0.5) (2,0) (1.5,0) (1,0) " + "(0.5,0) (0,0) (0,0.5) (0,0.5001) (0.5,0.5001) (1.4999,1.5) " + "(1.4999,2) (1.5,2) } } "); + } } + { - std::ostringstream ostr; - w.out(ostr, Inf); - // "Polygon #1" - EXPECT_EQ(ostr.str(), - "{ { (1.5,2) (2,2) (2,1.5) (2,1) (2,0.5) (2,0) (1.5,0) (1,0) " - "(0.5,0) (0,0) (0,0.5) (0,0.5001) (0.5,0.5001) (1.4999,1.5) " - "(1.4999,2) (1.5,2) } } "); + TestPolygonWriter w; + { + PolygonRingAppender appender(w); + const double levels[] = {155.0}; + FixedLevelRangeIterator levelGenerator( + levels, 1, -std::numeric_limits::infinity(), + std::numeric_limits::infinity()); + SegmentMerger, + FixedLevelRangeIterator> + writer(appender, levelGenerator, /* polygonize */ true); + writer.setSkipLevels({1}); + ContourGenerator cg( + 2, 2, false, NaN, writer, levelGenerator); + cg.feedLine(&data[0]); + cg.feedLine(&data[2]); + } + { + EXPECT_EQ(w.polygons_.size(), 2); + EXPECT_EQ(w.polygons_.find(Inf)->second.size(), 0); + } + + { + std::ostringstream ostr; + w.out(ostr, 155.0); + // "Polygon #0" + EXPECT_EQ( + ostr.str(), + "{ { (1.4999,2) (1.4999,1.5) (0.5,0.5001) (0,0.5001) (0,1) " + "(0,1.5) (0,2) (0.5,2) (1,2) (1.4999,2) } } "); + } } } diff --git a/autotest/utilities/test_gdal_contour.py b/autotest/utilities/test_gdal_contour.py index 538096b63b28..d5cb48287620 100755 --- a/autotest/utilities/test_gdal_contour.py +++ b/autotest/utilities/test_gdal_contour.py @@ -100,7 +100,8 @@ def test_gdal_contour_1(gdal_contour_path, testdata_tif, tmp_path): (_, err) = gdaltest.runexternal_out_and_err( gdal_contour_path + f" -a elev -i 10 {testdata_tif} {contour_shp}" ) - assert err is None or err == "", "got error/warning" + + assert err is None or err == "", "got error/warning %s" % err ds = ogr.Open(contour_shp) @@ -632,3 +633,40 @@ def test_gdal_contour_gt(gdal_contour_path, tmp_path, gt): expected_heights = [75, 76, 81, 112, 243, 441] assert lyr.GetFeatureCount() == len(expected_heights) + + +############################################################################### +# Test there MIN and MAX values are correctly computed with polygonize fixed +# levels + + +def test_gdal_contour_fl_and_i__polygonize(gdal_contour_path, testdata_tif, tmp_path): + + contour_shp = str(tmp_path / "contour.shp") + + _, err = gdaltest.runexternal_out_and_err( + gdal_contour_path + + f" -amin elev -amax elev2 -fl MIN 6 16 20 MAX -p {testdata_tif} {contour_shp}" + ) + + assert err is None or err == "", "got error/warning" + + ds = ogr.Open(contour_shp) + + with ds.ExecuteSQL("select elev, elev2 from contour order by elev asc") as lyr: + + expected_heights = [0, 6, 16, 20] + + assert lyr.GetFeatureCount() == len(expected_heights) + + i = 0 + feat = lyr.GetNextFeature() + while feat is not None: + assert feat.GetField("elev") == expected_heights[i] + assert ( + feat.GetField("elev2") == expected_heights[i + 1] + if i < len(expected_heights) - 2 + else expected_heights[i] + 5 + ) + i = i + 1 + feat = lyr.GetNextFeature() diff --git a/doc/source/programs/gdal_contour.rst b/doc/source/programs/gdal_contour.rst index 40bcfa8134c5..96bc286babd9 100644 --- a/doc/source/programs/gdal_contour.rst +++ b/doc/source/programs/gdal_contour.rst @@ -94,7 +94,7 @@ be on the right, i.e. a line string goes clockwise around a top. .. option:: -i Elevation interval between contours. - Must specify either -i or -fl or -e. + Must specify either :option:`-i` or :option:`-fl` or :option:`-e`. .. option:: -off @@ -105,7 +105,7 @@ be on the right, i.e. a line string goes clockwise around a top. .. option:: -fl - Name one or more "fixed levels" to extract. + Name one or more "fixed levels" to extract, in ascending order separated by spaces. .. option:: -e @@ -122,6 +122,15 @@ be on the right, i.e. a line string goes clockwise around a top. Generate contour polygons rather than contour lines. + When this mode is selected the polygons are created for values between + each level specified by :option:`-i` or :option:`-fl` or :option:`-e`, + in case :option:`-fl` is used alone at least two fixed levels must be specified. + + The minimum and maximum values from the raster are not automatically added to + the fixed levels list but the special values `MIN`` and `MAX`` (case insensitive) + can be used to include them. + + .. versionadded:: 2.4.0 .. option:: -gt @@ -146,12 +155,112 @@ Examples -------- .. example:: - :title: Creating contours from a DEM - .. code-block:: + :title: Creating contours from a DEM + + .. code-block:: + + gdal_contour -a elev dem.tif contour.shp -i 10.0 + + This would create 10-meter contours from the DEM data in :file:`dem.tif` and + produce a shapefile in :file:`contour.shp|shx|dbf` with the contour elevations + in the ``elev`` attribute. + +.. example:: + + :title: Creating polygonal contours from a DEM + + .. code-block:: bash + + $ cat test.asc + 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 + + This would create 10-meter polygonal contours from the DEM data in :file:`test.asc` + and produce a GeoJSON output with the contour min and max elevations in the ``min`` + and ``max`` attributes, including the minimum and maximum values from the raster. + + .. code-block:: bash + + { + "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.499999909090926, 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.500000047619043 ], [ 1.5, 0.500000047619043 ], [ 1.499999909090926, 0.5 ], [ 1.499999909090926, 0.0 ] ] ] ] } } + ] + } + +.. example:: + + :title: Creating contours from a DEM with fixed levels + + .. code-block:: bash + + $ cat test.asc + ncols 2 + nrows 2 + xllcorner 0 + yllcorner 0 + cellsize 1 + 4 15 + 25 36 + + $ gdal_contour test.asc -f GeoJSON /vsistdout/ -fl 10 20 -p -amin min -amax max + + This would create a single polygonal contour between 10 and 20 meters from the DEM data in :file:`test.asc` + and produce a GeoJSON output with the contour min and max elevations in the ``min`` and ``max`` attributes. + + + If the minimum and maximum values from the raster are desired, the special values `MIN`` and `MAX`` + (case insensitive) can be used: + + .. code-block:: bash + + $ cat test.asc + ncols 2 + nrows 2 + xllcorner 0 + yllcorner 0 + cellsize 1 + 4 15 + 25 36 + + $ gdal_contour test.asc -f GeoJSON /vsistdout/ -fl MIN 10 20 MAX -p -amin min -amax max + + This would create three polygonal contours from the DEM data in :file:`test.asc` and produce a GeoJSON output + with the contour min and max elevations in the ``min`` and ``max`` attributes, the values of these fields will + be: (4.0, 10.0), (10, 20.0) and (20.0, 36.0). + +.. example:: + + :title: Creating contours from a DEM specifying an interval and fixed levels at the same time + + .. code-block:: bash + + $ cat test.asc + ncols 2 + nrows 2 + xllcorner 0 + yllcorner 0 + cellsize 1 + 4 15 + 25 36 + + $ gdal_contour test.asc -f GeoJSON /vsistdout/ -i 10 -fl 15 -p -amin min -amax max + + Creates contours at regular 10 meter intervals and adds extra contour for a fixed 15 m level. + Finally turns areas between the contours into polygons with the contour min and max elevations + in the ``min`` and ``max`` attributes, the values of these fields will be: + (4.0, 10.0), (10, 15.0), (15, 20.0), (20.0, 30.0) and (30.0, 36.0). - gdal_contour -a elev dem.tif contour.shp -i 10.0 - This would create 10-meter contours from the DEM data in :file:`dem.tif` and - produce a shapefile in :file:`contour.shp|shx|dbf` with the contour elevations - in the ``elev`` attribute.