Skip to content

Commit

Permalink
[gdal_contour] Do not include min/max when polygonize (#11673)
Browse files Browse the repository at this point in the history 
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
  • Loading branch information
@elpaso
elpaso authored Jan 16, 2025
1 parent 25b46ab commit 621a86a
Show file tree
Hide file tree
Showing 9 changed files with 493 additions and 64 deletions.
48 changes: 44 additions & 4 deletions alg/contour.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
*
Expand Down Expand Up @@ -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<double>::lowest();
}
else if (EQUAL(aosLevels[i], "MAX"))
{
fixedLevels[i] = std::numeric_limits<double>::max();
}
else
{
fixedLevels[i] = CPLAtof(aosLevels[i]);
}
if (i > 0 && !(fixedLevels[i] >= fixedLevels[i - 1]))
{
CPLError(CE_Failure, CPLE_AppDefined,
Expand Down Expand Up @@ -712,13 +729,27 @@ 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<double>::lowest())
{
fixedLevels[0] = dfMinimum;
}
if (fixedLevels.back() == std::numeric_limits<double>::max())
{
fixedLevels.back() = dfMaximum;
}
}

try
{
if (polygonize)
{

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
Expand Down Expand Up @@ -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)
{
Expand All @@ -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())
Expand All @@ -814,6 +849,11 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
-std::numeric_limits<double>::infinity(), dfMaximum);
SegmentMerger<RingAppender, FixedLevelRangeIterator> writer(
appender, levels, /* polygonize */ true);
std::vector<int> aoiSkipLevels;
// Skip first and last levels (min/max) in polygonal case
aoiSkipLevels.push_back(0);
aoiSkipLevels.push_back(static_cast<int>(levels.levelsCount()));
writer.setSkipLevels(aoiSkipLevels);
ContourGeneratorFromRaster<decltype(writer),
FixedLevelRangeIterator>
cg(hBand, useNoData, noDataValue, writer, levels);
Expand Down
5 changes: 5 additions & 0 deletions alg/marching_squares/level_generator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,11 @@ class FixedLevelRangeIterator
return minLevel_;
}

size_t levelsCount() const
{
return count_;
}

private:
const double *levels_;
size_t count_;
Expand Down
6 changes: 5 additions & 1 deletion alg/marching_squares/polygon_ring_appender.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,10 +137,14 @@ template <typename PolygonWriter> 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<Ring *> queue;
std::transform(levelRings.begin(), levelRings.end(),
Expand Down
38 changes: 37 additions & 1 deletion alg/marching_squares/segment_merger.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#ifndef MARCHING_SQUARES_SEGMENT_MERGER_H
#define MARCHING_SQUARES_SEGMENT_MERGER_H

#include "cpl_error.h"
#include "point.h"

#include <list>
Expand All @@ -38,7 +39,7 @@ template <typename LineWriter, typename LevelGenerator> struct SegmentMerger
SegmentMerger(LineWriter &lineWriter, const LevelGenerator &levelGenerator,
bool polygonize_)
: polygonize(polygonize_), lineWriter_(lineWriter), lines_(),
levelGenerator_(levelGenerator)
levelGenerator_(levelGenerator), m_anSkipLevels()
{
}

Expand All @@ -56,6 +57,13 @@ template <typename LineWriter, typename LevelGenerator> 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),
Expand Down Expand Up @@ -125,6 +133,23 @@ template <typename LineWriter, typename LevelGenerator> struct SegmentMerger
SegmentMerger<LineWriter, LevelGenerator> &
operator=(const SegmentMerger<LineWriter, LevelGenerator> &) = 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<int> &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:
Expand All @@ -133,8 +158,12 @@ template <typename LineWriter, typename LevelGenerator> struct SegmentMerger
std::map<int, Lines> lines_;
const LevelGenerator &levelGenerator_;

// Store 0-indexed levels to skip when polygonize option is set
std::vector<int> m_anSkipLevels;

void addSegment_(int levelIdx, const Point &start, const Point &end)
{

Lines &lines = lines_[levelIdx];

if (start == end)
Expand Down Expand Up @@ -257,11 +286,18 @@ template <typename LineWriter, typename LevelGenerator> 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);
}
Expand Down
56 changes: 39 additions & 17 deletions apps/gdal_contour.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ struct GDALContourOptions
std::string osElevAttrib;
std::string osElevAttribMin;
std::string osElevAttribMax;
std::vector<double> adfFixedLevels;
std::vector<std::string> aosFixedLevels;
CPLStringList aosOpenOptions;
CPLStringList aosCreationOptions;
bool bQuiet = false;
Expand Down Expand Up @@ -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("<level>").help(
argParser->add_argument("-fl").metavar("<level>").help(
_("Name one or more \"fixed levels\" to extract."));

argParser->add_argument("-off")
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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());
Expand Down Expand Up @@ -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());
}
Expand Down
Loading

0 comments on commit 621a86a

Please sign in to comment.