Skip to content

Commit

Permalink
Merge pull request #10708 from rouault/warp_int16_nodata_32767
Browse files Browse the repository at this point in the history 
Fix issues related to warping and source nodata value
  • Loading branch information
@rouault
rouault committed Sep 7, 2024
2 parents 19a9093 + cccdd7c commit 99651ea
Show file tree
Hide file tree
Showing 8 changed files with 427 additions and 25 deletions.
26 changes: 26 additions & 0 deletions alg/gdalwarper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1536,6 +1536,32 @@ void CPL_STDCALL GDALWarpResolveWorkingDataType(GDALWarpOptions *psOptions)

psOptions->eWorkingDataType = GDT_Byte;

// If none of the provided input nodata values can be represented in the
// data type of the corresponding source band, ignore them.
if (psOptions->hSrcDS && psOptions->padfSrcNoDataReal)
{
int nCountInvalidSrcNoDataReal = 0;
for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
{
GDALRasterBandH hSrcBand = GDALGetRasterBand(
psOptions->hSrcDS, psOptions->panSrcBands[iBand]);

if (hSrcBand &&
!GDALIsValueExactAs(psOptions->padfSrcNoDataReal[iBand],
GDALGetRasterDataType(hSrcBand)))
{
nCountInvalidSrcNoDataReal++;
}
}
if (nCountInvalidSrcNoDataReal == psOptions->nBandCount)
{
CPLFree(psOptions->padfSrcNoDataReal);
psOptions->padfSrcNoDataReal = nullptr;
CPLFree(psOptions->padfSrcNoDataImag);
psOptions->padfSrcNoDataImag = nullptr;
}
}

for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
{
if (psOptions->hDstDS != nullptr)
Expand Down
29 changes: 29 additions & 0 deletions autotest/cpp/test_alg.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,35 @@ TEST_F(test_alg, GDALWarpResolveWorkingDataType_padfSrcNoDataReal)
GDALDestroyWarpOptions(psOptions);
}

// GDALWarpResolveWorkingDataType: effect of padfSrcNoDataReal
TEST_F(test_alg, GDALWarpResolveWorkingDataType_padfSrcNoDataReal_with_band)
{
GDALDatasetUniquePtr poDS(GDALDriver::FromHandle(GDALGetDriverByName("MEM"))
->Create("", 1, 1, 1, GDT_Byte, nullptr));
GDALWarpOptions *psOptions = GDALCreateWarpOptions();
psOptions->hSrcDS = GDALDataset::ToHandle(poDS.get());
psOptions->nBandCount = 1;
psOptions->panSrcBands =
static_cast<int *>(CPLMalloc(psOptions->nBandCount * sizeof(int)));
psOptions->panSrcBands[0] = 1;
psOptions->padfSrcNoDataReal =
static_cast<double *>(CPLMalloc(sizeof(double)));
psOptions->padfSrcNoDataReal[0] = 0.0;
GDALWarpResolveWorkingDataType(psOptions);
EXPECT_EQ(psOptions->eWorkingDataType, GDT_Byte);

psOptions->padfSrcNoDataReal[0] = -1.0;
GDALWarpResolveWorkingDataType(psOptions);
EXPECT_EQ(psOptions->eWorkingDataType, GDT_Byte);

psOptions->eWorkingDataType = GDT_Unknown;
psOptions->padfSrcNoDataReal[0] = 2.0;
GDALWarpResolveWorkingDataType(psOptions);
EXPECT_EQ(psOptions->eWorkingDataType, GDT_Byte);

GDALDestroyWarpOptions(psOptions);
}

// GDALWarpResolveWorkingDataType: effect of padfSrcNoDataImag
TEST_F(test_alg, GDALWarpResolveWorkingDataType_padfSrcNoDataImag)
{
Expand Down
172 changes: 167 additions & 5 deletions autotest/cpp/test_gdal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,12 +260,74 @@ TEST_F(test_gdal, GDALDataTypeUnion_special_cases)
false /* complex */),
GDT_Int64);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -128, 0), GDT_Int16);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32768, 0), GDT_Int16);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32769, 0), GDT_Int32);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999, 0), GDT_Float32);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999.9876, 0),
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -128, false), GDT_Int16);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32768, false), GDT_Int16);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32769, false), GDT_Int32);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int8, 127, false), GDT_Int8);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int8, 128, false), GDT_Int16);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int16, 32767, false), GDT_Int16);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int16, 32768, false), GDT_Int32);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt16, 65535, false), GDT_UInt16);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt16, 65536, false), GDT_UInt32);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int32, INT32_MAX, false),
GDT_Int32);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int32, INT32_MAX + 1.0, false),
GDT_Int64);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt32, UINT32_MAX, false),
GDT_UInt32);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt32, UINT32_MAX + 1.0, false),
GDT_UInt64);

