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