GTiff: add support to read RAT stored in ArcGIS style .tif.vat.dbf au…

…xiliary file

autotest/gcore/tiff_read.py

@@ -5357,3 +5357,108 @@ def test_tiff_read_ovr_dimap_pleiades(tmp_path):
     ds = None
     # Check that cleaning overviews did not suppress the DIMAP XML file
     assert os.path.exists(tmp_path / "bundle" / "DIM_foo.XML")
+
+
+###############################################################################
+# Test reading ArcGIS .tif.vat.dbf auxiliary file
+
+
+@pytest.mark.parametrize("GDAL_DISABLE_READDIR_ON_OPEN", ["NO", "YES", "EMPTY_DIR"])
+def test_tiff_read_vat_dbf(GDAL_DISABLE_READDIR_ON_OPEN):
+
+    with gdal.config_option(
+        "GDAL_DISABLE_READDIR_ON_OPEN", GDAL_DISABLE_READDIR_ON_OPEN
+    ):
+        with gdal.Open("data/gtiff/testrat.tif") as ds:
+            band = ds.GetRasterBand(1)
+            rat = band.GetDefaultRAT()
+
+            if GDAL_DISABLE_READDIR_ON_OPEN == "EMPTY_DIR":
+                assert rat is None
+                return
+
+            assert rat
+            assert rat.GetColumnCount() == 9
+            assert rat.GetRowCount() == 2
+            assert [rat.GetNameOfCol(i) for i in range(9)] == [
+                "VALUE",
+                "COUNT",
+                "CLASS",
+                "Red",
+                "Green",
+                "Blue",
+                "OtherInt",
+                "OtherReal",
+                "OtherStr",
+            ]
+            assert [rat.GetUsageOfCol(i) for i in range(9)] == [
+                gdal.GFU_MinMax,
+                gdal.GFU_PixelCount,
+                gdal.GFU_Name,
+                gdal.GFU_Red,
+                gdal.GFU_Green,
+                gdal.GFU_Blue,
+                gdal.GFU_Generic,
+                gdal.GFU_Generic,
+                gdal.GFU_Generic,
+            ]
+            assert [rat.GetTypeOfCol(i) for i in range(9)] == [
+                gdal.GFT_Integer,
+                gdal.GFT_Integer,
+                gdal.GFT_String,
+                gdal.GFT_Integer,
+                gdal.GFT_Integer,
+                gdal.GFT_Integer,
+                gdal.GFT_Integer,
+                gdal.GFT_Real,
+                gdal.GFT_String,
+            ]
+            assert rat.GetValueAsInt(0, 0) == 1
+            assert rat.GetValueAsInt(0, 1) == 10
+            assert rat.GetValueAsString(0, 2) == "my class"
+            assert rat.GetValueAsInt(0, 3) == 26
+            assert rat.GetValueAsInt(0, 4) == 51
+            assert rat.GetValueAsInt(0, 5) == 128
+            assert rat.GetValueAsInt(0, 6) == 2
+            assert rat.GetValueAsDouble(0, 7) == 2.5
+            assert rat.GetValueAsString(0, 8) == "foo"
+
+            assert rat.GetValueAsInt(1, 0) == 2
+            assert rat.GetValueAsString(1, 2) == "my class2"
+            assert rat.GetValueAsString(1, 8) == "foo2"
+
+            rat = band.GetDefaultRAT()
+            assert rat
+            assert rat.GetColumnCount() == 9
+
+
+###############################################################################
+# Test reading absent ArcGIS .tif.vat.dbf auxiliary file
+
+
+@pytest.mark.parametrize("GDAL_DISABLE_READDIR_ON_OPEN", ["NO", "YES", "EMPTY_DIR"])
+def test_tiff_read_no_vat_dbf(GDAL_DISABLE_READDIR_ON_OPEN):
+
+    with gdal.config_option(
+        "GDAL_DISABLE_READDIR_ON_OPEN", GDAL_DISABLE_READDIR_ON_OPEN
+    ):
+        with gdal.Open("data/byte.tif") as ds:
+            band = ds.GetRasterBand(1)
+            assert band.GetDefaultRAT() is None
+
+
+###############################################################################
+# Test reading corrupted ArcGIS .tif.vat.dbf auxiliary file
+
+
+def test_tiff_read_corrupted_vat_dbf(tmp_vsimem):
+
+    filename = str(tmp_vsimem / "test.tif")
+    gdal.GetDriverByName("GTiff").Create(filename, 1, 1)
+    vat_dbf_filename = filename + ".vat.dbf"
+    gdal.FileFromMemBuffer(vat_dbf_filename, "")
+
+    with gdal.Open(filename) as ds:
+        band = ds.GetRasterBand(1)
+        with pytest.raises(Exception):
+            band.GetDefaultRAT()

doc/source/drivers/raster/gtiff.rst

@@ -246,6 +246,13 @@ GDALGeoTIFF. Note that all bands must use the same nodata value. When
 BASELINE or GeoTIFF profile are used, the nodata value is stored into a
 PAM .aux.xml file.
 