// (1 << 63) - 1024
EXPECT_EQ(
GDALDataTypeUnionWithValue(GDT_Int64, 9223372036854774784.0, false),
GDT_Int64);
// (1 << 63) - 512
EXPECT_EQ(
GDALDataTypeUnionWithValue(GDT_Int64, 9223372036854775296.0, false),
GDT_Float64);

// (1 << 64) - 2048
EXPECT_EQ(
GDALDataTypeUnionWithValue(GDT_UInt64, 18446744073709549568.0, false),
GDT_UInt64);
// (1 << 64) + 4096
EXPECT_EQ(
GDALDataTypeUnionWithValue(GDT_UInt64, 18446744073709555712.0, false),
GDT_Float64);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999, false),
GDT_Float32);
EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999.9876, false),
GDT_Float64);
EXPECT_EQ(GDALDataTypeUnionWithValue(
GDT_Float32, std::numeric_limits<double>::quiet_NaN(), false),
GDT_Float32);
EXPECT_EQ(GDALDataTypeUnionWithValue(
GDT_Float32, -std::numeric_limits<double>::infinity(), false),
GDT_Float32);
EXPECT_EQ(GDALDataTypeUnionWithValue(
GDT_Float32, -std::numeric_limits<double>::infinity(), false),
GDT_Float32);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float64, -99999.9876, false),
GDT_Float64);
EXPECT_EQ(GDALDataTypeUnionWithValue(
GDT_Float64, std::numeric_limits<double>::quiet_NaN(), false),
GDT_Float64);
EXPECT_EQ(GDALDataTypeUnionWithValue(
GDT_Float64, -std::numeric_limits<double>::infinity(), false),
GDT_Float64);
EXPECT_EQ(GDALDataTypeUnionWithValue(
GDT_Float64, -std::numeric_limits<double>::infinity(), false),
GDT_Float64);

EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Unknown, 0, false), GDT_Byte);
}

// Test GDALAdjustValueToDataType()
Expand Down Expand Up @@ -807,6 +869,106 @@ TEST_F(test_gdal, GDALIsValueExactAs)
GDALIsValueExactAs<double>(std::numeric_limits<double>::quiet_NaN()));
}

// Test GDALIsValueExactAs()
TEST_F(test_gdal, GDALIsValueExactAs_C_func)
{
EXPECT_TRUE(GDALIsValueExactAs(0, GDT_Byte));
EXPECT_TRUE(GDALIsValueExactAs(255, GDT_Byte));
EXPECT_FALSE(GDALIsValueExactAs(-1, GDT_Byte));
EXPECT_FALSE(GDALIsValueExactAs(256, GDT_Byte));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Byte));

EXPECT_TRUE(GDALIsValueExactAs(-128, GDT_Int8));
EXPECT_TRUE(GDALIsValueExactAs(127, GDT_Int8));
EXPECT_FALSE(GDALIsValueExactAs(-129, GDT_Int8));
EXPECT_FALSE(GDALIsValueExactAs(128, GDT_Int8));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int8));

EXPECT_TRUE(GDALIsValueExactAs(0, GDT_UInt16));
EXPECT_TRUE(GDALIsValueExactAs(65535, GDT_UInt16));
EXPECT_FALSE(GDALIsValueExactAs(-1, GDT_UInt16));
EXPECT_FALSE(GDALIsValueExactAs(65536, GDT_UInt16));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt16));

EXPECT_TRUE(GDALIsValueExactAs(-32768, GDT_Int16));
EXPECT_TRUE(GDALIsValueExactAs(32767, GDT_Int16));
EXPECT_FALSE(GDALIsValueExactAs(-32769, GDT_Int16));
EXPECT_FALSE(GDALIsValueExactAs(32768, GDT_Int16));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int16));

EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<uint32_t>::lowest(),
GDT_UInt32));
EXPECT_TRUE(
GDALIsValueExactAs(std::numeric_limits<uint32_t>::max(), GDT_UInt32));
EXPECT_FALSE(GDALIsValueExactAs(
std::numeric_limits<uint32_t>::lowest() - 1.0, GDT_UInt32));
EXPECT_FALSE(GDALIsValueExactAs(std::numeric_limits<uint32_t>::max() + 1.0,
GDT_UInt32));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt32));

EXPECT_TRUE(
GDALIsValueExactAs(std::numeric_limits<int32_t>::lowest(), GDT_Int32));
EXPECT_TRUE(
GDALIsValueExactAs(std::numeric_limits<int32_t>::max(), GDT_Int32));
EXPECT_FALSE(GDALIsValueExactAs(
std::numeric_limits<int32_t>::lowest() - 1.0, GDT_Int32));
EXPECT_FALSE(GDALIsValueExactAs(std::numeric_limits<int32_t>::max() + 1.0,
GDT_Int32));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int32));

EXPECT_TRUE(GDALIsValueExactAs(
static_cast<double>(std::numeric_limits<uint64_t>::lowest()),
GDT_UInt64));
// (1 << 64) - 2048
EXPECT_TRUE(GDALIsValueExactAs(18446744073709549568.0, GDT_UInt64));
EXPECT_FALSE(GDALIsValueExactAs(
static_cast<double>(std::numeric_limits<uint64_t>::lowest()) - 1.0,
GDT_UInt64));
// (1 << 64)
EXPECT_FALSE(GDALIsValueExactAs(18446744073709551616.0, GDT_UInt64));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt64));

EXPECT_TRUE(GDALIsValueExactAs(
static_cast<double>(std::numeric_limits<int64_t>::lowest()),
GDT_Int64));
// (1 << 63) - 1024
EXPECT_TRUE(GDALIsValueExactAs(9223372036854774784.0, GDT_Int64));
EXPECT_FALSE(GDALIsValueExactAs(
static_cast<double>(std::numeric_limits<int64_t>::lowest()) - 2048.0,
GDT_Int64));
// (1 << 63) - 512
EXPECT_FALSE(GDALIsValueExactAs(9223372036854775296.0, GDT_Int64));
EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int64));

EXPECT_TRUE(
GDALIsValueExactAs(-std::numeric_limits<float>::max(), GDT_Float32));
EXPECT_TRUE(
GDALIsValueExactAs(std::numeric_limits<float>::max(), GDT_Float32));
EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits<float>::infinity(),
GDT_Float32));
EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<float>::infinity(),
GDT_Float32));
EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<double>::quiet_NaN(),
GDT_Float32));
EXPECT_TRUE(
!GDALIsValueExactAs(-std::numeric_limits<double>::max(), GDT_Float32));
EXPECT_TRUE(
!GDALIsValueExactAs(std::numeric_limits<double>::max(), GDT_Float32));

EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits<double>::infinity(),
GDT_Float64));
EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<double>::infinity(),
GDT_Float64));
EXPECT_TRUE(
GDALIsValueExactAs(-std::numeric_limits<double>::max(), GDT_Float64));
EXPECT_TRUE(
GDALIsValueExactAs(std::numeric_limits<double>::max(), GDT_Float64));
EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<double>::quiet_NaN(),
GDT_Float64));

EXPECT_TRUE(GDALIsValueExactAs(0, GDT_CInt16));
}

#ifdef _MSC_VER
#pragma warning(pop)
#endif
Expand Down
28 changes: 28 additions & 0 deletions autotest/gdrivers/vrtwarp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -732,3 +732,31 @@ def test_vrtwarp_irasterio_optim_window_splitting():
with gdaltest.config_option("GDAL_VRT_WARP_USE_DATASET_RASTERIO", "NO"):
expected_data = warped_vrt_ds.ReadRaster()
assert warped_vrt_ds.ReadRaster() == expected_data


###############################################################################
# Test gdal.AutoCreateWarpedVRT() on a Int16 band with nodata = 32767


def test_vrtwarp_autocreatewarpedvrt_int16_nodata_32767():

ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1, gdal.GDT_Int16)
ds.SetGeoTransform([0, 1, 0, 0, 0, -1])
ds.GetRasterBand(1).SetNoDataValue(32767)
vrt_ds = gdal.AutoCreateWarpedVRT(ds)
assert vrt_ds.GetRasterBand(1).DataType == gdal.GDT_Int16
assert vrt_ds.GetRasterBand(1).GetNoDataValue() == 32767


