Merge pull request #10127 from rouault/vrt_lut_nan

VRT: fix processing of LUT where the first source value is NaN

autotest/gdrivers/aaigrid.py

@@ -31,6 +31,7 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
+import math
 import os
 import struct
 
@@ -523,3 +524,14 @@ def test_aaigrid_starting_with_nan():
     ds = gdal.Open("data/aaigrid/starting_with_nan.asc")
     assert ds.GetRasterBand(1).DataType == gdal.GDT_Float32
     assert ds.GetRasterBand(1).Checksum() == 65300
+
+
+###############################################################################
+# Test reading a file starting with nan as nodata value
+
+
+def test_aaigrid_nodata_nan():
+
+    ds = gdal.Open("data/aaigrid/nodata_nan.asc")
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_Float32
+    assert math.isnan(ds.GetRasterBand(1).GetNoDataValue())

autotest/gdrivers/data/aaigrid/nodata_nan.asc

@@ -0,0 +1,8 @@
+ncols        3
+nrows        2
+xllcorner    0
+yllcorner    0
+cellsize     1
+NODATA_value nan
+nan -1 0
+1 2 3

autotest/gdrivers/data/vrt/lut_with_nan.vrt

@@ -0,0 +1,9 @@
+<VRTDataset rasterXSize="3" rasterYSize="2">
+  <VRTRasterBand dataType="Byte" band="1">
+    <ComplexSource>
+      <SourceFilename relativeToVRT="1">../aaigrid/nodata_nan.asc</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <LUT>nan:0,0:10,2:20</LUT>
+    </ComplexSource>
+  </VRTRasterBand>
+</VRTDataset>

autotest/gdrivers/vrtlut.py

@@ -28,8 +28,12 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
+import struct
 
 import gdaltest
+import pytest
+
+from osgeo import gdal
 
 ###############################################################################
 # Simple test
@@ -39,3 +43,13 @@ def test_vrtlut_1():
 
     tst = gdaltest.GDALTest("VRT", "vrt/byte_lut.vrt", 1, 4655)
     tst.testOpen()
+
+
+###############################################################################
+
+
+@pytest.mark.require_driver("AAIGRID")
+def test_vrtlut_with_nan():
+
+    ds = gdal.Open("data/vrt/lut_with_nan.vrt")
+    assert struct.unpack("B" * 2 * 3, ds.ReadRaster()) == (0, 10, 10, 15, 20, 20)

doc/source/drivers/raster/vrt.rst

@@ -424,8 +424,9 @@ the following form:
 
 The intermediary values are calculated using a linear interpolation
 between the bounding destination values of the corresponding range.
-Source values should be monotonically non-decreasing. Clamping is performed for
-input pixel values outside of the range specified by the LUT. That is, if an
+Source values should be listed in a monotonically non-decreasing order.
+If there is a Not-A-Number (NaN) source value, it should be the first one.
+Clamping is performed for input pixel values outside of the range specified by the LUT. That is, if an
 input pixel value is lower than the minimum source value, then the destination
 value corresponding to that minimum source value is used as the output pixel value.
 And similarly for an input pixel value that is greater than the maximum source value.

frmts/aaigrid/aaigriddataset.cpp

@@ -606,6 +606,7 @@ int AAIGDataset::ParseHeader(const char *pszHeader, const char *pszDataType)
             if (pszDataType == nullptr &&
                 (strchr(pszNoData, '.') != nullptr ||
                  strchr(pszNoData, ',') != nullptr ||
+                 std::isnan(dfNoDataValue) ||
                  std::numeric_limits<int>::min() > dfNoDataValue ||
                  dfNoDataValue > std::numeric_limits<int>::max()))
             {
@@ -718,7 +719,7 @@ int GRASSASCIIDataset::ParseHeader(const char *pszHeader,
         dfNoDataValue = CPLAtofM(pszNoData);
         if (pszDataType == nullptr &&
             (strchr(pszNoData, '.') != nullptr ||
-             strchr(pszNoData, ',') != nullptr ||
+             strchr(pszNoData, ',') != nullptr || std::isnan(dfNoDataValue) ||
              std::numeric_limits<int>::min() > dfNoDataValue ||
              dfNoDataValue > std::numeric_limits<int>::max()))
         {

frmts/vrt/vrtsources.cpp

@@ -2633,9 +2633,21 @@ VRTComplexSource::XMLInit(const CPLXMLNode *psSrc, const char *pszVRTPath,
 
             // Enforce the requirement that the LUT input array is
             // monotonically non-decreasing.
-            if (nIndex > 0 &&
-                m_adfLUTInputs[nIndex] < m_adfLUTInputs[nIndex - 1])
+            if (std::isnan(m_adfLUTInputs[nIndex]) && nIndex != 0)
             {
+                CPLError(CE_Failure, CPLE_AppDefined,
+                         "A Not-A-Number (NaN) source value should be the "
+                         "first one of the LUT.");
+                m_adfLUTInputs.clear();
+                m_adfLUTOutputs.clear();
+                return CE_Failure;
+            }
+            else if (nIndex > 0 &&
+                     m_adfLUTInputs[nIndex] < m_adfLUTInputs[nIndex - 1])
+            {
+                CPLError(CE_Failure, CPLE_AppDefined,
+                         "Source values of the LUT are not listed in a "
+                         "monotonically non-decreasing order");
                 m_adfLUTInputs.clear();
                 m_adfLUTOutputs.clear();
                 return CE_Failure;
@@ -2662,20 +2674,29 @@ VRTComplexSource::XMLInit(const CPLXMLNode *psSrc, const char *pszVRTPath,
 
 double VRTComplexSource::LookupValue(double dfInput)
 {
+    auto beginIter = m_adfLUTInputs.begin();
+    auto endIter = m_adfLUTInputs.end();
+    size_t offset = 0;
+    if (std::isnan(m_adfLUTInputs[0]))
+    {
+        if (std::isnan(dfInput) || m_adfLUTInputs.size() == 1)
+            return m_adfLUTOutputs[0];
+        ++beginIter;
+        offset = 1;
+    }
+
     // Find the index of the first element in the LUT input array that
     // is not smaller than the input value.
-    int i = static_cast<int>(
-        std::lower_bound(m_adfLUTInputs.data(),
-                         m_adfLUTInputs.data() + m_adfLUTInputs.size(),
-                         dfInput) -
-        m_adfLUTInputs.data());
+    const size_t i =
+        offset +
+        std::distance(beginIter, std::lower_bound(beginIter, endIter, dfInput));
 
-    if (i == 0)
-        return m_adfLUTOutputs[0];
+    if (i == offset)
+        return m_adfLUTOutputs[offset];
 
     // If the index is beyond the end of the LUT input array, the input
     // value is larger than all the values in the array.
-    if (i == static_cast<int>(m_adfLUTInputs.size()))
+    if (i == m_adfLUTInputs.size())
         return m_adfLUTOutputs.back();
 
     if (m_adfLUTInputs[i] == dfInput)