+Raster Attribute Table
+----------------------
+
+Starting with GDAL 3.11, Raster attribute tables stored in auxiliary
+.tif.vat.dbf files, as written by ArcGIS, can be read as GDAL Raster Attribute
+Table.
+
 Sparse files
 ------------
 
@@ -487,7 +494,7 @@ This driver supports the following creation options:
       * ``JXL`` is for JPEG-XL, and is only available when using internal libtiff and building GDAL against
         https://github.com/libjxl/libjxl . It is recommended to use JXL compression with the ``TILED=YES`` creation
         option and block size of 256x256, 512x512, or 1024x1024 pixels. Supported data types are ``Byte``,
-        ``UInt16`` and ``Float32`` only. For GDAL < 3.6.0, JXL compression may only be used alongside 
+        ``UInt16`` and ``Float32`` only. For GDAL < 3.6.0, JXL compression may only be used alongside
         ``INTERLEAVE=PIXEL`` (the default) on datasets with 4 bands or less.
 
       * ``NONE`` is the default.

frmts/gtiff/gtiffrasterband.h

@@ -15,6 +15,7 @@
 #define GTIFFRASTERBAND_H_INCLUDED
 
 #include "gdal_pam.h"
+#include "gdal_rat.h"
 
 #include "gtiff.h"
 
@@ -41,6 +42,7 @@ class GTiffRasterBand CPL_NON_FINAL : public GDALPamRasterBand
     GDALColorInterp m_eBandInterp = GCI_Undefined;
     std::set<GTiffRasterBand **> m_aSetPSelf{};
     bool m_bHaveOffsetScale = false;
+    std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
 
     int DirectIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
                  int nYSize, void *pData, int nBufXSize, int nBufYSize,

frmts/gtiff/gtiffrasterband_read.cpp

@@ -34,8 +34,39 @@
 GDALRasterAttributeTable *GTiffRasterBand::GetDefaultRAT()
 
 {
+    if (m_poRAT)
+        return m_poRAT.get();
+
     m_poGDS->LoadGeoreferencingAndPamIfNeeded();
-    return GDALPamRasterBand::GetDefaultRAT();
+    auto poRAT = GDALPamRasterBand::GetDefaultRAT();
+    if (poRAT)
+        return poRAT;
+
+    if (!GDALCanFileAcceptSidecarFile(m_poGDS->m_pszFilename))
+        return nullptr;
+    const std::string osVATDBF =
+        std::string(m_poGDS->m_pszFilename) + ".vat.dbf";
+    CSLConstList papszSiblingFiles = m_poGDS->GetSiblingFiles();
+    if (papszSiblingFiles &&
+        // cppcheck-suppress knownConditionTrueFalse
+        GDALCanReliablyUseSiblingFileList(osVATDBF.c_str()))
+    {
+        int iSibling =
+            CSLFindString(papszSiblingFiles, CPLGetFilename(osVATDBF.c_str()));
+        if (iSibling >= 0)
+        {
+            CPLString osFilename = m_poGDS->m_pszFilename;
+            osFilename.resize(strlen(m_poGDS->m_pszFilename) -
+                              strlen(CPLGetFilename(m_poGDS->m_pszFilename)));
+            osFilename += papszSiblingFiles[iSibling];
+            m_poRAT = GDALLoadVATDBF(osFilename.c_str());
+        }
+        return m_poRAT.get();
+    }
+    VSIStatBufL sStatBuf;
+    if (VSIStatL(osVATDBF.c_str(), &sStatBuf) == 0)
+        m_poRAT = GDALLoadVATDBF(osVATDBF.c_str());
+    return m_poRAT.get();
 }
 
 /************************************************************************/

gcore/CMakeLists.txt

@@ -17,6 +17,7 @@ add_library(
   gdaljp2box.cpp
   gdalmultidomainmetadata.cpp
   gdal_rat.cpp
+  gdal_rat_vat_dbf.cpp
   gdalpamproxydb.cpp
   gdalallvalidmaskband.cpp
   gdalnodatamaskband.cpp

gcore/gdal_rat.h

@@ -386,4 +386,7 @@ class CPL_DLL GDALDefaultRasterAttributeTable : public GDALRasterAttributeTable
     void RemoveStatistics() override;
 };
 
+std::unique_ptr<GDALRasterAttributeTable>
+    CPL_DLL GDALLoadVATDBF(const char *pszFilename);
+
 #endif /* ndef GDAL_RAT_H_INCLUDED */

gcore/gdal_rat_vat_dbf.cpp

