diff --git a/autotest/gcore/data/gtiff/testrat.tif b/autotest/gcore/data/gtiff/testrat.tif new file mode 100644 index 000000000000..6fddac128624 Binary files /dev/null and b/autotest/gcore/data/gtiff/testrat.tif differ diff --git a/autotest/gcore/data/gtiff/testrat.tif.vat.dbf b/autotest/gcore/data/gtiff/testrat.tif.vat.dbf new file mode 100644 index 000000000000..34472371c776 Binary files /dev/null and b/autotest/gcore/data/gtiff/testrat.tif.vat.dbf differ diff --git a/autotest/gcore/tiff_read.py b/autotest/gcore/tiff_read.py index 48c6910ae2d8..32b694bb35a7 100755 --- a/autotest/gcore/tiff_read.py +++ b/autotest/gcore/tiff_read.py @@ -5361,3 +5361,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() diff --git a/doc/source/drivers/raster/gtiff.rst b/doc/source/drivers/raster/gtiff.rst index f86a9fef8546..f6d4f722c42a 100644 --- a/doc/source/drivers/raster/gtiff.rst +++ b/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. diff --git a/frmts/gtiff/gtiffrasterband.h b/frmts/gtiff/gtiffrasterband.h index 4543b651ecb0..8aab99209a79 100644 --- a/frmts/gtiff/gtiffrasterband.h +++ b/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 m_aSetPSelf{}; bool m_bHaveOffsetScale = false; + std::unique_ptr m_poRAT{}; int DirectIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, void *pData, int nBufXSize, int nBufYSize, diff --git a/frmts/gtiff/gtiffrasterband_read.cpp b/frmts/gtiff/gtiffrasterband_read.cpp index b888dace5a2c..63b381671375 100644 --- a/frmts/gtiff/gtiffrasterband_read.cpp +++ b/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(); } /************************************************************************/ diff --git a/gcore/CMakeLists.txt b/gcore/CMakeLists.txt index 1625d408500b..7643c4a7f249 100644 --- a/gcore/CMakeLists.txt +++ b/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 diff --git a/gcore/gdal_rat.cpp b/gcore/gdal_rat.cpp index 5378a28ec824..fcc21fb553e1 100644 --- a/gcore/gdal_rat.cpp +++ b/gcore/gdal_rat.cpp @@ -648,11 +648,57 @@ CPLXMLNode *GDALRasterAttributeTable::Serialize() const snprintf(szValue, sizeof(szValue), "%d", static_cast(GetTypeOfCol(iCol))); - CPLCreateXMLElementAndValue(psCol, "Type", szValue); + CPLXMLNode *psType = + CPLCreateXMLElementAndValue(psCol, "Type", szValue); + const char *pszTypeStr = "String"; + switch (GetTypeOfCol(iCol)) + { + case GFT_Integer: + pszTypeStr = "Integer"; + break; + case GFT_Real: + pszTypeStr = "Real"; + break; + case GFT_String: + break; + } + CPLAddXMLAttributeAndValue(psType, "typeAsString", pszTypeStr); snprintf(szValue, sizeof(szValue), "%d", static_cast(GetUsageOfCol(iCol))); - CPLCreateXMLElementAndValue(psCol, "Usage", szValue); + CPLXMLNode *psUsage = + CPLCreateXMLElementAndValue(psCol, "Usage", szValue); + const char *pszUsageStr = ""; + +#define USAGE_STR(x) \ + case GFU_##x: \ + pszUsageStr = #x; \ + break + switch (GetUsageOfCol(iCol)) + { + USAGE_STR(Generic); + USAGE_STR(PixelCount); + USAGE_STR(Name); + USAGE_STR(Min); + USAGE_STR(Max); + USAGE_STR(MinMax); + USAGE_STR(Red); + USAGE_STR(Green); + USAGE_STR(Blue); + USAGE_STR(Alpha); + USAGE_STR(RedMin); + USAGE_STR(GreenMin); + USAGE_STR(BlueMin); + USAGE_STR(AlphaMin); + USAGE_STR(RedMax); + USAGE_STR(GreenMax); + USAGE_STR(BlueMax); + USAGE_STR(AlphaMax); + case GFU_MaxCount: + break; + } +#undef USAGE_STR + CPLAddXMLAttributeAndValue(psUsage, "usageAsString", pszUsageStr); } /* -------------------------------------------------------------------- */ diff --git a/gcore/gdal_rat.h b/gcore/gdal_rat.h index adb89455ea9e..e19449f94057 100644 --- a/gcore/gdal_rat.h +++ b/gcore/gdal_rat.h @@ -386,4 +386,7 @@ class CPL_DLL GDALDefaultRasterAttributeTable : public GDALRasterAttributeTable void RemoveStatistics() override; }; +std::unique_ptr + CPL_DLL GDALLoadVATDBF(const char *pszFilename); + #endif /* ndef GDAL_RAT_H_INCLUDED */ diff --git a/gcore/gdal_rat_vat_dbf.cpp b/gcore/gdal_rat_vat_dbf.cpp new file mode 100644 index 000000000000..e2bfd3090186 --- /dev/null +++ b/gcore/gdal_rat_vat_dbf.cpp @@ -0,0 +1,141 @@ +/****************************************************************************** + * + * Project: GDAL Core + * Purpose: Read ArcGIS .vat.dbf raster attribute table + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2024, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#include "gdal_priv.h" +#include "gdal_rat.h" +#include "ogrsf_frmts.h" + +#include + +/************************************************************************/ +/* GDALLoadVATDBF() */ +/************************************************************************/ + +/** + * \brief Load a ESRI .vat.dbf auxiliary file as a GDAL attribute table. + * + * @since GDAL 3.11 + */ +std::unique_ptr +GDALLoadVATDBF(const char *pszFilename) +{ + auto poDS = std::unique_ptr( + 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(); + + 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( + 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; +}