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