@@ -0,0 +1,141 @@
+/******************************************************************************
+ *
+ * Project:  GDAL Core
+ * Purpose:  Read ArcGIS .vat.dbf raster attribute table
+ * Author:   Even Rouault <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#include "gdal_priv.h"
+#include "gdal_rat.h"
+#include "ogrsf_frmts.h"
+
+#include <algorithm>
+
+/************************************************************************/
+/*                            GDALLoadVATDBF()                          */
+/************************************************************************/
+
+/**
+ * \brief Load a ESRI .vat.dbf auxiliary file as a GDAL attribute table.
+ *
+ * @since GDAL 3.11
+ */
+std::unique_ptr<GDALRasterAttributeTable>
+GDALLoadVATDBF(const char *pszFilename)
+{
+    auto poDS = std::unique_ptr<GDALDataset>(
+        GDALDataset::Open(pszFilename, GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR,
+                          nullptr, nullptr, nullptr));
+    auto poLayer = poDS ? poDS->GetLayer(0) : nullptr;
+    if (!poLayer)
+        return nullptr;
+    auto poRAT = std::make_unique<GDALDefaultRasterAttributeTable>();
+
+    const auto poFDefn = poLayer->GetLayerDefn();
+    int iRedIdxFloat = -1;
+    int iGreenIdxFloat = -1;
+    int iBlueIdxFloat = -1;
+    const int nFieldCount = poFDefn->GetFieldCount();
+    for (int i = 0; i < nFieldCount; ++i)
+    {
+        const auto poFieldDefn = poFDefn->GetFieldDefn(i);
+        const auto eFieldType = poFieldDefn->GetType();
+        const char *pszName = poFieldDefn->GetNameRef();
+        if (EQUAL(pszName, "VALUE"))
+        {
+            if (eFieldType == OFTReal)
+                poRAT->CreateColumn(pszName, GFT_Real, GFU_MinMax);
+            else
+                poRAT->CreateColumn(pszName, GFT_Integer, GFU_MinMax);
+        }
+        else if (EQUAL(pszName, "COUNT") &&
+                 (eFieldType == OFTInteger || eFieldType == OFTInteger64))
+        {
+            poRAT->CreateColumn(pszName, GFT_Integer, GFU_PixelCount);
+        }
+        else if ((STARTS_WITH_CI(pszName, "CLASS") || EQUAL(pszName, "NAME")) &&
+                 eFieldType == OFTString)
+        {
+            poRAT->CreateColumn(pszName, GFT_String, GFU_Name);
+        }
+        else if (EQUAL(pszName, "RED") && !strstr(pszName, "min") &&
+                 !strstr(pszName, "max") && eFieldType == OFTReal)
+        {
+            // Convert from [0,1] to [0,255]
+            iRedIdxFloat = i;
+            poRAT->CreateColumn(pszName, GFT_Integer, GFU_Red);
+        }
+        else if (EQUAL(pszName, "GREEN") && !strstr(pszName, "min") &&
+                 !strstr(pszName, "max") && eFieldType == OFTReal)
+        {
+            // Convert from [0,1] to [0,255]
+            iGreenIdxFloat = i;
+            poRAT->CreateColumn(pszName, GFT_Integer, GFU_Green);
+        }
+        else if (EQUAL(pszName, "BLUE") && !strstr(pszName, "min") &&
+                 !strstr(pszName, "max") && eFieldType == OFTReal)
+        {
+            // Convert from [0,1] to [0,255]
+            iBlueIdxFloat = i;
+            poRAT->CreateColumn(pszName, GFT_Integer, GFU_Blue);
+        }
+        else
+        {
+            poRAT->CreateColumn(
+                pszName,
+                eFieldType == OFTReal ? GFT_Real
+                : (eFieldType == OFTInteger || eFieldType == OFTInteger64)
+                    ? GFT_Integer
+                    : GFT_String,
+                GFU_Generic);
+        }
+    }
+
+    int iRow = 0;
+    for (auto &&poFeature : *poLayer)
+    {
+        for (int i = 0; i < nFieldCount; ++i)
+        {
+            if (i == iRedIdxFloat || i == iGreenIdxFloat || i == iBlueIdxFloat)
+            {
+                poRAT->SetValue(
+                    iRow, i,
+                    static_cast<int>(
+                        std::clamp(255.0 * poFeature->GetFieldAsDouble(i) + 0.5,
+                                   0.0, 255.0)));
+            }
+            else
+            {
+                switch (poRAT->GDALDefaultRasterAttributeTable::GetTypeOfCol(i))
+                {
+                    case GFT_Integer:
+                    {
+                        poRAT->GDALDefaultRasterAttributeTable::SetValue(
+                            iRow, i, poFeature->GetFieldAsInteger(i));
+                        break;
+                    }
+                    case GFT_Real:
+                    {
+                        poRAT->GDALDefaultRasterAttributeTable::SetValue(
+                            iRow, i, poFeature->GetFieldAsDouble(i));
+                        break;
+                    }
+                    case GFT_String:
+                    {
+                        poRAT->GDALDefaultRasterAttributeTable::SetValue(
+                            iRow, i, poFeature->GetFieldAsString(i));
+                        break;
+                    }
+                }
+            }
+        }
+        ++iRow;
+    }
+
+    return poRAT;
+}