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