Line data    Source code

       1             : // SPDX-License-Identifier: MIT
       2             : // Copyright 2024, Even Rouault <even.rouault at spatialys.com>
       3             : 
       4             : #include "cpl_port.h"
       5             : #include "cpl_minixml.h"
       6             : #include "cpl_vsi_virtual.h"
       7             : #include "gdal_pam.h"
       8             : #include "rawdataset.h"
       9             : 
      10             : #define LIBERTIFF_NS GDAL_libertiff
      11             : #include "libertiff.hpp"
      12             : 
      13             : #include <algorithm>
      14             : #include <cmath>
      15             : 
      16             : constexpr const char *SNAP_TIFF_PREFIX = "SNAP_TIFF:";
      17             : 
      18             : // Non-standard TIFF tag holding DIMAP XML for SNAP TIFF products
      19             : constexpr uint16_t DIMAP_TAG = 65000;
      20             : 
      21             : // Number of values per GCP in the GeoTIFFTiePoints tag
      22             : constexpr int VALUES_PER_GCP = 6;
      23             : 
      24             : constexpr int GCP_PIXEL = 0;
      25             : constexpr int GCP_LINE = 1;
      26             : // constexpr int GCP_DEPTH = 2;
      27             : constexpr int GCP_X = 3;
      28             : constexpr int GCP_Y = 4;
      29             : constexpr int GCP_Z = 5;
      30             : 
      31             : /************************************************************************/
      32             : /*                            SNAPTIFFDataset                           */
      33             : /************************************************************************/
      34             : 
      35             : class SNAPTIFFDataset final : public GDALPamDataset
      36             : {
      37             :   public:
      38           3 :     SNAPTIFFDataset() = default;
      39             : 
      40             :     static int Identify(GDALOpenInfo *poOpenInfo);
      41             :     static GDALDataset *Open(GDALOpenInfo *poOpenInfo);
      42             : 
      43             :     char **GetMetadataDomainList() override;
      44             :     char **GetMetadata(const char *pszDomain = "") override;
      45             :     const char *GetMetadataItem(const char *pszName,
      46             :                                 const char *pszDomain = "") override;
      47             : 
      48           1 :     const OGRSpatialReference *GetGCPSpatialRef() const override
      49             :     {
      50           1 :         return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
      51             :     }
      52             : 
      53           2 :     int GetGCPCount() override
      54             :     {
      55           2 :         return int(m_asGCPs.size());
      56             :     }
      57             : 
      58           1 :     const GDAL_GCP *GetGCPs() override
      59             :     {
      60           1 :         return gdal::GCP::c_ptr(m_asGCPs);
      61             :     }
      62             : 
      63             :   private:
      64             :     VSIVirtualHandleUniquePtr m_poFile{};
      65             :     std::unique_ptr<const LIBERTIFF_NS::Image> m_poImage{};
      66             : 
      67             :     //! Whether this dataset is actually the geolocation array
      68             :     bool m_bIsGeolocArray = false;
      69             : 
      70             :     //! Content of "xml:DIMAP" metadata domain
      71             :     CPLStringList m_aosDIMAPMetadata{};
      72             : 
      73             :     CPLStringList m_aosGEOLOCATION{};
      74             :     int m_nGeolocArrayWidth = 0;
      75             :     int m_nGeolocArrayHeight = 0;
      76             : 
      77             :     CPLStringList m_aosSUBDATASETS{};
      78             : 
      79             :     std::vector<gdal::GCP> m_asGCPs{};
      80             :     OGRSpatialReference m_oSRS{};
      81             : 
      82             :     bool GetGeolocationMetadata();
      83             : 
      84             :     void ReadSRS();
      85             : };
      86             : 
      87             : /************************************************************************/
      88             : /*                              Identify()                              */
      89             : /************************************************************************/
      90             : 
      91       84112 : int SNAPTIFFDataset::Identify(GDALOpenInfo *poOpenInfo)
      92             : {
      93       84112 :     if (STARTS_WITH(poOpenInfo->pszFilename, SNAP_TIFF_PREFIX))
      94          12 :         return true;
      95             : 
      96       84100 :     if (poOpenInfo->fpL == nullptr || poOpenInfo->nHeaderBytes < 16 ||
      97             :         // BigEndian classic TIFF
      98       29676 :         memcmp(poOpenInfo->pabyHeader, "\x4D\x4D\x00\x2A", 4) != 0)
      99             :     {
     100       83968 :         return false;
     101             :     }
     102             : 
     103             :     struct MyFileReader final : public LIBERTIFF_NS::FileReader
     104             :     {
     105             :         const GDALOpenInfo *m_poOpenInfo;
     106             : 
     107         132 :         explicit MyFileReader(const GDALOpenInfo *poOpenInfoIn)
     108         132 :             : m_poOpenInfo(poOpenInfoIn)
     109             :         {
     110         132 :         }
     111             : 
     112         136 :         uint64_t size() const override
     113             :         {
     114         136 :             return m_poOpenInfo->nHeaderBytes;
     115             :         }
     116             : 
     117        6875 :         size_t read(uint64_t offset, size_t count, void *buffer) const override
     118             :         {
     119        6875 :             if (offset + count >
     120        6875 :                 static_cast<uint64_t>(m_poOpenInfo->nHeaderBytes))
     121         178 :                 return 0;
     122        6697 :             memcpy(buffer,
     123        6697 :                    m_poOpenInfo->pabyHeader + static_cast<size_t>(offset),
     124             :                    count);
     125        6697 :             return count;
     126             :         }
     127             : 
     128             :         CPL_DISALLOW_COPY_ASSIGN(MyFileReader)
     129             :     };
     130             : 
     131         264 :     auto f = std::make_shared<const MyFileReader>(poOpenInfo);
     132             : #ifdef DEBUG
     133             :     // Just to increase coverage testing
     134         132 :     CPLAssert(f->size() == uint64_t(poOpenInfo->nHeaderBytes));
     135         132 :     char dummy = 0;
     136         132 :     CPLAssert(f->read(poOpenInfo->nHeaderBytes, 1, &dummy) == 0);
     137         132 :     CPL_IGNORE_RET_VAL(dummy);
     138             : #endif
     139         264 :     auto image = LIBERTIFF_NS::open</*acceptBigTIFF = */ false>(std::move(f));
     140             :     // Checks that it is a single-band Float32 uncompressed dataset, made
     141             :     // of a single strip.
     142         339 :     if (!image || image->nextImageOffset() != 0 ||
     143         185 :         image->compression() != LIBERTIFF_NS::Compression::None ||
     144          98 :         image->sampleFormat() != LIBERTIFF_NS::SampleFormat::IEEEFP ||
     145          23 :         image->samplesPerPixel() != 1 || image->bitsPerSample() != 32 ||
     146          18 :         image->isTiled() || image->strileCount() != 1 || image->width() == 0 ||
     147          12 :         image->width() > INT_MAX / sizeof(float) || image->height() == 0 ||
     148          12 :         image->height() > INT_MAX || image->rowsPerStrip() != image->height() ||
     149          12 :         !image->tag(LIBERTIFF_NS::TagCode::GeoTIFFPixelScale) ||
     150          12 :         !image->tag(LIBERTIFF_NS::TagCode::GeoTIFFTiePoints) ||
     151         249 :         !image->tag(LIBERTIFF_NS::TagCode::GeoTIFFGeoKeyDirectory) ||
     152           6 :         !image->tag(DIMAP_TAG))
     153             :     {
     154         128 :         return false;
     155             :     }
     156             : 
     157           4 :     return true;
     158             : }
     159             : 
     160             : /************************************************************************/
     161             : /*                               Open()                                 */
     162             : /************************************************************************/
     163             : 
     164           7 : GDALDataset *SNAPTIFFDataset::Open(GDALOpenInfo *poOpenInfo)
     165             : {
     166           7 :     if (poOpenInfo->eAccess == GA_Update || !Identify(poOpenInfo))
     167           0 :         return nullptr;
     168             : 
     169           7 :     bool bIsGeolocation = false;
     170             :     // Check if it is SNAP_TIFF:"filename":{subdataset_component} syntax
     171           7 :     if (STARTS_WITH(poOpenInfo->pszFilename, SNAP_TIFF_PREFIX))
     172             :     {
     173             :         const CPLStringList aosTokens(CSLTokenizeString2(
     174           6 :             poOpenInfo->pszFilename, ":", CSLT_HONOURSTRINGS));
     175           6 :         if (aosTokens.size() != 3)
     176           1 :             return nullptr;
     177           5 :         bIsGeolocation = EQUAL(aosTokens[2], "GEOLOCATION");
     178           5 :         if (!bIsGeolocation && !EQUAL(aosTokens[2], "MAIN"))
     179           1 :             return nullptr;
     180           4 :         GDALOpenInfo oOpenInfo(aosTokens[1], GA_ReadOnly);
     181           4 :         if (!Identify(&oOpenInfo))
     182           2 :             return nullptr;
     183           2 :         std::swap(poOpenInfo->fpL, oOpenInfo.fpL);
     184           2 :         if (!poOpenInfo->fpL)
     185           0 :             return nullptr;
     186             :     }
     187             : 
     188             :     struct MyFileReader final : public LIBERTIFF_NS::FileReader
     189             :     {
     190             :         VSILFILE *m_fp;
     191             :         mutable uint64_t m_nFileSize = 0;
     192             : 
     193           3 :         explicit MyFileReader(VSILFILE *fp) : m_fp(fp)
     194             :         {
     195           3 :         }
     196             : 
     197           9 :         uint64_t size() const override
     198             :         {
     199           9 :             if (m_nFileSize == 0)
     200             :             {
     201           3 :                 m_fp->Seek(0, SEEK_END);
     202           3 :                 m_nFileSize = m_fp->Tell();
     203             :             }
     204           9 :             return m_nFileSize;
     205             :         }
     206             : 
     207         261 :         size_t read(uint64_t offset, size_t count, void *buffer) const override
     208             :         {
     209         261 :             return m_fp->Seek(offset, SEEK_SET) == 0
     210         261 :                        ? m_fp->Read(buffer, 1, count)
     211         261 :                        : 0;
     212             :         }
     213             : 
     214             :         CPL_DISALLOW_COPY_ASSIGN(MyFileReader)
     215             :     };
     216             : 
     217           6 :     auto f = std::make_shared<const MyFileReader>(poOpenInfo->fpL);
     218             : #ifdef DEBUG
     219             :     // Just to increase coverage testing
     220           3 :     char dummy = 0;
     221           3 :     CPLAssert(f->read(f->size(), 1, &dummy) == 0);
     222           3 :     CPL_IGNORE_RET_VAL(dummy);
     223             : #endif
     224             : 
     225           6 :     auto poDS = std::make_unique<SNAPTIFFDataset>();
     226             : 
     227           3 :     poDS->m_poImage = LIBERTIFF_NS::open(std::move(f));
     228           3 :     if (!poDS->m_poImage)
     229           0 :         return nullptr;
     230           3 :     poDS->m_poFile.reset(poOpenInfo->fpL);
     231           3 :     poOpenInfo->fpL = nullptr;
     232             : 
     233           3 :     poDS->nRasterXSize = static_cast<int>(poDS->m_poImage->width());
     234           3 :     poDS->nRasterYSize = static_cast<int>(poDS->m_poImage->height());
     235           3 :     poDS->SetDescription(poOpenInfo->pszFilename);
     236             : 
     237           3 :     if (bIsGeolocation)
     238             :     {
     239           1 :         poDS->m_bIsGeolocArray = true;
     240           1 :         if (!poDS->GetGeolocationMetadata())
     241             :         {
     242           0 :             return nullptr;
     243             :         }
     244           1 :         poDS->nRasterXSize = poDS->m_nGeolocArrayWidth;
     245           1 :         poDS->nRasterYSize = poDS->m_nGeolocArrayHeight;
     246             : 
     247             :         const auto *psTag =
     248           1 :             poDS->m_poImage->tag(LIBERTIFF_NS::TagCode::GeoTIFFTiePoints);
     249             : 
     250           3 :         for (int iBand = 1; iBand <= 2; ++iBand)
     251             :         {
     252             :             auto poBand = std::make_unique<RawRasterBand>(
     253           2 :                 poDS->m_poFile.get(),
     254           2 :                 psTag->value_offset +
     255           2 :                     (iBand == 1 ? GCP_X : GCP_Y) * sizeof(double),
     256           0 :                 int(sizeof(double)) * VALUES_PER_GCP,
     257           2 :                 int(sizeof(double)) * VALUES_PER_GCP * poDS->nRasterXSize,
     258           2 :                 GDT_Float64, !poDS->m_poImage->mustByteSwap(),
     259           2 :                 poDS->nRasterXSize, poDS->nRasterYSize,
     260           4 :                 RawRasterBand::OwnFP::NO);
     261           2 :             if (!poBand->IsValid())
     262           0 :                 return nullptr;
     263           2 :             if (iBand == 1)
     264           1 :                 poBand->SetDescription("Longitude");
     265             :             else
     266           1 :                 poBand->SetDescription("Latitude");
     267           2 :             poDS->SetBand(iBand, std::move(poBand));
     268             :         }
     269             : 
     270           1 :         return poDS.release();
     271             :     }
     272             : 
     273           2 :     poDS->ReadSRS();
     274           2 :     poDS->GetGeolocationMetadata();
     275             : 
     276             :     {
     277           2 :         bool okStrileOffset = false;
     278             :         auto poBand = std::make_unique<RawRasterBand>(
     279           2 :             poDS->m_poFile.get(),
     280           2 :             poDS->m_poImage->strileOffset(0, okStrileOffset),
     281           2 :             int(sizeof(float)), int(sizeof(float)) * poDS->nRasterXSize,
     282           2 :             GDT_Float32, !poDS->m_poImage->mustByteSwap(), poDS->nRasterXSize,
     283           4 :             poDS->nRasterYSize, RawRasterBand::OwnFP::NO);
     284           2 :         if (!poBand->IsValid())
     285           0 :             return nullptr;
     286           2 :         poDS->SetBand(1, std::move(poBand));
     287             :     }
     288           2 :     auto poBand = poDS->papoBands[0];
     289             : 
     290             :     const auto &psImageDescription =
     291           2 :         poDS->m_poImage->tag(LIBERTIFF_NS::TagCode::ImageDescription);
     292           2 :     if (psImageDescription &&
     293           2 :         psImageDescription->type == LIBERTIFF_NS::TagType::ASCII &&
     294           2 :         !psImageDescription->invalid_value_offset &&
     295             :         // Sanity check
     296           2 :         psImageDescription->count < 100 * 1000)
     297             :     {
     298           2 :         bool ok = true;
     299             :         const std::string s =
     300           4 :             poDS->m_poImage->readTagAsString(*psImageDescription, ok);
     301           2 :         if (ok)
     302             :         {
     303           2 :             poDS->GDALDataset::SetMetadataItem("IMAGE_DESCRIPTION", s.c_str());
     304             :         }
     305             :     }
     306             : 
     307           2 :     const auto *psDimapTag = poDS->m_poImage->tag(DIMAP_TAG);
     308           2 :     if (psDimapTag && psDimapTag->type == LIBERTIFF_NS::TagType::ASCII &&
     309           2 :         !psDimapTag->invalid_value_offset)
     310             :     {
     311             :         try
     312             :         {
     313           2 :             bool ok = true;
     314             :             // Just read the first 10 kB max to fetch essential band metadata
     315           2 :             const std::string s = poDS->m_poImage->readContext()->readString(
     316           2 :                 psDimapTag->value_offset,
     317             :                 static_cast<size_t>(
     318           2 :                     std::min<uint64_t>(psDimapTag->count, 10000)),
     319           6 :                 ok);
     320           2 :             const auto posStart = s.find("<Spectral_Band_Info>");
     321           2 :             if (posStart != std::string::npos)
     322             :             {
     323           2 :                 const char *markerEnd = "</Spectral_Band_Info>";
     324           2 :                 const auto posEnd = s.find(markerEnd, posStart);
     325           2 :                 if (posEnd != std::string::npos)
     326             :                 {
     327             :                     const std::string osSubstr = s.substr(
     328           4 :                         posStart, posEnd - posStart + strlen(markerEnd));
     329           4 :                     CPLXMLTreeCloser oRoot(CPLParseXMLString(osSubstr.c_str()));
     330           2 :                     if (oRoot)
     331             :                     {
     332           2 :                         const char *pszNoDataValueUsed = CPLGetXMLValue(
     333           2 :                             oRoot.get(), "NO_DATA_VALUE_USED", nullptr);
     334           2 :                         const char *pszNoDataValue = CPLGetXMLValue(
     335           2 :                             oRoot.get(), "NO_DATA_VALUE", nullptr);
     336           4 :                         if (pszNoDataValueUsed && pszNoDataValue &&
     337           2 :                             CPLTestBool(pszNoDataValueUsed))
     338             :                         {
     339           2 :                             poBand->SetNoDataValue(CPLAtof(pszNoDataValue));
     340             :                         }
     341             : 
     342           2 :                         if (const char *pszScalingFactor = CPLGetXMLValue(
     343           2 :                                 oRoot.get(), "SCALING_FACTOR", nullptr))
     344             :                         {
     345           2 :                             poBand->SetScale(CPLAtof(pszScalingFactor));
     346             :                         }
     347             : 
     348           2 :                         if (const char *pszScalingOffset = CPLGetXMLValue(
     349           2 :                                 oRoot.get(), "SCALING_OFFSET", nullptr))
     350             :                         {
     351           2 :                             poBand->SetOffset(CPLAtof(pszScalingOffset));
     352             :                         }
     353             : 
     354           2 :                         if (const char *pszBandName = CPLGetXMLValue(
     355           2 :                                 oRoot.get(), "BAND_NAME", nullptr))
     356             :                         {
     357           2 :                             poBand->SetDescription(pszBandName);
     358             :                         }
     359             : 
     360           2 :                         if (const char *pszUnit = CPLGetXMLValue(
     361           2 :                                 oRoot.get(), "PHYSICAL_UNIT", nullptr))
     362             :                         {
     363           2 :                             poBand->SetUnitType(pszUnit);
     364             :                         }
     365             :                     }
     366             :                 }
     367             :             }
     368             :         }
     369           0 :         catch (const std::exception &)
     370             :         {
     371             :         }
     372             :     }
     373             : 
     374             :     // Initialize PAM
     375           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     376           2 :     poDS->TryLoadXML();
     377             : 
     378             :     // Check for overviews.
     379           2 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     380             : 
     381           2 :     return poDS.release();
     382             : }
     383             : 
     384             : /************************************************************************/
     385             : /*                      GetMetadataDomainList()                         */
     386             : /************************************************************************/
     387             : 
     388           1 : char **SNAPTIFFDataset::GetMetadataDomainList()
     389             : {
     390           1 :     return BuildMetadataDomainList(GDALPamDataset::GetMetadataDomainList(),
     391             :                                    TRUE, "GEOLOCATION", "SUBDATASETS",
     392           1 :                                    "xml:DIMAP", nullptr);
     393             : }
     394             : 
     395             : /************************************************************************/
     396             : /*                      GetGeolocationMetadata()                        */
     397             : /************************************************************************/
     398             : 
     399             : // (Partially) read the content of the GeoTIFFTiePoints tag to check if the
     400             : // tie points form a regular geolocation array, and extract the width, height,
     401             : // and spacing of that geolocation array. Also fills the m_aosGEOLOCATION
     402             : // metadata domain
     403           4 : bool SNAPTIFFDataset::GetGeolocationMetadata()
     404             : {
     405           4 :     const auto *psTag = m_poImage->tag(LIBERTIFF_NS::TagCode::GeoTIFFTiePoints);
     406           4 :     if (psTag && psTag->type == LIBERTIFF_NS::TagType::Double &&
     407           4 :         !psTag->invalid_value_offset && (psTag->count % VALUES_PER_GCP) == 0 &&
     408             :         // Sanity check
     409           4 :         psTag->count <= uint64_t(nRasterXSize) * nRasterYSize * VALUES_PER_GCP)
     410             :     {
     411           4 :         const int numGCPs = int(psTag->count / VALUES_PER_GCP);
     412             : 
     413             :         // Assume that the geolocation array has the same proportion as the
     414             :         // main array
     415             :         const double dfGCPArrayWidth =
     416           4 :             sqrt(double(nRasterXSize) * numGCPs / nRasterYSize);
     417             :         const double dfGCPArrayHeight =
     418           4 :             sqrt(double(nRasterYSize) * numGCPs / nRasterXSize);
     419           4 :         if (dfGCPArrayWidth > INT_MAX || dfGCPArrayHeight > INT_MAX)
     420             :         {
     421           0 :             return false;
     422             :         }
     423             : 
     424           4 :         const int nGCPArrayWidth = int(std::round(dfGCPArrayWidth));
     425           4 :         const int nGCPArrayHeight = int(std::round(dfGCPArrayHeight));
     426           4 :         constexpr int NUM_LINES = 3;
     427           4 :         if (nGCPArrayWidth * nGCPArrayHeight != numGCPs ||
     428             :             nGCPArrayHeight < NUM_LINES)
     429             :         {
     430           0 :             return false;
     431             :         }
     432             : 
     433           4 :         bool ok = true;
     434             :         // Just read the first 3 lines of the geolocation array
     435           4 :         const auto numValuesPerLine = nGCPArrayWidth * VALUES_PER_GCP;
     436           4 :         const auto values = m_poImage->readContext()->readArray<double>(
     437           8 :             psTag->value_offset, numValuesPerLine * NUM_LINES, ok);
     438           4 :         if (!ok || values.empty())
     439           0 :             return false;
     440             : 
     441           4 :         if (values[GCP_LINE] != 0.5 && values[GCP_PIXEL] != 0.5)
     442           0 :             return false;
     443             : 
     444           4 :         constexpr double RELATIVE_TOLERANCE = 1e-5;
     445           4 :         constexpr double PIXEL_TOLERANCE = 1e-3;
     446             : 
     447             :         // Check that pixel spacing is constant on all three first lines
     448             :         const double pixelSpacing =
     449           4 :             values[VALUES_PER_GCP + GCP_PIXEL] - values[GCP_PIXEL];
     450           4 :         if (!(pixelSpacing >= 1))
     451           0 :             return false;
     452           4 :         if (std::fabs(pixelSpacing * (nGCPArrayWidth - 1) -
     453           4 :                       (nRasterXSize - 1)) > PIXEL_TOLERANCE)
     454           0 :             return false;
     455             : 
     456             :         double adfY[NUM_LINES];
     457          16 :         for (int iLine = 0; iLine < NUM_LINES; ++iLine)
     458             :         {
     459          12 :             adfY[iLine] = values[iLine * numValuesPerLine + GCP_LINE];
     460          12 :             for (int i = iLine * numValuesPerLine + VALUES_PER_GCP;
     461       19140 :                  i < (iLine + 1) * numValuesPerLine; i += VALUES_PER_GCP)
     462             :             {
     463       38256 :                 if (values[i + GCP_LINE] !=
     464       19128 :                     values[i - VALUES_PER_GCP + GCP_LINE])
     465             :                 {
     466           0 :                     return false;
     467             :                 }
     468             :                 const double pixelSpacingNew =
     469       19128 :                     values[i + GCP_PIXEL] -
     470       19128 :                     values[i - VALUES_PER_GCP + GCP_PIXEL];
     471       19128 :                 if (std::fabs(pixelSpacingNew - pixelSpacing) >
     472       19128 :                     RELATIVE_TOLERANCE * std::fabs(pixelSpacing))
     473             :                 {
     474           0 :                     return false;
     475             :                 }
     476             :             }
     477             :         }
     478             : 
     479             :         // Check that line spacing is constant on the three first lines
     480           4 :         const double lineSpacing = adfY[1] - adfY[0];
     481           4 :         if (!(lineSpacing >= 1))
     482           0 :             return false;
     483           4 :         if (std::fabs(lineSpacing * (nGCPArrayHeight - 1) -
     484           4 :                       (nRasterYSize - 1)) > PIXEL_TOLERANCE)
     485           0 :             return false;
     486             : 
     487           8 :         for (int iLine = 1; iLine + 1 < NUM_LINES; ++iLine)
     488             :         {
     489           4 :             const double lineSpacingNew = adfY[iLine + 1] - adfY[iLine];
     490           4 :             if (std::fabs(lineSpacingNew - lineSpacing) >
     491           4 :                 RELATIVE_TOLERANCE * std::fabs(lineSpacing))
     492             :             {
     493           0 :                 return false;
     494             :             }
     495             :         }
     496             : 
     497             :         // Read last line
     498           4 :         const auto lastLineValue = m_poImage->readContext()->readArray<double>(
     499           4 :             psTag->value_offset + uint64_t(nGCPArrayHeight - 1) *
     500           4 :                                       numValuesPerLine * sizeof(double),
     501           8 :             numValuesPerLine, ok);
     502           4 :         if (!ok)
     503           0 :             return false;
     504             : 
     505           4 :         if (!m_bIsGeolocArray)
     506             :         {
     507             :             // Expose the 4 corner GCPs for rough georeferencing.
     508             : 
     509           3 :             const uint32_t nShift = numValuesPerLine - VALUES_PER_GCP;
     510           3 :             m_asGCPs.emplace_back("TL", "Top Left", values[GCP_PIXEL],
     511           3 :                                   values[GCP_LINE], values[GCP_X],
     512           6 :                                   values[GCP_Y], values[GCP_Z]);
     513             :             m_asGCPs.emplace_back(
     514           3 :                 "TR", "Top Right", values[nShift + GCP_PIXEL],
     515           3 :                 values[nShift + GCP_LINE], values[nShift + GCP_X],
     516           6 :                 values[nShift + GCP_Y], values[nShift + GCP_Z]);
     517           3 :             m_asGCPs.emplace_back("BL", "Bottom Left", lastLineValue[GCP_PIXEL],
     518           3 :                                   lastLineValue[GCP_LINE], lastLineValue[GCP_X],
     519           6 :                                   lastLineValue[GCP_Y], lastLineValue[GCP_Z]);
     520             :             m_asGCPs.emplace_back(
     521           3 :                 "BR", "Bottom Right", lastLineValue[nShift + GCP_PIXEL],
     522           3 :                 lastLineValue[nShift + GCP_LINE], lastLineValue[nShift + GCP_X],
     523           6 :                 lastLineValue[nShift + GCP_Y], lastLineValue[nShift + GCP_Z]);
     524             :         }
     525             : 
     526           4 :         m_nGeolocArrayWidth = nGCPArrayWidth;
     527           4 :         m_nGeolocArrayHeight = nGCPArrayHeight;
     528             : 
     529           4 :         if (!m_bIsGeolocArray)
     530             :         {
     531           3 :             if (!m_oSRS.IsEmpty())
     532             :             {
     533           3 :                 char *pszWKT = nullptr;
     534           3 :                 m_oSRS.exportToWkt(&pszWKT);
     535           3 :                 m_aosGEOLOCATION.SetNameValue("SRS", pszWKT);
     536           3 :                 CPLFree(pszWKT);
     537             :             }
     538             : 
     539             :             m_aosGEOLOCATION.SetNameValue(
     540             :                 "X_DATASET", CPLSPrintf("%s\"%s\":GEOLOCATION",
     541           3 :                                         SNAP_TIFF_PREFIX, GetDescription()));
     542           3 :             m_aosGEOLOCATION.SetNameValue("X_BAND", "1");
     543             :             m_aosGEOLOCATION.SetNameValue(
     544             :                 "Y_DATASET", CPLSPrintf("%s\"%s\":GEOLOCATION",
     545           3 :                                         SNAP_TIFF_PREFIX, GetDescription()));
     546           3 :             m_aosGEOLOCATION.SetNameValue("Y_BAND", "2");
     547           3 :             m_aosGEOLOCATION.SetNameValue("PIXEL_OFFSET", "0");
     548             :             m_aosGEOLOCATION.SetNameValue("PIXEL_STEP",
     549           3 :                                           CPLSPrintf("%.17g", pixelSpacing));
     550           3 :             m_aosGEOLOCATION.SetNameValue("LINE_OFFSET", "0");
     551             :             m_aosGEOLOCATION.SetNameValue("LINE_STEP",
     552           3 :                                           CPLSPrintf("%.17g", lineSpacing));
     553             :         }
     554             : 
     555           4 :         return true;
     556             :     }
     557           0 :     return false;
     558             : }
     559             : 
     560             : /************************************************************************/
     561             : /*                             ReadSRS()                                */
     562             : /************************************************************************/
     563             : 
     564             : // Simplified GeoTIFF SRS reader, assuming the SRS is encoded as a EPSG code
     565           2 : void SNAPTIFFDataset::ReadSRS()
     566             : {
     567             :     const auto &psGeoKeysTag =
     568           2 :         m_poImage->tag(LIBERTIFF_NS::TagCode::GeoTIFFGeoKeyDirectory);
     569           2 :     constexpr int VALUES_PER_GEOKEY = 4;
     570           2 :     if (psGeoKeysTag && psGeoKeysTag->type == LIBERTIFF_NS::TagType::Short &&
     571           2 :         !psGeoKeysTag->invalid_value_offset &&
     572           2 :         psGeoKeysTag->count >= VALUES_PER_GEOKEY &&
     573           2 :         (psGeoKeysTag->count % VALUES_PER_GEOKEY) == 0 &&
     574             :         // Sanity check
     575           2 :         psGeoKeysTag->count < 1000)
     576             :     {
     577           2 :         bool ok = true;
     578             :         const auto values =
     579           4 :             m_poImage->readTagAsVector<uint16_t>(*psGeoKeysTag, ok);
     580           2 :         if (values.size() >= 4)
     581             :         {
     582           2 :             const uint16_t geokeysCount = values[3];
     583           2 :             constexpr uint16_t GEOTIFF_KEY_DIRECTORY_VERSION_V1 = 1;
     584           2 :             constexpr uint16_t GEOTIFF_KEY_VERSION_MAJOR_V1 = 1;
     585           2 :             if (values[0] == GEOTIFF_KEY_DIRECTORY_VERSION_V1 &&
     586             :                 // GeoTIFF 1.x
     587           4 :                 values[1] == GEOTIFF_KEY_VERSION_MAJOR_V1 &&
     588           2 :                 geokeysCount == psGeoKeysTag->count / VALUES_PER_GEOKEY - 1)
     589             :             {
     590           2 :                 constexpr uint16_t GeoTIFFTypeShort = 0;
     591           2 :                 constexpr uint16_t GeodeticCRSGeoKey = 2048;
     592           2 :                 constexpr uint16_t ProjectedCRSGeoKey = 3072;
     593           2 :                 uint16_t nEPSGCode = 0;
     594           8 :                 for (uint32_t i = 1; i <= geokeysCount; ++i)
     595             :                 {
     596           6 :                     const auto geokey = values[VALUES_PER_GEOKEY * i];
     597           6 :                     const auto geokeyType = values[VALUES_PER_GEOKEY * i + 1];
     598           6 :                     const auto geokeyCount = values[VALUES_PER_GEOKEY * i + 2];
     599           6 :                     const auto geokeyValue = values[VALUES_PER_GEOKEY * i + 3];
     600           6 :                     if ((geokey == GeodeticCRSGeoKey ||
     601           2 :                          geokey == ProjectedCRSGeoKey) &&
     602           2 :                         geokeyType == GeoTIFFTypeShort && geokeyCount == 1 &&
     603             :                         geokeyValue > 0)
     604             :                     {
     605           2 :                         nEPSGCode = geokeyValue;
     606           2 :                         if (geokey == ProjectedCRSGeoKey)
     607           0 :                             break;
     608             :                     }
     609             :                 }
     610           2 :                 if (nEPSGCode > 0)
     611             :                 {
     612           2 :                     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     613           2 :                     m_oSRS.importFromEPSG(nEPSGCode);
     614             :                 }
     615             :             }
     616             :         }
     617             :     }
     618           2 : }
     619             : 
     620             : /************************************************************************/
     621             : /*                            GetMetadata()                             */
     622             : /************************************************************************/
     623             : 
     624          11 : char **SNAPTIFFDataset::GetMetadata(const char *pszDomain)
     625             : {
     626          11 :     if (!m_bIsGeolocArray)
     627             :     {
     628          10 :         if (pszDomain && EQUAL(pszDomain, "xml:DIMAP"))
     629             :         {
     630           2 :             if (m_aosDIMAPMetadata.empty())
     631             :             {
     632           1 :                 const auto *psDimapTag = m_poImage->tag(DIMAP_TAG);
     633           1 :                 if (psDimapTag &&
     634           1 :                     psDimapTag->type == LIBERTIFF_NS::TagType::ASCII &&
     635           1 :                     !psDimapTag->invalid_value_offset &&
     636             :                     // Sanity check
     637           1 :                     psDimapTag->count < 100 * 1000 * 1000)
     638             :                 {
     639             :                     try
     640             :                     {
     641           1 :                         bool ok = true;
     642             :                         const std::string s =
     643           2 :                             m_poImage->readTagAsString(*psDimapTag, ok);
     644           1 :                         if (ok)
     645             :                         {
     646           1 :                             m_aosDIMAPMetadata.AddString(s.c_str());
     647             :                         }
     648             :                     }
     649           0 :                     catch (const std::exception &)
     650             :                     {
     651           0 :                         CPLError(CE_Failure, CPLE_OutOfMemory,
     652             :                                  "Out of memory in GetMetadata()");
     653             :                     }
     654             :                 }
     655             :             }
     656           2 :             return m_aosDIMAPMetadata.List();
     657             :         }
     658           8 :         else if (pszDomain && EQUAL(pszDomain, "GEOLOCATION"))
     659             :         {
     660           3 :             return m_aosGEOLOCATION.List();
     661             :         }
     662           5 :         else if (pszDomain && EQUAL(pszDomain, "SUBDATASETS"))
     663             :         {
     664           4 :             if (m_aosSUBDATASETS.empty() && GetGeolocationMetadata())
     665             :             {
     666             :                 m_aosSUBDATASETS.SetNameValue("SUBDATASET_1_NAME",
     667             :                                               CPLSPrintf("%s\"%s\":MAIN",
     668             :                                                          SNAP_TIFF_PREFIX,
     669           1 :                                                          GetDescription()));
     670             :                 m_aosSUBDATASETS.SetNameValue("SUBDATASET_1_DESC",
     671           2 :                                               std::string("Main content of ")
     672           1 :                                                   .append(GetDescription())
     673           2 :                                                   .c_str());
     674             : 
     675             :                 m_aosSUBDATASETS.SetNameValue(
     676             :                     "SUBDATASET_2_NAME",
     677           1 :                     m_aosGEOLOCATION.FetchNameValue("X_DATASET"));
     678             :                 m_aosSUBDATASETS.SetNameValue(
     679           2 :                     "SUBDATASET_2_DESC", std::string("Geolocation array of ")
     680           1 :                                              .append(GetDescription())
     681           2 :                                              .c_str());
     682             :             }
     683           4 :             return m_aosSUBDATASETS.List();
     684             :         }
     685             :     }
     686             : 
     687           2 :     return GDALPamDataset::GetMetadata(pszDomain);
     688             : }
     689             : 
     690             : /************************************************************************/
     691             : /*                         GetMetadataItem()                            */
     692             : /************************************************************************/
     693             : 
     694           4 : const char *SNAPTIFFDataset::GetMetadataItem(const char *pszName,
     695             :                                              const char *pszDomain)
     696             : {
     697           4 :     if (!m_bIsGeolocArray && pszDomain &&
     698           4 :         (EQUAL(pszDomain, "GEOLOCATION") || EQUAL(pszDomain, "SUBDATASETS")))
     699             :     {
     700           3 :         return CSLFetchNameValue(GetMetadata(pszDomain), pszName);
     701             :     }
     702           1 :     return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
     703             : }
     704             : 
     705             : /************************************************************************/
     706             : /*                     GDALRegister_SNAP_TIFF()                         */
     707             : /************************************************************************/
     708             : 
     709        1889 : void GDALRegister_SNAP_TIFF()
     710             : {
     711        1889 :     if (GDALGetDriverByName("SNAP_TIFF") != nullptr)
     712         282 :         return;
     713             : 
     714        1607 :     GDALDriver *poDriver = new GDALDriver();
     715        1607 :     poDriver->SetDescription("SNAP_TIFF");
     716        1607 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     717        1607 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     718        1607 :                               "Sentinel Application Processing GeoTIFF");
     719        1607 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
     720        1607 :                               "drivers/raster/snap_tiff.html");
     721             :     // Declaring the tif extension confuses QGIS
     722             :     // Cf https://github.com/qgis/QGIS/issues/59112
     723             :     // This driver is of too marginal usage to justify causing chaos downstream.
     724             :     // poDriver->SetMetadataItem(GDAL_DMD_MIMETYPE, "image/tiff");
     725             :     // poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "tif tiff");
     726        1607 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     727             : 
     728        1607 :     poDriver->pfnOpen = SNAPTIFFDataset::Open;
     729        1607 :     poDriver->pfnIdentify = SNAPTIFFDataset::Identify;
     730             : 
     731        1607 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     732             : }