Skip to content

Commit

Permalink
Merge pull request #10601 from OSGeo/backport-10596-to-release/3.9
Browse files Browse the repository at this point in the history 
[Backport release/3.9] gdaldem color-relief: fix issues with entry at 0 and -exact_color_entry mode, and other issues
  • Loading branch information
@rouault
rouault authored Aug 25, 2024
2 parents 4f2dfee + f18254e commit fa98dc4
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 44 deletions.
104 changes: 61 additions & 43 deletions apps/gdaldem_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1435,7 +1435,8 @@ static int GDALColorReliefSortColors(const ColorAssociation &pA,
static void
GDALColorReliefProcessColors(ColorAssociation **ppasColorAssociation,
int *pnColorAssociation, int bSrcHasNoData,
double dfSrcNoDataValue)
double dfSrcNoDataValue,
ColorSelectionMode eColorSelectionMode)
{
ColorAssociation *pasColorAssociation = *ppasColorAssociation;
int nColorAssociation = *pnColorAssociation;
Expand All @@ -1453,12 +1454,13 @@ GDALColorReliefProcessColors(ColorAssociation **ppasColorAssociation,
ColorAssociation *pCurrent = &pasColorAssociation[i];

// NaN comparison is always false, so it handles itself
if (bSrcHasNoData && pCurrent->dfVal == dfSrcNoDataValue)
if (eColorSelectionMode != COLOR_SELECTION_EXACT_ENTRY &&
bSrcHasNoData && pCurrent->dfVal == dfSrcNoDataValue)
{
// Check if there is enough distance between the nodata value and
// its predecessor.
const double dfNewValue =
pCurrent->dfVal - std::abs(pCurrent->dfVal) * DBL_EPSILON;
const double dfNewValue = std::nextafter(
pCurrent->dfVal, -std::numeric_limits<double>::infinity());
if (dfNewValue > pPrevious->dfVal)
{
// add one just below the nodata value
Expand All @@ -1474,12 +1476,13 @@ GDALColorReliefProcessColors(ColorAssociation **ppasColorAssociation,
dfNewValue;
}
}
else if (bSrcHasNoData && pPrevious->dfVal == dfSrcNoDataValue)
else if (eColorSelectionMode != COLOR_SELECTION_EXACT_ENTRY &&
bSrcHasNoData && pPrevious->dfVal == dfSrcNoDataValue)
{
// Check if there is enough distance between the nodata value and
// its successor.
const double dfNewValue =
pPrevious->dfVal + std::abs(pPrevious->dfVal) * DBL_EPSILON;
const double dfNewValue = std::nextafter(
pPrevious->dfVal, std::numeric_limits<double>::infinity());
if (dfNewValue < pCurrent->dfVal)
{
// add one just above the nodata value
Expand Down Expand Up @@ -1650,8 +1653,25 @@ static bool GDALColorReliefGetRGBA(ColorAssociation *pasColorAssociation,
}
else
{
if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
pasColorAssociation[i - 1].dfVal != dfVal)
if (pasColorAssociation[i - 1].dfVal == dfVal)
{
*pnR = pasColorAssociation[i - 1].nR;
*pnG = pasColorAssociation[i - 1].nG;
*pnB = pasColorAssociation[i - 1].nB;
*pnA = pasColorAssociation[i - 1].nA;
return true;
}

if (pasColorAssociation[i].dfVal == dfVal)
{
*pnR = pasColorAssociation[i].nR;
*pnG = pasColorAssociation[i].nG;
*pnB = pasColorAssociation[i].nB;
*pnA = pasColorAssociation[i].nA;
return true;
}

if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
{
*pnR = 0;
*pnG = 0;
Expand All @@ -1677,15 +1697,6 @@ static bool GDALColorReliefGetRGBA(ColorAssociation *pasColorAssociation,
return true;
}

if (pasColorAssociation[i - 1].dfVal == dfVal)
{
*pnR = pasColorAssociation[i - 1].nR;
*pnG = pasColorAssociation[i - 1].nG;
*pnB = pasColorAssociation[i - 1].nB;
*pnA = pasColorAssociation[i - 1].nA;
return true;
}

if (CPLIsNan(pasColorAssociation[i - 1].dfVal))
{
*pnR = pasColorAssociation[i].nR;
Expand Down Expand Up @@ -1789,7 +1800,8 @@ static bool GDALColorReliefFindNamedColor(const char *pszColorName, int *pnR,

static ColorAssociation *
GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
const char *pszColorFilename, int *pnColors)
const char *pszColorFilename, int *pnColors,
ColorSelectionMode eColorSelectionMode)
{
VSILFILE *fpColorFile = VSIFOpenL(pszColorFilename, "rt");
if (fpColorFile == nullptr)
Expand Down Expand Up @@ -1970,7 +1982,8 @@ GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
}

GDALColorReliefProcessColors(&pasColorAssociation, &nColorAssociation,
bSrcHasNoData, dfSrcNoDataValue);
bSrcHasNoData, dfSrcNoDataValue,
eColorSelectionMode);

*pnColors = nColorAssociation;
return pasColorAssociation;
Expand Down Expand Up @@ -2080,7 +2093,7 @@ GDALColorReliefDataset::GDALColorReliefDataset(
panSourceBuf(nullptr), nCurBlockXOff(-1), nCurBlockYOff(-1)
{
pasColorAssociation = GDALColorReliefParseColorFile(
hSrcBand, pszColorFilename, &nColorAssociation);
hSrcBand, pszColorFilename, &nColorAssociation, eColorSelectionMode);

nRasterXSize = GDALGetRasterXSize(hSrcDS);
nRasterYSize = GDALGetRasterYSize(hSrcDS);
Expand Down Expand Up @@ -2220,7 +2233,7 @@ GDALColorRelief(GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand1,

int nColorAssociation = 0;
ColorAssociation *pasColorAssociation = GDALColorReliefParseColorFile(
hSrcBand, pszColorFilename, &nColorAssociation);
hSrcBand, pszColorFilename, &nColorAssociation, eColorSelectionMode);
if (pasColorAssociation == nullptr)
return CE_Failure;

Expand Down Expand Up @@ -2425,7 +2438,7 @@ static CPLErr GDALGenerateVRTColorRelief(const char *pszDstFilename,
{
int nColorAssociation = 0;
ColorAssociation *pasColorAssociation = GDALColorReliefParseColorFile(
hSrcBand, pszColorFilename, &nColorAssociation);
hSrcBand, pszColorFilename, &nColorAssociation, eColorSelectionMode);
if (pasColorAssociation == nullptr)
return CE_Failure;

Expand Down Expand Up @@ -2511,29 +2524,23 @@ static CPLErr GDALGenerateVRTColorRelief(const char *pszDstFilename,

for (int iColor = 0; iColor < nColorAssociation; iColor++)
{
if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
{
if (iColor > 1)
bOK &= VSIFPrintfL(fp, ",") > 0;
}
else if (iColor > 0)
bOK &= VSIFPrintfL(fp, ",") > 0;

const double dfVal = pasColorAssociation[iColor].dfVal;
if (iColor > 0)
bOK &= VSIFPrintfL(fp, ",") > 0;

if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
{
bOK &= VSIFPrintfL(fp, "%.18g:0,",
dfVal - fabs(dfVal) * DBL_EPSILON) > 0;
}
else if (iColor > 0 &&
eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY)
if (iColor > 0 &&
eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
dfVal !=
std::nextafter(pasColorAssociation[iColor - 1].dfVal,
std::numeric_limits<double>::infinity()))
{
const double dfMidVal =
(dfVal + pasColorAssociation[iColor - 1].dfVal) / 2.0;
bOK &=
VSIFPrintfL(
fp, "%.18g:%d", dfMidVal - fabs(dfMidVal) * DBL_EPSILON,
fp, "%.18g:%d",
std::nextafter(
dfMidVal, -std::numeric_limits<double>::infinity()),
(iBand == 0) ? pasColorAssociation[iColor - 1].nR
: (iBand == 1) ? pasColorAssociation[iColor - 1].nG
: (iBand == 2)
Expand All @@ -2546,9 +2553,17 @@ static CPLErr GDALGenerateVRTColorRelief(const char *pszDstFilename,
: (iBand == 2) ? pasColorAssociation[iColor].nB
: pasColorAssociation[iColor].nA) > 0;
}

if (eColorSelectionMode != COLOR_SELECTION_NEAREST_ENTRY)
else
{
if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
{
bOK &=
VSIFPrintfL(
fp, "%.18g:0,",
std::nextafter(
dfVal,
-std::numeric_limits<double>::infinity())) > 0;
}
if (dfVal != static_cast<double>(static_cast<int>(dfVal)))
bOK &= VSIFPrintfL(fp, "%.18g", dfVal) > 0;
else
Expand All @@ -2563,8 +2578,11 @@ static CPLErr GDALGenerateVRTColorRelief(const char *pszDstFilename,

if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY)
{
bOK &= VSIFPrintfL(fp, ",%.18g:0",
dfVal + fabs(dfVal) * DBL_EPSILON) > 0;
bOK &= VSIFPrintfL(
fp, ",%.18g:0",
std::nextafter(
dfVal,
std::numeric_limits<double>::infinity())) > 0;
}
}
bOK &= VSIFPrintfL(fp, "</LUT>\n") > 0;
Expand Down
78 changes: 77 additions & 1 deletion autotest/utilities/test_gdaldem_lib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,7 +469,7 @@ def test_gdaldem_lib_color_relief():
colorFilename="data/color_file.txt",
colorSelection="exact_color_entry",
)
assert ds.GetRasterBand(1).Checksum() == 0
assert ds.GetRasterBand(1).Checksum() == 8073

ds = gdal.DEMProcessing(
"",
Expand Down Expand Up @@ -526,6 +526,82 @@ def test_gdaldem_lib_color_relief_nodata_value(tmp_vsimem):
gdal.Unlink(colorFilename)


@pytest.mark.parametrize(
"colorSelection",
["nearest_color_entry", "exact_color_entry", "linear_interpolation"],
)
@pytest.mark.parametrize("format", ["MEM", "VRT"])
def test_gdaldem_lib_color_relief_synthetic(tmp_path, colorSelection, format):

src_filename = "" if format == "MEM" else str(tmp_path / "in.tif")
src_ds = gdal.GetDriverByName("MEM" if format == "MEM" else "GTiff").Create(
src_filename, 4, 1
)
src_ds.GetRasterBand(1).SetNoDataValue(0)
src_ds.GetRasterBand(1).WriteRaster(0, 0, 4, 1, b"\x00\x01\x02\x03")
if format != "MEM":
src_ds.Close()
src_ds = gdal.Open(src_filename)
colorFilename = tmp_path / "color_file.txt"
with open(colorFilename, "wb") as f:
f.write(b"""0 0 0 0\n1 10 11 12\n2 20 21 22\n3 30 31 32\n""")

out_filename = "" if format == "MEM" else str(tmp_path / "out.vrt")
ds = gdal.DEMProcessing(
out_filename,
src_ds,
"color-relief",
format=format,
colorFilename=colorFilename,
colorSelection=colorSelection,
)
if format != "MEM":
ds.Close()
ds = gdal.Open(out_filename)
assert struct.unpack("B" * 4, ds.GetRasterBand(1).ReadRaster()) == (0, 10, 20, 30)
assert struct.unpack("B" * 4, ds.GetRasterBand(2).ReadRaster()) == (0, 11, 21, 31)
assert struct.unpack("B" * 4, ds.GetRasterBand(3).ReadRaster()) == (0, 12, 22, 32)


@pytest.mark.parametrize(
"colorSelection",
["nearest_color_entry", "exact_color_entry", "linear_interpolation"],
)
@pytest.mark.parametrize("format", ["MEM", "VRT"])
def test_gdaldem_lib_color_relief_synthetic_nodata_255(
tmp_path, colorSelection, format
):

src_filename = "" if format == "MEM" else str(tmp_path / "in.tif")
src_ds = gdal.GetDriverByName("MEM" if format == "MEM" else "GTiff").Create(
src_filename, 4, 1
)
src_ds.GetRasterBand(1).SetNoDataValue(255)
src_ds.GetRasterBand(1).WriteRaster(0, 0, 4, 1, b"\x00\x01\x02\xFF")
if format != "MEM":
src_ds.Close()
src_ds = gdal.Open(src_filename)
colorFilename = tmp_path / "color_file.txt"
with open(colorFilename, "wb") as f:
f.write(b"""0 0 1 2\n1 10 11 12\n2 20 21 22\n255 255 255 255\n""")

out_filename = "" if format == "MEM" else str(tmp_path / "out.vrt")
ds = gdal.DEMProcessing(
out_filename,
src_ds,
"color-relief",
format=format,
colorFilename=colorFilename,
colorSelection=colorSelection,
)
if format != "MEM":
ds.Close()
ds = gdal.Open(out_filename)
assert struct.unpack("B" * 4, ds.GetRasterBand(1).ReadRaster()) == (0, 10, 20, 255)
assert struct.unpack("B" * 4, ds.GetRasterBand(2).ReadRaster()) == (1, 11, 21, 255)
assert struct.unpack("B" * 4, ds.GetRasterBand(3).ReadRaster()) == (2, 12, 22, 255)


###############################################################################
# Test gdaldem tpi

Expand Down

0 comments on commit fa98dc4

Please sign in to comment.