Permalink
Browse files

VRT: deal with serialized nodata value that is slightly outside Float…

…32 validity range (fixes #1071)
  • Loading branch information...
rouault committed Nov 2, 2018
1 parent 3c003cb commit 35bffe8063631f82da77398a37900c784b2b414f
@@ -0,0 +1,11 @@
<VRTDataset rasterXSize="4" rasterYSize="1">
<VRTRasterBand dataType="Float32" band="1">
<NoDataValue>-3.402823466385289e+38</NoDataValue>
<ColorInterp>Gray</ColorInterp>
<ComplexSource>
<SourceFilename relativeToVRT="1">minfloat.tif</SourceFilename>
<SourceBand>1</SourceBand>
<NODATA>-3.402823466385289e+38</NODATA>
</ComplexSource>
</VRTRasterBand>
</VRTDataset>
View
@@ -1391,6 +1391,47 @@ def vrt_read_31():
return 'success'
###############################################################################
# Test reading a VRT where the NODATA & NoDataValue are slighly below the
# minimum float value (https://github.com/OSGeo/gdal/issues/1071)
def vrt_float32_with_nodata_slightly_below_float_min():
shutil.copyfile('data/minfloat.tif', 'tmp/minfloat.tif')
shutil.copyfile('data/minfloat_nodata_slightly_out_of_float.vrt',
'tmp/minfloat_nodata_slightly_out_of_float.vrt')
gdal.Unlink('tmp/minfloat_nodata_slightly_out_of_float.vrt.aux.xml')
ds = gdal.Open('tmp/minfloat_nodata_slightly_out_of_float.vrt')
nodata = ds.GetRasterBand(1).GetNoDataValue()
stats = ds.GetRasterBand(1).ComputeStatistics(False)
ds = None
vrt_content = open('tmp/minfloat_nodata_slightly_out_of_float.vrt', 'rt').read()
gdal.Unlink('tmp/minfloat.tif')
gdal.Unlink('tmp/minfloat_nodata_slightly_out_of_float.vrt')
# Check that the values were 'normalized' when regenerating the VRT
if vrt_content.find('-3.402823466385289') >= 0:
gdaltest.post_reason('did not get expected nodata in rewritten VRT')
print(vrt_content)
return 'fail'
if nodata != -3.4028234663852886e+38:
gdaltest.post_reason('did not get expected nodata')
print("%.18g" % nodata)
return 'fail'
if stats != [-3.0, 5.0, 1.0, 4.0]:
gdaltest.post_reason('did not get expected stats')
print(stats)
return 'fail'
return 'success'
for item in init_list:
ut = gdaltest.GDALTest('VRT', item[0], item[1], item[2])
if ut is None:
@@ -1429,6 +1470,7 @@ def vrt_read_31():
gdaltest_list.append(vrt_read_29)
gdaltest_list.append(vrt_read_30)
gdaltest_list.append(vrt_read_31)
gdaltest_list.append(vrt_float32_with_nodata_slightly_below_float_min)
if __name__ == '__main__':
@@ -33,6 +33,7 @@
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <limits>
#include <memory>
#include <vector>
@@ -620,10 +621,26 @@ CPLXMLNode *VRTRasterBand::SerializeToXML( const char *pszVRTPath )
if( m_bNoDataValueSet )
{
if (CPLIsNan(m_dfNoDataValue))
{
CPLSetXMLValue( psTree, "NoDataValue", "nan");
}
else if( eDataType == GDT_Float32 &&
m_dfNoDataValue == -std::numeric_limits<float>::max() )
{
// To avoid rounding out of the range of float
CPLSetXMLValue( psTree, "NoDataValue", "-3.4028234663852886e+38");
}
else if( eDataType == GDT_Float32 &&
m_dfNoDataValue == std::numeric_limits<float>::max() )
{
// To avoid rounding out of the range of float
CPLSetXMLValue( psTree, "NoDataValue", "3.4028234663852886e+38");
}
else
{
CPLSetXMLValue( psTree, "NoDataValue",
CPLSPrintf( "%.16g", m_dfNoDataValue ) );
}
}
if( m_bHideNoDataValue )
@@ -775,6 +792,11 @@ CPLXMLNode *VRTRasterBand::SerializeToXML( const char *pszVRTPath )
CPLErr VRTRasterBand::SetNoDataValue( double dfNewValue )
{
if( eDataType == GDT_Float32 )
{
dfNewValue = GDALAdjustNoDataCloseToFloatMax(dfNewValue);
}
m_bNoDataValueSet = TRUE;
m_dfNoDataValue = dfNewValue;
@@ -38,6 +38,7 @@
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <limits>
#include <string>
#include "cpl_conv.h"
@@ -2091,9 +2092,23 @@ CPLXMLNode *VRTComplexSource::SerializeToXML( const char *pszVRTPath )
{
if( CPLIsNan(m_dfNoDataValue) )
CPLSetXMLValue( psSrc, "NODATA", "nan");
else if( m_poRasterBand->GetRasterDataType() == GDT_Float32 &&
m_dfNoDataValue == -std::numeric_limits<float>::max() )
{
// To avoid rounding out of the range of float
CPLSetXMLValue( psSrc, "NODATA", "-3.4028234663852886e+38");
}
else if( m_poRasterBand->GetRasterDataType() == GDT_Float32 &&
m_dfNoDataValue == std::numeric_limits<float>::max() )
{
// To avoid rounding out of the range of float
CPLSetXMLValue( psSrc, "NODATA", "3.4028234663852886e+38");
}
else
{
CPLSetXMLValue( psSrc, "NODATA",
CPLSPrintf("%.16g", m_dfNoDataValue) );
}
}
switch( m_eScalingType )
@@ -2229,6 +2244,10 @@ CPLErr VRTComplexSource::XMLInit( CPLXMLNode *psSrc, const char *pszVRTPath,
{
m_bNoDataSet = TRUE;
m_dfNoDataValue = CPLAtofM( CPLGetXMLValue(psSrc, "NODATA", "0") );
if( m_poRasterBand->GetRasterDataType() == GDT_Float32 )
{
m_dfNoDataValue = GDALAdjustNoDataCloseToFloatMax(m_dfNoDataValue);
}
}
if( CPLGetXMLValue(psSrc, "LUT", nullptr) != nullptr )
@@ -2500,7 +2519,9 @@ CPLErr VRTComplexSource::RasterIOInternal( int nReqXOff, int nReqYOff,
const bool bIsComplex = CPL_TO_BOOL( GDALDataTypeIsComplex(eBufType) );
const int nWordSize = GDALGetDataTypeSizeBytes(eWrkDataType);
const bool bNoDataSetIsNan = m_bNoDataSet && CPLIsNan(m_dfNoDataValue);
const bool bNoDataSetAndNotNan = m_bNoDataSet && !CPLIsNan(m_dfNoDataValue);
const bool bNoDataSetAndNotNan = m_bNoDataSet && !CPLIsNan(m_dfNoDataValue) &&
GDALIsValueInRange<WorkingDT>(m_dfNoDataValue);
const auto fWorkingDataTypeNoData = static_cast<WorkingDT>(m_dfNoDataValue);
WorkingDT *pafData = nullptr;
if( m_eScalingType == VRT_SCALING_LINEAR &&
@@ -2582,8 +2603,7 @@ CPLErr VRTComplexSource::RasterIOInternal( int nReqXOff, int nReqYOff,
if( bNoDataSetIsNan && CPLIsNan(fResult) )
continue;
if( bNoDataSetAndNotNan &&
GDALIsValueInRange<WorkingDT>(m_dfNoDataValue) &&
ARE_REAL_EQUAL(fResult, static_cast<WorkingDT>(m_dfNoDataValue)) )
ARE_REAL_EQUAL(fResult, fWorkingDataTypeNoData) )
continue;
if( m_nColorTableComponent )
View
@@ -4019,3 +4019,17 @@ int GDALCanFileAcceptSidecarFile(const char* pszFilename)
return FALSE;
return TRUE;
}
/************************************************************************/
/* GDALAdjustNoDataCloseToFloatMax() */
/************************************************************************/
double GDALAdjustNoDataCloseToFloatMax(double dfVal)
{
const auto kMaxFloat = std::numeric_limits<float>::max();
if( std::fabs(dfVal - -kMaxFloat) < 1e-10 * kMaxFloat )
return -kMaxFloat;
if( std::fabs(dfVal - kMaxFloat) < 1e-10 * kMaxFloat )
return kMaxFloat;
return dfVal;
}
View
@@ -1811,6 +1811,8 @@ template<class T> inline bool ARE_REAL_EQUAL(T fVal1, T fVal2, int ulp = 2)
std::abs(fVal1 - fVal2) < std::numeric_limits<float>::epsilon() * std::abs(fVal1+fVal2) * ulp;
}
double GDALAdjustNoDataCloseToFloatMax(double dfVal);
#define DIV_ROUND_UP(a, b) ( ((a) % (b)) == 0 ? ((a) / (b)) : (((a) / (b)) + 1) )
// Number of data samples that will be used to compute approximate statistics

0 comments on commit 35bffe8

Please sign in to comment.