###############################################################################
# Test gdal.AutoCreateWarpedVRT() on a source nodata value that does not fit
# the source band type


def test_vrtwarp_autocreatewarpedvrt_invalid_nodata():

ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1, gdal.GDT_Byte)
ds.SetGeoTransform([0, 1, 0, 0, 0, -1])
ds.GetRasterBand(1).SetNoDataValue(-9999)
vrt_ds = gdal.AutoCreateWarpedVRT(ds)
assert vrt_ds.GetRasterBand(1).DataType == gdal.GDT_Byte
1 change: 1 addition & 0 deletions doc/source/spelling_wordlist.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -661,6 +661,7 @@ dfTargetHeight
dfTimeout
dfTolerance
dfURY
dfValue
dfVisibleVal
dfXOff
dfXSize
Expand Down
49 changes: 36 additions & 13 deletions frmts/vrt/vrtwarped.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,9 @@ GDALDatasetH CPL_STDCALL GDALAutoCreateWarpedVRTEx(
if (psWO->padfSrcNoDataReal == nullptr &&
psWO->padfDstNoDataReal == nullptr && psWO->nSrcAlphaBand == 0)
{
// If none of the provided input nodata values can be represented in the
// data type of the corresponding source band, ignore them.
int nCountInvalidSrcNoDataReal = 0;
for (int i = 0; i < psWO->nBandCount; i++)
{
GDALRasterBandH rasterBand =
Expand All @@ -189,22 +192,42 @@ GDALDatasetH CPL_STDCALL GDALAutoCreateWarpedVRTEx(
double noDataValue =
GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);

if (hasNoDataValue)
if (hasNoDataValue &&
!GDALIsValueExactAs(noDataValue,
GDALGetRasterDataType(rasterBand)))
{
// Check if the nodata value is out of range
int bClamped = FALSE;
int bRounded = FALSE;
CPL_IGNORE_RET_VAL(GDALAdjustValueToDataType(
GDALGetRasterDataType(rasterBand), noDataValue, &bClamped,
&bRounded));
if (!bClamped)
nCountInvalidSrcNoDataReal++;
}
}

if (nCountInvalidSrcNoDataReal != psWO->nBandCount)
{
for (int i = 0; i < psWO->nBandCount; i++)
{
GDALRasterBandH rasterBand =
GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[i]);

int hasNoDataValue;
double noDataValue =
GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);

if (hasNoDataValue)
{
GDALWarpInitNoDataReal(psWO, -1e10);
if (psWO->padfSrcNoDataReal != nullptr &&
psWO->padfDstNoDataReal != nullptr)
// Check if the nodata value is out of range
int bClamped = FALSE;
int bRounded = FALSE;
CPL_IGNORE_RET_VAL(GDALAdjustValueToDataType(
GDALGetRasterDataType(rasterBand), noDataValue,
&bClamped, &bRounded));
if (!bClamped)
{
psWO->padfSrcNoDataReal[i] = noDataValue;
psWO->padfDstNoDataReal[i] = noDataValue;
GDALWarpInitNoDataReal(psWO, -1e10);
if (psWO->padfSrcNoDataReal != nullptr &&
psWO->padfDstNoDataReal != nullptr)
{
psWO->padfSrcNoDataReal[i] = noDataValue;
psWO->padfDstNoDataReal[i] = noDataValue;
}
}
}
}
Expand Down
1 change: 1 addition & 0 deletions gcore/gdal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,7 @@ GDALDataType CPL_DLL CPL_STDCALL GDALFindDataTypeForValue(double dValue,
int bComplex);
double CPL_DLL GDALAdjustValueToDataType(GDALDataType eDT, double dfValue,
int *pbClamped, int *pbRounded);
bool CPL_DLL GDALIsValueExactAs(double dfValue, GDALDataType eDT);
GDALDataType CPL_DLL CPL_STDCALL GDALGetNonComplexDataType(GDALDataType);
int CPL_DLL CPL_STDCALL GDALDataTypeIsConversionLossy(GDALDataType eTypeFrom,
GDALDataType eTypeTo);
Expand Down
Loading

0 comments on commit 99651ea

Please sign in to comment.