LCOV - code coverage report
Current view: top level - frmts/pds - pds4dataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 2087 2355 88.6 %
Date: 2025-07-09 17:50:03 Functions: 75 77 97.4 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  PDS 4 Driver; Planetary Data System Format
       4             :  * Purpose:  Implementation of PDS4Dataset
       5             :  * Author:   Even Rouault, even.rouault at spatialys.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2017, Hobu Inc
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_vsi_error.h"
      14             : #include "gdal_proxy.h"
      15             : #include "rawdataset.h"
      16             : #include "vrtdataset.h"
      17             : #include "ogrsf_frmts.h"
      18             : #include "ogr_spatialref.h"
      19             : #include "gdal_priv_templates.hpp"
      20             : #include "ogreditablelayer.h"
      21             : #include "pds4dataset.h"
      22             : #include "pdsdrivercore.h"
      23             : 
      24             : #ifdef EMBED_RESOURCE_FILES
      25             : #include "embedded_resources.h"
      26             : #endif
      27             : 
      28             : #include <cstdlib>
      29             : #include <vector>
      30             : #include <algorithm>
      31             : 
      32             : #define TIFF_GEOTIFF_STRING "TIFF 6.0"
      33             : #define BIGTIFF_GEOTIFF_STRING "TIFF 6.0"
      34             : #define PREEXISTING_BINARY_FILE                                                \
      35             :     "Binary file pre-existing PDS4 label. This comment is used by GDAL to "    \
      36             :     "avoid deleting the binary file when the label is deleted. Keep it to "    \
      37             :     "preserve this behavior."
      38             : 
      39             : #define CURRENT_CART_VERSION "1G00_1950"
      40             : 
      41             : extern "C" void GDALRegister_PDS4();
      42             : 
      43             : /************************************************************************/
      44             : /*                        PDS4WrapperRasterBand()                      */
      45             : /************************************************************************/
      46             : 
      47          17 : PDS4WrapperRasterBand::PDS4WrapperRasterBand(GDALRasterBand *poBaseBandIn)
      48             :     : m_poBaseBand(poBaseBandIn), m_bHasOffset(false), m_bHasScale(false),
      49          17 :       m_bHasNoData(false), m_dfOffset(0.0), m_dfScale(1.0), m_dfNoData(0.0)
      50             : {
      51          17 :     eDataType = m_poBaseBand->GetRasterDataType();
      52          17 :     m_poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
      53          17 : }
      54             : 
      55             : /************************************************************************/
      56             : /*                             SetMaskBand()                            */
      57             : /************************************************************************/
      58             : 
      59           0 : void PDS4WrapperRasterBand::SetMaskBand(
      60             :     std::unique_ptr<GDALRasterBand> poMaskBand)
      61             : {
      62           0 :     poMask.reset(std::move(poMaskBand));
      63           0 :     nMaskFlags = 0;
      64           0 : }
      65             : 
      66             : /************************************************************************/
      67             : /*                              GetOffset()                             */
      68             : /************************************************************************/
      69             : 
      70          11 : double PDS4WrapperRasterBand::GetOffset(int *pbSuccess)
      71             : {
      72          11 :     if (pbSuccess)
      73          11 :         *pbSuccess = m_bHasOffset;
      74          11 :     return m_dfOffset;
      75             : }
      76             : 
      77             : /************************************************************************/
      78             : /*                              GetScale()                              */
      79             : /************************************************************************/
      80             : 
      81          11 : double PDS4WrapperRasterBand::GetScale(int *pbSuccess)
      82             : {
      83          11 :     if (pbSuccess)
      84          11 :         *pbSuccess = m_bHasScale;
      85          11 :     return m_dfScale;
      86             : }
      87             : 
      88             : /************************************************************************/
      89             : /*                              SetOffset()                             */
      90             : /************************************************************************/
      91             : 
      92           1 : CPLErr PDS4WrapperRasterBand::SetOffset(double dfNewOffset)
      93             : {
      94           1 :     m_dfOffset = dfNewOffset;
      95           1 :     m_bHasOffset = true;
      96             : 
      97           1 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
      98           1 :     if (poGDS->m_poExternalDS && eAccess == GA_Update)
      99           1 :         poGDS->m_poExternalDS->GetRasterBand(nBand)->SetOffset(dfNewOffset);
     100             : 
     101           1 :     return CE_None;
     102             : }
     103             : 
     104             : /************************************************************************/
     105             : /*                              SetScale()                              */
     106             : /************************************************************************/
     107             : 
     108           1 : CPLErr PDS4WrapperRasterBand::SetScale(double dfNewScale)
     109             : {
     110           1 :     m_dfScale = dfNewScale;
     111           1 :     m_bHasScale = true;
     112             : 
     113           1 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     114           1 :     if (poGDS->m_poExternalDS && eAccess == GA_Update)
     115           1 :         poGDS->m_poExternalDS->GetRasterBand(nBand)->SetScale(dfNewScale);
     116             : 
     117           1 :     return CE_None;
     118             : }
     119             : 
     120             : /************************************************************************/
     121             : /*                           GetNoDataValue()                           */
     122             : /************************************************************************/
     123             : 
     124          22 : double PDS4WrapperRasterBand::GetNoDataValue(int *pbSuccess)
     125             : {
     126          22 :     if (pbSuccess)
     127          22 :         *pbSuccess = m_bHasNoData;
     128          22 :     return m_dfNoData;
     129             : }
     130             : 
     131             : /************************************************************************/
     132             : /*                           SetNoDataValue()                           */
     133             : /************************************************************************/
     134             : 
     135           4 : CPLErr PDS4WrapperRasterBand::SetNoDataValue(double dfNewNoData)
     136             : {
     137           4 :     m_dfNoData = dfNewNoData;
     138           4 :     m_bHasNoData = true;
     139             : 
     140           4 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     141           4 :     if (poGDS->m_poExternalDS && eAccess == GA_Update)
     142           4 :         poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValue(
     143           4 :             dfNewNoData);
     144             : 
     145           4 :     return CE_None;
     146             : }
     147             : 
     148             : /************************************************************************/
     149             : /*                               Fill()                                 */
     150             : /************************************************************************/
     151             : 
     152           2 : CPLErr PDS4WrapperRasterBand::Fill(double dfRealValue, double dfImaginaryValue)
     153             : {
     154           2 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     155           2 :     if (poGDS->m_bMustInitImageFile)
     156             :     {
     157           2 :         if (!poGDS->InitImageFile())
     158           0 :             return CE_Failure;
     159             :     }
     160           2 :     return GDALProxyRasterBand::Fill(dfRealValue, dfImaginaryValue);
     161             : }
     162             : 
     163             : /************************************************************************/
     164             : /*                             IWriteBlock()                             */
     165             : /************************************************************************/
     166             : 
     167           0 : CPLErr PDS4WrapperRasterBand::IWriteBlock(int nXBlock, int nYBlock,
     168             :                                           void *pImage)
     169             : 
     170             : {
     171           0 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     172           0 :     if (poGDS->m_bMustInitImageFile)
     173             :     {
     174           0 :         if (!poGDS->InitImageFile())
     175           0 :             return CE_Failure;
     176             :     }
     177           0 :     return GDALProxyRasterBand::IWriteBlock(nXBlock, nYBlock, pImage);
     178             : }
     179             : 
     180             : /************************************************************************/
     181             : /*                             IRasterIO()                              */
     182             : /************************************************************************/
     183             : 
     184         113 : CPLErr PDS4WrapperRasterBand::IRasterIO(
     185             :     GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
     186             :     void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
     187             :     GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
     188             : 
     189             : {
     190         113 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     191         113 :     if (eRWFlag == GF_Write && poGDS->m_bMustInitImageFile)
     192             :     {
     193           7 :         if (!poGDS->InitImageFile())
     194           0 :             return CE_Failure;
     195             :     }
     196         113 :     return GDALProxyRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
     197             :                                           pData, nBufXSize, nBufYSize, eBufType,
     198         113 :                                           nPixelSpace, nLineSpace, psExtraArg);
     199             : }
     200             : 
     201             : /************************************************************************/
     202             : /*                       PDS4RawRasterBand()                            */
     203             : /************************************************************************/
     204             : 
     205         422 : PDS4RawRasterBand::PDS4RawRasterBand(GDALDataset *l_poDS, int l_nBand,
     206             :                                      VSILFILE *l_fpRaw,
     207             :                                      vsi_l_offset l_nImgOffset,
     208             :                                      int l_nPixelOffset, int l_nLineOffset,
     209             :                                      GDALDataType l_eDataType,
     210         422 :                                      RawRasterBand::ByteOrder eByteOrderIn)
     211             :     : RawRasterBand(l_poDS, l_nBand, l_fpRaw, l_nImgOffset, l_nPixelOffset,
     212             :                     l_nLineOffset, l_eDataType, eByteOrderIn,
     213             :                     RawRasterBand::OwnFP::NO),
     214             :       m_bHasOffset(false), m_bHasScale(false), m_bHasNoData(false),
     215         422 :       m_dfOffset(0.0), m_dfScale(1.0), m_dfNoData(0.0)
     216             : {
     217         422 : }
     218             : 
     219             : /************************************************************************/
     220             : /*                             SetMaskBand()                            */
     221             : /************************************************************************/
     222             : 
     223          39 : void PDS4RawRasterBand::SetMaskBand(std::unique_ptr<GDALRasterBand> poMaskBand)
     224             : {
     225          39 :     poMask.reset(std::move(poMaskBand));
     226          39 :     nMaskFlags = 0;
     227          39 : }
     228             : 
     229             : /************************************************************************/
     230             : /*                              GetOffset()                             */
     231             : /************************************************************************/
     232             : 
     233         118 : double PDS4RawRasterBand::GetOffset(int *pbSuccess)
     234             : {
     235         118 :     if (pbSuccess)
     236         110 :         *pbSuccess = m_bHasOffset;
     237         118 :     return m_dfOffset;
     238             : }
     239             : 
     240             : /************************************************************************/
     241             : /*                              GetScale()                              */
     242             : /************************************************************************/
     243             : 
     244         118 : double PDS4RawRasterBand::GetScale(int *pbSuccess)
     245             : {
     246         118 :     if (pbSuccess)
     247         110 :         *pbSuccess = m_bHasScale;
     248         118 :     return m_dfScale;
     249             : }
     250             : 
     251             : /************************************************************************/
     252             : /*                              SetOffset()                             */
     253             : /************************************************************************/
     254             : 
     255         264 : CPLErr PDS4RawRasterBand::SetOffset(double dfNewOffset)
     256             : {
     257         264 :     m_dfOffset = dfNewOffset;
     258         264 :     m_bHasOffset = true;
     259         264 :     return CE_None;
     260             : }
     261             : 
     262             : /************************************************************************/
     263             : /*                              SetScale()                              */
     264             : /************************************************************************/
     265             : 
     266         264 : CPLErr PDS4RawRasterBand::SetScale(double dfNewScale)
     267             : {
     268         264 :     m_dfScale = dfNewScale;
     269         264 :     m_bHasScale = true;
     270         264 :     return CE_None;
     271             : }
     272             : 
     273             : /************************************************************************/
     274             : /*                           GetNoDataValue()                           */
     275             : /************************************************************************/
     276             : 
     277         212 : double PDS4RawRasterBand::GetNoDataValue(int *pbSuccess)
     278             : {
     279         212 :     if (pbSuccess)
     280         212 :         *pbSuccess = m_bHasNoData;
     281         212 :     return m_dfNoData;
     282             : }
     283             : 
     284             : /************************************************************************/
     285             : /*                           SetNoDataValue()                           */
     286             : /************************************************************************/
     287             : 
     288          76 : CPLErr PDS4RawRasterBand::SetNoDataValue(double dfNewNoData)
     289             : {
     290          76 :     m_dfNoData = dfNewNoData;
     291          76 :     m_bHasNoData = true;
     292          76 :     return CE_None;
     293             : }
     294             : 
     295             : /************************************************************************/
     296             : /*                             IReadBlock()                             */
     297             : /************************************************************************/
     298             : 
     299         132 : CPLErr PDS4RawRasterBand::IWriteBlock(int nXBlock, int nYBlock, void *pImage)
     300             : 
     301             : {
     302         132 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     303         132 :     if (poGDS->m_bMustInitImageFile)
     304             :     {
     305           0 :         if (!poGDS->InitImageFile())
     306           0 :             return CE_Failure;
     307             :     }
     308             : 
     309         132 :     return RawRasterBand::IWriteBlock(nXBlock, nYBlock, pImage);
     310             : }
     311             : 
     312             : /************************************************************************/
     313             : /*                             IRasterIO()                              */
     314             : /************************************************************************/
     315             : 
     316        1098 : CPLErr PDS4RawRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
     317             :                                     int nXSize, int nYSize, void *pData,
     318             :                                     int nBufXSize, int nBufYSize,
     319             :                                     GDALDataType eBufType, GSpacing nPixelSpace,
     320             :                                     GSpacing nLineSpace,
     321             :                                     GDALRasterIOExtraArg *psExtraArg)
     322             : 
     323             : {
     324        1098 :     PDS4Dataset *poGDS = reinterpret_cast<PDS4Dataset *>(poDS);
     325        1098 :     if (eRWFlag == GF_Write && poGDS->m_bMustInitImageFile)
     326             :     {
     327          10 :         if (!poGDS->InitImageFile())
     328           0 :             return CE_Failure;
     329             :     }
     330             : 
     331        1098 :     return RawRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
     332             :                                     pData, nBufXSize, nBufYSize, eBufType,
     333        1098 :                                     nPixelSpace, nLineSpace, psExtraArg);
     334             : }
     335             : 
     336             : /************************************************************************/
     337             : /*                            PDS4MaskBand()                            */
     338             : /************************************************************************/
     339             : 
     340          39 : PDS4MaskBand::PDS4MaskBand(GDALRasterBand *poBaseBand,
     341          39 :                            const std::vector<double> &adfConstants)
     342          39 :     : m_poBaseBand(poBaseBand), m_pBuffer(nullptr), m_adfConstants(adfConstants)
     343             : {
     344          39 :     eDataType = GDT_Byte;
     345          39 :     poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
     346          39 :     nRasterXSize = poBaseBand->GetXSize();
     347          39 :     nRasterYSize = poBaseBand->GetYSize();
     348          39 : }
     349             : 
     350             : /************************************************************************/
     351             : /*                           ~PDS4MaskBand()                            */
     352             : /************************************************************************/
     353             : 
     354          78 : PDS4MaskBand::~PDS4MaskBand()
     355             : {
     356          39 :     VSIFree(m_pBuffer);
     357          78 : }
     358             : 
     359             : /************************************************************************/
     360             : /*                             FillMask()                               */
     361             : /************************************************************************/
     362             : 
     363             : template <class T>
     364          88 : static void FillMask(void *pvBuffer, GByte *pabyDst, int nReqXSize,
     365             :                      int nReqYSize, int nBlockXSize,
     366             :                      const std::vector<double> &adfConstants)
     367             : {
     368          88 :     const T *pSrc = static_cast<T *>(pvBuffer);
     369         176 :     std::vector<T> aConstants;
     370         248 :     for (size_t i = 0; i < adfConstants.size(); i++)
     371             :     {
     372             :         T cst;
     373         160 :         GDALCopyWord(adfConstants[i], cst);
     374         160 :         aConstants.push_back(cst);
     375             :     }
     376             : 
     377         176 :     for (int y = 0; y < nReqYSize; y++)
     378             :     {
     379        1696 :         for (int x = 0; x < nReqXSize; x++)
     380             :         {
     381        1608 :             const T nSrc = pSrc[y * nBlockXSize + x];
     382        1608 :             if (std::find(aConstants.begin(), aConstants.end(), nSrc) !=
     383             :                 aConstants.end())
     384             :             {
     385           6 :                 pabyDst[y * nBlockXSize + x] = 0;
     386             :             }
     387             :             else
     388             :             {
     389        1602 :                 pabyDst[y * nBlockXSize + x] = 255;
     390             :             }
     391             :         }
     392             :     }
     393          88 : }
     394             : 
     395             : /************************************************************************/
     396             : /*                           IReadBlock()                               */
     397             : /************************************************************************/
     398             : 
     399          88 : CPLErr PDS4MaskBand::IReadBlock(int nXBlock, int nYBlock, void *pImage)
     400             : 
     401             : {
     402          88 :     const GDALDataType eSrcDT = m_poBaseBand->GetRasterDataType();
     403          88 :     const int nSrcDTSize = GDALGetDataTypeSizeBytes(eSrcDT);
     404          88 :     if (m_pBuffer == nullptr)
     405             :     {
     406          12 :         m_pBuffer = VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, nSrcDTSize);
     407          12 :         if (m_pBuffer == nullptr)
     408           0 :             return CE_Failure;
     409             :     }
     410             : 
     411          88 :     int nXOff = nXBlock * nBlockXSize;
     412          88 :     int nReqXSize = nBlockXSize;
     413          88 :     if (nXOff + nReqXSize > nRasterXSize)
     414           0 :         nReqXSize = nRasterXSize - nXOff;
     415          88 :     int nYOff = nYBlock * nBlockYSize;
     416          88 :     int nReqYSize = nBlockYSize;
     417          88 :     if (nYOff + nReqYSize > nRasterYSize)
     418           0 :         nReqYSize = nRasterYSize - nYOff;
     419             : 
     420         176 :     if (m_poBaseBand->RasterIO(GF_Read, nXOff, nYOff, nReqXSize, nReqYSize,
     421             :                                m_pBuffer, nReqXSize, nReqYSize, eSrcDT,
     422             :                                nSrcDTSize,
     423          88 :                                static_cast<GSpacing>(nSrcDTSize) * nBlockXSize,
     424          88 :                                nullptr) != CE_None)
     425             :     {
     426           0 :         return CE_Failure;
     427             :     }
     428             : 
     429          88 :     GByte *pabyDst = static_cast<GByte *>(pImage);
     430          88 :     if (eSrcDT == GDT_Byte)
     431             :     {
     432          81 :         FillMask<GByte>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     433          81 :                         m_adfConstants);
     434             :     }
     435           7 :     else if (eSrcDT == GDT_Int8)
     436             :     {
     437           1 :         FillMask<GInt8>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     438           1 :                         m_adfConstants);
     439             :     }
     440           6 :     else if (eSrcDT == GDT_UInt16)
     441             :     {
     442           1 :         FillMask<GUInt16>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     443           1 :                           m_adfConstants);
     444             :     }
     445           5 :     else if (eSrcDT == GDT_Int16)
     446             :     {
     447           1 :         FillMask<GInt16>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     448           1 :                          m_adfConstants);
     449             :     }
     450           4 :     else if (eSrcDT == GDT_UInt32)
     451             :     {
     452           1 :         FillMask<GUInt32>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     453           1 :                           m_adfConstants);
     454             :     }
     455           3 :     else if (eSrcDT == GDT_Int32)
     456             :     {
     457           1 :         FillMask<GInt32>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     458           1 :                          m_adfConstants);
     459             :     }
     460           2 :     else if (eSrcDT == GDT_Float32)
     461             :     {
     462           1 :         FillMask<float>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     463           1 :                         m_adfConstants);
     464             :     }
     465           1 :     else if (eSrcDT == GDT_Float64)
     466             :     {
     467           1 :         FillMask<double>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
     468           1 :                          m_adfConstants);
     469             :     }
     470             : 
     471          88 :     return CE_None;
     472             : }
     473             : 
     474             : /************************************************************************/
     475             : /*                            PDS4Dataset()                             */
     476             : /************************************************************************/
     477             : 
     478         453 : PDS4Dataset::PDS4Dataset()
     479             : {
     480         453 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     481         453 : }
     482             : 
     483             : /************************************************************************/
     484             : /*                           ~PDS4Dataset()                             */
     485             : /************************************************************************/
     486             : 
     487         906 : PDS4Dataset::~PDS4Dataset()
     488             : {
     489         453 :     PDS4Dataset::Close();
     490         906 : }
     491             : 
     492             : /************************************************************************/
     493             : /*                              Close()                                 */
     494             : /************************************************************************/
     495             : 
     496         772 : CPLErr PDS4Dataset::Close()
     497             : {
     498         772 :     CPLErr eErr = CE_None;
     499         772 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     500             :     {
     501         453 :         if (m_bMustInitImageFile)
     502             :         {
     503          58 :             if (!InitImageFile())
     504           0 :                 eErr = CE_Failure;
     505             :         }
     506             : 
     507         453 :         if (PDS4Dataset::FlushCache(true) != CE_None)
     508           0 :             eErr = CE_Failure;
     509             : 
     510         453 :         if (m_bCreateHeader || m_bDirtyHeader)
     511         177 :             WriteHeader();
     512         453 :         if (m_fpImage)
     513         309 :             VSIFCloseL(m_fpImage);
     514         453 :         CSLDestroy(m_papszCreationOptions);
     515         453 :         PDS4Dataset::CloseDependentDatasets();
     516             : 
     517         453 :         if (GDALPamDataset::Close() != CE_None)
     518           0 :             eErr = CE_Failure;
     519             :     }
     520         772 :     return eErr;
     521             : }
     522             : 
     523             : /************************************************************************/
     524             : /*                        GetRawBinaryLayout()                          */
     525             : /************************************************************************/
     526             : 
     527           1 : bool PDS4Dataset::GetRawBinaryLayout(GDALDataset::RawBinaryLayout &sLayout)
     528             : {
     529           1 :     if (!RawDataset::GetRawBinaryLayout(sLayout))
     530           0 :         return false;
     531           1 :     sLayout.osRawFilename = m_osImageFilename;
     532           1 :     return true;
     533             : }
     534             : 
     535             : /************************************************************************/
     536             : /*                        CloseDependentDatasets()                      */
     537             : /************************************************************************/
     538             : 
     539         453 : int PDS4Dataset::CloseDependentDatasets()
     540             : {
     541         453 :     int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
     542             : 
     543         453 :     if (m_poExternalDS)
     544             :     {
     545          11 :         bHasDroppedRef = FALSE;
     546          11 :         delete m_poExternalDS;
     547          11 :         m_poExternalDS = nullptr;
     548             : 
     549          28 :         for (int iBand = 0; iBand < nBands; iBand++)
     550             :         {
     551          17 :             delete papoBands[iBand];
     552          17 :             papoBands[iBand] = nullptr;
     553             :         }
     554          11 :         nBands = 0;
     555             :     }
     556             : 
     557         453 :     return bHasDroppedRef;
     558             : }
     559             : 
     560             : /************************************************************************/
     561             : /*                         GetSpatialRef()                              */
     562             : /************************************************************************/
     563             : 
     564          44 : const OGRSpatialReference *PDS4Dataset::GetSpatialRef() const
     565             : {
     566          44 :     if (!m_oSRS.IsEmpty())
     567          40 :         return &m_oSRS;
     568           4 :     return GDALPamDataset::GetSpatialRef();
     569             : }
     570             : 
     571             : /************************************************************************/
     572             : /*                           SetSpatialRef()                            */
     573             : /************************************************************************/
     574             : 
     575          95 : CPLErr PDS4Dataset::SetSpatialRef(const OGRSpatialReference *poSRS)
     576             : 
     577             : {
     578          95 :     if (eAccess == GA_ReadOnly)
     579           0 :         return CE_Failure;
     580          95 :     m_oSRS.Clear();
     581          95 :     if (poSRS)
     582          95 :         m_oSRS = *poSRS;
     583          95 :     if (m_poExternalDS)
     584           8 :         m_poExternalDS->SetSpatialRef(poSRS);
     585          95 :     return CE_None;
     586             : }
     587             : 
     588             : /************************************************************************/
     589             : /*                          GetGeoTransform()                           */
     590             : /************************************************************************/
     591             : 
     592          44 : CPLErr PDS4Dataset::GetGeoTransform(GDALGeoTransform &gt) const
     593             : 
     594             : {
     595          44 :     if (m_bGotTransform)
     596             :     {
     597          42 :         gt = m_gt;
     598          42 :         return CE_None;
     599             :     }
     600             : 
     601           2 :     return GDALPamDataset::GetGeoTransform(gt);
     602             : }
     603             : 
     604             : /************************************************************************/
     605             : /*                          SetGeoTransform()                           */
     606             : /************************************************************************/
     607             : 
     608          95 : CPLErr PDS4Dataset::SetGeoTransform(const GDALGeoTransform &gt)
     609             : 
     610             : {
     611          97 :     if (!((gt[1] > 0.0 && gt[2] == 0.0 && gt[4] == 0.0 && gt[5] < 0.0) ||
     612           2 :           (gt[1] == 0.0 && gt[2] > 0.0 && gt[4] > 0.0 && gt[5] == 0.0)))
     613             :     {
     614           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     615             :                  "Only north-up geotransform or map_projection_rotation=90 "
     616             :                  "supported");
     617           0 :         return CE_Failure;
     618             :     }
     619          95 :     m_gt = gt;
     620          95 :     m_bGotTransform = true;
     621          95 :     if (m_poExternalDS)
     622           8 :         m_poExternalDS->SetGeoTransform(m_gt);
     623          95 :     return CE_None;
     624             : }
     625             : 
     626             : /************************************************************************/
     627             : /*                             SetMetadata()                            */
     628             : /************************************************************************/
     629             : 
     630           7 : CPLErr PDS4Dataset::SetMetadata(char **papszMD, const char *pszDomain)
     631             : {
     632           7 :     if (m_bUseSrcLabel && eAccess == GA_Update && pszDomain != nullptr &&
     633           7 :         EQUAL(pszDomain, "xml:PDS4"))
     634             :     {
     635           6 :         if (papszMD != nullptr && papszMD[0] != nullptr)
     636             :         {
     637           6 :             m_osXMLPDS4 = papszMD[0];
     638             :         }
     639           6 :         return CE_None;
     640             :     }
     641           1 :     return GDALPamDataset::SetMetadata(papszMD, pszDomain);
     642             : }
     643             : 
     644             : /************************************************************************/
     645             : /*                            GetFileList()                             */
     646             : /************************************************************************/
     647             : 
     648         135 : char **PDS4Dataset::GetFileList()
     649             : {
     650         135 :     char **papszFileList = GDALPamDataset::GetFileList();
     651         270 :     if (!m_osXMLFilename.empty() &&
     652         135 :         CSLFindString(papszFileList, m_osXMLFilename) < 0)
     653             :     {
     654           0 :         papszFileList = CSLAddString(papszFileList, m_osXMLFilename);
     655             :     }
     656         135 :     if (!m_osImageFilename.empty())
     657             :     {
     658          86 :         papszFileList = CSLAddString(papszFileList, m_osImageFilename);
     659             :     }
     660         199 :     for (const auto &poLayer : m_apoLayers)
     661             :     {
     662          64 :         auto papszTemp = poLayer->GetFileList();
     663          64 :         papszFileList = CSLInsertStrings(papszFileList, -1, papszTemp);
     664          64 :         CSLDestroy(papszTemp);
     665             :     }
     666         135 :     return papszFileList;
     667             : }
     668             : 
     669             : /************************************************************************/
     670             : /*                            GetLinearValue()                          */
     671             : /************************************************************************/
     672             : 
     673             : static const struct
     674             : {
     675             :     const char *pszUnit;
     676             :     double dfToMeter;
     677             : } apsLinearUnits[] = {
     678             :     {"AU", 149597870700.0}, {"Angstrom", 1e-10}, {"cm", 1e-2}, {"km", 1e3},
     679             :     {"micrometer", 1e-6},   {"mm", 1e-3},        {"nm", 1e-9}};
     680             : 
     681         762 : static double GetLinearValue(const CPLXMLNode *psParent,
     682             :                              const char *pszElementName)
     683             : {
     684         762 :     const CPLXMLNode *psNode = CPLGetXMLNode(psParent, pszElementName);
     685         762 :     if (psNode == nullptr)
     686           0 :         return 0.0;
     687         762 :     double dfVal = CPLAtof(CPLGetXMLValue(psNode, nullptr, ""));
     688         762 :     const char *pszUnit = CPLGetXMLValue(psNode, "unit", nullptr);
     689         762 :     if (pszUnit && !EQUAL(pszUnit, "m"))
     690             :     {
     691          19 :         bool bFound = false;
     692          76 :         for (size_t i = 0; i < CPL_ARRAYSIZE(apsLinearUnits); i++)
     693             :         {
     694          76 :             if (EQUAL(pszUnit, apsLinearUnits[i].pszUnit))
     695             :             {
     696          19 :                 dfVal *= apsLinearUnits[i].dfToMeter;
     697          19 :                 bFound = true;
     698          19 :                 break;
     699             :             }
     700             :         }
     701          19 :         if (!bFound)
     702             :         {
     703           0 :             CPLError(CE_Warning, CPLE_AppDefined, "Unknown unit '%s' for '%s'",
     704             :                      pszUnit, pszElementName);
     705             :         }
     706             :     }
     707         762 :     return dfVal;
     708             : }
     709             : 
     710             : /************************************************************************/
     711             : /*                          GetResolutionValue()                        */
     712             : /************************************************************************/
     713             : 
     714             : static const struct
     715             : {
     716             :     const char *pszUnit;
     717             :     double dfToMeter;
     718             : } apsResolutionUnits[] = {
     719             :     {"km/pixel", 1e3},
     720             :     {"mm/pixel", 1e-3},
     721             : };
     722             : 
     723         300 : static double GetResolutionValue(CPLXMLNode *psParent,
     724             :                                  const char *pszElementName)
     725             : {
     726         300 :     CPLXMLNode *psNode = CPLGetXMLNode(psParent, pszElementName);
     727         300 :     if (psNode == nullptr)
     728           0 :         return 0.0;
     729         300 :     double dfVal = CPLAtof(CPLGetXMLValue(psNode, nullptr, ""));
     730         300 :     const char *pszUnit = CPLGetXMLValue(psNode, "unit", nullptr);
     731         300 :     if (pszUnit && !EQUAL(pszUnit, "m/pixel"))
     732             :     {
     733          19 :         bool bFound = false;
     734          19 :         for (size_t i = 0; i < CPL_ARRAYSIZE(apsResolutionUnits); i++)
     735             :         {
     736          19 :             if (EQUAL(pszUnit, apsResolutionUnits[i].pszUnit))
     737             :             {
     738          19 :                 dfVal *= apsResolutionUnits[i].dfToMeter;
     739          19 :                 bFound = true;
     740          19 :                 break;
     741             :             }
     742             :         }
     743          19 :         if (!bFound)
     744             :         {
     745           0 :             CPLError(CE_Warning, CPLE_AppDefined, "Unknown unit '%s' for '%s'",
     746             :                      pszUnit, pszElementName);
     747             :         }
     748             :     }
     749         300 :     return dfVal;
     750             : }
     751             : 
     752             : /************************************************************************/
     753             : /*                            GetAngularValue()                         */
     754             : /************************************************************************/
     755             : 
     756             : static const struct
     757             : {
     758             :     const char *pszUnit;
     759             :     double dfToDeg;
     760             : } apsAngularUnits[] = {{"arcmin", 1. / 60.},
     761             :                        {"arcsec", 1. / 3600},
     762             :                        {"hr", 15.0},
     763             :                        {"mrad", 180.0 / M_PI / 1000.},
     764             :                        {"rad", 180.0 / M_PI}};
     765             : 
     766         779 : static double GetAngularValue(CPLXMLNode *psParent, const char *pszElementName,
     767             :                               bool *pbGotVal = nullptr)
     768             : {
     769         779 :     CPLXMLNode *psNode = CPLGetXMLNode(psParent, pszElementName);
     770         779 :     if (psNode == nullptr)
     771             :     {
     772         406 :         if (pbGotVal)
     773         255 :             *pbGotVal = false;
     774         406 :         return 0.0;
     775             :     }
     776         373 :     double dfVal = CPLAtof(CPLGetXMLValue(psNode, nullptr, ""));
     777         373 :     const char *pszUnit = CPLGetXMLValue(psNode, "unit", nullptr);
     778         373 :     if (pszUnit && !EQUAL(pszUnit, "deg"))
     779             :     {
     780           0 :         bool bFound = false;
     781           0 :         for (size_t i = 0; i < CPL_ARRAYSIZE(apsAngularUnits); i++)
     782             :         {
     783           0 :             if (EQUAL(pszUnit, apsAngularUnits[i].pszUnit))
     784             :             {
     785           0 :                 dfVal *= apsAngularUnits[i].dfToDeg;
     786           0 :                 bFound = true;
     787           0 :                 break;
     788             :             }
     789             :         }
     790           0 :         if (!bFound)
     791             :         {
     792           0 :             CPLError(CE_Warning, CPLE_AppDefined, "Unknown unit '%s' for '%s'",
     793             :                      pszUnit, pszElementName);
     794             :         }
     795             :     }
     796         373 :     if (pbGotVal)
     797         206 :         *pbGotVal = true;
     798         373 :     return dfVal;
     799             : }
     800             : 
     801             : /************************************************************************/
     802             : /*                          ReadGeoreferencing()                       */
     803             : /************************************************************************/
     804             : 
     805             : // See https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1G00_1950.xsd, (GDAL 3.4)
     806             : //     https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1D00_1933.xsd,
     807             : //     https://raw.githubusercontent.com/nasa-pds-data-dictionaries/ldd-cart/master/build/1.B.0.0/PDS4_CART_1B00.xsd,
     808             : //     https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1700.xsd
     809             : // and the corresponding .sch files
     810         276 : void PDS4Dataset::ReadGeoreferencing(CPLXMLNode *psProduct)
     811             : {
     812         276 :     CPLXMLNode *psCart = CPLGetXMLNode(
     813             :         psProduct, "Observation_Area.Discipline_Area.Cartography");
     814         276 :     if (psCart == nullptr)
     815             :     {
     816         122 :         CPLDebug("PDS4",
     817             :                  "Did not find Observation_Area.Discipline_Area.Cartography");
     818         122 :         return;
     819             :     }
     820             : 
     821             :     // Bounding box: informative only
     822             :     CPLXMLNode *psBounding =
     823         154 :         CPLGetXMLNode(psCart, "Spatial_Domain.Bounding_Coordinates");
     824         154 :     if (psBounding)
     825             :     {
     826             :         const char *pszWest =
     827         154 :             CPLGetXMLValue(psBounding, "west_bounding_coordinate", nullptr);
     828             :         const char *pszEast =
     829         154 :             CPLGetXMLValue(psBounding, "east_bounding_coordinate", nullptr);
     830             :         const char *pszNorth =
     831         154 :             CPLGetXMLValue(psBounding, "north_bounding_coordinate", nullptr);
     832             :         const char *pszSouth =
     833         154 :             CPLGetXMLValue(psBounding, "south_bounding_coordinate", nullptr);
     834         154 :         if (pszWest)
     835         154 :             CPLDebug("PDS4", "West: %s", pszWest);
     836         154 :         if (pszEast)
     837         154 :             CPLDebug("PDS4", "East: %s", pszEast);
     838         154 :         if (pszNorth)
     839         154 :             CPLDebug("PDS4", "North: %s", pszNorth);
     840         154 :         if (pszSouth)
     841         154 :             CPLDebug("PDS4", "South: %s", pszSouth);
     842             :     }
     843             : 
     844             :     CPLXMLNode *psSR =
     845         154 :         CPLGetXMLNode(psCart, "Spatial_Reference_Information.Horizontal_"
     846             :                               "Coordinate_System_Definition");
     847         154 :     if (psSR == nullptr)
     848             :     {
     849           0 :         CPLDebug("PDS4", "Did not find Spatial_Reference_Information."
     850             :                          "Horizontal_Coordinate_System_Definition");
     851           0 :         return;
     852             :     }
     853             : 
     854         154 :     double dfLongitudeMultiplier = 1;
     855         154 :     const CPLXMLNode *psGeodeticModel = CPLGetXMLNode(psSR, "Geodetic_Model");
     856         154 :     if (psGeodeticModel != nullptr)
     857             :     {
     858         154 :         if (EQUAL(CPLGetXMLValue(psGeodeticModel, "longitude_direction", ""),
     859             :                   "Positive West"))
     860             :         {
     861           3 :             dfLongitudeMultiplier = -1;
     862             :         }
     863             :     }
     864             : 
     865         308 :     OGRSpatialReference oSRS;
     866             :     CPLXMLNode *psGridCoordinateSystem =
     867         154 :         CPLGetXMLNode(psSR, "Planar.Grid_Coordinate_System");
     868         154 :     CPLXMLNode *psMapProjection = CPLGetXMLNode(psSR, "Planar.Map_Projection");
     869         308 :     CPLString osProjName;
     870         154 :     double dfCenterLon = 0.0;
     871         154 :     double dfCenterLat = 0.0;
     872         154 :     double dfStdParallel1 = 0.0;
     873         154 :     double dfStdParallel2 = 0.0;
     874         154 :     double dfScale = 1.0;
     875         154 :     double dfMapProjectionRotation = 0.0;
     876         154 :     if (psGridCoordinateSystem != nullptr)
     877             :     {
     878             :         osProjName = CPLGetXMLValue(psGridCoordinateSystem,
     879           0 :                                     "grid_coordinate_system_name", "");
     880           0 :         if (!osProjName.empty())
     881             :         {
     882           0 :             if (osProjName == "Universal Transverse Mercator")
     883             :             {
     884           0 :                 CPLXMLNode *psUTMZoneNumber = CPLGetXMLNode(
     885             :                     psGridCoordinateSystem,
     886             :                     "Universal_Transverse_Mercator.utm_zone_number");
     887           0 :                 if (psUTMZoneNumber)
     888             :                 {
     889             :                     int nZone =
     890           0 :                         atoi(CPLGetXMLValue(psUTMZoneNumber, nullptr, ""));
     891           0 :                     oSRS.SetUTM(std::abs(nZone), nZone >= 0);
     892             :                 }
     893             :             }
     894           0 :             else if (osProjName == "Universal Polar Stereographic")
     895             :             {
     896           0 :                 CPLXMLNode *psProjParamNode = CPLGetXMLNode(
     897             :                     psGridCoordinateSystem,
     898             :                     "Universal_Polar_Stereographic.Polar_Stereographic");
     899           0 :                 if (psProjParamNode)
     900             :                 {
     901           0 :                     dfCenterLon =
     902           0 :                         GetAngularValue(psProjParamNode,
     903             :                                         "longitude_of_central_meridian") *
     904             :                         dfLongitudeMultiplier;
     905           0 :                     dfCenterLat = GetAngularValue(
     906             :                         psProjParamNode, "latitude_of_projection_origin");
     907           0 :                     dfScale = CPLAtof(CPLGetXMLValue(
     908             :                         psProjParamNode, "scale_factor_at_projection_origin",
     909             :                         "1"));
     910           0 :                     oSRS.SetPS(dfCenterLat, dfCenterLon, dfScale, 0, 0);
     911             :                 }
     912             :             }
     913             :             else
     914             :             {
     915           0 :                 CPLError(CE_Warning, CPLE_NotSupported,
     916             :                          "grid_coordinate_system_name = %s not supported",
     917             :                          osProjName.c_str());
     918             :             }
     919             :         }
     920             :     }
     921         154 :     else if (psMapProjection != nullptr)
     922             :     {
     923         153 :         osProjName = CPLGetXMLValue(psMapProjection, "map_projection_name", "");
     924         153 :         if (!osProjName.empty())
     925             :         {
     926         153 :             CPLXMLNode *psProjParamNode = CPLGetXMLNode(
     927             :                 psMapProjection,
     928         306 :                 CPLString(osProjName).replaceAll(' ', '_').c_str());
     929         153 :             if (psProjParamNode == nullptr &&
     930             :                 // typo in https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1700.sch
     931           0 :                 EQUAL(osProjName, "Orothographic"))
     932             :             {
     933             :                 psProjParamNode =
     934           0 :                     CPLGetXMLNode(psMapProjection, "Orthographic");
     935             :             }
     936         153 :             bool bGotStdParallel1 = false;
     937         153 :             bool bGotStdParallel2 = false;
     938         153 :             bool bGotScale = false;
     939         153 :             if (psProjParamNode)
     940             :             {
     941         153 :                 bool bGotCenterLon = false;
     942         153 :                 dfCenterLon = GetAngularValue(psProjParamNode,
     943             :                                               "longitude_of_central_meridian",
     944             :                                               &bGotCenterLon) *
     945             :                               dfLongitudeMultiplier;
     946         153 :                 if (!bGotCenterLon)
     947             :                 {
     948           2 :                     dfCenterLon =
     949           2 :                         GetAngularValue(psProjParamNode,
     950             :                                         "straight_vertical_longitude_from_pole",
     951             :                                         &bGotCenterLon) *
     952             :                         dfLongitudeMultiplier;
     953             :                 }
     954         153 :                 dfCenterLat = GetAngularValue(psProjParamNode,
     955             :                                               "latitude_of_projection_origin");
     956         153 :                 dfStdParallel1 = GetAngularValue(
     957             :                     psProjParamNode, "standard_parallel_1", &bGotStdParallel1);
     958         153 :                 dfStdParallel2 = GetAngularValue(
     959             :                     psProjParamNode, "standard_parallel_2", &bGotStdParallel2);
     960             :                 const char *pszScaleParam =
     961         153 :                     (osProjName == "Transverse Mercator")
     962         153 :                         ? "scale_factor_at_central_meridian"
     963         153 :                         : "scale_factor_at_projection_origin";
     964             :                 const char *pszScaleVal =
     965         153 :                     CPLGetXMLValue(psProjParamNode, pszScaleParam, nullptr);
     966         153 :                 bGotScale = pszScaleVal != nullptr;
     967         153 :                 dfScale = (pszScaleVal) ? CPLAtof(pszScaleVal) : 1.0;
     968             : 
     969             :                 dfMapProjectionRotation =
     970         153 :                     GetAngularValue(psProjParamNode, "map_projection_rotation");
     971             :             }
     972             : 
     973             :             CPLXMLNode *psObliqueAzimuth =
     974         153 :                 CPLGetXMLNode(psProjParamNode, "Oblique_Line_Azimuth");
     975             :             CPLXMLNode *psObliquePoint =
     976         153 :                 CPLGetXMLNode(psProjParamNode, "Oblique_Line_Point");
     977             : 
     978         153 :             if (EQUAL(osProjName, "Equirectangular"))
     979             :             {
     980          49 :                 oSRS.SetEquirectangular2(dfCenterLat, dfCenterLon,
     981             :                                          dfStdParallel1, 0, 0);
     982             :             }
     983         104 :             else if (EQUAL(osProjName, "Lambert Conformal Conic"))
     984             :             {
     985           4 :                 if (bGotScale)
     986             :                 {
     987           2 :                     if ((bGotStdParallel1 && dfStdParallel1 != dfCenterLat) ||
     988           0 :                         (bGotStdParallel2 && dfStdParallel2 != dfCenterLat))
     989             :                     {
     990           0 :                         CPLError(
     991             :                             CE_Warning, CPLE_AppDefined,
     992             :                             "Ignoring standard_parallel_1 and/or "
     993             :                             "standard_parallel_2 with LCC_1SP formulation");
     994             :                     }
     995           2 :                     oSRS.SetLCC1SP(dfCenterLat, dfCenterLon, dfScale, 0, 0);
     996             :                 }
     997             :                 else
     998             :                 {
     999           2 :                     oSRS.SetLCC(dfStdParallel1, dfStdParallel2, dfCenterLat,
    1000             :                                 dfCenterLon, 0, 0);
    1001             :                 }
    1002             :             }
    1003         100 :             else if (EQUAL(osProjName, "Mercator"))
    1004             :             {
    1005           6 :                 if (bGotScale)
    1006             :                 {
    1007           4 :                     oSRS.SetMercator(dfCenterLat,  // should be 0 normally
    1008             :                                      dfCenterLon, dfScale, 0.0, 0.0);
    1009             :                 }
    1010             :                 else
    1011             :                 {
    1012           2 :                     oSRS.SetMercator2SP(dfStdParallel1,
    1013             :                                         dfCenterLat,  // should be 0 normally
    1014             :                                         dfCenterLon, 0.0, 0.0);
    1015             :                 }
    1016             :             }
    1017          94 :             else if (EQUAL(osProjName, "Orthographic"))
    1018             :             {
    1019           2 :                 oSRS.SetOrthographic(dfCenterLat, dfCenterLon, 0.0, 0.0);
    1020             :             }
    1021          94 :             else if (EQUAL(osProjName, "Oblique Mercator") &&
    1022           2 :                      (psObliqueAzimuth != nullptr || psObliquePoint != nullptr))
    1023             :             {
    1024           4 :                 if (psObliqueAzimuth)
    1025             :                 {
    1026             :                     // Not sure of this
    1027           2 :                     dfCenterLon = CPLAtof(
    1028             :                         CPLGetXMLValue(psObliqueAzimuth,
    1029             :                                        "azimuth_measure_point_longitude", "0"));
    1030             : 
    1031           2 :                     double dfAzimuth = CPLAtof(CPLGetXMLValue(
    1032             :                         psObliqueAzimuth, "azimuthal_angle", "0"));
    1033           2 :                     oSRS.SetProjection(
    1034             :                         SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER);
    1035           2 :                     oSRS.SetNormProjParm(SRS_PP_LATITUDE_OF_CENTER,
    1036             :                                          dfCenterLat);
    1037           2 :                     oSRS.SetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER,
    1038             :                                          dfCenterLon);
    1039           2 :                     oSRS.SetNormProjParm(SRS_PP_AZIMUTH, dfAzimuth);
    1040             :                     // SetNormProjParm( SRS_PP_RECTIFIED_GRID_ANGLE,
    1041             :                     // dfRectToSkew );
    1042           2 :                     oSRS.SetNormProjParm(SRS_PP_SCALE_FACTOR, dfScale);
    1043           2 :                     oSRS.SetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
    1044           2 :                     oSRS.SetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
    1045             :                 }
    1046             :                 else
    1047             :                 {
    1048           2 :                     double dfLat1 = 0.0;
    1049           2 :                     double dfLong1 = 0.0;
    1050           2 :                     double dfLat2 = 0.0;
    1051           2 :                     double dfLong2 = 0.0;
    1052           2 :                     CPLXMLNode *psPoint = CPLGetXMLNode(
    1053             :                         psObliquePoint, "Oblique_Line_Point_Group");
    1054           2 :                     if (psPoint)
    1055             :                     {
    1056           2 :                         dfLat1 = CPLAtof(CPLGetXMLValue(
    1057             :                             psPoint, "oblique_line_latitude", "0.0"));
    1058           2 :                         dfLong1 = CPLAtof(CPLGetXMLValue(
    1059             :                             psPoint, "oblique_line_longitude", "0.0"));
    1060           2 :                         psPoint = psPoint->psNext;
    1061           2 :                         if (psPoint && psPoint->eType == CXT_Element &&
    1062           2 :                             EQUAL(psPoint->pszValue,
    1063             :                                   "Oblique_Line_Point_Group"))
    1064             :                         {
    1065           2 :                             dfLat2 = CPLAtof(CPLGetXMLValue(
    1066             :                                 psPoint, "oblique_line_latitude", "0.0"));
    1067           2 :                             dfLong2 = CPLAtof(CPLGetXMLValue(
    1068             :                                 psPoint, "oblique_line_longitude", "0.0"));
    1069             :                         }
    1070             :                     }
    1071           2 :                     oSRS.SetHOM2PNO(dfCenterLat, dfLat1, dfLong1, dfLat2,
    1072             :                                     dfLong2, dfScale, 0.0, 0.0);
    1073             :                 }
    1074             :             }
    1075          88 :             else if (EQUAL(osProjName, "Polar Stereographic"))
    1076             :             {
    1077           2 :                 oSRS.SetPS(dfCenterLat, dfCenterLon, dfScale, 0, 0);
    1078             :             }
    1079          86 :             else if (EQUAL(osProjName, "Polyconic"))
    1080             :             {
    1081           2 :                 oSRS.SetPolyconic(dfCenterLat, dfCenterLon, 0, 0);
    1082             :             }
    1083          84 :             else if (EQUAL(osProjName, "Sinusoidal"))
    1084             :             {
    1085           4 :                 oSRS.SetSinusoidal(dfCenterLon, 0, 0);
    1086             :             }
    1087          80 :             else if (EQUAL(osProjName, "Transverse Mercator"))
    1088             :             {
    1089          74 :                 oSRS.SetTM(dfCenterLat, dfCenterLon, dfScale, 0, 0);
    1090             :             }
    1091             : 
    1092             :             // Below values are valid map_projection_name according to
    1093             :             // the schematron but they don't have a dedicated element to
    1094             :             // hold the projection parameter. Assumed the schema is extended
    1095             :             // similarly to the existing for a few obvious ones
    1096           6 :             else if (EQUAL(osProjName, "Albers Conical Equal Area"))
    1097             :             {
    1098           0 :                 oSRS.SetACEA(dfStdParallel1, dfStdParallel2, dfCenterLat,
    1099             :                              dfCenterLon, 0.0, 0.0);
    1100             :             }
    1101           6 :             else if (EQUAL(osProjName, "Azimuthal Equidistant"))
    1102             :             {
    1103           0 :                 oSRS.SetAE(dfCenterLat, dfCenterLon, 0, 0);
    1104             :             }
    1105           6 :             else if (EQUAL(osProjName, "Equidistant Conic"))
    1106             :             {
    1107           0 :                 oSRS.SetEC(dfStdParallel1, dfStdParallel2, dfCenterLat,
    1108             :                            dfCenterLon, 0.0, 0.0);
    1109             :             }
    1110             :             // Unhandled: General Vertical Near-sided Projection
    1111           6 :             else if (EQUAL(osProjName, "Gnomonic"))
    1112             :             {
    1113           0 :                 oSRS.SetGnomonic(dfCenterLat, dfCenterLon, 0, 0);
    1114             :             }
    1115           6 :             else if (EQUAL(osProjName, "Lambert Azimuthal Equal Area"))
    1116             :             {
    1117           2 :                 oSRS.SetLAEA(dfCenterLat, dfCenterLon, 0, 0);
    1118             :             }
    1119           4 :             else if (EQUAL(osProjName, "Miller Cylindrical"))
    1120             :             {
    1121           0 :                 oSRS.SetMC(dfCenterLat, dfCenterLon, 0, 0);
    1122             :             }
    1123           4 :             else if (EQUAL(osProjName, "Orothographic")  // typo
    1124           4 :                      || EQUAL(osProjName, "Orthographic"))
    1125             :             {
    1126           0 :                 osProjName = "Orthographic";
    1127           0 :                 oSRS.SetOrthographic(dfCenterLat, dfCenterLon, 0, 0);
    1128             :             }
    1129           4 :             else if (EQUAL(osProjName, "Robinson"))
    1130             :             {
    1131           0 :                 oSRS.SetRobinson(dfCenterLon, 0, 0);
    1132             :             }
    1133             :             // Unhandled: Space Oblique Mercator
    1134           4 :             else if (EQUAL(osProjName, "Stereographic"))
    1135             :             {
    1136           0 :                 oSRS.SetStereographic(dfCenterLat, dfCenterLon, dfScale, 0, 0);
    1137             :             }
    1138           4 :             else if (EQUAL(osProjName, "van der Grinten"))
    1139             :             {
    1140           0 :                 oSRS.SetVDG(dfCenterLon, 0, 0);
    1141             :             }
    1142           4 :             else if (EQUAL(osProjName, "Oblique Cylindrical"))
    1143             :             {
    1144           4 :                 const double poleLatitude = GetAngularValue(
    1145             :                     psProjParamNode, "oblique_proj_pole_latitude");
    1146             :                 const double poleLongitude =
    1147           4 :                     GetAngularValue(psProjParamNode,
    1148             :                                     "oblique_proj_pole_longitude") *
    1149           4 :                     dfLongitudeMultiplier;
    1150           4 :                 const double poleRotation = GetAngularValue(
    1151             :                     psProjParamNode, "oblique_proj_pole_rotation");
    1152             : 
    1153           8 :                 CPLString oProj4String;
    1154             :                 // Cf isis3dataset.cpp comments for ObliqueCylindrical
    1155             :                 oProj4String.Printf("+proj=ob_tran +o_proj=eqc +o_lon_p=%.17g "
    1156             :                                     "+o_lat_p=%.17g +lon_0=%.17g",
    1157             :                                     -poleRotation, 180 - poleLatitude,
    1158           4 :                                     poleLongitude);
    1159           4 :                 oSRS.SetFromUserInput(oProj4String);
    1160             :             }
    1161             :             else
    1162             :             {
    1163           0 :                 CPLError(CE_Warning, CPLE_NotSupported,
    1164             :                          "map_projection_name = %s not supported",
    1165             :                          osProjName.c_str());
    1166             :             }
    1167             :         }
    1168             :     }
    1169             :     else
    1170             :     {
    1171           1 :         CPLXMLNode *psGeographic = CPLGetXMLNode(psSR, "Geographic");
    1172           1 :         if (GetLayerCount() && psGeographic)
    1173             :         {
    1174             :             // do nothing
    1175             :         }
    1176             :         else
    1177             :         {
    1178           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    1179             :                      "Planar.Map_Projection not found");
    1180             :         }
    1181             :     }
    1182             : 
    1183         154 :     if (oSRS.IsProjected())
    1184             :     {
    1185         153 :         oSRS.SetLinearUnits("Metre", 1.0);
    1186             :     }
    1187             : 
    1188         154 :     if (psGeodeticModel != nullptr)
    1189             :     {
    1190             :         const char *pszLatitudeType =
    1191         154 :             CPLGetXMLValue(psGeodeticModel, "latitude_type", "");
    1192         154 :         bool bIsOgraphic = EQUAL(pszLatitudeType, "Planetographic");
    1193             : 
    1194             :         const bool bUseLDD1930RadiusNames =
    1195         154 :             CPLGetXMLNode(psGeodeticModel, "a_axis_radius") != nullptr;
    1196             : 
    1197             :         // Before PDS CART schema pre-1.B.10.0 (pre LDD version 1.9.3.0),
    1198             :         // the confusing semi_major_radius, semi_minor_radius and polar_radius
    1199             :         // were used but did not follow the recommended
    1200             :         // FGDC names. Using both "semi" and "radius" in the same keyword,
    1201             :         // which both mean half, does not make sense.
    1202         154 :         const char *pszAAxis =
    1203         154 :             bUseLDD1930RadiusNames ? "a_axis_radius" : "semi_major_radius";
    1204         154 :         const char *pszBAxis =
    1205         154 :             bUseLDD1930RadiusNames ? "b_axis_radius" : "semi_minor_radius";
    1206         154 :         const char *pszCAxis =
    1207         154 :             bUseLDD1930RadiusNames ? "c_axis_radius" : "polar_radius";
    1208             : 
    1209         154 :         const double dfSemiMajor = GetLinearValue(psGeodeticModel, pszAAxis);
    1210             : 
    1211             :         // a_axis_radius and b_axis_radius should be the same in most cases
    1212             :         // unless a triaxial body is being defined. This should be extremely
    1213             :         // rare (and not used) since the IAU generally defines a best-fit sphere
    1214             :         // for triaxial bodies: https://astrogeology.usgs.gov/groups/IAU-WGCCRE
    1215         154 :         const double dfBValue = GetLinearValue(psGeodeticModel, pszBAxis);
    1216         154 :         if (dfSemiMajor != dfBValue)
    1217             :         {
    1218           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    1219             :                      "%s = %f m, different from "
    1220             :                      "%s = %f, will be ignored",
    1221             :                      pszBAxis, dfBValue, pszAAxis, dfSemiMajor);
    1222             :         }
    1223             : 
    1224         154 :         const double dfPolarRadius = GetLinearValue(psGeodeticModel, pszCAxis);
    1225             :         // Use the polar_radius as the actual semi minor
    1226         154 :         const double dfSemiMinor = dfPolarRadius;
    1227             : 
    1228             :         // Compulsory
    1229         154 :         const char *pszTargetName = CPLGetXMLValue(
    1230             :             psProduct, "Observation_Area.Target_Identification.name",
    1231             :             "unknown");
    1232             : 
    1233         154 :         if (oSRS.IsProjected())
    1234             :         {
    1235         459 :             CPLString osProjTargetName = osProjName + " " + pszTargetName;
    1236         153 :             oSRS.SetProjCS(osProjTargetName);
    1237             :         }
    1238             : 
    1239         462 :         CPLString osGeogName = CPLString("GCS_") + pszTargetName;
    1240             : 
    1241             :         CPLString osSphereName =
    1242         308 :             CPLGetXMLValue(psGeodeticModel, "spheroid_name", pszTargetName);
    1243         308 :         CPLString osDatumName = "D_" + osSphereName;
    1244             : 
    1245             :         // calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
    1246         154 :         double dfInvFlattening = 0;
    1247         154 :         if ((dfSemiMajor - dfSemiMinor) >= 0.00000001)
    1248             :         {
    1249          76 :             dfInvFlattening = dfSemiMajor / (dfSemiMajor - dfSemiMinor);
    1250             :         }
    1251             : 
    1252             :         //(if stereographic with center lat ==90) or (polar stereographic )
    1253         308 :         if ((EQUAL(osProjName, "STEREOGRAPHIC") && fabs(dfCenterLat) == 90) ||
    1254         154 :             (EQUAL(osProjName, "POLAR STEREOGRAPHIC")))
    1255             :         {
    1256           2 :             if (bIsOgraphic)
    1257             :             {
    1258           0 :                 oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
    1259             :                                dfSemiMajor, dfInvFlattening,
    1260             :                                "Reference_Meridian", 0.0);
    1261             :             }
    1262             :             else
    1263             :             {
    1264           2 :                 osSphereName += "_polarRadius";
    1265           2 :                 oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
    1266             :                                dfPolarRadius, 0.0, "Reference_Meridian", 0.0);
    1267             :             }
    1268             :         }
    1269         152 :         else if ((EQUAL(osProjName, "EQUIRECTANGULAR")) ||
    1270         103 :                  (EQUAL(osProjName, "ORTHOGRAPHIC")) ||
    1271         356 :                  (EQUAL(osProjName, "STEREOGRAPHIC")) ||
    1272         101 :                  (EQUAL(osProjName, "SINUSOIDAL")))
    1273             :         {
    1274          55 :             oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName, dfSemiMajor,
    1275             :                            0.0, "Reference_Meridian", 0.0);
    1276             :         }
    1277             :         else
    1278             :         {
    1279          97 :             if (bIsOgraphic)
    1280             :             {
    1281           2 :                 oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
    1282             :                                dfSemiMajor, dfInvFlattening,
    1283             :                                "Reference_Meridian", 0.0);
    1284             :             }
    1285             :             else
    1286             :             {
    1287          95 :                 oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
    1288             :                                dfSemiMajor, 0.0, "Reference_Meridian", 0.0);
    1289             :             }
    1290             :         }
    1291             :     }
    1292             : 
    1293             :     CPLXMLNode *psPCI =
    1294         154 :         CPLGetXMLNode(psSR, "Planar.Planar_Coordinate_Information");
    1295         154 :     CPLXMLNode *psGT = CPLGetXMLNode(psSR, "Planar.Geo_Transformation");
    1296         154 :     if (psPCI && psGT)
    1297             :     {
    1298             :         const char *pszPCIEncoding =
    1299         150 :             CPLGetXMLValue(psPCI, "planar_coordinate_encoding_method", "");
    1300         150 :         CPLXMLNode *psCR = CPLGetXMLNode(psPCI, "Coordinate_Representation");
    1301         150 :         if (!EQUAL(pszPCIEncoding, "Coordinate Pair"))
    1302             :         {
    1303           0 :             CPLError(CE_Warning, CPLE_NotSupported,
    1304             :                      "planar_coordinate_encoding_method = %s not supported",
    1305             :                      pszPCIEncoding);
    1306             :         }
    1307         150 :         else if (psCR != nullptr)
    1308             :         {
    1309         150 :             double dfXRes = GetResolutionValue(psCR, "pixel_resolution_x");
    1310         150 :             double dfYRes = GetResolutionValue(psCR, "pixel_resolution_y");
    1311         150 :             double dfULX = GetLinearValue(psGT, "upperleft_corner_x");
    1312         150 :             double dfULY = GetLinearValue(psGT, "upperleft_corner_y");
    1313             : 
    1314             :             // The PDS4 specification is not really clear about the
    1315             :             // origin convention, but it appears from
    1316             :             // https://github.com/OSGeo/gdal/issues/735 that it matches GDAL
    1317             :             // top-left corner of top-left pixel
    1318         150 :             m_gt[0] = dfULX;
    1319         150 :             m_gt[1] = dfXRes;
    1320         150 :             m_gt[2] = 0.0;
    1321         150 :             m_gt[3] = dfULY;
    1322         150 :             m_gt[4] = 0.0;
    1323         150 :             m_gt[5] = -dfYRes;
    1324         150 :             m_bGotTransform = true;
    1325             : 
    1326         150 :             if (dfMapProjectionRotation != 0)
    1327             :             {
    1328           4 :                 const double sin_rot =
    1329             :                     dfMapProjectionRotation == 90
    1330           4 :                         ? 1.0
    1331           0 :                         : sin(dfMapProjectionRotation / 180 * M_PI);
    1332           4 :                 const double cos_rot =
    1333             :                     dfMapProjectionRotation == 90
    1334           4 :                         ? 0.0
    1335           0 :                         : cos(dfMapProjectionRotation / 180 * M_PI);
    1336           4 :                 const double gt_1 = cos_rot * m_gt[1] - sin_rot * m_gt[4];
    1337           4 :                 const double gt_2 = cos_rot * m_gt[2] - sin_rot * m_gt[5];
    1338           4 :                 const double gt_0 = cos_rot * m_gt[0] - sin_rot * m_gt[3];
    1339           4 :                 const double gt_4 = sin_rot * m_gt[1] + cos_rot * m_gt[4];
    1340           4 :                 const double gt_5 = sin_rot * m_gt[2] + cos_rot * m_gt[5];
    1341           4 :                 const double gt_3 = sin_rot * m_gt[0] + cos_rot * m_gt[3];
    1342           4 :                 m_gt[1] = gt_1;
    1343           4 :                 m_gt[2] = gt_2;
    1344           4 :                 m_gt[0] = gt_0;
    1345           4 :                 m_gt[4] = gt_4;
    1346           4 :                 m_gt[5] = gt_5;
    1347           4 :                 m_gt[3] = gt_3;
    1348             :             }
    1349             :         }
    1350             :     }
    1351             : 
    1352         154 :     if (!oSRS.IsEmpty())
    1353             :     {
    1354         154 :         if (GetRasterCount())
    1355             :         {
    1356         150 :             m_oSRS = std::move(oSRS);
    1357             :         }
    1358           4 :         else if (GetLayerCount())
    1359             :         {
    1360           8 :             for (auto &poLayer : m_apoLayers)
    1361             :             {
    1362           4 :                 if (poLayer->GetGeomType() != wkbNone)
    1363             :                 {
    1364           4 :                     auto poSRSClone = oSRS.Clone();
    1365           4 :                     poLayer->SetSpatialRef(poSRSClone);
    1366           4 :                     poSRSClone->Release();
    1367             :                 }
    1368             :             }
    1369             :         }
    1370             :     }
    1371             : }
    1372             : 
    1373             : /************************************************************************/
    1374             : /*                              GetLayer()                              */
    1375             : /************************************************************************/
    1376             : 
    1377         321 : OGRLayer *PDS4Dataset::GetLayer(int nIndex)
    1378             : {
    1379         321 :     if (nIndex < 0 || nIndex >= GetLayerCount())
    1380           8 :         return nullptr;
    1381         313 :     return m_apoLayers[nIndex].get();
    1382             : }
    1383             : 
    1384             : /************************************************************************/
    1385             : /*                       FixupTableFilename()                           */
    1386             : /************************************************************************/
    1387             : 
    1388          97 : static std::string FixupTableFilename(const std::string &osFilename)
    1389             : {
    1390             :     VSIStatBufL sStat;
    1391          97 :     if (VSIStatL(osFilename.c_str(), &sStat) == 0)
    1392             :     {
    1393          97 :         return osFilename;
    1394             :     }
    1395           0 :     const std::string osExt = CPLGetExtensionSafe(osFilename.c_str());
    1396           0 :     if (!osExt.empty())
    1397             :     {
    1398           0 :         std::string osTry(osFilename);
    1399           0 :         if (osExt[0] >= 'a' && osExt[0] <= 'z')
    1400             :         {
    1401           0 :             osTry = CPLResetExtensionSafe(osFilename.c_str(),
    1402           0 :                                           CPLString(osExt).toupper());
    1403             :         }
    1404             :         else
    1405             :         {
    1406           0 :             osTry = CPLResetExtensionSafe(osFilename.c_str(),
    1407           0 :                                           CPLString(osExt).tolower());
    1408             :         }
    1409           0 :         if (VSIStatL(osTry.c_str(), &sStat) == 0)
    1410             :         {
    1411           0 :             return osTry;
    1412             :         }
    1413             :     }
    1414           0 :     return osFilename;
    1415             : }
    1416             : 
    1417             : /************************************************************************/
    1418             : /*                       OpenTableCharacter()                           */
    1419             : /************************************************************************/
    1420             : 
    1421          22 : bool PDS4Dataset::OpenTableCharacter(const char *pszFilename,
    1422             :                                      const CPLXMLNode *psTable)
    1423             : {
    1424          44 :     const std::string osLayerName(CPLGetBasenameSafe(pszFilename));
    1425          44 :     const std::string osFullFilename = FixupTableFilename(CPLFormFilenameSafe(
    1426          66 :         CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(), pszFilename, nullptr));
    1427             :     std::unique_ptr<PDS4TableCharacter> poLayer(new PDS4TableCharacter(
    1428          44 :         this, osLayerName.c_str(), osFullFilename.c_str()));
    1429          22 :     if (!poLayer->ReadTableDef(psTable))
    1430             :     {
    1431           0 :         return false;
    1432             :     }
    1433             :     std::unique_ptr<PDS4EditableLayer> poEditableLayer(
    1434          22 :         new PDS4EditableLayer(poLayer.release()));
    1435          22 :     m_apoLayers.push_back(std::move(poEditableLayer));
    1436          22 :     return true;
    1437             : }
    1438             : 
    1439             : /************************************************************************/
    1440             : /*                       OpenTableBinary()                              */
    1441             : /************************************************************************/
    1442             : 
    1443          11 : bool PDS4Dataset::OpenTableBinary(const char *pszFilename,
    1444             :                                   const CPLXMLNode *psTable)
    1445             : {
    1446          22 :     const std::string osLayerName(CPLGetBasenameSafe(pszFilename));
    1447          22 :     const std::string osFullFilename = FixupTableFilename(CPLFormFilenameSafe(
    1448          33 :         CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(), pszFilename, nullptr));
    1449             :     std::unique_ptr<PDS4TableBinary> poLayer(
    1450          22 :         new PDS4TableBinary(this, osLayerName.c_str(), osFullFilename.c_str()));
    1451          11 :     if (!poLayer->ReadTableDef(psTable))
    1452             :     {
    1453           0 :         return false;
    1454             :     }
    1455             :     std::unique_ptr<PDS4EditableLayer> poEditableLayer(
    1456          11 :         new PDS4EditableLayer(poLayer.release()));
    1457          11 :     m_apoLayers.push_back(std::move(poEditableLayer));
    1458          11 :     return true;
    1459             : }
    1460             : 
    1461             : /************************************************************************/
    1462             : /*                      OpenTableDelimited()                            */
    1463             : /************************************************************************/
    1464             : 
    1465          64 : bool PDS4Dataset::OpenTableDelimited(const char *pszFilename,
    1466             :                                      const CPLXMLNode *psTable)
    1467             : {
    1468         128 :     const std::string osLayerName(CPLGetBasenameSafe(pszFilename));
    1469         128 :     const std::string osFullFilename = FixupTableFilename(CPLFormFilenameSafe(
    1470         192 :         CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(), pszFilename, nullptr));
    1471             :     std::unique_ptr<PDS4DelimitedTable> poLayer(new PDS4DelimitedTable(
    1472         128 :         this, osLayerName.c_str(), osFullFilename.c_str()));
    1473          64 :     if (!poLayer->ReadTableDef(psTable))
    1474             :     {
    1475           0 :         return false;
    1476             :     }
    1477             :     std::unique_ptr<PDS4EditableLayer> poEditableLayer(
    1478          64 :         new PDS4EditableLayer(poLayer.release()));
    1479          64 :     m_apoLayers.push_back(std::move(poEditableLayer));
    1480          64 :     return true;
    1481             : }
    1482             : 
    1483             : /************************************************************************/
    1484             : /*                                Open()                                */
    1485             : /************************************************************************/
    1486             : 
    1487             : // See https://pds.nasa.gov/pds4/pds/v1/PDS4_PDS_1800.xsd
    1488             : // and https://pds.nasa.gov/pds4/pds/v1/PDS4_PDS_1800.sch
    1489         287 : PDS4Dataset *PDS4Dataset::OpenInternal(GDALOpenInfo *poOpenInfo)
    1490             : {
    1491         287 :     if (!PDS4DriverIdentify(poOpenInfo))
    1492           0 :         return nullptr;
    1493             : 
    1494         574 :     CPLString osXMLFilename(poOpenInfo->pszFilename);
    1495         287 :     int nFAOIdxLookup = -1;
    1496         287 :     int nArrayIdxLookup = -1;
    1497         287 :     if (STARTS_WITH_CI(poOpenInfo->pszFilename, "PDS4:"))
    1498             :     {
    1499             :         char **papszTokens =
    1500          15 :             CSLTokenizeString2(poOpenInfo->pszFilename, ":", 0);
    1501          15 :         int nCount = CSLCount(papszTokens);
    1502          15 :         if (nCount == 5 && strlen(papszTokens[1]) == 1 &&
    1503           1 :             (papszTokens[2][0] == '\\' || papszTokens[2][0] == '/'))
    1504             :         {
    1505           1 :             osXMLFilename = CPLString(papszTokens[1]) + ":" + papszTokens[2];
    1506           1 :             nFAOIdxLookup = atoi(papszTokens[3]);
    1507           1 :             nArrayIdxLookup = atoi(papszTokens[4]);
    1508             :         }
    1509          14 :         else if (nCount == 5 && (EQUAL(papszTokens[1], "/vsicurl/http") ||
    1510           0 :                                  EQUAL(papszTokens[1], "/vsicurl/https")))
    1511             :         {
    1512           0 :             osXMLFilename = CPLString(papszTokens[1]) + ":" + papszTokens[2];
    1513           0 :             nFAOIdxLookup = atoi(papszTokens[3]);
    1514           0 :             nArrayIdxLookup = atoi(papszTokens[4]);
    1515             :         }
    1516          14 :         else if (nCount == 4)
    1517             :         {
    1518          13 :             osXMLFilename = papszTokens[1];
    1519          13 :             nFAOIdxLookup = atoi(papszTokens[2]);
    1520          13 :             nArrayIdxLookup = atoi(papszTokens[3]);
    1521             :         }
    1522             :         else
    1523             :         {
    1524           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1525             :                      "Invalid syntax for PDS4 subdataset name");
    1526           1 :             CSLDestroy(papszTokens);
    1527           1 :             return nullptr;
    1528             :         }
    1529          14 :         CSLDestroy(papszTokens);
    1530             :     }
    1531             : 
    1532         572 :     CPLXMLTreeCloser oCloser(CPLParseXMLFile(osXMLFilename));
    1533         286 :     CPLXMLNode *psRoot = oCloser.get();
    1534         286 :     CPLStripXMLNamespace(psRoot, nullptr, TRUE);
    1535             : 
    1536         572 :     GDALAccess eAccess = STARTS_WITH_CI(poOpenInfo->pszFilename, "PDS4:")
    1537         286 :                              ? GA_ReadOnly
    1538             :                              : poOpenInfo->eAccess;
    1539             : 
    1540         286 :     CPLXMLNode *psProduct = CPLGetXMLNode(psRoot, "=Product_Observational");
    1541         286 :     if (psProduct == nullptr)
    1542             :     {
    1543           4 :         eAccess = GA_ReadOnly;
    1544           4 :         psProduct = CPLGetXMLNode(psRoot, "=Product_Ancillary");
    1545           4 :         if (psProduct == nullptr)
    1546             :         {
    1547           4 :             psProduct = CPLGetXMLNode(psRoot, "=Product_Collection");
    1548             :         }
    1549             :     }
    1550         286 :     if (psProduct == nullptr)
    1551             :     {
    1552           3 :         return nullptr;
    1553             :     }
    1554             : 
    1555             :     // Test case:
    1556             :     // https://starbase.jpl.nasa.gov/pds4/1700/dph_example_products/test_Images_DisplaySettings/TestPattern_Image/TestPattern.xml
    1557         283 :     const char *pszVertDir = CPLGetXMLValue(
    1558             :         psProduct,
    1559             :         "Observation_Area.Discipline_Area.Display_Settings.Display_Direction."
    1560             :         "vertical_display_direction",
    1561             :         "");
    1562         283 :     const bool bBottomToTop = EQUAL(pszVertDir, "Bottom to Top");
    1563             : 
    1564         283 :     const char *pszHorizDir = CPLGetXMLValue(
    1565             :         psProduct,
    1566             :         "Observation_Area.Discipline_Area.Display_Settings.Display_Direction."
    1567             :         "horizontal_display_direction",
    1568             :         "");
    1569         283 :     const bool bRightToLeft = EQUAL(pszHorizDir, "Right to Left");
    1570             : 
    1571         566 :     auto poDS = std::make_unique<PDS4Dataset>();
    1572         283 :     poDS->m_osXMLFilename = osXMLFilename;
    1573         283 :     poDS->eAccess = eAccess;
    1574         283 :     poDS->papszOpenOptions = CSLDuplicate(poOpenInfo->papszOpenOptions);
    1575             : 
    1576         566 :     CPLStringList aosSubdatasets;
    1577         283 :     int nFAOIdx = 0;
    1578        2765 :     for (CPLXMLNode *psIter = psProduct->psChild; psIter != nullptr;
    1579        2482 :          psIter = psIter->psNext)
    1580             :     {
    1581        2484 :         if (psIter->eType != CXT_Element ||
    1582         878 :             (strcmp(psIter->pszValue, "File_Area_Observational") != 0 &&
    1583         566 :              strcmp(psIter->pszValue, "File_Area_Ancillary") != 0 &&
    1584         566 :              strcmp(psIter->pszValue, "File_Area_Inventory") != 0))
    1585             :         {
    1586        2171 :             continue;
    1587             :         }
    1588             : 
    1589         313 :         nFAOIdx++;
    1590         313 :         CPLXMLNode *psFile = CPLGetXMLNode(psIter, "File");
    1591         313 :         if (psFile == nullptr)
    1592             :         {
    1593           1 :             continue;
    1594             :         }
    1595         312 :         const char *pszFilename = CPLGetXMLValue(psFile, "file_name", nullptr);
    1596         312 :         if (pszFilename == nullptr)
    1597             :         {
    1598           1 :             continue;
    1599             :         }
    1600             : 
    1601         655 :         for (CPLXMLNode *psSubIter = psFile->psChild; psSubIter != nullptr;
    1602         344 :              psSubIter = psSubIter->psNext)
    1603             :         {
    1604         344 :             if (psSubIter->eType == CXT_Comment &&
    1605          14 :                 EQUAL(psSubIter->pszValue, PREEXISTING_BINARY_FILE))
    1606             :             {
    1607          14 :                 poDS->m_bCreatedFromExistingBinaryFile = true;
    1608             :             }
    1609             :         }
    1610             : 
    1611         311 :         int nArrayIdx = 0;
    1612         311 :         for (CPLXMLNode *psSubIter = psIter->psChild;
    1613         986 :              (nFAOIdxLookup < 0 || nFAOIdxLookup == nFAOIdx) &&
    1614             :              psSubIter != nullptr;
    1615         675 :              psSubIter = psSubIter->psNext)
    1616             :         {
    1617         677 :             if (psSubIter->eType != CXT_Element)
    1618             :             {
    1619         478 :                 continue;
    1620             :             }
    1621         677 :             int nDIM = 0;
    1622         677 :             if (STARTS_WITH(psSubIter->pszValue, "Array_1D"))
    1623             :             {
    1624           0 :                 nDIM = 1;
    1625             :             }
    1626         677 :             else if (STARTS_WITH(psSubIter->pszValue, "Array_2D"))
    1627             :             {
    1628           6 :                 nDIM = 2;
    1629             :             }
    1630         671 :             else if (STARTS_WITH(psSubIter->pszValue, "Array_3D"))
    1631             :             {
    1632         222 :                 nDIM = 3;
    1633             :             }
    1634         449 :             else if (strcmp(psSubIter->pszValue, "Array") == 0)
    1635             :             {
    1636           3 :                 nDIM = atoi(CPLGetXMLValue(psSubIter, "axes", "0"));
    1637             :             }
    1638         446 :             else if (strcmp(psSubIter->pszValue, "Table_Character") == 0)
    1639             :             {
    1640          22 :                 poDS->OpenTableCharacter(pszFilename, psSubIter);
    1641          22 :                 continue;
    1642             :             }
    1643         424 :             else if (strcmp(psSubIter->pszValue, "Table_Binary") == 0)
    1644             :             {
    1645          11 :                 poDS->OpenTableBinary(pszFilename, psSubIter);
    1646          11 :                 continue;
    1647             :             }
    1648         413 :             else if (strcmp(psSubIter->pszValue, "Table_Delimited") == 0 ||
    1649         350 :                      strcmp(psSubIter->pszValue, "Inventory") == 0)
    1650             :             {
    1651          64 :                 poDS->OpenTableDelimited(pszFilename, psSubIter);
    1652          64 :                 continue;
    1653             :             }
    1654         580 :             if (nDIM == 0)
    1655             :             {
    1656         349 :                 continue;
    1657             :             }
    1658             : 
    1659         231 :             nArrayIdx++;
    1660             :             // Does it match a selected subdataset ?
    1661         231 :             if (nArrayIdxLookup > 0 && nArrayIdx != nArrayIdxLookup)
    1662             :             {
    1663          13 :                 continue;
    1664             :             }
    1665             : 
    1666             :             const char *pszArrayName =
    1667         218 :                 CPLGetXMLValue(psSubIter, "name", nullptr);
    1668             :             const char *pszArrayId =
    1669         218 :                 CPLGetXMLValue(psSubIter, "local_identifier", nullptr);
    1670             :             vsi_l_offset nOffset = static_cast<vsi_l_offset>(
    1671         218 :                 CPLAtoGIntBig(CPLGetXMLValue(psSubIter, "offset", "0")));
    1672             : 
    1673             :             const char *pszAxisIndexOrder =
    1674         218 :                 CPLGetXMLValue(psSubIter, "axis_index_order", "");
    1675         218 :             if (!EQUAL(pszAxisIndexOrder, "Last Index Fastest"))
    1676             :             {
    1677           1 :                 CPLError(CE_Warning, CPLE_NotSupported,
    1678             :                          "axis_index_order = '%s' unhandled",
    1679             :                          pszAxisIndexOrder);
    1680           1 :                 continue;
    1681             :             }
    1682             : 
    1683             :             // Figure out data type
    1684             :             const char *pszDataType =
    1685         217 :                 CPLGetXMLValue(psSubIter, "Element_Array.data_type", "");
    1686         217 :             GDALDataType eDT = GDT_Byte;
    1687         217 :             bool bLSBOrder = strstr(pszDataType, "LSB") != nullptr;
    1688             : 
    1689             :             // ComplexLSB16', 'ComplexLSB8', 'ComplexMSB16', 'ComplexMSB8',
    1690             :             // 'IEEE754LSBDouble', 'IEEE754LSBSingle', 'IEEE754MSBDouble',
    1691             :             // 'IEEE754MSBSingle', 'SignedBitString', 'SignedByte',
    1692             :             // 'SignedLSB2', 'SignedLSB4', 'SignedLSB8', 'SignedMSB2',
    1693             :             // 'SignedMSB4', 'SignedMSB8', 'UnsignedBitString', 'UnsignedByte',
    1694             :             // 'UnsignedLSB2', 'UnsignedLSB4', 'UnsignedLSB8', 'UnsignedMSB2',
    1695             :             // 'UnsignedMSB4', 'UnsignedMSB8'
    1696         217 :             if (EQUAL(pszDataType, "ComplexLSB16") ||
    1697         213 :                 EQUAL(pszDataType, "ComplexMSB16"))
    1698             :             {
    1699           6 :                 eDT = GDT_CFloat64;
    1700             :             }
    1701         211 :             else if (EQUAL(pszDataType, "ComplexLSB8") ||
    1702         207 :                      EQUAL(pszDataType, "ComplexMSB8"))
    1703             :             {
    1704           4 :                 eDT = GDT_CFloat32;
    1705             :             }
    1706         207 :             else if (EQUAL(pszDataType, "IEEE754LSBDouble") ||
    1707         203 :                      EQUAL(pszDataType, "IEEE754MSBDouble"))
    1708             :             {
    1709           4 :                 eDT = GDT_Float64;
    1710             :             }
    1711         203 :             else if (EQUAL(pszDataType, "IEEE754LSBSingle") ||
    1712         195 :                      EQUAL(pszDataType, "IEEE754MSBSingle"))
    1713             :             {
    1714           8 :                 eDT = GDT_Float32;
    1715             :             }
    1716             :             // SignedBitString unhandled
    1717         195 :             else if (EQUAL(pszDataType, "SignedByte"))
    1718             :             {
    1719           6 :                 eDT = GDT_Int8;
    1720             :             }
    1721         189 :             else if (EQUAL(pszDataType, "SignedLSB2") ||
    1722         185 :                      EQUAL(pszDataType, "SignedMSB2"))
    1723             :             {
    1724           6 :                 eDT = GDT_Int16;
    1725             :             }
    1726         183 :             else if (EQUAL(pszDataType, "SignedLSB4") ||
    1727         179 :                      EQUAL(pszDataType, "SignedMSB4"))
    1728             :             {
    1729           4 :                 eDT = GDT_Int32;
    1730             :             }
    1731             :             // SignedLSB8 and SignedMSB8 unhandled
    1732         179 :             else if (EQUAL(pszDataType, "UnsignedByte"))
    1733             :             {
    1734         168 :                 eDT = GDT_Byte;
    1735             :             }
    1736          11 :             else if (EQUAL(pszDataType, "UnsignedLSB2") ||
    1737           5 :                      EQUAL(pszDataType, "UnsignedMSB2"))
    1738             :             {
    1739           6 :                 eDT = GDT_UInt16;
    1740             :             }
    1741           5 :             else if (EQUAL(pszDataType, "UnsignedLSB4") ||
    1742           1 :                      EQUAL(pszDataType, "UnsignedMSB4"))
    1743             :             {
    1744           4 :                 eDT = GDT_UInt32;
    1745             :             }
    1746             :             // UnsignedLSB8 and UnsignedMSB8 unhandled
    1747             :             else
    1748             :             {
    1749           1 :                 CPLDebug("PDS4", "data_type = '%s' unhandled", pszDataType);
    1750           1 :                 continue;
    1751             :             }
    1752             : 
    1753         216 :             poDS->m_osUnits =
    1754         432 :                 CPLGetXMLValue(psSubIter, "Element_Array.unit", "");
    1755             : 
    1756         216 :             double dfValueOffset = CPLAtof(
    1757             :                 CPLGetXMLValue(psSubIter, "Element_Array.value_offset", "0"));
    1758         216 :             double dfValueScale = CPLAtof(
    1759             :                 CPLGetXMLValue(psSubIter, "Element_Array.scaling_factor", "1"));
    1760             : 
    1761             :             // Parse Axis_Array elements
    1762         216 :             char szOrder[4] = {0};
    1763         216 :             int l_nBands = 1;
    1764         216 :             int nLines = 0;
    1765         216 :             int nSamples = 0;
    1766         216 :             int nAxisFound = 0;
    1767         216 :             int anElements[3] = {0};
    1768         216 :             for (CPLXMLNode *psAxisIter = psSubIter->psChild;
    1769        1980 :                  psAxisIter != nullptr; psAxisIter = psAxisIter->psNext)
    1770             :             {
    1771        1764 :                 if (psAxisIter->eType != CXT_Element ||
    1772        1764 :                     strcmp(psAxisIter->pszValue, "Axis_Array") != 0)
    1773             :                 {
    1774        1118 :                     continue;
    1775             :                 }
    1776             :                 const char *pszAxisName =
    1777         646 :                     CPLGetXMLValue(psAxisIter, "axis_name", nullptr);
    1778             :                 const char *pszElements =
    1779         646 :                     CPLGetXMLValue(psAxisIter, "elements", nullptr);
    1780             :                 const char *pszSequenceNumber =
    1781         646 :                     CPLGetXMLValue(psAxisIter, "sequence_number", nullptr);
    1782         646 :                 if (pszAxisName == nullptr || pszElements == nullptr ||
    1783             :                     pszSequenceNumber == nullptr)
    1784             :                 {
    1785           1 :                     continue;
    1786             :                 }
    1787         645 :                 int nSeqNumber = atoi(pszSequenceNumber);
    1788         645 :                 if (nSeqNumber < 1 || nSeqNumber > nDIM)
    1789             :                 {
    1790           2 :                     CPLError(CE_Warning, CPLE_AppDefined,
    1791             :                              "Invalid sequence_number = %s", pszSequenceNumber);
    1792           2 :                     continue;
    1793             :                 }
    1794         643 :                 int nElements = atoi(pszElements);
    1795         643 :                 if (nElements <= 0)
    1796             :                 {
    1797           1 :                     CPLError(CE_Warning, CPLE_AppDefined,
    1798             :                              "Invalid elements = %s", pszElements);
    1799           1 :                     continue;
    1800             :                 }
    1801         642 :                 nSeqNumber--;
    1802         642 :                 if (szOrder[nSeqNumber] != '\0')
    1803             :                 {
    1804           1 :                     CPLError(CE_Warning, CPLE_AppDefined,
    1805             :                              "Invalid sequence_number = %s", pszSequenceNumber);
    1806           1 :                     continue;
    1807             :                 }
    1808         641 :                 if (EQUAL(pszAxisName, "Band") && nDIM == 3)
    1809             :                 {
    1810         210 :                     szOrder[nSeqNumber] = 'B';
    1811         210 :                     l_nBands = nElements;
    1812         210 :                     anElements[nSeqNumber] = nElements;
    1813         210 :                     nAxisFound++;
    1814             :                 }
    1815         431 :                 else if (EQUAL(pszAxisName, "Line"))
    1816             :                 {
    1817         215 :                     szOrder[nSeqNumber] = 'L';
    1818         215 :                     nLines = nElements;
    1819         215 :                     anElements[nSeqNumber] = nElements;
    1820         215 :                     nAxisFound++;
    1821             :                 }
    1822         216 :                 else if (EQUAL(pszAxisName, "Sample"))
    1823             :                 {
    1824         215 :                     szOrder[nSeqNumber] = 'S';
    1825         215 :                     nSamples = nElements;
    1826         215 :                     anElements[nSeqNumber] = nElements;
    1827         215 :                     nAxisFound++;
    1828             :                 }
    1829             :                 else
    1830             :                 {
    1831           1 :                     CPLError(CE_Warning, CPLE_NotSupported,
    1832             :                              "Unsupported axis_name = %s", pszAxisName);
    1833           1 :                     continue;
    1834             :                 }
    1835             :             }
    1836         216 :             if (nAxisFound != nDIM)
    1837             :             {
    1838           1 :                 CPLError(CE_Warning, CPLE_AppDefined,
    1839             :                          "Found only %d Axis_Array elements. %d expected",
    1840             :                          nAxisFound, nDIM);
    1841           1 :                 continue;
    1842             :             }
    1843             : 
    1844         430 :             if (!GDALCheckDatasetDimensions(nSamples, nLines) ||
    1845         215 :                 !GDALCheckBandCount(l_nBands, FALSE))
    1846             :             {
    1847           1 :                 continue;
    1848             :             }
    1849             : 
    1850             :             // Compute pixel, line and band spacing
    1851         214 :             vsi_l_offset nSpacing = GDALGetDataTypeSizeBytes(eDT);
    1852         214 :             int nPixelOffset = 0;
    1853         214 :             int nLineOffset = 0;
    1854         214 :             vsi_l_offset nBandOffset = 0;
    1855         214 :             int nCountPreviousDim = 1;
    1856         848 :             for (int i = nDIM - 1; i >= 0; i--)
    1857             :             {
    1858         636 :                 if (szOrder[i] == 'S')
    1859             :                 {
    1860         214 :                     if (nSpacing >
    1861         214 :                         static_cast<vsi_l_offset>(INT_MAX / nCountPreviousDim))
    1862             :                     {
    1863           1 :                         CPLError(CE_Failure, CPLE_NotSupported,
    1864             :                                  "Integer overflow");
    1865           2 :                         return nullptr;
    1866             :                     }
    1867         213 :                     nPixelOffset =
    1868         213 :                         static_cast<int>(nSpacing * nCountPreviousDim);
    1869         213 :                     nSpacing = nPixelOffset;
    1870             :                 }
    1871         422 :                 else if (szOrder[i] == 'L')
    1872             :                 {
    1873         214 :                     if (nSpacing >
    1874         214 :                         static_cast<vsi_l_offset>(INT_MAX / nCountPreviousDim))
    1875             :                     {
    1876           1 :                         CPLError(CE_Failure, CPLE_NotSupported,
    1877             :                                  "Integer overflow");
    1878           1 :                         return nullptr;
    1879             :                     }
    1880         213 :                     nLineOffset =
    1881         213 :                         static_cast<int>(nSpacing * nCountPreviousDim);
    1882         213 :                     nSpacing = nLineOffset;
    1883             :                 }
    1884             :                 else
    1885             :                 {
    1886         208 :                     nBandOffset = nSpacing * nCountPreviousDim;
    1887         208 :                     nSpacing = nBandOffset;
    1888             :                 }
    1889         634 :                 nCountPreviousDim = anElements[i];
    1890             :             }
    1891             : 
    1892             :             // Retrieve no data value
    1893         212 :             bool bNoDataSet = false;
    1894         212 :             double dfNoData = 0.0;
    1895         212 :             std::vector<double> adfConstants;
    1896         212 :             CPLXMLNode *psSC = CPLGetXMLNode(psSubIter, "Special_Constants");
    1897         212 :             if (psSC)
    1898             :             {
    1899             :                 const char *pszMC =
    1900          60 :                     CPLGetXMLValue(psSC, "missing_constant", nullptr);
    1901          60 :                 if (pszMC)
    1902             :                 {
    1903          60 :                     bNoDataSet = true;
    1904          60 :                     dfNoData = CPLAtof(pszMC);
    1905             :                 }
    1906             : 
    1907          60 :                 const char *apszConstantNames[] = {
    1908             :                     "saturated_constant",
    1909             :                     "missing_constant",
    1910             :                     "error_constant",
    1911             :                     "invalid_constant",
    1912             :                     "unknown_constant",
    1913             :                     "not_applicable_constant",
    1914             :                     "high_instrument_saturation",
    1915             :                     "high_representation_saturation",
    1916             :                     "low_instrument_saturation",
    1917             :                     "low_representation_saturation"};
    1918         660 :                 for (size_t i = 0; i < CPL_ARRAYSIZE(apszConstantNames); i++)
    1919             :                 {
    1920             :                     const char *pszConstant =
    1921         600 :                         CPLGetXMLValue(psSC, apszConstantNames[i], nullptr);
    1922         600 :                     if (pszConstant)
    1923          94 :                         adfConstants.push_back(CPLAtof(pszConstant));
    1924             :                 }
    1925             :             }
    1926             : 
    1927             :             // Add subdatasets
    1928         212 :             const int nSDSIdx = 1 + aosSubdatasets.size() / 2;
    1929             :             aosSubdatasets.SetNameValue(
    1930             :                 CPLSPrintf("SUBDATASET_%d_NAME", nSDSIdx),
    1931             :                 CPLSPrintf("PDS4:%s:%d:%d", osXMLFilename.c_str(), nFAOIdx,
    1932         212 :                            nArrayIdx));
    1933             :             aosSubdatasets.SetNameValue(
    1934             :                 CPLSPrintf("SUBDATASET_%d_DESC", nSDSIdx),
    1935             :                 CPLSPrintf("Image file %s, array %s", pszFilename,
    1936             :                            pszArrayName ? pszArrayName
    1937         212 :                            : pszArrayId ? pszArrayId
    1938         424 :                                         : CPLSPrintf("%d", nArrayIdx)));
    1939             : 
    1940         212 :             if (poDS->nBands != 0)
    1941          14 :                 continue;
    1942             : 
    1943             :             const std::string osImageFullFilename = CPLFormFilenameSafe(
    1944         198 :                 CPLGetPathSafe(osXMLFilename.c_str()).c_str(), pszFilename,
    1945         198 :                 nullptr);
    1946         198 :             VSILFILE *fp = VSIFOpenExL(
    1947             :                 osImageFullFilename.c_str(),
    1948         198 :                 (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb", true);
    1949         198 :             if (fp == nullptr)
    1950             :             {
    1951           1 :                 CPLError(CE_Warning, CPLE_FileIO, "Cannt open %s: %s",
    1952             :                          osImageFullFilename.c_str(), VSIGetLastErrorMsg());
    1953           1 :                 continue;
    1954             :             }
    1955         197 :             poDS->nRasterXSize = nSamples;
    1956         197 :             poDS->nRasterYSize = nLines;
    1957         197 :             poDS->m_osImageFilename = osImageFullFilename;
    1958         197 :             poDS->m_fpImage = fp;
    1959         197 :             poDS->m_bIsLSB = bLSBOrder;
    1960             : 
    1961         197 :             if (memcmp(szOrder, "BLS", 3) == 0)
    1962             :             {
    1963         174 :                 poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "BAND",
    1964             :                                                    "IMAGE_STRUCTURE");
    1965             :             }
    1966          23 :             else if (memcmp(szOrder, "LSB", 3) == 0)
    1967             :             {
    1968          18 :                 poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
    1969             :                                                    "IMAGE_STRUCTURE");
    1970             :             }
    1971             : 
    1972         197 :             CPLXMLNode *psOS = CPLGetXMLNode(psSubIter, "Object_Statistics");
    1973         197 :             const char *pszMin = nullptr;
    1974         197 :             const char *pszMax = nullptr;
    1975         197 :             const char *pszMean = nullptr;
    1976         197 :             const char *pszStdDev = nullptr;
    1977         197 :             if (psOS)
    1978             :             {
    1979          17 :                 pszMin = CPLGetXMLValue(psOS, "minimum", nullptr);
    1980          17 :                 pszMax = CPLGetXMLValue(psOS, "maximum", nullptr);
    1981          17 :                 pszMean = CPLGetXMLValue(psOS, "mean", nullptr);
    1982          17 :                 pszStdDev = CPLGetXMLValue(psOS, "standard_deviation", nullptr);
    1983             :             }
    1984             : 
    1985         457 :             for (int i = 0; i < l_nBands; i++)
    1986             :             {
    1987         260 :                 vsi_l_offset nThisBandOffset = nOffset + nBandOffset * i;
    1988         260 :                 if (bBottomToTop)
    1989             :                 {
    1990           2 :                     nThisBandOffset +=
    1991           2 :                         static_cast<vsi_l_offset>(nLines - 1) * nLineOffset;
    1992             :                 }
    1993         260 :                 if (bRightToLeft)
    1994             :                 {
    1995           1 :                     nThisBandOffset +=
    1996           1 :                         static_cast<vsi_l_offset>(nSamples - 1) * nPixelOffset;
    1997             :                 }
    1998             :                 auto poBand = std::make_unique<PDS4RawRasterBand>(
    1999         260 :                     poDS.get(), i + 1, poDS->m_fpImage, nThisBandOffset,
    2000         260 :                     bRightToLeft ? -nPixelOffset : nPixelOffset,
    2001         260 :                     bBottomToTop ? -nLineOffset : nLineOffset, eDT,
    2002         260 :                     bLSBOrder ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
    2003         260 :                               : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN);
    2004         260 :                 if (!poBand->IsValid())
    2005           0 :                     return nullptr;
    2006         260 :                 if (bNoDataSet)
    2007             :                 {
    2008          62 :                     poBand->SetNoDataValue(dfNoData);
    2009             :                 }
    2010         260 :                 poBand->SetOffset(dfValueOffset);
    2011         260 :                 poBand->SetScale(dfValueScale);
    2012             : 
    2013         260 :                 if (l_nBands == 1)
    2014             :                 {
    2015         165 :                     if (pszMin)
    2016             :                     {
    2017          17 :                         poBand->GDALRasterBand::SetMetadataItem(
    2018             :                             "STATISTICS_MINIMUM", pszMin);
    2019             :                     }
    2020         165 :                     if (pszMax)
    2021             :                     {
    2022          17 :                         poBand->GDALRasterBand::SetMetadataItem(
    2023             :                             "STATISTICS_MAXIMUM", pszMax);
    2024             :                     }
    2025         165 :                     if (pszMean)
    2026             :                     {
    2027          17 :                         poBand->GDALRasterBand::SetMetadataItem(
    2028             :                             "STATISTICS_MEAN", pszMean);
    2029             :                     }
    2030         165 :                     if (pszStdDev)
    2031             :                     {
    2032          17 :                         poBand->GDALRasterBand::SetMetadataItem(
    2033             :                             "STATISTICS_STDDEV", pszStdDev);
    2034             :                     }
    2035             :                 }
    2036             : 
    2037             :                 // Only instantiate explicit mask band if we have at least one
    2038             :                 // special constant (that is not the missing_constant,
    2039             :                 // already exposed as nodata value)
    2040         508 :                 if (!GDALDataTypeIsComplex(eDT) &&
    2041         488 :                     (CPLTestBool(CPLGetConfigOption("PDS4_FORCE_MASK", "NO")) ||
    2042         449 :                      adfConstants.size() >= 2 ||
    2043         240 :                      (adfConstants.size() == 1 && !bNoDataSet)))
    2044             :                 {
    2045          78 :                     poBand->SetMaskBand(std::make_unique<PDS4MaskBand>(
    2046          78 :                         poBand.get(), adfConstants));
    2047             :                 }
    2048             : 
    2049         260 :                 poDS->SetBand(i + 1, std::move(poBand));
    2050             :             }
    2051             :         }
    2052             :     }
    2053             : 
    2054         281 :     if (nFAOIdxLookup < 0 && aosSubdatasets.size() > 2)
    2055             :     {
    2056          10 :         poDS->GDALDataset::SetMetadata(aosSubdatasets.List(), "SUBDATASETS");
    2057             :     }
    2058         271 :     else if (poDS->nBands == 0 &&
    2059         281 :              (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0 &&
    2060          10 :              (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) == 0)
    2061             :     {
    2062           5 :         return nullptr;
    2063             :     }
    2064         266 :     else if (poDS->m_apoLayers.empty() &&
    2065         268 :              (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) != 0 &&
    2066           2 :              (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0)
    2067             :     {
    2068           0 :         return nullptr;
    2069             :     }
    2070             : 
    2071             :     // Expose XML content in xml:PDS4 metadata domain
    2072         276 :     GByte *pabyRet = nullptr;
    2073         276 :     CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, osXMLFilename, &pabyRet, nullptr,
    2074             :                                      10 * 1024 * 1024));
    2075         276 :     if (pabyRet)
    2076             :     {
    2077         276 :         char *apszXML[2] = {reinterpret_cast<char *>(pabyRet), nullptr};
    2078         276 :         poDS->GDALDataset::SetMetadata(apszXML, "xml:PDS4");
    2079             :     }
    2080         276 :     VSIFree(pabyRet);
    2081             : 
    2082             :     /*--------------------------------------------------------------------------*/
    2083             :     /*  Parse georeferencing info */
    2084             :     /*--------------------------------------------------------------------------*/
    2085         276 :     poDS->ReadGeoreferencing(psProduct);
    2086             : 
    2087             :     /*--------------------------------------------------------------------------*/
    2088             :     /*  Check for overviews */
    2089             :     /*--------------------------------------------------------------------------*/
    2090         276 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
    2091             : 
    2092             :     /*--------------------------------------------------------------------------*/
    2093             :     /*  Initialize any PAM information */
    2094             :     /*--------------------------------------------------------------------------*/
    2095         276 :     poDS->SetDescription(poOpenInfo->pszFilename);
    2096         276 :     poDS->TryLoadXML();
    2097             : 
    2098         276 :     return poDS.release();
    2099             : }
    2100             : 
    2101             : /************************************************************************/
    2102             : /*                         IsCARTVersionGTE()                           */
    2103             : /************************************************************************/
    2104             : 
    2105             : // Returns true is pszCur >= pszRef
    2106             : // Must be things like 1900, 1B00, 1D00_1933 ...
    2107         368 : static bool IsCARTVersionGTE(const char *pszCur, const char *pszRef)
    2108             : {
    2109         368 :     return strcmp(pszCur, pszRef) >= 0;
    2110             : }
    2111             : 
    2112             : /************************************************************************/
    2113             : /*                         WriteGeoreferencing()                        */
    2114             : /************************************************************************/
    2115             : 
    2116          92 : void PDS4Dataset::WriteGeoreferencing(CPLXMLNode *psCart,
    2117             :                                       const char *pszCARTVersion)
    2118             : {
    2119          92 :     bool bHasBoundingBox = false;
    2120          92 :     double adfX[4] = {0};
    2121          92 :     double adfY[4] = {0};
    2122          92 :     CPLString osPrefix;
    2123          92 :     const char *pszColon = strchr(psCart->pszValue, ':');
    2124          92 :     if (pszColon)
    2125          92 :         osPrefix.assign(psCart->pszValue, pszColon - psCart->pszValue + 1);
    2126             : 
    2127          92 :     if (m_bGotTransform)
    2128             :     {
    2129          90 :         bHasBoundingBox = true;
    2130             : 
    2131             :         // upper left
    2132          90 :         adfX[0] = m_gt[0];
    2133          90 :         adfY[0] = m_gt[3];
    2134             : 
    2135             :         // upper right
    2136          90 :         adfX[1] = m_gt[0] + m_gt[1] * nRasterXSize;
    2137          90 :         adfY[1] = m_gt[3];
    2138             : 
    2139             :         // lower left
    2140          90 :         adfX[2] = m_gt[0];
    2141          90 :         adfY[2] = m_gt[3] + m_gt[5] * nRasterYSize;
    2142             : 
    2143             :         // lower right
    2144          90 :         adfX[3] = m_gt[0] + m_gt[1] * nRasterXSize;
    2145          90 :         adfY[3] = m_gt[3] + m_gt[5] * nRasterYSize;
    2146             :     }
    2147             :     else
    2148             :     {
    2149           2 :         OGRLayer *poLayer = GetLayer(0);
    2150           2 :         OGREnvelope sEnvelope;
    2151           2 :         if (poLayer->GetExtent(&sEnvelope) == OGRERR_NONE)
    2152             :         {
    2153           1 :             bHasBoundingBox = true;
    2154             : 
    2155           1 :             adfX[0] = sEnvelope.MinX;
    2156           1 :             adfY[0] = sEnvelope.MaxY;
    2157             : 
    2158           1 :             adfX[1] = sEnvelope.MaxX;
    2159           1 :             adfY[1] = sEnvelope.MaxY;
    2160             : 
    2161           1 :             adfX[2] = sEnvelope.MinX;
    2162           1 :             adfY[2] = sEnvelope.MinY;
    2163             : 
    2164           1 :             adfX[3] = sEnvelope.MaxX;
    2165           1 :             adfY[3] = sEnvelope.MinY;
    2166             :         }
    2167             :     }
    2168             : 
    2169          92 :     if (bHasBoundingBox && !m_oSRS.IsGeographic())
    2170             :     {
    2171          31 :         bHasBoundingBox = false;
    2172          31 :         OGRSpatialReference *poSRSLongLat = m_oSRS.CloneGeogCS();
    2173          31 :         if (poSRSLongLat)
    2174             :         {
    2175          31 :             poSRSLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2176             :             OGRCoordinateTransformation *poCT =
    2177          31 :                 OGRCreateCoordinateTransformation(&m_oSRS, poSRSLongLat);
    2178          31 :             if (poCT)
    2179             :             {
    2180          31 :                 if (poCT->Transform(4, adfX, adfY))
    2181             :                 {
    2182          31 :                     bHasBoundingBox = true;
    2183             :                 }
    2184          31 :                 delete poCT;
    2185             :             }
    2186          31 :             delete poSRSLongLat;
    2187             :         }
    2188             :     }
    2189             : 
    2190          92 :     if (!bHasBoundingBox)
    2191             :     {
    2192             :         // Write dummy values
    2193           1 :         adfX[0] = -180.0;
    2194           1 :         adfY[0] = 90.0;
    2195           1 :         adfX[1] = 180.0;
    2196           1 :         adfY[1] = 90.0;
    2197           1 :         adfX[2] = -180.0;
    2198           1 :         adfY[2] = -90.0;
    2199           1 :         adfX[3] = 180.0;
    2200           1 :         adfY[3] = -90.0;
    2201             :     }
    2202             : 
    2203         184 :     const char *pszLongitudeDirection = CSLFetchNameValueDef(
    2204          92 :         m_papszCreationOptions, "LONGITUDE_DIRECTION", "Positive East");
    2205          92 :     const double dfLongitudeMultiplier =
    2206          92 :         EQUAL(pszLongitudeDirection, "Positive West") ? -1 : 1;
    2207         200 :     const auto FixLong = [dfLongitudeMultiplier](double dfLon)
    2208         200 :     { return dfLon * dfLongitudeMultiplier; };
    2209             : 
    2210             :     // Note: starting with CART 1900, Spatial_Domain is actually optional
    2211          92 :     CPLXMLNode *psSD = CPLCreateXMLNode(psCart, CXT_Element,
    2212         184 :                                         (osPrefix + "Spatial_Domain").c_str());
    2213          92 :     CPLXMLNode *psBC = CPLCreateXMLNode(
    2214         184 :         psSD, CXT_Element, (osPrefix + "Bounding_Coordinates").c_str());
    2215             : 
    2216             :     const char *pszBoundingDegrees =
    2217          92 :         CSLFetchNameValue(m_papszCreationOptions, "BOUNDING_DEGREES");
    2218          92 :     double dfWest = FixLong(
    2219          92 :         std::min(std::min(adfX[0], adfX[1]), std::min(adfX[2], adfX[3])));
    2220          92 :     double dfEast = FixLong(
    2221          92 :         std::max(std::max(adfX[0], adfX[1]), std::max(adfX[2], adfX[3])));
    2222             :     double dfNorth =
    2223          92 :         std::max(std::max(adfY[0], adfY[1]), std::max(adfY[2], adfY[3]));
    2224             :     double dfSouth =
    2225          92 :         std::min(std::min(adfY[0], adfY[1]), std::min(adfY[2], adfY[3]));
    2226          92 :     if (pszBoundingDegrees)
    2227             :     {
    2228           1 :         char **papszTokens = CSLTokenizeString2(pszBoundingDegrees, ",", 0);
    2229           1 :         if (CSLCount(papszTokens) == 4)
    2230             :         {
    2231           1 :             dfWest = CPLAtof(papszTokens[0]);
    2232           1 :             dfSouth = CPLAtof(papszTokens[1]);
    2233           1 :             dfEast = CPLAtof(papszTokens[2]);
    2234           1 :             dfNorth = CPLAtof(papszTokens[3]);
    2235             :         }
    2236           1 :         CSLDestroy(papszTokens);
    2237             :     }
    2238             : 
    2239         184 :     CPLAddXMLAttributeAndValue(
    2240             :         CPLCreateXMLElementAndValue(
    2241         184 :             psBC, (osPrefix + "west_bounding_coordinate").c_str(),
    2242             :             CPLSPrintf("%.17g", dfWest)),
    2243             :         "unit", "deg");
    2244         184 :     CPLAddXMLAttributeAndValue(
    2245             :         CPLCreateXMLElementAndValue(
    2246         184 :             psBC, (osPrefix + "east_bounding_coordinate").c_str(),
    2247             :             CPLSPrintf("%.17g", dfEast)),
    2248             :         "unit", "deg");
    2249         184 :     CPLAddXMLAttributeAndValue(
    2250             :         CPLCreateXMLElementAndValue(
    2251         184 :             psBC, (osPrefix + "north_bounding_coordinate").c_str(),
    2252             :             CPLSPrintf("%.17g", dfNorth)),
    2253             :         "unit", "deg");
    2254         184 :     CPLAddXMLAttributeAndValue(
    2255             :         CPLCreateXMLElementAndValue(
    2256         184 :             psBC, (osPrefix + "south_bounding_coordinate").c_str(),
    2257             :             CPLSPrintf("%.17g", dfSouth)),
    2258             :         "unit", "deg");
    2259             : 
    2260             :     CPLXMLNode *psSRI =
    2261          92 :         CPLCreateXMLNode(psCart, CXT_Element,
    2262         184 :                          (osPrefix + "Spatial_Reference_Information").c_str());
    2263          92 :     CPLXMLNode *psHCSD = CPLCreateXMLNode(
    2264             :         psSRI, CXT_Element,
    2265         184 :         (osPrefix + "Horizontal_Coordinate_System_Definition").c_str());
    2266             : 
    2267          92 :     double dfUnrotatedULX = m_gt[0];
    2268          92 :     double dfUnrotatedULY = m_gt[3];
    2269          92 :     double dfUnrotatedResX = m_gt[1];
    2270          92 :     double dfUnrotatedResY = m_gt[5];
    2271          92 :     double dfMapProjectionRotation = 0.0;
    2272          92 :     if (m_gt[1] == 0.0 && m_gt[2] > 0.0 && m_gt[4] > 0.0 && m_gt[5] == 0.0)
    2273             :     {
    2274           1 :         dfUnrotatedULX = m_gt[3];
    2275           1 :         dfUnrotatedULY = -m_gt[0];
    2276           1 :         dfUnrotatedResX = m_gt[4];
    2277           1 :         dfUnrotatedResY = -m_gt[2];
    2278           1 :         dfMapProjectionRotation = 90.0;
    2279             :     }
    2280             : 
    2281          92 :     if (GetRasterCount() || m_oSRS.IsProjected())
    2282             :     {
    2283          91 :         CPLXMLNode *psPlanar = CPLCreateXMLNode(psHCSD, CXT_Element,
    2284         182 :                                                 (osPrefix + "Planar").c_str());
    2285          91 :         CPLXMLNode *psMP = CPLCreateXMLNode(
    2286         182 :             psPlanar, CXT_Element, (osPrefix + "Map_Projection").c_str());
    2287          91 :         const char *pszProjection = m_oSRS.GetAttrValue("PROJECTION");
    2288         182 :         CPLString pszPDS4ProjectionName = "";
    2289             :         typedef std::pair<const char *, double> ProjParam;
    2290         182 :         std::vector<ProjParam> aoProjParams;
    2291             : 
    2292             :         const bool bUse_CART_1933_Or_Later =
    2293          91 :             IsCARTVersionGTE(pszCARTVersion, "1D00_1933");
    2294             : 
    2295             :         const bool bUse_CART_1950_Or_Later =
    2296          91 :             IsCARTVersionGTE(pszCARTVersion, "1G00_1950");
    2297             : 
    2298          91 :         if (pszProjection == nullptr)
    2299             :         {
    2300          59 :             pszPDS4ProjectionName = "Equirectangular";
    2301          59 :             if (bUse_CART_1933_Or_Later)
    2302             :             {
    2303          57 :                 aoProjParams.push_back(
    2304          57 :                     ProjParam("latitude_of_projection_origin", 0.0));
    2305          57 :                 aoProjParams.push_back(ProjParam("standard_parallel_1", 0.0));
    2306          57 :                 aoProjParams.push_back(
    2307         114 :                     ProjParam("longitude_of_central_meridian", 0.0));
    2308             :             }
    2309             :             else
    2310             :             {
    2311           2 :                 aoProjParams.push_back(ProjParam("standard_parallel_1", 0.0));
    2312           2 :                 aoProjParams.push_back(
    2313           2 :                     ProjParam("longitude_of_central_meridian", 0.0));
    2314           2 :                 aoProjParams.push_back(
    2315           4 :                     ProjParam("latitude_of_projection_origin", 0.0));
    2316             :             }
    2317             :         }
    2318             : 
    2319          32 :         else if (EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR))
    2320             :         {
    2321           2 :             pszPDS4ProjectionName = "Equirectangular";
    2322           2 :             if (bUse_CART_1933_Or_Later)
    2323             :             {
    2324           2 :                 aoProjParams.push_back(ProjParam(
    2325           0 :                     "latitude_of_projection_origin",
    2326           2 :                     m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2327           2 :                 aoProjParams.push_back(ProjParam(
    2328           0 :                     "standard_parallel_1",
    2329           2 :                     m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 1.0)));
    2330           2 :                 aoProjParams.push_back(
    2331           0 :                     ProjParam("longitude_of_central_meridian",
    2332           4 :                               FixLong(m_oSRS.GetNormProjParm(
    2333             :                                   SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2334             :             }
    2335             :             else
    2336             :             {
    2337           0 :                 aoProjParams.push_back(ProjParam(
    2338           0 :                     "standard_parallel_1",
    2339           0 :                     m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 1.0)));
    2340           0 :                 aoProjParams.push_back(
    2341           0 :                     ProjParam("longitude_of_central_meridian",
    2342           0 :                               FixLong(m_oSRS.GetNormProjParm(
    2343             :                                   SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2344           0 :                 aoProjParams.push_back(ProjParam(
    2345           0 :                     "latitude_of_projection_origin",
    2346           0 :                     m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2347             :             }
    2348             :         }
    2349             : 
    2350          30 :         else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
    2351             :         {
    2352           1 :             pszPDS4ProjectionName = "Lambert Conformal Conic";
    2353           1 :             if (bUse_CART_1933_Or_Later)
    2354             :             {
    2355           1 :                 aoProjParams.push_back(
    2356           0 :                     ProjParam("longitude_of_central_meridian",
    2357           1 :                               FixLong(m_oSRS.GetNormProjParm(
    2358             :                                   SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2359           1 :                 aoProjParams.push_back(ProjParam(
    2360           0 :                     "latitude_of_projection_origin",
    2361           1 :                     m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2362           1 :                 aoProjParams.push_back(ProjParam(
    2363           0 :                     "scale_factor_at_projection_origin",
    2364           2 :                     m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
    2365             :             }
    2366             :             else
    2367             :             {
    2368           0 :                 aoProjParams.push_back(ProjParam(
    2369           0 :                     "scale_factor_at_projection_origin",
    2370           0 :                     m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
    2371           0 :                 aoProjParams.push_back(
    2372           0 :                     ProjParam("longitude_of_central_meridian",
    2373           0 :                               FixLong(m_oSRS.GetNormProjParm(
    2374             :                                   SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2375           0 :                 aoProjParams.push_back(ProjParam(
    2376           0 :                     "latitude_of_projection_origin",
    2377           0 :                     m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2378             :             }
    2379             :         }
    2380             : 
    2381          29 :         else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
    2382             :         {
    2383           1 :             pszPDS4ProjectionName = "Lambert Conformal Conic";
    2384           1 :             aoProjParams.push_back(ProjParam(
    2385           0 :                 "standard_parallel_1",
    2386           1 :                 m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0)));
    2387           1 :             aoProjParams.push_back(ProjParam(
    2388           0 :                 "standard_parallel_2",
    2389           1 :                 m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0)));
    2390           1 :             aoProjParams.push_back(ProjParam(
    2391           0 :                 "longitude_of_central_meridian",
    2392           1 :                 FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2393           1 :             aoProjParams.push_back(ProjParam(
    2394           0 :                 "latitude_of_projection_origin",
    2395           2 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2396             :         }
    2397             : 
    2398          28 :         else if (EQUAL(pszProjection,
    2399             :                        SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
    2400             :         {
    2401           1 :             pszPDS4ProjectionName = "Oblique Mercator";
    2402             :             // Proj params defined later
    2403             :         }
    2404             : 
    2405          27 :         else if (EQUAL(pszProjection,
    2406             :                        SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN))
    2407             :         {
    2408           1 :             pszPDS4ProjectionName = "Oblique Mercator";
    2409             :             // Proj params defined later
    2410             :         }
    2411             : 
    2412          26 :         else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
    2413             :         {
    2414           1 :             pszPDS4ProjectionName = "Polar Stereographic";
    2415           1 :             if (bUse_CART_1950_Or_Later)
    2416             :             {
    2417           1 :                 aoProjParams.push_back(
    2418           0 :                     ProjParam("longitude_of_central_meridian",
    2419           1 :                               FixLong(m_oSRS.GetNormProjParm(
    2420             :                                   SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2421           1 :                 aoProjParams.push_back(ProjParam(
    2422           0 :                     "latitude_of_projection_origin",
    2423           1 :                     m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2424           1 :                 aoProjParams.push_back(ProjParam(
    2425           0 :                     "scale_factor_at_projection_origin",
    2426           2 :                     m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
    2427             :             }
    2428             :             else
    2429             :             {
    2430           0 :                 aoProjParams.push_back(
    2431           0 :                     ProjParam(bUse_CART_1933_Or_Later
    2432           0 :                                   ? "longitude_of_central_meridian"
    2433             :                                   : "straight_vertical_longitude_from_pole",
    2434           0 :                               FixLong(m_oSRS.GetNormProjParm(
    2435             :                                   SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2436           0 :                 aoProjParams.push_back(ProjParam(
    2437           0 :                     "scale_factor_at_projection_origin",
    2438           0 :                     m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
    2439           0 :                 aoProjParams.push_back(ProjParam(
    2440           0 :                     "latitude_of_projection_origin",
    2441           0 :                     m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2442             :             }
    2443             :         }
    2444             : 
    2445          25 :         else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
    2446             :         {
    2447           1 :             pszPDS4ProjectionName = "Polyconic";
    2448           1 :             aoProjParams.push_back(ProjParam(
    2449           0 :                 "longitude_of_central_meridian",
    2450           1 :                 m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)));
    2451           1 :             aoProjParams.push_back(ProjParam(
    2452           0 :                 "latitude_of_projection_origin",
    2453           2 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2454             :         }
    2455          24 :         else if (EQUAL(pszProjection, SRS_PT_SINUSOIDAL))
    2456             :         {
    2457           2 :             pszPDS4ProjectionName = "Sinusoidal";
    2458           2 :             aoProjParams.push_back(ProjParam(
    2459           0 :                 "longitude_of_central_meridian",
    2460           2 :                 FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2461           2 :             aoProjParams.push_back(ProjParam(
    2462           0 :                 "latitude_of_projection_origin",
    2463           4 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2464             :         }
    2465             : 
    2466          22 :         else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
    2467             :         {
    2468          16 :             pszPDS4ProjectionName = "Transverse Mercator";
    2469          16 :             aoProjParams.push_back(
    2470           0 :                 ProjParam("scale_factor_at_central_meridian",
    2471          16 :                           m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
    2472          16 :             aoProjParams.push_back(ProjParam(
    2473           0 :                 "longitude_of_central_meridian",
    2474          16 :                 m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)));
    2475          16 :             aoProjParams.push_back(ProjParam(
    2476           0 :                 "latitude_of_projection_origin",
    2477          32 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2478             :         }
    2479           6 :         else if (EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC))
    2480             :         {
    2481           1 :             pszPDS4ProjectionName = "Orthographic";
    2482           1 :             aoProjParams.push_back(ProjParam(
    2483           0 :                 "longitude_of_central_meridian",
    2484           1 :                 FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2485           1 :             aoProjParams.push_back(ProjParam(
    2486           0 :                 "latitude_of_projection_origin",
    2487           2 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2488             :         }
    2489             : 
    2490           5 :         else if (EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
    2491             :         {
    2492           2 :             pszPDS4ProjectionName = "Mercator";
    2493           2 :             aoProjParams.push_back(ProjParam(
    2494           0 :                 "longitude_of_central_meridian",
    2495           2 :                 FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2496           2 :             aoProjParams.push_back(ProjParam(
    2497           0 :                 "latitude_of_projection_origin",
    2498           2 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2499           2 :             aoProjParams.push_back(
    2500           0 :                 ProjParam("scale_factor_at_projection_origin",
    2501           4 :                           m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
    2502             :         }
    2503             : 
    2504           3 :         else if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
    2505             :         {
    2506           1 :             pszPDS4ProjectionName = "Mercator";
    2507           1 :             aoProjParams.push_back(ProjParam(
    2508           0 :                 "standard_parallel_1",
    2509           1 :                 m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0)));
    2510           1 :             aoProjParams.push_back(ProjParam(
    2511           0 :                 "longitude_of_central_meridian",
    2512           1 :                 FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2513           1 :             aoProjParams.push_back(ProjParam(
    2514           0 :                 "latitude_of_projection_origin",
    2515           2 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2516             :         }
    2517             : 
    2518           2 :         else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
    2519             :         {
    2520           1 :             pszPDS4ProjectionName = "Lambert Azimuthal Equal Area";
    2521           1 :             aoProjParams.push_back(ProjParam(
    2522           0 :                 "longitude_of_central_meridian",
    2523           1 :                 FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
    2524           1 :             aoProjParams.push_back(ProjParam(
    2525           0 :                 "latitude_of_projection_origin",
    2526           2 :                 m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
    2527             :         }
    2528             : 
    2529           1 :         else if (EQUAL(pszProjection, "custom_proj4"))
    2530             :         {
    2531             :             const char *pszProj4 =
    2532           1 :                 m_oSRS.GetExtension("PROJCS", "PROJ4", nullptr);
    2533           1 :             if (pszProj4 && strstr(pszProj4, "+proj=ob_tran") &&
    2534           1 :                 strstr(pszProj4, "+o_proj=eqc"))
    2535             :             {
    2536           1 :                 pszPDS4ProjectionName = "Oblique Cylindrical";
    2537             :                 const auto FetchParam =
    2538           3 :                     [](const char *pszProj4Str, const char *pszKey)
    2539             :                 {
    2540           6 :                     CPLString needle;
    2541           3 :                     needle.Printf("+%s=", pszKey);
    2542           3 :                     const char *pszVal = strstr(pszProj4Str, needle.c_str());
    2543           3 :                     if (pszVal)
    2544           3 :                         return CPLAtof(pszVal + needle.size());
    2545           0 :                     return 0.0;
    2546             :                 };
    2547             : 
    2548           1 :                 double dfLonP = FetchParam(pszProj4, "o_lon_p");
    2549           1 :                 double dfLatP = FetchParam(pszProj4, "o_lat_p");
    2550           1 :                 double dfLon0 = FetchParam(pszProj4, "lon_0");
    2551           1 :                 double dfPoleRotation = -dfLonP;
    2552           1 :                 double dfPoleLatitude = 180 - dfLatP;
    2553           1 :                 double dfPoleLongitude = dfLon0;
    2554             : 
    2555           1 :                 aoProjParams.push_back(ProjParam("map_projection_rotation",
    2556             :                                                  dfMapProjectionRotation));
    2557           1 :                 aoProjParams.push_back(
    2558           1 :                     ProjParam("oblique_proj_pole_latitude", dfPoleLatitude));
    2559           1 :                 aoProjParams.push_back(ProjParam("oblique_proj_pole_longitude",
    2560           1 :                                                  FixLong(dfPoleLongitude)));
    2561           1 :                 aoProjParams.push_back(
    2562           2 :                     ProjParam("oblique_proj_pole_rotation", dfPoleRotation));
    2563             :             }
    2564             :             else
    2565             :             {
    2566           0 :                 CPLError(CE_Warning, CPLE_NotSupported,
    2567             :                          "Projection %s not supported", pszProjection);
    2568             :             }
    2569             :         }
    2570             :         else
    2571             :         {
    2572           0 :             CPLError(CE_Warning, CPLE_NotSupported,
    2573             :                      "Projection %s not supported", pszProjection);
    2574             :         }
    2575         182 :         CPLCreateXMLElementAndValue(psMP,
    2576         182 :                                     (osPrefix + "map_projection_name").c_str(),
    2577             :                                     pszPDS4ProjectionName);
    2578          91 :         CPLXMLNode *psProj = CPLCreateXMLNode(
    2579             :             psMP, CXT_Element,
    2580         182 :             CPLString(osPrefix + pszPDS4ProjectionName).replaceAll(' ', '_'));
    2581         355 :         for (size_t i = 0; i < aoProjParams.size(); i++)
    2582             :         {
    2583         528 :             CPLXMLNode *psParam = CPLCreateXMLElementAndValue(
    2584         528 :                 psProj, (osPrefix + aoProjParams[i].first).c_str(),
    2585         264 :                 CPLSPrintf("%.17g", aoProjParams[i].second));
    2586         264 :             if (!STARTS_WITH(aoProjParams[i].first, "scale_factor"))
    2587             :             {
    2588         244 :                 CPLAddXMLAttributeAndValue(psParam, "unit", "deg");
    2589             :             }
    2590             :         }
    2591             : 
    2592          91 :         if (pszProjection &&
    2593          32 :             EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
    2594             :         {
    2595             :             CPLXMLNode *psOLA =
    2596           1 :                 CPLCreateXMLNode(nullptr, CXT_Element,
    2597           2 :                                  (osPrefix + "Oblique_Line_Azimuth").c_str());
    2598           2 :             CPLAddXMLAttributeAndValue(
    2599             :                 CPLCreateXMLElementAndValue(
    2600           2 :                     psOLA, (osPrefix + "azimuthal_angle").c_str(),
    2601             :                     CPLSPrintf("%.17g",
    2602             :                                m_oSRS.GetNormProjParm(SRS_PP_AZIMUTH, 0.0))),
    2603             :                 "unit", "deg");
    2604             :             ;
    2605             :             // Not completely sure of this
    2606           2 :             CPLAddXMLAttributeAndValue(
    2607             :                 CPLCreateXMLElementAndValue(
    2608             :                     psOLA,
    2609           2 :                     (osPrefix + "azimuth_measure_point_longitude").c_str(),
    2610             :                     CPLSPrintf("%.17g", FixLong(m_oSRS.GetNormProjParm(
    2611             :                                             SRS_PP_CENTRAL_MERIDIAN, 0.0)))),
    2612             :                 "unit", "deg");
    2613             : 
    2614           1 :             if (bUse_CART_1933_Or_Later)
    2615             :             {
    2616           1 :                 CPLAddXMLChild(psProj, psOLA);
    2617             : 
    2618           1 :                 CPLAddXMLAttributeAndValue(
    2619             :                     CPLCreateXMLElementAndValue(
    2620             :                         psProj,
    2621           2 :                         (osPrefix + "longitude_of_central_meridian").c_str(),
    2622             :                         "0"),
    2623             :                     "unit", "deg");
    2624             : 
    2625             :                 const double dfScaleFactor =
    2626           1 :                     m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0);
    2627           1 :                 if (dfScaleFactor != 1.0)
    2628             :                 {
    2629           0 :                     CPLError(CE_Warning, CPLE_NotSupported,
    2630             :                              "Scale factor on initial support = %.17g cannot "
    2631             :                              "be encoded in PDS4",
    2632             :                              dfScaleFactor);
    2633             :                 }
    2634             :             }
    2635             :             else
    2636             :             {
    2637           0 :                 CPLCreateXMLElementAndValue(
    2638             :                     psProj,
    2639           0 :                     (osPrefix + "scale_factor_at_projection_origin").c_str(),
    2640             :                     CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
    2641             :                                             SRS_PP_SCALE_FACTOR, 0.0)));
    2642             : 
    2643           0 :                 CPLAddXMLChild(psProj, psOLA);
    2644             :             }
    2645             : 
    2646           2 :             CPLAddXMLAttributeAndValue(
    2647             :                 CPLCreateXMLElementAndValue(
    2648             :                     psProj,
    2649           2 :                     (osPrefix + "latitude_of_projection_origin").c_str(),
    2650             :                     CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
    2651             :                                             SRS_PP_LATITUDE_OF_ORIGIN, 0.0))),
    2652           1 :                 "unit", "deg");
    2653             :         }
    2654          90 :         else if (pszProjection &&
    2655          31 :                  EQUAL(pszProjection,
    2656             :                        SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN))
    2657             :         {
    2658           1 :             if (bUse_CART_1933_Or_Later)
    2659             :             {
    2660             :                 const double dfScaleFactor =
    2661           1 :                     m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0);
    2662           1 :                 if (dfScaleFactor != 1.0)
    2663             :                 {
    2664           0 :                     CPLError(CE_Warning, CPLE_NotSupported,
    2665             :                              "Scale factor on initial support = %.17g cannot "
    2666             :                              "be encoded in PDS4",
    2667             :                              dfScaleFactor);
    2668             :                 }
    2669             :             }
    2670             :             else
    2671             :             {
    2672           0 :                 CPLCreateXMLElementAndValue(
    2673             :                     psProj,
    2674           0 :                     (osPrefix + "scale_factor_at_projection_origin").c_str(),
    2675             :                     CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
    2676             :                                             SRS_PP_SCALE_FACTOR, 0.0)));
    2677             :             }
    2678             : 
    2679           1 :             CPLXMLNode *psOLP = CPLCreateXMLNode(
    2680           2 :                 psProj, CXT_Element, (osPrefix + "Oblique_Line_Point").c_str());
    2681           1 :             CPLXMLNode *psOLPG1 = CPLCreateXMLNode(
    2682             :                 psOLP, CXT_Element,
    2683           2 :                 (osPrefix + "Oblique_Line_Point_Group").c_str());
    2684           2 :             CPLAddXMLAttributeAndValue(
    2685             :                 CPLCreateXMLElementAndValue(
    2686           2 :                     psOLPG1, (osPrefix + "oblique_line_latitude").c_str(),
    2687             :                     CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
    2688             :                                             SRS_PP_LATITUDE_OF_POINT_1, 0.0))),
    2689             :                 "unit", "deg");
    2690           2 :             CPLAddXMLAttributeAndValue(
    2691             :                 CPLCreateXMLElementAndValue(
    2692           2 :                     psOLPG1, (osPrefix + "oblique_line_longitude").c_str(),
    2693             :                     CPLSPrintf("%.17g",
    2694             :                                FixLong(m_oSRS.GetNormProjParm(
    2695             :                                    SRS_PP_LONGITUDE_OF_POINT_1, 0.0)))),
    2696             :                 "unit", "deg");
    2697           1 :             CPLXMLNode *psOLPG2 = CPLCreateXMLNode(
    2698             :                 psOLP, CXT_Element,
    2699           2 :                 (osPrefix + "Oblique_Line_Point_Group").c_str());
    2700           2 :             CPLAddXMLAttributeAndValue(
    2701             :                 CPLCreateXMLElementAndValue(
    2702           2 :                     psOLPG2, (osPrefix + "oblique_line_latitude").c_str(),
    2703             :                     CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
    2704             :                                             SRS_PP_LATITUDE_OF_POINT_2, 0.0))),
    2705             :                 "unit", "deg");
    2706           2 :             CPLAddXMLAttributeAndValue(
    2707             :                 CPLCreateXMLElementAndValue(
    2708           2 :                     psOLPG2, (osPrefix + "oblique_line_longitude").c_str(),
    2709             :                     CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
    2710             :                                             SRS_PP_LONGITUDE_OF_POINT_2, 0.0))),
    2711             :                 "unit", "deg");
    2712             : 
    2713           1 :             if (bUse_CART_1933_Or_Later)
    2714             :             {
    2715           1 :                 CPLAddXMLAttributeAndValue(
    2716             :                     CPLCreateXMLElementAndValue(
    2717             :                         psProj,
    2718           2 :                         (osPrefix + "longitude_of_central_meridian").c_str(),
    2719             :                         "0"),
    2720             :                     "unit", "deg");
    2721             :             }
    2722             : 
    2723           2 :             CPLAddXMLAttributeAndValue(
    2724             :                 CPLCreateXMLElementAndValue(
    2725             :                     psProj,
    2726           2 :                     (osPrefix + "latitude_of_projection_origin").c_str(),
    2727             :                     CPLSPrintf("%.17g", FixLong(m_oSRS.GetNormProjParm(
    2728             :                                             SRS_PP_LATITUDE_OF_ORIGIN, 0.0)))),
    2729             :                 "unit", "deg");
    2730             :         }
    2731             : 
    2732          91 :         CPLXMLNode *psCR = nullptr;
    2733          91 :         if (m_bGotTransform || !IsCARTVersionGTE(pszCARTVersion, "1B00"))
    2734             :         {
    2735          90 :             CPLXMLNode *psPCI = CPLCreateXMLNode(
    2736             :                 psPlanar, CXT_Element,
    2737         180 :                 (osPrefix + "Planar_Coordinate_Information").c_str());
    2738          90 :             CPLCreateXMLElementAndValue(
    2739         180 :                 psPCI, (osPrefix + "planar_coordinate_encoding_method").c_str(),
    2740             :                 "Coordinate Pair");
    2741          90 :             psCR = CPLCreateXMLNode(
    2742             :                 psPCI, CXT_Element,
    2743         180 :                 (osPrefix + "Coordinate_Representation").c_str());
    2744             :         }
    2745          91 :         const double dfLinearUnits = m_oSRS.GetLinearUnits();
    2746          91 :         const double dfDegToMeter = m_oSRS.GetSemiMajor() * M_PI / 180.0;
    2747             : 
    2748          91 :         if (psCR == nullptr)
    2749             :         {
    2750             :             // do nothing
    2751             :         }
    2752          90 :         else if (!m_bGotTransform)
    2753             :         {
    2754           0 :             CPLAddXMLAttributeAndValue(
    2755             :                 CPLCreateXMLElementAndValue(
    2756           0 :                     psCR, (osPrefix + "pixel_resolution_x").c_str(), "0"),
    2757             :                 "unit", "m/pixel");
    2758           0 :             CPLAddXMLAttributeAndValue(
    2759             :                 CPLCreateXMLElementAndValue(
    2760           0 :                     psCR, (osPrefix + "pixel_resolution_y").c_str(), "0"),
    2761             :                 "unit", "m/pixel");
    2762           0 :             CPLAddXMLAttributeAndValue(
    2763             :                 CPLCreateXMLElementAndValue(
    2764           0 :                     psCR, (osPrefix + "pixel_scale_x").c_str(), "0"),
    2765             :                 "unit", "pixel/deg");
    2766           0 :             CPLAddXMLAttributeAndValue(
    2767             :                 CPLCreateXMLElementAndValue(
    2768           0 :                     psCR, (osPrefix + "pixel_scale_y").c_str(), "0"),
    2769             :                 "unit", "pixel/deg");
    2770             :         }
    2771          90 :         else if (m_oSRS.IsGeographic())
    2772             :         {
    2773         118 :             CPLAddXMLAttributeAndValue(
    2774             :                 CPLCreateXMLElementAndValue(
    2775         118 :                     psCR, (osPrefix + "pixel_resolution_x").c_str(),
    2776             :                     CPLSPrintf("%.17g", dfUnrotatedResX * dfDegToMeter)),
    2777             :                 "unit", "m/pixel");
    2778          59 :             CPLAddXMLAttributeAndValue(
    2779             :                 CPLCreateXMLElementAndValue(
    2780         118 :                     psCR, (osPrefix + "pixel_resolution_y").c_str(),
    2781          59 :                     CPLSPrintf("%.17g", -dfUnrotatedResY * dfDegToMeter)),
    2782             :                 "unit", "m/pixel");
    2783         118 :             CPLAddXMLAttributeAndValue(
    2784             :                 CPLCreateXMLElementAndValue(
    2785         118 :                     psCR, (osPrefix + "pixel_scale_x").c_str(),
    2786             :                     CPLSPrintf("%.17g", 1.0 / (dfUnrotatedResX))),
    2787             :                 "unit", "pixel/deg");
    2788         118 :             CPLAddXMLAttributeAndValue(
    2789             :                 CPLCreateXMLElementAndValue(
    2790         118 :                     psCR, (osPrefix + "pixel_scale_y").c_str(),
    2791             :                     CPLSPrintf("%.17g", 1.0 / (-dfUnrotatedResY))),
    2792             :                 "unit", "pixel/deg");
    2793             :         }
    2794          31 :         else if (m_oSRS.IsProjected())
    2795             :         {
    2796          62 :             CPLAddXMLAttributeAndValue(
    2797             :                 CPLCreateXMLElementAndValue(
    2798          62 :                     psCR, (osPrefix + "pixel_resolution_x").c_str(),
    2799             :                     CPLSPrintf("%.17g", dfUnrotatedResX * dfLinearUnits)),
    2800             :                 "unit", "m/pixel");
    2801          31 :             CPLAddXMLAttributeAndValue(
    2802             :                 CPLCreateXMLElementAndValue(
    2803          62 :                     psCR, (osPrefix + "pixel_resolution_y").c_str(),
    2804          31 :                     CPLSPrintf("%.17g", -dfUnrotatedResY * dfLinearUnits)),
    2805             :                 "unit", "m/pixel");
    2806          31 :             CPLAddXMLAttributeAndValue(
    2807             :                 CPLCreateXMLElementAndValue(
    2808          62 :                     psCR, (osPrefix + "pixel_scale_x").c_str(),
    2809             :                     CPLSPrintf("%.17g", dfDegToMeter /
    2810          31 :                                             (dfUnrotatedResX * dfLinearUnits))),
    2811             :                 "unit", "pixel/deg");
    2812          31 :             CPLAddXMLAttributeAndValue(
    2813             :                 CPLCreateXMLElementAndValue(
    2814          62 :                     psCR, (osPrefix + "pixel_scale_y").c_str(),
    2815          31 :                     CPLSPrintf("%.17g", dfDegToMeter / (-dfUnrotatedResY *
    2816             :                                                         dfLinearUnits))),
    2817             :                 "unit", "pixel/deg");
    2818             :         }
    2819             : 
    2820          91 :         if (m_bGotTransform)
    2821             :         {
    2822             :             CPLXMLNode *psGT =
    2823          90 :                 CPLCreateXMLNode(psPlanar, CXT_Element,
    2824         180 :                                  (osPrefix + "Geo_Transformation").c_str());
    2825             :             const double dfFalseEasting =
    2826          90 :                 m_oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
    2827             :             const double dfFalseNorthing =
    2828          90 :                 m_oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
    2829          90 :             const double dfULX = -dfFalseEasting + dfUnrotatedULX;
    2830          90 :             const double dfULY = -dfFalseNorthing + dfUnrotatedULY;
    2831          90 :             if (m_oSRS.IsGeographic())
    2832             :             {
    2833         118 :                 CPLAddXMLAttributeAndValue(
    2834             :                     CPLCreateXMLElementAndValue(
    2835         118 :                         psGT, (osPrefix + "upperleft_corner_x").c_str(),
    2836             :                         CPLSPrintf("%.17g", dfULX * dfDegToMeter)),
    2837             :                     "unit", "m");
    2838         118 :                 CPLAddXMLAttributeAndValue(
    2839             :                     CPLCreateXMLElementAndValue(
    2840         118 :                         psGT, (osPrefix + "upperleft_corner_y").c_str(),
    2841             :                         CPLSPrintf("%.17g", dfULY * dfDegToMeter)),
    2842             :                     "unit", "m");
    2843             :             }
    2844          31 :             else if (m_oSRS.IsProjected())
    2845             :             {
    2846          62 :                 CPLAddXMLAttributeAndValue(
    2847             :                     CPLCreateXMLElementAndValue(
    2848          62 :                         psGT, (osPrefix + "upperleft_corner_x").c_str(),
    2849             :                         CPLSPrintf("%.17g", dfULX * dfLinearUnits)),
    2850             :                     "unit", "m");
    2851          62 :                 CPLAddXMLAttributeAndValue(
    2852             :                     CPLCreateXMLElementAndValue(
    2853          62 :                         psGT, (osPrefix + "upperleft_corner_y").c_str(),
    2854             :                         CPLSPrintf("%.17g", dfULY * dfLinearUnits)),
    2855             :                     "unit", "m");
    2856             :             }
    2857             :         }
    2858             :     }
    2859             :     else
    2860             :     {
    2861           1 :         CPLXMLNode *psGeographic = CPLCreateXMLNode(
    2862           2 :             psHCSD, CXT_Element, (osPrefix + "Geographic").c_str());
    2863           1 :         if (!IsCARTVersionGTE(pszCARTVersion, "1B00"))
    2864             :         {
    2865           0 :             CPLAddXMLAttributeAndValue(
    2866             :                 CPLCreateXMLElementAndValue(
    2867           0 :                     psGeographic, (osPrefix + "latitude_resolution").c_str(),
    2868             :                     "0"),
    2869             :                 "unit", "deg");
    2870           0 :             CPLAddXMLAttributeAndValue(
    2871             :                 CPLCreateXMLElementAndValue(
    2872           0 :                     psGeographic, (osPrefix + "longitude_resolution").c_str(),
    2873             :                     "0"),
    2874             :                 "unit", "deg");
    2875             :         }
    2876             :     }
    2877             : 
    2878          92 :     CPLXMLNode *psGM = CPLCreateXMLNode(psHCSD, CXT_Element,
    2879         184 :                                         (osPrefix + "Geodetic_Model").c_str());
    2880         184 :     const char *pszLatitudeType = CSLFetchNameValueDef(
    2881          92 :         m_papszCreationOptions, "LATITUDE_TYPE", "Planetocentric");
    2882             :     // Fix case
    2883          92 :     if (EQUAL(pszLatitudeType, "Planetocentric"))
    2884          91 :         pszLatitudeType = "Planetocentric";
    2885           1 :     else if (EQUAL(pszLatitudeType, "Planetographic"))
    2886           1 :         pszLatitudeType = "Planetographic";
    2887          92 :     CPLCreateXMLElementAndValue(psGM, (osPrefix + "latitude_type").c_str(),
    2888             :                                 pszLatitudeType);
    2889             : 
    2890          92 :     const char *pszDatum = m_oSRS.GetAttrValue("DATUM");
    2891          92 :     if (pszDatum && STARTS_WITH(pszDatum, "D_"))
    2892             :     {
    2893           4 :         CPLCreateXMLElementAndValue(psGM, (osPrefix + "spheroid_name").c_str(),
    2894             :                                     pszDatum + 2);
    2895             :     }
    2896          88 :     else if (pszDatum)
    2897             :     {
    2898          88 :         CPLCreateXMLElementAndValue(psGM, (osPrefix + "spheroid_name").c_str(),
    2899             :                                     pszDatum);
    2900             :     }
    2901             : 
    2902          92 :     double dfSemiMajor = m_oSRS.GetSemiMajor();
    2903          92 :     double dfSemiMinor = m_oSRS.GetSemiMinor();
    2904          92 :     const char *pszRadii = CSLFetchNameValue(m_papszCreationOptions, "RADII");
    2905          92 :     if (pszRadii)
    2906             :     {
    2907           1 :         char **papszTokens = CSLTokenizeString2(pszRadii, " ,", 0);
    2908           1 :         if (CSLCount(papszTokens) == 2)
    2909             :         {
    2910           1 :             dfSemiMajor = CPLAtof(papszTokens[0]);
    2911           1 :             dfSemiMinor = CPLAtof(papszTokens[1]);
    2912             :         }
    2913           1 :         CSLDestroy(papszTokens);
    2914             :     }
    2915             : 
    2916             :     const bool bUseLDD1930RadiusNames =
    2917          92 :         IsCARTVersionGTE(pszCARTVersion, "1B10_1930");
    2918             : 
    2919         184 :     CPLAddXMLAttributeAndValue(
    2920             :         CPLCreateXMLElementAndValue(
    2921             :             psGM,
    2922         184 :             (osPrefix +
    2923             :              (bUseLDD1930RadiusNames ? "a_axis_radius" : "semi_major_radius"))
    2924             :                 .c_str(),
    2925             :             CPLSPrintf("%.17g", dfSemiMajor)),
    2926             :         "unit", "m");
    2927             :     // No, this is not a bug. The PDS4  b_axis_radius/semi_minor_radius is the
    2928             :     // minor radius on the equatorial plane. Which in WKT doesn't really exist,
    2929             :     // so reuse the WKT semi major
    2930         184 :     CPLAddXMLAttributeAndValue(
    2931             :         CPLCreateXMLElementAndValue(
    2932             :             psGM,
    2933         184 :             (osPrefix +
    2934             :              (bUseLDD1930RadiusNames ? "b_axis_radius" : "semi_minor_radius"))
    2935             :                 .c_str(),
    2936             :             CPLSPrintf("%.17g", dfSemiMajor)),
    2937             :         "unit", "m");
    2938         184 :     CPLAddXMLAttributeAndValue(
    2939             :         CPLCreateXMLElementAndValue(
    2940             :             psGM,
    2941         184 :             (osPrefix +
    2942             :              (bUseLDD1930RadiusNames ? "c_axis_radius" : "polar_radius"))
    2943             :                 .c_str(),
    2944             :             CPLSPrintf("%.17g", dfSemiMinor)),
    2945             :         "unit", "m");
    2946             : 
    2947             :     // Fix case
    2948          92 :     if (EQUAL(pszLongitudeDirection, "Positive East"))
    2949          91 :         pszLongitudeDirection = "Positive East";
    2950           1 :     else if (EQUAL(pszLongitudeDirection, "Positive West"))
    2951           1 :         pszLongitudeDirection = "Positive West";
    2952          92 :     CPLCreateXMLElementAndValue(psGM,
    2953         184 :                                 (osPrefix + "longitude_direction").c_str(),
    2954             :                                 pszLongitudeDirection);
    2955          92 : }
    2956             : 
    2957             : /************************************************************************/
    2958             : /*                         SubstituteVariables()                        */
    2959             : /************************************************************************/
    2960             : 
    2961       15242 : void PDS4Dataset::SubstituteVariables(CPLXMLNode *psNode, char **papszDict)
    2962             : {
    2963       15242 :     if (psNode->eType == CXT_Text && psNode->pszValue &&
    2964        6172 :         strstr(psNode->pszValue, "${"))
    2965             :     {
    2966        1254 :         CPLString osVal(psNode->pszValue);
    2967             : 
    2968        1412 :         if (strstr(psNode->pszValue, "${TITLE}") != nullptr &&
    2969         158 :             CSLFetchNameValue(papszDict, "VAR_TITLE") == nullptr)
    2970             :         {
    2971         144 :             const CPLString osTitle(CPLGetFilename(GetDescription()));
    2972         144 :             CPLError(CE_Warning, CPLE_AppDefined,
    2973             :                      "VAR_TITLE not defined. Using %s by default",
    2974             :                      osTitle.c_str());
    2975         144 :             osVal.replaceAll("${TITLE}", osTitle);
    2976             :         }
    2977             : 
    2978        4311 :         for (char **papszIter = papszDict; papszIter && *papszIter; papszIter++)
    2979             :         {
    2980        3057 :             if (STARTS_WITH_CI(*papszIter, "VAR_"))
    2981             :             {
    2982        2650 :                 char *pszKey = nullptr;
    2983        2650 :                 const char *pszValue = CPLParseNameValue(*papszIter, &pszKey);
    2984        2650 :                 if (pszKey && pszValue)
    2985             :                 {
    2986        2650 :                     const char *pszVarName = pszKey + strlen("VAR_");
    2987             :                     osVal.replaceAll(
    2988        2650 :                         (CPLString("${") + pszVarName + "}").c_str(), pszValue);
    2989             :                     osVal.replaceAll(
    2990        5300 :                         CPLString(CPLString("${") + pszVarName + "}")
    2991        2650 :                             .tolower()
    2992             :                             .c_str(),
    2993        7950 :                         CPLString(pszValue).tolower());
    2994        2650 :                     CPLFree(pszKey);
    2995             :                 }
    2996             :             }
    2997             :         }
    2998        1254 :         if (osVal.find("${") != std::string::npos)
    2999             :         {
    3000         684 :             CPLError(CE_Warning, CPLE_AppDefined, "%s could not be substituted",
    3001             :                      osVal.c_str());
    3002             :         }
    3003        1254 :         CPLFree(psNode->pszValue);
    3004        1254 :         psNode->pszValue = CPLStrdup(osVal);
    3005             :     }
    3006             : 
    3007       30319 :     for (CPLXMLNode *psIter = psNode->psChild; psIter; psIter = psIter->psNext)
    3008             :     {
    3009       15077 :         SubstituteVariables(psIter, papszDict);
    3010             :     }
    3011       15242 : }
    3012             : 
    3013             : /************************************************************************/
    3014             : /*                         InitImageFile()                             */
    3015             : /************************************************************************/
    3016             : 
    3017          77 : bool PDS4Dataset::InitImageFile()
    3018             : {
    3019          77 :     m_bMustInitImageFile = false;
    3020             : 
    3021          77 :     if (m_poExternalDS)
    3022             :     {
    3023             :         int nBlockXSize, nBlockYSize;
    3024          11 :         GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
    3025          11 :         const GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
    3026          11 :         const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
    3027          11 :         const int nBlockSizeBytes = nBlockXSize * nBlockYSize * nDTSize;
    3028          11 :         const int l_nBlocksPerColumn = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
    3029             : 
    3030          11 :         int bHasNoData = FALSE;
    3031          11 :         double dfNoData = GetRasterBand(1)->GetNoDataValue(&bHasNoData);
    3032          11 :         if (!bHasNoData)
    3033           7 :             dfNoData = 0;
    3034             : 
    3035          11 :         if (nBands == 1 || EQUAL(m_osInterleave, "BSQ"))
    3036             :         {
    3037             :             // We need to make sure that blocks are written in the right order
    3038          24 :             for (int i = 0; i < nBands; i++)
    3039             :             {
    3040          14 :                 if (m_poExternalDS->GetRasterBand(i + 1)->Fill(dfNoData) !=
    3041             :                     CE_None)
    3042             :                 {
    3043           0 :                     return false;
    3044             :                 }
    3045             :             }
    3046          10 :             m_poExternalDS->FlushCache(false);
    3047             : 
    3048             :             // Check that blocks are effectively written in expected order.
    3049          10 :             GIntBig nLastOffset = 0;
    3050          24 :             for (int i = 0; i < nBands; i++)
    3051             :             {
    3052         319 :                 for (int y = 0; y < l_nBlocksPerColumn; y++)
    3053             :                 {
    3054             :                     const char *pszBlockOffset =
    3055         610 :                         m_poExternalDS->GetRasterBand(i + 1)->GetMetadataItem(
    3056         305 :                             CPLSPrintf("BLOCK_OFFSET_%d_%d", 0, y), "TIFF");
    3057         305 :                     if (pszBlockOffset)
    3058             :                     {
    3059         305 :                         GIntBig nOffset = CPLAtoGIntBig(pszBlockOffset);
    3060         305 :                         if (i != 0 || y != 0)
    3061             :                         {
    3062         295 :                             if (nOffset != nLastOffset + nBlockSizeBytes)
    3063             :                             {
    3064           0 :                                 CPLError(CE_Warning, CPLE_AppDefined,
    3065             :                                          "Block %d,%d band %d not at expected "
    3066             :                                          "offset",
    3067             :                                          0, y, i + 1);
    3068           0 :                                 return false;
    3069             :                             }
    3070             :                         }
    3071         305 :                         nLastOffset = nOffset;
    3072             :                     }
    3073             :                     else
    3074             :                     {
    3075           0 :                         CPLError(CE_Warning, CPLE_AppDefined,
    3076             :                                  "Block %d,%d band %d not at expected "
    3077             :                                  "offset",
    3078             :                                  0, y, i + 1);
    3079           0 :                         return false;
    3080             :                     }
    3081             :                 }
    3082             :             }
    3083             :         }
    3084             :         else
    3085             :         {
    3086           1 :             void *pBlockData = VSI_MALLOC_VERBOSE(nBlockSizeBytes);
    3087           1 :             if (pBlockData == nullptr)
    3088           0 :                 return false;
    3089           1 :             GDALCopyWords(&dfNoData, GDT_Float64, 0, pBlockData, eDT, nDTSize,
    3090             :                           nBlockXSize * nBlockYSize);
    3091           2 :             for (int y = 0; y < l_nBlocksPerColumn; y++)
    3092             :             {
    3093           4 :                 for (int i = 0; i < nBands; i++)
    3094             :                 {
    3095           3 :                     if (m_poExternalDS->GetRasterBand(i + 1)->WriteBlock(
    3096           3 :                             0, y, pBlockData) != CE_None)
    3097             :                     {
    3098           0 :                         VSIFree(pBlockData);
    3099           0 :                         return false;
    3100             :                     }
    3101             :                 }
    3102             :             }
    3103           1 :             VSIFree(pBlockData);
    3104           1 :             m_poExternalDS->FlushCache(false);
    3105             : 
    3106             :             // Check that blocks are effectively written in expected order.
    3107           1 :             GIntBig nLastOffset = 0;
    3108           2 :             for (int y = 0; y < l_nBlocksPerColumn; y++)
    3109             :             {
    3110             :                 const char *pszBlockOffset =
    3111           2 :                     m_poExternalDS->GetRasterBand(1)->GetMetadataItem(
    3112           1 :                         CPLSPrintf("BLOCK_OFFSET_%d_%d", 0, y), "TIFF");
    3113           1 :                 if (pszBlockOffset)
    3114             :                 {
    3115           1 :                     GIntBig nOffset = CPLAtoGIntBig(pszBlockOffset);
    3116           1 :                     if (y != 0)
    3117             :                     {
    3118           0 :                         if (nOffset !=
    3119           0 :                             nLastOffset +
    3120           0 :                                 static_cast<GIntBig>(nBlockSizeBytes) * nBands)
    3121             :                         {
    3122           0 :                             CPLError(CE_Warning, CPLE_AppDefined,
    3123             :                                      "Block %d,%d not at expected "
    3124             :                                      "offset",
    3125             :                                      0, y);
    3126           0 :                             return false;
    3127             :                         }
    3128             :                     }
    3129           1 :                     nLastOffset = nOffset;
    3130             :                 }
    3131             :                 else
    3132             :                 {
    3133           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
    3134             :                              "Block %d,%d not at expected "
    3135             :                              "offset",
    3136             :                              0, y);
    3137           0 :                     return false;
    3138             :                 }
    3139             :             }
    3140             :         }
    3141             : 
    3142          11 :         return true;
    3143             :     }
    3144             : 
    3145          66 :     int bHasNoData = FALSE;
    3146          66 :     const double dfNoData = GetRasterBand(1)->GetNoDataValue(&bHasNoData);
    3147          66 :     const GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
    3148          66 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
    3149          66 :     const vsi_l_offset nFileSize = static_cast<vsi_l_offset>(nRasterXSize) *
    3150          66 :                                    nRasterYSize * nBands * nDTSize;
    3151          66 :     if (dfNoData == 0 || !bHasNoData)
    3152             :     {
    3153          62 :         if (VSIFTruncateL(m_fpImage, nFileSize) != 0)
    3154             :         {
    3155           0 :             CPLError(CE_Failure, CPLE_FileIO,
    3156             :                      "Cannot create file of size " CPL_FRMT_GUIB " bytes",
    3157             :                      nFileSize);
    3158           0 :             return false;
    3159             :         }
    3160             :     }
    3161             :     else
    3162             :     {
    3163           4 :         size_t nLineSize = static_cast<size_t>(nRasterXSize) * nDTSize;
    3164           4 :         void *pData = VSI_MALLOC_VERBOSE(nLineSize);
    3165           4 :         if (pData == nullptr)
    3166           0 :             return false;
    3167           4 :         GDALCopyWords(&dfNoData, GDT_Float64, 0, pData, eDT, nDTSize,
    3168             :                       nRasterXSize);
    3169             : #ifdef CPL_MSB
    3170             :         if (GDALDataTypeIsComplex(eDT))
    3171             :         {
    3172             :             GDALSwapWords(pData, nDTSize / 2, nRasterXSize * 2, nDTSize / 2);
    3173             :         }
    3174             :         else
    3175             :         {
    3176             :             GDALSwapWords(pData, nDTSize, nRasterXSize, nDTSize);
    3177             :         }
    3178             : #endif
    3179           8 :         for (vsi_l_offset i = 0;
    3180           8 :              i < static_cast<vsi_l_offset>(nRasterYSize) * nBands; i++)
    3181             :         {
    3182           4 :             size_t nBytesWritten = VSIFWriteL(pData, 1, nLineSize, m_fpImage);
    3183           4 :             if (nBytesWritten != nLineSize)
    3184             :             {
    3185           0 :                 CPLError(CE_Failure, CPLE_FileIO,
    3186             :                          "Cannot create file of size " CPL_FRMT_GUIB " bytes",
    3187             :                          nFileSize);
    3188           0 :                 VSIFree(pData);
    3189           0 :                 return false;
    3190             :             }
    3191             :         }
    3192           4 :         VSIFree(pData);
    3193             :     }
    3194          66 :     return true;
    3195             : }
    3196             : 
    3197             : /************************************************************************/
    3198             : /*                          GetSpecialConstants()                       */
    3199             : /************************************************************************/
    3200             : 
    3201          12 : static CPLXMLNode *GetSpecialConstants(const CPLString &osPrefix,
    3202             :                                        CPLXMLNode *psFileAreaObservational)
    3203             : {
    3204          27 :     for (CPLXMLNode *psIter = psFileAreaObservational->psChild; psIter;
    3205          15 :          psIter = psIter->psNext)
    3206             :     {
    3207          48 :         if (psIter->eType == CXT_Element &&
    3208          48 :             STARTS_WITH(psIter->pszValue, (osPrefix + "Array").c_str()))
    3209             :         {
    3210             :             CPLXMLNode *psSC =
    3211          11 :                 CPLGetXMLNode(psIter, (osPrefix + "Special_Constants").c_str());
    3212          11 :             if (psSC)
    3213             :             {
    3214           9 :                 CPLXMLNode *psNext = psSC->psNext;
    3215           9 :                 psSC->psNext = nullptr;
    3216           9 :                 CPLXMLNode *psRet = CPLCloneXMLTree(psSC);
    3217           9 :                 psSC->psNext = psNext;
    3218           9 :                 return psRet;
    3219             :             }
    3220             :         }
    3221             :     }
    3222           3 :     return nullptr;
    3223             : }
    3224             : 
    3225             : /************************************************************************/
    3226             : /*                          WriteHeaderAppendCase()                     */
    3227             : /************************************************************************/
    3228             : 
    3229           4 : void PDS4Dataset::WriteHeaderAppendCase()
    3230             : {
    3231           4 :     CPLXMLTreeCloser oCloser(CPLParseXMLFile(GetDescription()));
    3232           4 :     CPLXMLNode *psRoot = oCloser.get();
    3233           4 :     if (psRoot == nullptr)
    3234           0 :         return;
    3235           4 :     CPLString osPrefix;
    3236           4 :     CPLXMLNode *psProduct = CPLGetXMLNode(psRoot, "=Product_Observational");
    3237           4 :     if (psProduct == nullptr)
    3238             :     {
    3239           0 :         psProduct = CPLGetXMLNode(psRoot, "=pds:Product_Observational");
    3240           0 :         if (psProduct)
    3241           0 :             osPrefix = "pds:";
    3242             :     }
    3243           4 :     if (psProduct == nullptr)
    3244             :     {
    3245           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    3246             :                  "Cannot find Product_Observational element");
    3247           0 :         return;
    3248             :     }
    3249           4 :     CPLXMLNode *psFAO = CPLGetXMLNode(
    3250           8 :         psProduct, (osPrefix + "File_Area_Observational").c_str());
    3251           4 :     if (psFAO == nullptr)
    3252             :     {
    3253           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    3254             :                  "Cannot find File_Area_Observational element");
    3255           0 :         return;
    3256             :     }
    3257             : 
    3258           4 :     WriteArray(osPrefix, psFAO, nullptr, nullptr);
    3259             : 
    3260           4 :     CPLSerializeXMLTreeToFile(psRoot, GetDescription());
    3261             : }
    3262             : 
    3263             : /************************************************************************/
    3264             : /*                              WriteArray()                            */
    3265             : /************************************************************************/
    3266             : 
    3267         119 : void PDS4Dataset::WriteArray(const CPLString &osPrefix, CPLXMLNode *psFAO,
    3268             :                              const char *pszLocalIdentifierDefault,
    3269             :                              CPLXMLNode *psTemplateSpecialConstants)
    3270             : {
    3271         238 :     const char *pszArrayType = CSLFetchNameValueDef(
    3272         119 :         m_papszCreationOptions, "ARRAY_TYPE", "Array_3D_Image");
    3273         119 :     const bool bIsArray2D = STARTS_WITH(pszArrayType, "Array_2D");
    3274             :     CPLXMLNode *psArray =
    3275         119 :         CPLCreateXMLNode(psFAO, CXT_Element, (osPrefix + pszArrayType).c_str());
    3276             : 
    3277         238 :     const char *pszLocalIdentifier = CSLFetchNameValueDef(
    3278         119 :         m_papszCreationOptions, "ARRAY_IDENTIFIER", pszLocalIdentifierDefault);
    3279         119 :     if (pszLocalIdentifier)
    3280             :     {
    3281         115 :         CPLCreateXMLElementAndValue(psArray,
    3282         230 :                                     (osPrefix + "local_identifier").c_str(),
    3283             :                                     pszLocalIdentifier);
    3284             :     }
    3285             : 
    3286         119 :     GUIntBig nOffset = m_nBaseOffset;
    3287         119 :     if (m_poExternalDS)
    3288             :     {
    3289             :         const char *pszOffset =
    3290          11 :             m_poExternalDS->GetRasterBand(1)->GetMetadataItem(
    3291          11 :                 "BLOCK_OFFSET_0_0", "TIFF");
    3292          11 :         if (pszOffset)
    3293          11 :             nOffset = CPLAtoGIntBig(pszOffset);
    3294             :     }
    3295         238 :     CPLAddXMLAttributeAndValue(
    3296         238 :         CPLCreateXMLElementAndValue(psArray, (osPrefix + "offset").c_str(),
    3297             :                                     CPLSPrintf(CPL_FRMT_GUIB, nOffset)),
    3298             :         "unit", "byte");
    3299         119 :     CPLCreateXMLElementAndValue(psArray, (osPrefix + "axes").c_str(),
    3300             :                                 (bIsArray2D) ? "2" : "3");
    3301         119 :     CPLCreateXMLElementAndValue(
    3302         238 :         psArray, (osPrefix + "axis_index_order").c_str(), "Last Index Fastest");
    3303         119 :     CPLXMLNode *psElementArray = CPLCreateXMLNode(
    3304         238 :         psArray, CXT_Element, (osPrefix + "Element_Array").c_str());
    3305         119 :     GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
    3306         119 :     const char *pszDataType =
    3307         157 :         (eDT == GDT_Byte)     ? "UnsignedByte"
    3308          73 :         : (eDT == GDT_Int8)   ? "SignedByte"
    3309          65 :         : (eDT == GDT_UInt16) ? "UnsignedLSB2"
    3310          55 :         : (eDT == GDT_Int16)  ? (m_bIsLSB ? "SignedLSB2" : "SignedMSB2")
    3311          46 :         : (eDT == GDT_UInt32) ? (m_bIsLSB ? "UnsignedLSB4" : "UnsignedMSB4")
    3312          38 :         : (eDT == GDT_Int32)  ? (m_bIsLSB ? "SignedLSB4" : "SignedMSB4")
    3313             :         : (eDT == GDT_Float32)
    3314          29 :             ? (m_bIsLSB ? "IEEE754LSBSingle" : "IEEE754MSBSingle")
    3315             :         : (eDT == GDT_Float64)
    3316          20 :             ? (m_bIsLSB ? "IEEE754LSBDouble" : "IEEE754MSBDouble")
    3317          12 :         : (eDT == GDT_CFloat32) ? (m_bIsLSB ? "ComplexLSB8" : "ComplexMSB8")
    3318           4 :         : (eDT == GDT_CFloat64) ? (m_bIsLSB ? "ComplexLSB16" : "ComplexMSB16")
    3319             :                                 : "should not happen";
    3320         119 :     CPLCreateXMLElementAndValue(psElementArray,
    3321         238 :                                 (osPrefix + "data_type").c_str(), pszDataType);
    3322             : 
    3323         119 :     const char *pszUnits = GetRasterBand(1)->GetUnitType();
    3324         119 :     const char *pszUnitsCO = CSLFetchNameValue(m_papszCreationOptions, "UNIT");
    3325         119 :     if (pszUnitsCO)
    3326             :     {
    3327           0 :         pszUnits = pszUnitsCO;
    3328             :     }
    3329         119 :     if (pszUnits && pszUnits[0] != 0)
    3330             :     {
    3331           1 :         CPLCreateXMLElementAndValue(psElementArray, (osPrefix + "unit").c_str(),
    3332             :                                     pszUnits);
    3333             :     }
    3334             : 
    3335         119 :     int bHasScale = FALSE;
    3336         119 :     double dfScale = GetRasterBand(1)->GetScale(&bHasScale);
    3337         119 :     if (bHasScale && dfScale != 1.0)
    3338             :     {
    3339          10 :         CPLCreateXMLElementAndValue(psElementArray,
    3340          10 :                                     (osPrefix + "scaling_factor").c_str(),
    3341             :                                     CPLSPrintf("%.17g", dfScale));
    3342             :     }
    3343             : 
    3344         119 :     int bHasOffset = FALSE;
    3345         119 :     double dfOffset = GetRasterBand(1)->GetOffset(&bHasOffset);
    3346         119 :     if (bHasOffset && dfOffset != 1.0)
    3347             :     {
    3348          10 :         CPLCreateXMLElementAndValue(psElementArray,
    3349          10 :                                     (osPrefix + "value_offset").c_str(),
    3350             :                                     CPLSPrintf("%.17g", dfOffset));
    3351             :     }
    3352             : 
    3353             :     // Axis definitions
    3354             :     {
    3355         119 :         CPLXMLNode *psAxis = CPLCreateXMLNode(
    3356         238 :             psArray, CXT_Element, (osPrefix + "Axis_Array").c_str());
    3357         238 :         CPLCreateXMLElementAndValue(
    3358         238 :             psAxis, (osPrefix + "axis_name").c_str(),
    3359         119 :             EQUAL(m_osInterleave, "BSQ")
    3360             :                 ? "Band"
    3361             :                 :
    3362             :                 /* EQUAL(m_osInterleave, "BIL") ? "Line" : */
    3363             :                 "Line");
    3364         238 :         CPLCreateXMLElementAndValue(
    3365         238 :             psAxis, (osPrefix + "elements").c_str(),
    3366             :             CPLSPrintf("%d",
    3367         119 :                        EQUAL(m_osInterleave, "BSQ")
    3368             :                            ? nBands
    3369             :                            :
    3370             :                            /* EQUAL(m_osInterleave, "BIL") ? nRasterYSize : */
    3371             :                            nRasterYSize));
    3372         119 :         CPLCreateXMLElementAndValue(
    3373         238 :             psAxis, (osPrefix + "sequence_number").c_str(), "1");
    3374             :     }
    3375             :     {
    3376         119 :         CPLXMLNode *psAxis = CPLCreateXMLNode(
    3377         238 :             psArray, CXT_Element, (osPrefix + "Axis_Array").c_str());
    3378         119 :         CPLCreateXMLElementAndValue(psAxis, (osPrefix + "axis_name").c_str(),
    3379         119 :                                     EQUAL(m_osInterleave, "BSQ")   ? "Line"
    3380          11 :                                     : EQUAL(m_osInterleave, "BIL") ? "Band"
    3381             :                                                                    : "Sample");
    3382         238 :         CPLCreateXMLElementAndValue(
    3383         238 :             psAxis, (osPrefix + "elements").c_str(),
    3384         119 :             CPLSPrintf("%d", EQUAL(m_osInterleave, "BSQ")   ? nRasterYSize
    3385          11 :                              : EQUAL(m_osInterleave, "BIL") ? nBands
    3386             :                                                             : nRasterXSize));
    3387         119 :         CPLCreateXMLElementAndValue(
    3388         238 :             psAxis, (osPrefix + "sequence_number").c_str(), "2");
    3389             :     }
    3390         119 :     if (!bIsArray2D)
    3391             :     {
    3392         118 :         CPLXMLNode *psAxis = CPLCreateXMLNode(
    3393         236 :             psArray, CXT_Element, (osPrefix + "Axis_Array").c_str());
    3394         118 :         CPLCreateXMLElementAndValue(psAxis, (osPrefix + "axis_name").c_str(),
    3395         118 :                                     EQUAL(m_osInterleave, "BSQ")   ? "Sample"
    3396          10 :                                     : EQUAL(m_osInterleave, "BIL") ? "Sample"
    3397             :                                                                    : "Band");
    3398         236 :         CPLCreateXMLElementAndValue(
    3399         236 :             psAxis, (osPrefix + "elements").c_str(),
    3400         118 :             CPLSPrintf("%d", EQUAL(m_osInterleave, "BSQ")   ? nRasterXSize
    3401          10 :                              : EQUAL(m_osInterleave, "BIL") ? nRasterXSize
    3402             :                                                             : nBands));
    3403         118 :         CPLCreateXMLElementAndValue(
    3404         236 :             psAxis, (osPrefix + "sequence_number").c_str(), "3");
    3405             :     }
    3406             : 
    3407         119 :     int bHasNoData = FALSE;
    3408         119 :     double dfNoData = GetRasterBand(1)->GetNoDataValue(&bHasNoData);
    3409         119 :     if (psTemplateSpecialConstants)
    3410             :     {
    3411           9 :         CPLAddXMLChild(psArray, psTemplateSpecialConstants);
    3412           9 :         if (bHasNoData)
    3413             :         {
    3414             :             CPLXMLNode *psMC =
    3415           6 :                 CPLGetXMLNode(psTemplateSpecialConstants,
    3416          12 :                               (osPrefix + "missing_constant").c_str());
    3417           6 :             if (psMC != nullptr)
    3418             :             {
    3419           4 :                 if (psMC->psChild && psMC->psChild->eType == CXT_Text)
    3420             :                 {
    3421           4 :                     CPLFree(psMC->psChild->pszValue);
    3422           8 :                     psMC->psChild->pszValue =
    3423           4 :                         CPLStrdup(CPLSPrintf("%.17g", dfNoData));
    3424             :                 }
    3425             :             }
    3426             :             else
    3427             :             {
    3428             :                 CPLXMLNode *psSaturatedConstant =
    3429           2 :                     CPLGetXMLNode(psTemplateSpecialConstants,
    3430           4 :                                   (osPrefix + "saturated_constant").c_str());
    3431           4 :                 psMC = CPLCreateXMLElementAndValue(
    3432           4 :                     nullptr, (osPrefix + "missing_constant").c_str(),
    3433             :                     CPLSPrintf("%.17g", dfNoData));
    3434             :                 CPLXMLNode *psNext;
    3435           2 :                 if (psSaturatedConstant)
    3436             :                 {
    3437           1 :                     psNext = psSaturatedConstant->psNext;
    3438           1 :                     psSaturatedConstant->psNext = psMC;
    3439             :                 }
    3440             :                 else
    3441             :                 {
    3442           1 :                     psNext = psTemplateSpecialConstants->psChild;
    3443           1 :                     psTemplateSpecialConstants->psChild = psMC;
    3444             :                 }
    3445           2 :                 psMC->psNext = psNext;
    3446             :             }
    3447             :         }
    3448             :     }
    3449         110 :     else if (bHasNoData)
    3450             :     {
    3451          11 :         CPLXMLNode *psSC = CPLCreateXMLNode(
    3452          22 :             psArray, CXT_Element, (osPrefix + "Special_Constants").c_str());
    3453          22 :         CPLCreateXMLElementAndValue(psSC,
    3454          22 :                                     (osPrefix + "missing_constant").c_str(),
    3455             :                                     CPLSPrintf("%.17g", dfNoData));
    3456             :     }
    3457         119 : }
    3458             : 
    3459             : /************************************************************************/
    3460             : /*                          WriteVectorLayers()                         */
    3461             : /************************************************************************/
    3462             : 
    3463         172 : void PDS4Dataset::WriteVectorLayers(CPLXMLNode *psProduct)
    3464             : {
    3465         344 :     CPLString osPrefix;
    3466         172 :     if (STARTS_WITH(psProduct->pszValue, "pds:"))
    3467           0 :         osPrefix = "pds:";
    3468             : 
    3469         244 :     for (auto &poLayer : m_apoLayers)
    3470             :     {
    3471          72 :         if (!poLayer->IsDirtyHeader())
    3472           3 :             continue;
    3473             : 
    3474          69 :         if (poLayer->GetFeatureCount(false) == 0)
    3475             :         {
    3476          16 :             CPLError(CE_Warning, CPLE_AppDefined,
    3477             :                      "Writing header for layer %s which has 0 features. "
    3478             :                      "This is not legal in PDS4",
    3479          16 :                      poLayer->GetName());
    3480             :         }
    3481             : 
    3482          69 :         if (poLayer->GetRawFieldCount() == 0)
    3483             :         {
    3484          16 :             CPLError(CE_Warning, CPLE_AppDefined,
    3485             :                      "Writing header for layer %s which has 0 fields. "
    3486             :                      "This is not legal in PDS4",
    3487          16 :                      poLayer->GetName());
    3488             :         }
    3489             : 
    3490             :         const std::string osRelativePath(
    3491         138 :             CPLExtractRelativePath(CPLGetPathSafe(m_osXMLFilename).c_str(),
    3492         207 :                                    poLayer->GetFileName(), nullptr));
    3493             : 
    3494          69 :         bool bFound = false;
    3495         634 :         for (CPLXMLNode *psIter = psProduct->psChild; psIter != nullptr;
    3496         565 :              psIter = psIter->psNext)
    3497             :         {
    3498         733 :             if (psIter->eType == CXT_Element &&
    3499         164 :                 strcmp(psIter->pszValue,
    3500         733 :                        (osPrefix + "File_Area_Observational").c_str()) == 0)
    3501             :             {
    3502          23 :                 const char *pszFilename = CPLGetXMLValue(
    3503             :                     psIter,
    3504          46 :                     (osPrefix + "File." + osPrefix + "file_name").c_str(), "");
    3505          23 :                 if (strcmp(pszFilename, osRelativePath.c_str()) == 0)
    3506             :                 {
    3507           4 :                     poLayer->RefreshFileAreaObservational(psIter);
    3508           4 :                     bFound = true;
    3509           4 :                     break;
    3510             :                 }
    3511             :             }
    3512             :         }
    3513          69 :         if (!bFound)
    3514             :         {
    3515          65 :             CPLXMLNode *psFAO = CPLCreateXMLNode(
    3516             :                 psProduct, CXT_Element,
    3517         130 :                 (osPrefix + "File_Area_Observational").c_str());
    3518          65 :             CPLXMLNode *psFile = CPLCreateXMLNode(psFAO, CXT_Element,
    3519         130 :                                                   (osPrefix + "File").c_str());
    3520         130 :             CPLCreateXMLElementAndValue(psFile,
    3521         130 :                                         (osPrefix + "file_name").c_str(),
    3522             :                                         osRelativePath.c_str());
    3523          65 :             poLayer->RefreshFileAreaObservational(psFAO);
    3524             :         }
    3525             :     }
    3526         172 : }
    3527             : 
    3528             : /************************************************************************/
    3529             : /*                            CreateHeader()                            */
    3530             : /************************************************************************/
    3531             : 
    3532         165 : void PDS4Dataset::CreateHeader(CPLXMLNode *psProduct,
    3533             :                                const char *pszCARTVersion)
    3534             : {
    3535         165 :     CPLString osPrefix;
    3536         165 :     if (STARTS_WITH(psProduct->pszValue, "pds:"))
    3537           0 :         osPrefix = "pds:";
    3538             : 
    3539         165 :     OGREnvelope sExtent;
    3540         211 :     if (m_oSRS.IsEmpty() && GetLayerCount() >= 1 &&
    3541          46 :         GetLayer(0)->GetSpatialRef() != nullptr)
    3542             :     {
    3543           2 :         const auto poSRS = GetLayer(0)->GetSpatialRef();
    3544           2 :         if (poSRS)
    3545           2 :             m_oSRS = *poSRS;
    3546             :     }
    3547             : 
    3548         258 :     if (!m_oSRS.IsEmpty() &&
    3549          93 :         CSLFetchNameValue(m_papszCreationOptions, "VAR_TARGET") == nullptr)
    3550             :     {
    3551          92 :         const char *pszTarget = nullptr;
    3552          92 :         if (fabs(m_oSRS.GetSemiMajor() - 6378137) < 0.001 * 6378137)
    3553             :         {
    3554          71 :             pszTarget = "Earth";
    3555          71 :             m_papszCreationOptions = CSLSetNameValue(
    3556             :                 m_papszCreationOptions, "VAR_TARGET_TYPE", "Planet");
    3557             :         }
    3558             :         else
    3559             :         {
    3560          21 :             const char *pszDatum = m_oSRS.GetAttrValue("DATUM");
    3561          21 :             if (pszDatum && STARTS_WITH(pszDatum, "D_"))
    3562             :             {
    3563           3 :                 pszTarget = pszDatum + 2;
    3564             :             }
    3565          18 :             else if (pszDatum)
    3566             :             {
    3567          18 :                 pszTarget = pszDatum;
    3568             :             }
    3569             :         }
    3570          92 :         if (pszTarget)
    3571             :         {
    3572          92 :             m_papszCreationOptions = CSLSetNameValue(m_papszCreationOptions,
    3573             :                                                      "VAR_TARGET", pszTarget);
    3574             :         }
    3575             :     }
    3576         165 :     SubstituteVariables(psProduct, m_papszCreationOptions);
    3577             : 
    3578             :     // Remove <Discipline_Area>/<disp:Display_Settings> if there is no raster
    3579         165 :     if (GetRasterCount() == 0)
    3580             :     {
    3581             :         CPLXMLNode *psDisciplineArea =
    3582         141 :             CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area." +
    3583          94 :                                       osPrefix + "Discipline_Area")
    3584             :                                          .c_str());
    3585          47 :         if (psDisciplineArea)
    3586             :         {
    3587             :             CPLXMLNode *psDisplaySettings =
    3588          47 :                 CPLGetXMLNode(psDisciplineArea, "disp:Display_Settings");
    3589          47 :             if (psDisplaySettings)
    3590             :             {
    3591          47 :                 CPLRemoveXMLChild(psDisciplineArea, psDisplaySettings);
    3592          47 :                 CPLDestroyXMLNode(psDisplaySettings);
    3593             :             }
    3594             :         }
    3595             :     }
    3596             : 
    3597             :     // Depending on the version of the DISP schema, Local_Internal_Reference
    3598             :     // may be in the disp: namespace or the default one.
    3599             :     const auto GetLocalIdentifierReferenceFromDisciplineArea =
    3600         202 :         [](const CPLXMLNode *psDisciplineArea, const char *pszDefault)
    3601             :     {
    3602         202 :         return CPLGetXMLValue(
    3603             :             psDisciplineArea,
    3604             :             "disp:Display_Settings.Local_Internal_Reference."
    3605             :             "local_identifier_reference",
    3606             :             CPLGetXMLValue(
    3607             :                 psDisciplineArea,
    3608             :                 "disp:Display_Settings.disp:Local_Internal_Reference."
    3609             :                 "local_identifier_reference",
    3610         202 :                 pszDefault));
    3611             :     };
    3612             : 
    3613         165 :     if (GetRasterCount() || !m_oSRS.IsEmpty())
    3614             :     {
    3615             :         CPLXMLNode *psDisciplineArea =
    3616         360 :             CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area." +
    3617         240 :                                       osPrefix + "Discipline_Area")
    3618             :                                          .c_str());
    3619         120 :         if (GetRasterCount() && !(m_bGotTransform && !m_oSRS.IsEmpty()))
    3620             :         {
    3621             :             // if we have no georeferencing, strip any existing georeferencing
    3622             :             // from the template
    3623          27 :             if (psDisciplineArea)
    3624             :             {
    3625             :                 CPLXMLNode *psCart =
    3626          23 :                     CPLGetXMLNode(psDisciplineArea, "cart:Cartography");
    3627          23 :                 if (psCart == nullptr)
    3628          22 :                     psCart = CPLGetXMLNode(psDisciplineArea, "Cartography");
    3629          23 :                 if (psCart)
    3630             :                 {
    3631           1 :                     CPLRemoveXMLChild(psDisciplineArea, psCart);
    3632           1 :                     CPLDestroyXMLNode(psCart);
    3633             :                 }
    3634             : 
    3635          23 :                 if (CPLGetXMLNode(psDisciplineArea,
    3636          23 :                                   "sp:Spectral_Characteristics"))
    3637             :                 {
    3638             :                     const char *pszArrayType =
    3639           1 :                         CSLFetchNameValue(m_papszCreationOptions, "ARRAY_TYPE");
    3640             :                     // The schematron PDS4_SP_1100.sch requires that
    3641             :                     // sp:local_identifier_reference is used by
    3642             :                     // Array_[2D|3D]_Spectrum/pds:local_identifier
    3643           1 :                     if (pszArrayType == nullptr)
    3644             :                     {
    3645           1 :                         m_papszCreationOptions =
    3646           1 :                             CSLSetNameValue(m_papszCreationOptions,
    3647             :                                             "ARRAY_TYPE", "Array_3D_Spectrum");
    3648             :                     }
    3649           0 :                     else if (!EQUAL(pszArrayType, "Array_2D_Spectrum") &&
    3650           0 :                              !EQUAL(pszArrayType, "Array_3D_Spectrum"))
    3651             :                     {
    3652           0 :                         CPLError(CE_Warning, CPLE_AppDefined,
    3653             :                                  "PDS4_SP_xxxx.sch schematron requires the "
    3654             :                                  "use of ARRAY_TYPE=Array_2D_Spectrum or "
    3655             :                                  "Array_3D_Spectrum");
    3656             :                     }
    3657             :                 }
    3658             :             }
    3659             :         }
    3660             :         else
    3661             :         {
    3662          93 :             if (psDisciplineArea == nullptr)
    3663             :             {
    3664           2 :                 CPLXMLNode *psTI = CPLGetXMLNode(
    3665           4 :                     psProduct, (osPrefix + "Observation_Area." + osPrefix +
    3666             :                                 "Target_Identification")
    3667             :                                    .c_str());
    3668           2 :                 if (psTI == nullptr)
    3669             :                 {
    3670           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    3671             :                              "Cannot find Target_Identification element in "
    3672             :                              "template");
    3673           1 :                     return;
    3674             :                 }
    3675             :                 psDisciplineArea =
    3676           1 :                     CPLCreateXMLNode(nullptr, CXT_Element,
    3677           2 :                                      (osPrefix + "Discipline_Area").c_str());
    3678           1 :                 if (psTI->psNext)
    3679           0 :                     psDisciplineArea->psNext = psTI->psNext;
    3680           1 :                 psTI->psNext = psDisciplineArea;
    3681             :             }
    3682             :             CPLXMLNode *psCart =
    3683          92 :                 CPLGetXMLNode(psDisciplineArea, "cart:Cartography");
    3684          92 :             if (psCart == nullptr)
    3685          87 :                 psCart = CPLGetXMLNode(psDisciplineArea, "Cartography");
    3686          92 :             if (psCart == nullptr)
    3687             :             {
    3688          87 :                 psCart = CPLCreateXMLNode(psDisciplineArea, CXT_Element,
    3689             :                                           "cart:Cartography");
    3690          87 :                 if (CPLGetXMLNode(psProduct, "xmlns:cart") == nullptr)
    3691             :                 {
    3692             :                     CPLXMLNode *psNS =
    3693           1 :                         CPLCreateXMLNode(nullptr, CXT_Attribute, "xmlns:cart");
    3694           1 :                     CPLCreateXMLNode(psNS, CXT_Text,
    3695             :                                      "http://pds.nasa.gov/pds4/cart/v1");
    3696           1 :                     CPLAddXMLChild(psProduct, psNS);
    3697             :                     CPLXMLNode *psSchemaLoc =
    3698           1 :                         CPLGetXMLNode(psProduct, "xsi:schemaLocation");
    3699           1 :                     if (psSchemaLoc != nullptr &&
    3700           1 :                         psSchemaLoc->psChild != nullptr &&
    3701           1 :                         psSchemaLoc->psChild->pszValue != nullptr)
    3702             :                     {
    3703           2 :                         CPLString osCartSchema;
    3704           1 :                         if (strstr(psSchemaLoc->psChild->pszValue,
    3705             :                                    "PDS4_PDS_1800.xsd"))
    3706             :                         {
    3707             :                             // GDAL 2.4
    3708             :                             osCartSchema = "https://pds.nasa.gov/pds4/cart/v1/"
    3709           1 :                                            "PDS4_CART_1700.xsd";
    3710           1 :                             pszCARTVersion = "1700";
    3711             :                         }
    3712           0 :                         else if (strstr(psSchemaLoc->psChild->pszValue,
    3713             :                                         "PDS4_PDS_1B00.xsd"))
    3714             :                         {
    3715             :                             // GDAL 3.0
    3716             :                             osCartSchema =
    3717             :                                 "https://raw.githubusercontent.com/"
    3718             :                                 "nasa-pds-data-dictionaries/ldd-cart/master/"
    3719           0 :                                 "build/1.B.0.0/PDS4_CART_1B00.xsd";
    3720           0 :                             pszCARTVersion = "1B00";
    3721             :                         }
    3722           0 :                         else if (strstr(psSchemaLoc->psChild->pszValue,
    3723             :                                         "PDS4_PDS_1D00.xsd"))
    3724             :                         {
    3725             :                             // GDAL 3.1
    3726             :                             osCartSchema = "https://pds.nasa.gov/pds4/cart/v1/"
    3727           0 :                                            "PDS4_CART_1D00_1933.xsd";
    3728           0 :                             pszCARTVersion = "1D00_1933";
    3729             :                         }
    3730             :                         else
    3731             :                         {
    3732             :                             // GDAL 3.4
    3733             :                             osCartSchema =
    3734             :                                 "https://pds.nasa.gov/pds4/cart/v1/"
    3735           0 :                                 "PDS4_CART_" CURRENT_CART_VERSION ".xsd";
    3736           0 :                             pszCARTVersion = CURRENT_CART_VERSION;
    3737             :                         }
    3738           1 :                         CPLString osNewVal(psSchemaLoc->psChild->pszValue);
    3739             :                         osNewVal +=
    3740           1 :                             " http://pds.nasa.gov/pds4/cart/v1 " + osCartSchema;
    3741           1 :                         CPLFree(psSchemaLoc->psChild->pszValue);
    3742           1 :                         psSchemaLoc->psChild->pszValue = CPLStrdup(osNewVal);
    3743             :                     }
    3744             :                 }
    3745             :             }
    3746             :             else
    3747             :             {
    3748           5 :                 if (psCart->psChild)
    3749             :                 {
    3750           5 :                     CPLDestroyXMLNode(psCart->psChild);
    3751           5 :                     psCart->psChild = nullptr;
    3752             :                 }
    3753             :             }
    3754             : 
    3755          92 :             if (IsCARTVersionGTE(pszCARTVersion, "1900"))
    3756             :             {
    3757             :                 const char *pszLocalIdentifier =
    3758          87 :                     GetLocalIdentifierReferenceFromDisciplineArea(
    3759             :                         psDisciplineArea,
    3760          87 :                         GetRasterCount() == 0 && GetLayerCount() > 0
    3761           2 :                             ? GetLayer(0)->GetName()
    3762             :                             : "image");
    3763          87 :                 CPLXMLNode *psLIR = CPLCreateXMLNode(
    3764             :                     psCart, CXT_Element,
    3765         174 :                     (osPrefix + "Local_Internal_Reference").c_str());
    3766          87 :                 CPLCreateXMLElementAndValue(
    3767         174 :                     psLIR, (osPrefix + "local_identifier_reference").c_str(),
    3768             :                     pszLocalIdentifier);
    3769          87 :                 CPLCreateXMLElementAndValue(
    3770         174 :                     psLIR, (osPrefix + "local_reference_type").c_str(),
    3771             :                     "cartography_parameters_to_image_object");
    3772             :             }
    3773             : 
    3774          92 :             WriteGeoreferencing(psCart, pszCARTVersion);
    3775             :         }
    3776             : 
    3777         238 :         const char *pszVertDir = CSLFetchNameValue(
    3778         119 :             m_papszCreationOptions, "VAR_VERTICAL_DISPLAY_DIRECTION");
    3779         119 :         if (pszVertDir)
    3780             :         {
    3781             :             CPLXMLNode *psVertDirNode =
    3782           1 :                 CPLGetXMLNode(psDisciplineArea,
    3783             :                               "disp:Display_Settings.disp:Display_Direction."
    3784             :                               "disp:vertical_display_direction");
    3785           1 :             if (psVertDirNode == nullptr)
    3786             :             {
    3787           0 :                 CPLError(
    3788             :                     CE_Warning, CPLE_AppDefined,
    3789             :                     "PDS4 template lacks a disp:vertical_display_direction "
    3790             :                     "element where to write %s",
    3791             :                     pszVertDir);
    3792             :             }
    3793             :             else
    3794             :             {
    3795           1 :                 CPLDestroyXMLNode(psVertDirNode->psChild);
    3796           1 :                 psVertDirNode->psChild =
    3797           1 :                     CPLCreateXMLNode(nullptr, CXT_Text, pszVertDir);
    3798             :             }
    3799             :         }
    3800             :     }
    3801             :     else
    3802             :     {
    3803             :         // Remove Observation_Area.Discipline_Area if it contains only
    3804             :         // <disp:Display_Settings> or is empty
    3805             :         CPLXMLNode *psObservationArea =
    3806          45 :             CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area").c_str());
    3807          45 :         if (psObservationArea)
    3808             :         {
    3809          45 :             CPLXMLNode *psDisciplineArea = CPLGetXMLNode(
    3810          90 :                 psObservationArea, (osPrefix + "Discipline_Area").c_str());
    3811          45 :             if (psDisciplineArea &&
    3812          45 :                 (psDisciplineArea->psChild == nullptr ||
    3813           0 :                  (psDisciplineArea->psChild->eType == CXT_Element &&
    3814           0 :                   psDisciplineArea->psChild->psNext == nullptr &&
    3815           0 :                   strcmp(psDisciplineArea->psChild->pszValue,
    3816             :                          "disp:Display_Settings") == 0)))
    3817             :             {
    3818          45 :                 CPLRemoveXMLChild(psObservationArea, psDisciplineArea);
    3819          45 :                 CPLDestroyXMLNode(psDisciplineArea);
    3820             :             }
    3821             :         }
    3822             :     }
    3823             : 
    3824         164 :     if (m_bStripFileAreaObservationalFromTemplate)
    3825             :     {
    3826         164 :         m_bStripFileAreaObservationalFromTemplate = false;
    3827         164 :         CPLXMLNode *psObservationArea = nullptr;
    3828         164 :         CPLXMLNode *psPrev = nullptr;
    3829         164 :         CPLXMLNode *psTemplateSpecialConstants = nullptr;
    3830        1474 :         for (CPLXMLNode *psIter = psProduct->psChild; psIter != nullptr;)
    3831             :         {
    3832        1649 :             if (psIter->eType == CXT_Element &&
    3833        1649 :                 psIter->pszValue == osPrefix + "Observation_Area")
    3834             :             {
    3835         163 :                 psObservationArea = psIter;
    3836         163 :                 psPrev = psIter;
    3837         163 :                 psIter = psIter->psNext;
    3838             :             }
    3839        1487 :             else if (psIter->eType == CXT_Element &&
    3840         176 :                      (psIter->pszValue ==
    3841        1487 :                           osPrefix + "File_Area_Observational" ||
    3842         164 :                       psIter->pszValue ==
    3843        1311 :                           osPrefix + "File_Area_Observational_Supplemental"))
    3844             :             {
    3845          12 :                 if (psIter->pszValue == osPrefix + "File_Area_Observational")
    3846             :                 {
    3847             :                     psTemplateSpecialConstants =
    3848          12 :                         GetSpecialConstants(osPrefix, psIter);
    3849             :                 }
    3850          12 :                 if (psPrev)
    3851          12 :                     psPrev->psNext = psIter->psNext;
    3852             :                 else
    3853             :                 {
    3854           0 :                     CPLAssert(psProduct->psChild == psIter);
    3855           0 :                     psProduct->psChild = psIter->psNext;
    3856             :                 }
    3857          12 :                 CPLXMLNode *psNext = psIter->psNext;
    3858          12 :                 psIter->psNext = nullptr;
    3859          12 :                 CPLDestroyXMLNode(psIter);
    3860          12 :                 psIter = psNext;
    3861             :             }
    3862             :             else
    3863             :             {
    3864        1135 :                 psPrev = psIter;
    3865        1135 :                 psIter = psIter->psNext;
    3866             :             }
    3867             :         }
    3868         164 :         if (psObservationArea == nullptr)
    3869             :         {
    3870           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    3871             :                      "Cannot find Observation_Area in template");
    3872           1 :             CPLDestroyXMLNode(psTemplateSpecialConstants);
    3873           1 :             return;
    3874             :         }
    3875             : 
    3876         163 :         if (GetRasterCount())
    3877             :         {
    3878         116 :             CPLXMLNode *psFAOPrev = psObservationArea;
    3879         118 :             while (psFAOPrev->psNext != nullptr &&
    3880           4 :                    psFAOPrev->psNext->eType == CXT_Comment)
    3881             :             {
    3882           2 :                 psFAOPrev = psFAOPrev->psNext;
    3883             :             }
    3884         116 :             if (psFAOPrev->psNext != nullptr)
    3885             :             {
    3886             :                 // There may be an optional Reference_List element between
    3887             :                 // Observation_Area and File_Area_Observational
    3888           4 :                 if (!(psFAOPrev->psNext->eType == CXT_Element &&
    3889           2 :                       psFAOPrev->psNext->pszValue ==
    3890           4 :                           osPrefix + "Reference_List"))
    3891             :                 {
    3892           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    3893             :                              "Unexpected content found after Observation_Area "
    3894             :                              "in template");
    3895           1 :                     CPLDestroyXMLNode(psTemplateSpecialConstants);
    3896           1 :                     return;
    3897             :                 }
    3898           1 :                 psFAOPrev = psFAOPrev->psNext;
    3899           2 :                 while (psFAOPrev->psNext != nullptr &&
    3900           1 :                        psFAOPrev->psNext->eType == CXT_Comment)
    3901             :                 {
    3902           1 :                     psFAOPrev = psFAOPrev->psNext;
    3903             :                 }
    3904           1 :                 if (psFAOPrev->psNext != nullptr)
    3905             :                 {
    3906           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    3907             :                              "Unexpected content found after Reference_List in "
    3908             :                              "template");
    3909           0 :                     CPLDestroyXMLNode(psTemplateSpecialConstants);
    3910           0 :                     return;
    3911             :                 }
    3912             :             }
    3913             : 
    3914         115 :             CPLXMLNode *psFAO = CPLCreateXMLNode(
    3915             :                 nullptr, CXT_Element,
    3916         230 :                 (osPrefix + "File_Area_Observational").c_str());
    3917         115 :             psFAOPrev->psNext = psFAO;
    3918             : 
    3919         115 :             CPLXMLNode *psFile = CPLCreateXMLNode(psFAO, CXT_Element,
    3920         230 :                                                   (osPrefix + "File").c_str());
    3921         230 :             CPLCreateXMLElementAndValue(psFile,
    3922         230 :                                         (osPrefix + "file_name").c_str(),
    3923             :                                         CPLGetFilename(m_osImageFilename));
    3924         115 :             if (m_bCreatedFromExistingBinaryFile)
    3925             :             {
    3926           7 :                 CPLCreateXMLNode(psFile, CXT_Comment, PREEXISTING_BINARY_FILE);
    3927             :             }
    3928             :             CPLXMLNode *psDisciplineArea =
    3929         345 :                 CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area." +
    3930         230 :                                           osPrefix + "Discipline_Area")
    3931             :                                              .c_str());
    3932             :             const char *pszLocalIdentifier =
    3933         115 :                 GetLocalIdentifierReferenceFromDisciplineArea(psDisciplineArea,
    3934             :                                                               "image");
    3935             : 
    3936         124 :             if (m_poExternalDS && m_poExternalDS->GetDriver() &&
    3937           9 :                 EQUAL(m_poExternalDS->GetDriver()->GetDescription(), "GTiff"))
    3938             :             {
    3939             :                 VSILFILE *fpTemp =
    3940           9 :                     VSIFOpenL(m_poExternalDS->GetDescription(), "rb");
    3941           9 :                 if (fpTemp)
    3942             :                 {
    3943           9 :                     GByte abySignature[4] = {0};
    3944           9 :                     VSIFReadL(abySignature, 1, 4, fpTemp);
    3945           9 :                     VSIFCloseL(fpTemp);
    3946           9 :                     const bool bBigTIFF =
    3947           9 :                         abySignature[2] == 43 || abySignature[3] == 43;
    3948             :                     m_osHeaderParsingStandard =
    3949           9 :                         bBigTIFF ? BIGTIFF_GEOTIFF_STRING : TIFF_GEOTIFF_STRING;
    3950             :                     const char *pszOffset =
    3951           9 :                         m_poExternalDS->GetRasterBand(1)->GetMetadataItem(
    3952           9 :                             "BLOCK_OFFSET_0_0", "TIFF");
    3953           9 :                     if (pszOffset)
    3954           9 :                         m_nBaseOffset = CPLAtoGIntBig(pszOffset);
    3955             :                 }
    3956             :             }
    3957             : 
    3958         115 :             if (!m_osHeaderParsingStandard.empty() && m_nBaseOffset > 0)
    3959             :             {
    3960          15 :                 CPLXMLNode *psHeader = CPLCreateXMLNode(
    3961          30 :                     psFAO, CXT_Element, (osPrefix + "Header").c_str());
    3962          15 :                 CPLAddXMLAttributeAndValue(
    3963             :                     CPLCreateXMLElementAndValue(
    3964          30 :                         psHeader, (osPrefix + "offset").c_str(), "0"),
    3965             :                     "unit", "byte");
    3966          15 :                 CPLAddXMLAttributeAndValue(
    3967             :                     CPLCreateXMLElementAndValue(
    3968          30 :                         psHeader, (osPrefix + "object_length").c_str(),
    3969             :                         CPLSPrintf(CPL_FRMT_GUIB,
    3970          15 :                                    static_cast<GUIntBig>(m_nBaseOffset))),
    3971             :                     "unit", "byte");
    3972          30 :                 CPLCreateXMLElementAndValue(
    3973          30 :                     psHeader, (osPrefix + "parsing_standard_id").c_str(),
    3974             :                     m_osHeaderParsingStandard.c_str());
    3975          15 :                 if (m_osHeaderParsingStandard == TIFF_GEOTIFF_STRING)
    3976             :                 {
    3977          11 :                     CPLCreateXMLElementAndValue(
    3978          22 :                         psHeader, (osPrefix + "description").c_str(),
    3979             :                         "TIFF/GeoTIFF header. The TIFF/GeoTIFF format is used "
    3980             :                         "throughout the geospatial and science communities "
    3981             :                         "to share geographic image data. ");
    3982             :                 }
    3983           4 :                 else if (m_osHeaderParsingStandard == BIGTIFF_GEOTIFF_STRING)
    3984             :                 {
    3985           0 :                     CPLCreateXMLElementAndValue(
    3986           0 :                         psHeader, (osPrefix + "description").c_str(),
    3987             :                         "BigTIFF/GeoTIFF header. The BigTIFF/GeoTIFF format is "
    3988             :                         "used "
    3989             :                         "throughout the geospatial and science communities "
    3990             :                         "to share geographic image data. ");
    3991             :                 }
    3992             :             }
    3993             : 
    3994         115 :             WriteArray(osPrefix, psFAO, pszLocalIdentifier,
    3995             :                        psTemplateSpecialConstants);
    3996             :         }
    3997             :     }
    3998             : }
    3999             : 
    4000             : /************************************************************************/
    4001             : /*                             WriteHeader()                            */
    4002             : /************************************************************************/
    4003             : 
    4004         177 : void PDS4Dataset::WriteHeader()
    4005             : {
    4006             :     const bool bAppend =
    4007         177 :         CPLFetchBool(m_papszCreationOptions, "APPEND_SUBDATASET", false);
    4008         177 :     if (bAppend)
    4009             :     {
    4010           4 :         WriteHeaderAppendCase();
    4011           5 :         return;
    4012             :     }
    4013             : 
    4014             :     CPLXMLNode *psRoot;
    4015         173 :     if (m_bCreateHeader)
    4016             :     {
    4017             :         CPLString osTemplateFilename =
    4018         166 :             CSLFetchNameValueDef(m_papszCreationOptions, "TEMPLATE", "");
    4019         166 :         if (!osTemplateFilename.empty())
    4020             :         {
    4021          20 :             if (STARTS_WITH(osTemplateFilename, "http://") ||
    4022          10 :                 STARTS_WITH(osTemplateFilename, "https://"))
    4023             :             {
    4024           0 :                 osTemplateFilename = "/vsicurl_streaming/" + osTemplateFilename;
    4025             :             }
    4026          10 :             psRoot = CPLParseXMLFile(osTemplateFilename);
    4027             :         }
    4028         156 :         else if (!m_osXMLPDS4.empty())
    4029           6 :             psRoot = CPLParseXMLString(m_osXMLPDS4);
    4030             :         else
    4031             :         {
    4032             : #ifndef USE_ONLY_EMBEDDED_RESOURCE_FILES
    4033             : #ifdef EMBED_RESOURCE_FILES
    4034             :             CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
    4035             : #endif
    4036             :             const char *pszDefaultTemplateFilename =
    4037         150 :                 CPLFindFile("gdal", "pds4_template.xml");
    4038         150 :             if (pszDefaultTemplateFilename)
    4039             :             {
    4040         150 :                 psRoot = CPLParseXMLFile(pszDefaultTemplateFilename);
    4041             :             }
    4042             :             else
    4043             : #endif
    4044             :             {
    4045             : #ifdef EMBED_RESOURCE_FILES
    4046             :                 static const bool bOnce [[maybe_unused]] = []()
    4047             :                 {
    4048             :                     CPLDebug("PDS4", "Using embedded pds4_template.xml");
    4049             :                     return true;
    4050             :                 }();
    4051             :                 psRoot = CPLParseXMLString(PDS4GetEmbeddedTemplate());
    4052             : #else
    4053           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    4054             :                          "Cannot find pds4_template.xml and TEMPLATE "
    4055             :                          "creation option not specified");
    4056           0 :                 return;
    4057             : #endif
    4058             :             }
    4059             :         }
    4060             :     }
    4061             :     else
    4062             :     {
    4063           7 :         psRoot = CPLParseXMLFile(m_osXMLFilename);
    4064             :     }
    4065         173 :     CPLXMLTreeCloser oCloser(psRoot);
    4066         173 :     psRoot = oCloser.get();
    4067         173 :     if (psRoot == nullptr)
    4068           0 :         return;
    4069         173 :     CPLXMLNode *psProduct = CPLGetXMLNode(psRoot, "=Product_Observational");
    4070         173 :     if (psProduct == nullptr)
    4071             :     {
    4072           1 :         psProduct = CPLGetXMLNode(psRoot, "=pds:Product_Observational");
    4073             :     }
    4074         173 :     if (psProduct == nullptr)
    4075             :     {
    4076           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    4077             :                  "Cannot find Product_Observational element in template");
    4078           1 :         return;
    4079             :     }
    4080             : 
    4081         172 :     if (m_bCreateHeader)
    4082             :     {
    4083         330 :         CPLString osCARTVersion(CURRENT_CART_VERSION);
    4084         165 :         char *pszXML = CPLSerializeXMLTree(psRoot);
    4085         165 :         if (pszXML)
    4086             :         {
    4087         165 :             const char *pszIter = pszXML;
    4088             :             while (true)
    4089             :             {
    4090         323 :                 const char *pszCartSchema = strstr(pszIter, "PDS4_CART_");
    4091         323 :                 if (pszCartSchema)
    4092             :                 {
    4093         316 :                     const char *pszXSDExtension = strstr(pszCartSchema, ".xsd");
    4094         316 :                     if (pszXSDExtension &&
    4095         316 :                         pszXSDExtension - pszCartSchema <= 20)
    4096             :                     {
    4097         158 :                         osCARTVersion = pszCartSchema + strlen("PDS4_CART_");
    4098         158 :                         osCARTVersion.resize(pszXSDExtension - pszCartSchema -
    4099             :                                              strlen("PDS4_CART_"));
    4100         158 :                         break;
    4101             :                     }
    4102             :                     else
    4103             :                     {
    4104         158 :                         pszIter = pszCartSchema + 1;
    4105             :                     }
    4106             :                 }
    4107             :                 else
    4108             :                 {
    4109           7 :                     break;
    4110             :                 }
    4111         158 :             }
    4112             : 
    4113         165 :             CPLFree(pszXML);
    4114             :         }
    4115             : 
    4116         165 :         CreateHeader(psProduct, osCARTVersion.c_str());
    4117             :     }
    4118             : 
    4119         172 :     WriteVectorLayers(psProduct);
    4120             : 
    4121         172 :     CPLSerializeXMLTreeToFile(psRoot, GetDescription());
    4122             : }
    4123             : 
    4124             : /************************************************************************/
    4125             : /*                            ICreateLayer()                            */
    4126             : /************************************************************************/
    4127             : 
    4128          66 : OGRLayer *PDS4Dataset::ICreateLayer(const char *pszName,
    4129             :                                     const OGRGeomFieldDefn *poGeomFieldDefn,
    4130             :                                     CSLConstList papszOptions)
    4131             : {
    4132             :     const char *pszTableType =
    4133          66 :         CSLFetchNameValueDef(papszOptions, "TABLE_TYPE", "DELIMITED");
    4134          66 :     if (!EQUAL(pszTableType, "CHARACTER") && !EQUAL(pszTableType, "BINARY") &&
    4135          55 :         !EQUAL(pszTableType, "DELIMITED"))
    4136             :     {
    4137           0 :         return nullptr;
    4138             :     }
    4139             : 
    4140          66 :     const auto eGType = poGeomFieldDefn ? poGeomFieldDefn->GetType() : wkbNone;
    4141             :     const auto poSpatialRef =
    4142          66 :         poGeomFieldDefn ? poGeomFieldDefn->GetSpatialRef() : nullptr;
    4143             : 
    4144         126 :     const char *pszExt = EQUAL(pszTableType, "CHARACTER") ? "dat"
    4145          60 :                          : EQUAL(pszTableType, "BINARY")  ? "bin"
    4146             :                                                           : "csv";
    4147             : 
    4148             :     bool bSameDirectory =
    4149          66 :         CPLTestBool(CSLFetchNameValueDef(papszOptions, "SAME_DIRECTORY", "NO"));
    4150             : 
    4151         132 :     std::string osBasename(pszName);
    4152         629 :     for (char &ch : osBasename)
    4153             :     {
    4154         563 :         if (!isalnum(static_cast<unsigned char>(ch)) &&
    4155          53 :             static_cast<unsigned>(ch) <= 127)
    4156          53 :             ch = '_';
    4157             :     }
    4158             : 
    4159         132 :     CPLString osFullFilename;
    4160          66 :     if (bSameDirectory)
    4161             :     {
    4162             :         osFullFilename =
    4163           2 :             CPLFormFilenameSafe(CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(),
    4164           1 :                                 osBasename.c_str(), pszExt);
    4165             :         VSIStatBufL sStat;
    4166           1 :         if (VSIStatL(osFullFilename, &sStat) == 0)
    4167             :         {
    4168           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    4169             :                      "%s already exists. Please delete it before, or "
    4170             :                      "rename the layer",
    4171             :                      osFullFilename.c_str());
    4172           0 :             return nullptr;
    4173             :         }
    4174             :     }
    4175             :     else
    4176             :     {
    4177         195 :         CPLString osDirectory = CPLFormFilenameSafe(
    4178         130 :             CPLGetPathSafe(m_osXMLFilename).c_str(),
    4179         130 :             CPLGetBasenameSafe(m_osXMLFilename).c_str(), nullptr);
    4180             :         VSIStatBufL sStat;
    4181         108 :         if (VSIStatL(osDirectory, &sStat) != 0 &&
    4182          43 :             VSIMkdir(osDirectory, 0755) != 0)
    4183             :         {
    4184           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot create directory %s",
    4185             :                      osDirectory.c_str());
    4186           1 :             return nullptr;
    4187             :         }
    4188             :         osFullFilename =
    4189          64 :             CPLFormFilenameSafe(osDirectory, osBasename.c_str(), pszExt);
    4190             :     }
    4191             : 
    4192          65 :     if (EQUAL(pszTableType, "DELIMITED"))
    4193             :     {
    4194             :         std::unique_ptr<PDS4DelimitedTable> poLayer(
    4195          54 :             new PDS4DelimitedTable(this, pszName, osFullFilename));
    4196          54 :         if (!poLayer->InitializeNewLayer(poSpatialRef, false, eGType,
    4197             :                                          papszOptions))
    4198             :         {
    4199           0 :             return nullptr;
    4200             :         }
    4201             :         std::unique_ptr<PDS4EditableLayer> poEditableLayer(
    4202         108 :             new PDS4EditableLayer(poLayer.release()));
    4203          54 :         m_apoLayers.push_back(std::move(poEditableLayer));
    4204             :     }
    4205             :     else
    4206             :     {
    4207             :         std::unique_ptr<PDS4FixedWidthTable> poLayer(
    4208          11 :             EQUAL(pszTableType, "CHARACTER")
    4209             :                 ? static_cast<PDS4FixedWidthTable *>(
    4210           6 :                       new PDS4TableCharacter(this, pszName, osFullFilename))
    4211             :                 : static_cast<PDS4FixedWidthTable *>(
    4212          17 :                       new PDS4TableBinary(this, pszName, osFullFilename)));
    4213          11 :         if (!poLayer->InitializeNewLayer(poSpatialRef, false, eGType,
    4214             :                                          papszOptions))
    4215             :         {
    4216           0 :             return nullptr;
    4217             :         }
    4218             :         std::unique_ptr<PDS4EditableLayer> poEditableLayer(
    4219          22 :             new PDS4EditableLayer(poLayer.release()));
    4220          11 :         m_apoLayers.push_back(std::move(poEditableLayer));
    4221             :     }
    4222          65 :     return m_apoLayers.back().get();
    4223             : }
    4224             : 
    4225             : /************************************************************************/
    4226             : /*                           TestCapability()                           */
    4227             : /************************************************************************/
    4228             : 
    4229          69 : int PDS4Dataset::TestCapability(const char *pszCap)
    4230             : {
    4231          69 :     if (EQUAL(pszCap, ODsCCreateLayer))
    4232          35 :         return eAccess == GA_Update;
    4233          34 :     else if (EQUAL(pszCap, ODsCZGeometries))
    4234           6 :         return TRUE;
    4235             :     else
    4236          28 :         return FALSE;
    4237             : }
    4238             : 
    4239             : /************************************************************************/
    4240             : /*                             Create()                                 */
    4241             : /************************************************************************/
    4242             : 
    4243         130 : GDALDataset *PDS4Dataset::Create(const char *pszFilename, int nXSize,
    4244             :                                  int nYSize, int nBandsIn, GDALDataType eType,
    4245             :                                  char **papszOptions)
    4246             : {
    4247         130 :     return CreateInternal(pszFilename, nullptr, nXSize, nYSize, nBandsIn, eType,
    4248         130 :                           papszOptions);
    4249             : }
    4250             : 
    4251             : /************************************************************************/
    4252             : /*                           CreateInternal()                           */
    4253             : /************************************************************************/
    4254             : 
    4255         189 : PDS4Dataset *PDS4Dataset::CreateInternal(const char *pszFilename,
    4256             :                                          GDALDataset *poSrcDS, int nXSize,
    4257             :                                          int nYSize, int nBandsIn,
    4258             :                                          GDALDataType eType,
    4259             :                                          const char *const *papszOptionsIn)
    4260             : {
    4261         378 :     CPLStringList aosOptions(papszOptionsIn);
    4262             : 
    4263         189 :     if (nXSize == 0 && nYSize == 0 && nBandsIn == 0 && eType == GDT_Unknown)
    4264             :     {
    4265             :         // Vector file creation
    4266          47 :         PDS4Dataset *poDS = new PDS4Dataset();
    4267          47 :         poDS->SetDescription(pszFilename);
    4268          47 :         poDS->nRasterXSize = 0;
    4269          47 :         poDS->nRasterYSize = 0;
    4270          47 :         poDS->eAccess = GA_Update;
    4271          47 :         poDS->m_osXMLFilename = pszFilename;
    4272          47 :         poDS->m_bCreateHeader = true;
    4273          47 :         poDS->m_bStripFileAreaObservationalFromTemplate = true;
    4274          47 :         poDS->m_papszCreationOptions = CSLDuplicate(aosOptions.List());
    4275          47 :         poDS->m_bUseSrcLabel = aosOptions.FetchBool("USE_SRC_LABEL", true);
    4276          47 :         return poDS;
    4277             :     }
    4278             : 
    4279         142 :     if (nXSize == 0)
    4280           0 :         return nullptr;
    4281             : 
    4282         142 :     if (!(eType == GDT_Byte || eType == GDT_Int8 || eType == GDT_Int16 ||
    4283          40 :           eType == GDT_UInt16 || eType == GDT_Int32 || eType == GDT_UInt32 ||
    4284          27 :           eType == GDT_Float32 || eType == GDT_Float64 ||
    4285          18 :           eType == GDT_CFloat32 || eType == GDT_CFloat64))
    4286             :     {
    4287          10 :         CPLError(
    4288             :             CE_Failure, CPLE_NotSupported,
    4289             :             "The PDS4 driver does not supporting creating files of type %s.",
    4290             :             GDALGetDataTypeName(eType));
    4291          10 :         return nullptr;
    4292             :     }
    4293             : 
    4294         132 :     if (nBandsIn == 0)
    4295             :     {
    4296           1 :         CPLError(CE_Failure, CPLE_NotSupported, "Invalid number of bands");
    4297           1 :         return nullptr;
    4298             :     }
    4299             : 
    4300             :     const char *pszArrayType =
    4301         131 :         aosOptions.FetchNameValueDef("ARRAY_TYPE", "Array_3D_Image");
    4302         131 :     const bool bIsArray2D = STARTS_WITH(pszArrayType, "Array_2D");
    4303         131 :     if (nBandsIn > 1 && bIsArray2D)
    4304             :     {
    4305           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    4306             :                  "ARRAY_TYPE=%s is not supported for a multi-band raster",
    4307             :                  pszArrayType);
    4308           1 :         return nullptr;
    4309             :     }
    4310             : 
    4311             :     /* -------------------------------------------------------------------- */
    4312             :     /*      Compute pixel, line and band offsets                            */
    4313             :     /* -------------------------------------------------------------------- */
    4314         130 :     const int nItemSize = GDALGetDataTypeSizeBytes(eType);
    4315             :     int nLineOffset, nPixelOffset;
    4316             :     vsi_l_offset nBandOffset;
    4317             : 
    4318             :     const char *pszInterleave =
    4319         130 :         aosOptions.FetchNameValueDef("INTERLEAVE", "BSQ");
    4320         130 :     if (bIsArray2D)
    4321           1 :         pszInterleave = "BIP";
    4322             : 
    4323         130 :     if (EQUAL(pszInterleave, "BIP"))
    4324             :     {
    4325           4 :         nPixelOffset = nItemSize * nBandsIn;
    4326           4 :         if (nPixelOffset > INT_MAX / nBandsIn)
    4327             :         {
    4328           0 :             return nullptr;
    4329             :         }
    4330           4 :         nLineOffset = nPixelOffset * nXSize;
    4331           4 :         nBandOffset = nItemSize;
    4332             :     }
    4333         126 :     else if (EQUAL(pszInterleave, "BSQ"))
    4334             :     {
    4335         122 :         nPixelOffset = nItemSize;
    4336         122 :         if (nPixelOffset > INT_MAX / nXSize)
    4337             :         {
    4338           0 :             return nullptr;
    4339             :         }
    4340         122 :         nLineOffset = nPixelOffset * nXSize;
    4341         122 :         nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nYSize;
    4342             :     }
    4343           4 :     else if (EQUAL(pszInterleave, "BIL"))
    4344             :     {
    4345           3 :         nPixelOffset = nItemSize;
    4346           3 :         if (nPixelOffset > INT_MAX / nBandsIn ||
    4347           3 :             nPixelOffset * nBandsIn > INT_MAX / nXSize)
    4348             :         {
    4349           0 :             return nullptr;
    4350             :         }
    4351           3 :         nLineOffset = nItemSize * nBandsIn * nXSize;
    4352           3 :         nBandOffset = static_cast<vsi_l_offset>(nItemSize) * nXSize;
    4353             :     }
    4354             :     else
    4355             :     {
    4356           1 :         CPLError(CE_Failure, CPLE_NotSupported, "Invalid value for INTERLEAVE");
    4357           1 :         return nullptr;
    4358             :     }
    4359             : 
    4360             :     const char *pszImageFormat =
    4361         129 :         aosOptions.FetchNameValueDef("IMAGE_FORMAT", "RAW");
    4362         129 :     const char *pszImageExtension = aosOptions.FetchNameValueDef(
    4363         129 :         "IMAGE_EXTENSION", EQUAL(pszImageFormat, "RAW") ? "img" : "tif");
    4364             :     CPLString osImageFilename(aosOptions.FetchNameValueDef(
    4365             :         "IMAGE_FILENAME",
    4366         258 :         CPLResetExtensionSafe(pszFilename, pszImageExtension).c_str()));
    4367             : 
    4368         129 :     const bool bAppend = aosOptions.FetchBool("APPEND_SUBDATASET", false);
    4369         129 :     if (bAppend)
    4370             :     {
    4371           4 :         GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
    4372           4 :         auto poExistingPDS4 = static_cast<PDS4Dataset *>(Open(&oOpenInfo));
    4373           4 :         if (!poExistingPDS4)
    4374             :         {
    4375           0 :             return nullptr;
    4376             :         }
    4377           4 :         osImageFilename = poExistingPDS4->m_osImageFilename;
    4378           4 :         delete poExistingPDS4;
    4379             : 
    4380           4 :         auto poImageDS = GDALDataset::FromHandle(GDALOpenEx(
    4381             :             osImageFilename, GDAL_OF_RASTER, nullptr, nullptr, nullptr));
    4382           6 :         if (poImageDS && poImageDS->GetDriver() &&
    4383           2 :             EQUAL(poImageDS->GetDriver()->GetDescription(), "GTiff"))
    4384             :         {
    4385           2 :             pszImageFormat = "GEOTIFF";
    4386             :         }
    4387           4 :         delete poImageDS;
    4388             :     }
    4389             : 
    4390         129 :     GDALDataset *poExternalDS = nullptr;
    4391         129 :     VSILFILE *fpImage = nullptr;
    4392         129 :     vsi_l_offset nBaseOffset = 0;
    4393         129 :     bool bIsLSB = true;
    4394         258 :     CPLString osHeaderParsingStandard;
    4395             :     const bool bCreateLabelOnly =
    4396         129 :         aosOptions.FetchBool("CREATE_LABEL_ONLY", false);
    4397         129 :     if (bCreateLabelOnly)
    4398             :     {
    4399           8 :         if (poSrcDS == nullptr)
    4400             :         {
    4401           0 :             CPLError(
    4402             :                 CE_Failure, CPLE_AppDefined,
    4403             :                 "CREATE_LABEL_ONLY is only compatible of CreateCopy() mode");
    4404           1 :             return nullptr;
    4405             :         }
    4406           8 :         RawBinaryLayout sLayout;
    4407           8 :         if (!poSrcDS->GetRawBinaryLayout(sLayout))
    4408             :         {
    4409           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    4410             :                      "Source dataset is not compatible of a raw binary format");
    4411           1 :             return nullptr;
    4412             :         }
    4413           7 :         if ((nBandsIn > 1 &&
    4414           7 :              sLayout.eInterleaving == RawBinaryLayout::Interleaving::UNKNOWN) ||
    4415           6 :             (nBandsIn == 1 &&
    4416           6 :              !(sLayout.nPixelOffset == nItemSize &&
    4417           6 :                sLayout.nLineOffset == sLayout.nPixelOffset * nXSize)))
    4418             :         {
    4419           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    4420             :                      "Source dataset has an interleaving not handled in PDS4");
    4421           0 :             return nullptr;
    4422             :         }
    4423           7 :         fpImage = VSIFOpenL(sLayout.osRawFilename.c_str(), "rb");
    4424           7 :         if (fpImage == nullptr)
    4425             :         {
    4426           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot open raw image %s",
    4427             :                      sLayout.osRawFilename.c_str());
    4428           0 :             return nullptr;
    4429             :         }
    4430           7 :         osImageFilename = sLayout.osRawFilename;
    4431           7 :         if (nBandsIn == 1 ||
    4432           1 :             sLayout.eInterleaving == RawBinaryLayout::Interleaving::BIP)
    4433           7 :             pszInterleave = "BIP";
    4434           0 :         else if (sLayout.eInterleaving == RawBinaryLayout::Interleaving::BIL)
    4435           0 :             pszInterleave = "BIL";
    4436             :         else
    4437           0 :             pszInterleave = "BSQ";
    4438           7 :         nBaseOffset = sLayout.nImageOffset;
    4439           7 :         nPixelOffset = static_cast<int>(sLayout.nPixelOffset);
    4440           7 :         nLineOffset = static_cast<int>(sLayout.nLineOffset);
    4441           7 :         nBandOffset = static_cast<vsi_l_offset>(sLayout.nBandOffset);
    4442           7 :         bIsLSB = sLayout.bLittleEndianOrder;
    4443           7 :         auto poSrcDriver = poSrcDS->GetDriver();
    4444           7 :         if (poSrcDriver)
    4445             :         {
    4446           7 :             auto pszDriverName = poSrcDriver->GetDescription();
    4447           7 :             if (EQUAL(pszDriverName, "GTiff"))
    4448             :             {
    4449           2 :                 GByte abySignature[4] = {0};
    4450           2 :                 VSIFReadL(abySignature, 1, 4, fpImage);
    4451           2 :                 const bool bBigTIFF =
    4452           2 :                     abySignature[2] == 43 || abySignature[3] == 43;
    4453             :                 osHeaderParsingStandard =
    4454           2 :                     bBigTIFF ? BIGTIFF_GEOTIFF_STRING : TIFF_GEOTIFF_STRING;
    4455             :             }
    4456           5 :             else if (EQUAL(pszDriverName, "ISIS3"))
    4457             :             {
    4458           1 :                 osHeaderParsingStandard = "ISIS3";
    4459             :             }
    4460           4 :             else if (EQUAL(pszDriverName, "VICAR"))
    4461             :             {
    4462           1 :                 osHeaderParsingStandard = "VICAR2";
    4463             :             }
    4464           3 :             else if (EQUAL(pszDriverName, "PDS"))
    4465             :             {
    4466           1 :                 osHeaderParsingStandard = "PDS3";
    4467             :             }
    4468           2 :             else if (EQUAL(pszDriverName, "FITS"))
    4469             :             {
    4470           1 :                 osHeaderParsingStandard = "FITS 3.0";
    4471             :                 aosOptions.SetNameValue("VAR_VERTICAL_DISPLAY_DIRECTION",
    4472           1 :                                         "Bottom to Top");
    4473             :             }
    4474             :         }
    4475             :     }
    4476         121 :     else if (EQUAL(pszImageFormat, "GEOTIFF"))
    4477             :     {
    4478          13 :         if (EQUAL(pszInterleave, "BIL"))
    4479             :         {
    4480           2 :             if (aosOptions.FetchBool("@INTERLEAVE_ADDED_AUTOMATICALLY", false))
    4481             :             {
    4482           1 :                 pszInterleave = "BSQ";
    4483             :             }
    4484             :             else
    4485             :             {
    4486           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
    4487             :                          "INTERLEAVE=BIL not supported for GeoTIFF in PDS4");
    4488           1 :                 return nullptr;
    4489             :             }
    4490             :         }
    4491             :         GDALDriver *poDrv =
    4492          12 :             static_cast<GDALDriver *>(GDALGetDriverByName("GTiff"));
    4493          12 :         if (poDrv == nullptr)
    4494             :         {
    4495           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot find GTiff driver");
    4496           0 :             return nullptr;
    4497             :         }
    4498          12 :         char **papszGTiffOptions = nullptr;
    4499             : #ifdef notdef
    4500             :         // In practice I can't see which option we can really use
    4501             :         const char *pszGTiffOptions =
    4502             :             CSLFetchNameValueDef(papszOptions, "GEOTIFF_OPTIONS", "");
    4503             :         char **papszTokens = CSLTokenizeString2(pszGTiffOptions, ",", 0);
    4504             :         if (CPLFetchBool(papszTokens, "TILED", false))
    4505             :         {
    4506             :             CSLDestroy(papszTokens);
    4507             :             CPLError(CE_Failure, CPLE_AppDefined,
    4508             :                      "Tiled GeoTIFF is not supported for PDS4");
    4509             :             return NULL;
    4510             :         }
    4511             :         if (!EQUAL(CSLFetchNameValueDef(papszTokens, "COMPRESS", "NONE"),
    4512             :                    "NONE"))
    4513             :         {
    4514             :             CSLDestroy(papszTokens);
    4515             :             CPLError(CE_Failure, CPLE_AppDefined,
    4516             :                      "Compressed GeoTIFF is not supported for PDS4");
    4517             :             return NULL;
    4518             :         }
    4519             :         papszGTiffOptions =
    4520             :             CSLSetNameValue(papszGTiffOptions, "ENDIANNESS", "LITTLE");
    4521             :         for (int i = 0; papszTokens[i] != NULL; i++)
    4522             :         {
    4523             :             papszGTiffOptions = CSLAddString(papszGTiffOptions, papszTokens[i]);
    4524             :         }
    4525             :         CSLDestroy(papszTokens);
    4526             : #endif
    4527             : 
    4528             :         papszGTiffOptions =
    4529          12 :             CSLSetNameValue(papszGTiffOptions, "INTERLEAVE",
    4530          12 :                             EQUAL(pszInterleave, "BSQ") ? "BAND" : "PIXEL");
    4531             :         // Will make sure that our blocks at nodata are not optimized
    4532             :         // away but indeed well written
    4533          12 :         papszGTiffOptions = CSLSetNameValue(
    4534             :             papszGTiffOptions, "@WRITE_EMPTY_TILES_SYNCHRONOUSLY", "YES");
    4535          12 :         if (nBandsIn > 1 && EQUAL(pszInterleave, "BSQ"))
    4536             :         {
    4537             :             papszGTiffOptions =
    4538           2 :                 CSLSetNameValue(papszGTiffOptions, "BLOCKYSIZE", "1");
    4539             :         }
    4540             : 
    4541          12 :         if (bAppend)
    4542             :         {
    4543             :             papszGTiffOptions =
    4544           2 :                 CSLAddString(papszGTiffOptions, "APPEND_SUBDATASET=YES");
    4545             :         }
    4546             : 
    4547          12 :         poExternalDS = poDrv->Create(osImageFilename, nXSize, nYSize, nBandsIn,
    4548             :                                      eType, papszGTiffOptions);
    4549          12 :         CSLDestroy(papszGTiffOptions);
    4550          12 :         if (poExternalDS == nullptr)
    4551             :         {
    4552           1 :             CPLError(CE_Failure, CPLE_FileIO, "Cannot create %s",
    4553             :                      osImageFilename.c_str());
    4554           1 :             return nullptr;
    4555             :         }
    4556             :     }
    4557             :     else
    4558             :     {
    4559         214 :         fpImage = VSIFOpenL(
    4560             :             osImageFilename,
    4561             :             bAppend                                                 ? "rb+"
    4562         106 :             : VSISupportsRandomWrite(osImageFilename.c_str(), true) ? "wb+"
    4563             :                                                                     : "wb");
    4564         108 :         if (fpImage == nullptr)
    4565             :         {
    4566           3 :             CPLError(CE_Failure, CPLE_FileIO, "Cannot create %s",
    4567             :                      osImageFilename.c_str());
    4568           3 :             return nullptr;
    4569             :         }
    4570         105 :         if (bAppend)
    4571             :         {
    4572           2 :             VSIFSeekL(fpImage, 0, SEEK_END);
    4573           2 :             nBaseOffset = VSIFTellL(fpImage);
    4574             :         }
    4575             :     }
    4576             : 
    4577         246 :     auto poDS = std::make_unique<PDS4Dataset>();
    4578         123 :     poDS->SetDescription(pszFilename);
    4579         123 :     poDS->m_bMustInitImageFile = true;
    4580         123 :     poDS->m_fpImage = fpImage;
    4581         123 :     poDS->m_nBaseOffset = nBaseOffset;
    4582         123 :     poDS->m_poExternalDS = poExternalDS;
    4583         123 :     poDS->nRasterXSize = nXSize;
    4584         123 :     poDS->nRasterYSize = nYSize;
    4585         123 :     poDS->eAccess = GA_Update;
    4586         123 :     poDS->m_osImageFilename = std::move(osImageFilename);
    4587         123 :     poDS->m_bCreateHeader = true;
    4588         123 :     poDS->m_bStripFileAreaObservationalFromTemplate = true;
    4589         123 :     poDS->m_osInterleave = pszInterleave;
    4590         123 :     poDS->m_papszCreationOptions = CSLDuplicate(aosOptions.List());
    4591         123 :     poDS->m_bUseSrcLabel = aosOptions.FetchBool("USE_SRC_LABEL", true);
    4592         123 :     poDS->m_bIsLSB = bIsLSB;
    4593         123 :     poDS->m_osHeaderParsingStandard = std::move(osHeaderParsingStandard);
    4594         123 :     poDS->m_bCreatedFromExistingBinaryFile = bCreateLabelOnly;
    4595             : 
    4596         123 :     if (EQUAL(pszInterleave, "BIP"))
    4597             :     {
    4598          10 :         poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
    4599             :                                            "IMAGE_STRUCTURE");
    4600             :     }
    4601         113 :     else if (EQUAL(pszInterleave, "BSQ"))
    4602             :     {
    4603         112 :         poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "BAND",
    4604             :                                            "IMAGE_STRUCTURE");
    4605             :     }
    4606             : 
    4607         302 :     for (int i = 0; i < nBandsIn; i++)
    4608             :     {
    4609         179 :         if (poDS->m_poExternalDS != nullptr)
    4610             :         {
    4611             :             auto poBand = std::make_unique<PDS4WrapperRasterBand>(
    4612          17 :                 poDS->m_poExternalDS->GetRasterBand(i + 1));
    4613          17 :             poDS->SetBand(i + 1, std::move(poBand));
    4614             :         }
    4615             :         else
    4616             :         {
    4617             :             auto poBand = std::make_unique<PDS4RawRasterBand>(
    4618         162 :                 poDS.get(), i + 1, poDS->m_fpImage,
    4619         162 :                 poDS->m_nBaseOffset + nBandOffset * i, nPixelOffset,
    4620             :                 nLineOffset, eType,
    4621         162 :                 bIsLSB ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
    4622         162 :                        : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN);
    4623         162 :             poDS->SetBand(i + 1, std::move(poBand));
    4624             :         }
    4625             :     }
    4626             : 
    4627         123 :     return poDS.release();
    4628             : }
    4629             : 
    4630             : /************************************************************************/
    4631             : /*                      PDS4GetUnderlyingDataset()                      */
    4632             : /************************************************************************/
    4633             : 
    4634          63 : static GDALDataset *PDS4GetUnderlyingDataset(GDALDataset *poSrcDS)
    4635             : {
    4636         126 :     if (poSrcDS->GetDriver() != nullptr &&
    4637          63 :         poSrcDS->GetDriver() == GDALGetDriverByName("VRT"))
    4638             :     {
    4639           4 :         VRTDataset *poVRTDS = reinterpret_cast<VRTDataset *>(poSrcDS);
    4640           4 :         poSrcDS = poVRTDS->GetSingleSimpleSource();
    4641             :     }
    4642             : 
    4643          63 :     return poSrcDS;
    4644             : }
    4645             : 
    4646             : /************************************************************************/
    4647             : /*                            CreateCopy()                              */
    4648             : /************************************************************************/
    4649             : 
    4650          63 : GDALDataset *PDS4Dataset::CreateCopy(const char *pszFilename,
    4651             :                                      GDALDataset *poSrcDS, int bStrict,
    4652             :                                      char **papszOptions,
    4653             :                                      GDALProgressFunc pfnProgress,
    4654             :                                      void *pProgressData)
    4655             : {
    4656             :     const char *pszImageFormat =
    4657          63 :         CSLFetchNameValueDef(papszOptions, "IMAGE_FORMAT", "RAW");
    4658          63 :     GDALDataset *poSrcUnderlyingDS = PDS4GetUnderlyingDataset(poSrcDS);
    4659          63 :     if (poSrcUnderlyingDS == nullptr)
    4660           1 :         poSrcUnderlyingDS = poSrcDS;
    4661          69 :     if (EQUAL(pszImageFormat, "GEOTIFF") &&
    4662           6 :         strcmp(poSrcUnderlyingDS->GetDescription(),
    4663             :                CSLFetchNameValueDef(
    4664             :                    papszOptions, "IMAGE_FILENAME",
    4665          69 :                    CPLResetExtensionSafe(pszFilename, "tif").c_str())) == 0)
    4666             :     {
    4667           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    4668             :                  "Output file has same name as input file");
    4669           1 :         return nullptr;
    4670             :     }
    4671          62 :     if (poSrcDS->GetRasterCount() == 0)
    4672             :     {
    4673           1 :         CPLError(CE_Failure, CPLE_NotSupported, "Unsupported band count");
    4674           1 :         return nullptr;
    4675             :     }
    4676             : 
    4677          61 :     const bool bAppend = CPLFetchBool(papszOptions, "APPEND_SUBDATASET", false);
    4678          61 :     if (bAppend)
    4679             :     {
    4680           6 :         GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
    4681           6 :         GDALDataset *poExistingDS = Open(&oOpenInfo);
    4682           6 :         if (poExistingDS)
    4683             :         {
    4684           6 :             GDALGeoTransform existingGT;
    4685             :             const bool bExistingHasGT =
    4686           6 :                 poExistingDS->GetGeoTransform(existingGT) == CE_None;
    4687           6 :             GDALGeoTransform gt;
    4688           6 :             const bool bSrcHasGT = poSrcDS->GetGeoTransform(gt) == CE_None;
    4689             : 
    4690           6 :             OGRSpatialReference oExistingSRS;
    4691           6 :             OGRSpatialReference oSrcSRS;
    4692           6 :             const char *pszExistingSRS = poExistingDS->GetProjectionRef();
    4693           6 :             const char *pszSrcSRS = poSrcDS->GetProjectionRef();
    4694           6 :             CPLString osExistingProj4;
    4695           6 :             if (pszExistingSRS && pszExistingSRS[0])
    4696             :             {
    4697           6 :                 oExistingSRS.SetFromUserInput(
    4698             :                     pszExistingSRS,
    4699             :                     OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
    4700           6 :                 char *pszExistingProj4 = nullptr;
    4701           6 :                 oExistingSRS.exportToProj4(&pszExistingProj4);
    4702           6 :                 if (pszExistingProj4)
    4703           6 :                     osExistingProj4 = pszExistingProj4;
    4704           6 :                 CPLFree(pszExistingProj4);
    4705             :             }
    4706           6 :             CPLString osSrcProj4;
    4707           6 :             if (pszSrcSRS && pszSrcSRS[0])
    4708             :             {
    4709           6 :                 oSrcSRS.SetFromUserInput(
    4710             :                     pszSrcSRS,
    4711             :                     OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
    4712           6 :                 char *pszSrcProj4 = nullptr;
    4713           6 :                 oSrcSRS.exportToProj4(&pszSrcProj4);
    4714           6 :                 if (pszSrcProj4)
    4715           6 :                     osSrcProj4 = pszSrcProj4;
    4716           6 :                 CPLFree(pszSrcProj4);
    4717             :             }
    4718             : 
    4719           6 :             delete poExistingDS;
    4720             : 
    4721             :             const auto maxRelErrorGT =
    4722           6 :                 [](const GDALGeoTransform &gt1, const GDALGeoTransform &gt2)
    4723             :             {
    4724           6 :                 double maxRelError = 0.0;
    4725          42 :                 for (int i = 0; i < 6; i++)
    4726             :                 {
    4727          36 :                     if (gt1[i] == 0.0)
    4728             :                     {
    4729          12 :                         maxRelError = std::max(maxRelError, std::abs(gt2[i]));
    4730             :                     }
    4731             :                     else
    4732             :                     {
    4733          24 :                         maxRelError =
    4734          48 :                             std::max(maxRelError, std::abs(gt2[i] - gt1[i]) /
    4735          24 :                                                       std::abs(gt1[i]));
    4736             :                     }
    4737             :                 }
    4738           6 :                 return maxRelError;
    4739             :             };
    4740             : 
    4741           6 :             if ((bExistingHasGT && !bSrcHasGT) ||
    4742          18 :                 (!bExistingHasGT && bSrcHasGT) ||
    4743           6 :                 (bExistingHasGT && bSrcHasGT &&
    4744           6 :                  maxRelErrorGT(existingGT, gt) > 1e-10))
    4745             :             {
    4746           1 :                 CPLError(bStrict ? CE_Failure : CE_Warning, CPLE_NotSupported,
    4747             :                          "Appending to a dataset with a different "
    4748             :                          "geotransform is not supported");
    4749           1 :                 if (bStrict)
    4750           1 :                     return nullptr;
    4751             :             }
    4752             :             // Do proj string comparison, as it is unlikely that
    4753             :             // OGRSpatialReference::IsSame() will lead to identical reasons due
    4754             :             // to PDS changing CRS names, etc...
    4755           5 :             if (osExistingProj4 != osSrcProj4)
    4756             :             {
    4757           1 :                 CPLError(bStrict ? CE_Failure : CE_Warning, CPLE_NotSupported,
    4758             :                          "Appending to a dataset with a different "
    4759             :                          "coordinate reference system is not supported");
    4760           1 :                 if (bStrict)
    4761           1 :                     return nullptr;
    4762             :             }
    4763             :         }
    4764             :     }
    4765             : 
    4766          59 :     const int nXSize = poSrcDS->GetRasterXSize();
    4767          59 :     const int nYSize = poSrcDS->GetRasterYSize();
    4768          59 :     const int nBands = poSrcDS->GetRasterCount();
    4769          59 :     GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
    4770          59 :     PDS4Dataset *poDS = CreateInternal(pszFilename, poSrcDS, nXSize, nYSize,
    4771             :                                        nBands, eType, papszOptions);
    4772          59 :     if (poDS == nullptr)
    4773           6 :         return nullptr;
    4774             : 
    4775          53 :     GDALGeoTransform gt;
    4776          53 :     if (poSrcDS->GetGeoTransform(gt) == CE_None && gt != GDALGeoTransform())
    4777             :     {
    4778          50 :         poDS->SetGeoTransform(gt);
    4779             :     }
    4780             : 
    4781         106 :     if (poSrcDS->GetProjectionRef() != nullptr &&
    4782          53 :         strlen(poSrcDS->GetProjectionRef()) > 0)
    4783             :     {
    4784          50 :         poDS->SetProjection(poSrcDS->GetProjectionRef());
    4785             :     }
    4786             : 
    4787         133 :     for (int i = 1; i <= nBands; i++)
    4788             :     {
    4789          80 :         int bHasNoData = false;
    4790             :         const double dfNoData =
    4791          80 :             poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasNoData);
    4792          80 :         if (bHasNoData)
    4793          12 :             poDS->GetRasterBand(i)->SetNoDataValue(dfNoData);
    4794             : 
    4795          80 :         const double dfOffset = poSrcDS->GetRasterBand(i)->GetOffset();
    4796          80 :         if (dfOffset != 0.0)
    4797           3 :             poDS->GetRasterBand(i)->SetOffset(dfOffset);
    4798             : 
    4799          80 :         const double dfScale = poSrcDS->GetRasterBand(i)->GetScale();
    4800          80 :         if (dfScale != 1.0)
    4801           3 :             poDS->GetRasterBand(i)->SetScale(dfScale);
    4802             : 
    4803         160 :         poDS->GetRasterBand(i)->SetUnitType(
    4804          80 :             poSrcDS->GetRasterBand(i)->GetUnitType());
    4805             :     }
    4806             : 
    4807          53 :     if (poDS->m_bUseSrcLabel)
    4808             :     {
    4809          51 :         char **papszMD_PDS4 = poSrcDS->GetMetadata("xml:PDS4");
    4810          51 :         if (papszMD_PDS4 != nullptr)
    4811             :         {
    4812           6 :             poDS->SetMetadata(papszMD_PDS4, "xml:PDS4");
    4813             :         }
    4814             :     }
    4815             : 
    4816          53 :     if (poDS->m_poExternalDS == nullptr)
    4817             :     {
    4818             :         // We don't need to initialize the imagery as we are going to copy it
    4819             :         // completely
    4820          46 :         poDS->m_bMustInitImageFile = false;
    4821             :     }
    4822             : 
    4823          53 :     if (!CPLFetchBool(papszOptions, "CREATE_LABEL_ONLY", false))
    4824             :     {
    4825          46 :         CPLErr eErr = GDALDatasetCopyWholeRaster(poSrcDS, poDS, nullptr,
    4826             :                                                  pfnProgress, pProgressData);
    4827          46 :         poDS->FlushCache(false);
    4828          46 :         if (eErr != CE_None)
    4829             :         {
    4830           1 :             delete poDS;
    4831           1 :             return nullptr;
    4832             :         }
    4833             : 
    4834          45 :         char **papszISIS3MD = poSrcDS->GetMetadata("json:ISIS3");
    4835          45 :         if (papszISIS3MD)
    4836             :         {
    4837           1 :             poDS->SetMetadata(papszISIS3MD, "json:ISIS3");
    4838             :         }
    4839             :     }
    4840             : 
    4841          52 :     return poDS;
    4842             : }
    4843             : 
    4844             : /************************************************************************/
    4845             : /*                             Delete()                                 */
    4846             : /************************************************************************/
    4847             : 
    4848         108 : CPLErr PDS4Dataset::Delete(const char *pszFilename)
    4849             : 
    4850             : {
    4851             :     /* -------------------------------------------------------------------- */
    4852             :     /*      Collect file list.                                              */
    4853             :     /* -------------------------------------------------------------------- */
    4854         216 :     GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
    4855             :     auto poDS =
    4856         216 :         std::unique_ptr<PDS4Dataset>(PDS4Dataset::OpenInternal(&oOpenInfo));
    4857         108 :     if (poDS == nullptr)
    4858             :     {
    4859           0 :         if (CPLGetLastErrorNo() == 0)
    4860           0 :             CPLError(CE_Failure, CPLE_OpenFailed,
    4861             :                      "Unable to open %s to obtain file list.", pszFilename);
    4862             : 
    4863           0 :         return CE_Failure;
    4864             :     }
    4865             : 
    4866         108 :     char **papszFileList = poDS->GetFileList();
    4867         216 :     CPLString osImageFilename = poDS->m_osImageFilename;
    4868             :     bool bCreatedFromExistingBinaryFile =
    4869         108 :         poDS->m_bCreatedFromExistingBinaryFile;
    4870             : 
    4871         108 :     poDS.reset();
    4872             : 
    4873         108 :     if (CSLCount(papszFileList) == 0)
    4874             :     {
    4875           0 :         CPLError(CE_Failure, CPLE_NotSupported,
    4876             :                  "Unable to determine files associated with %s, "
    4877             :                  "delete fails.",
    4878             :                  pszFilename);
    4879           0 :         CSLDestroy(papszFileList);
    4880           0 :         return CE_Failure;
    4881             :     }
    4882             : 
    4883             :     /* -------------------------------------------------------------------- */
    4884             :     /*      Delete all files.                                               */
    4885             :     /* -------------------------------------------------------------------- */
    4886         108 :     CPLErr eErr = CE_None;
    4887         388 :     for (int i = 0; papszFileList[i] != nullptr; ++i)
    4888             :     {
    4889         294 :         if (bCreatedFromExistingBinaryFile &&
    4890          14 :             EQUAL(papszFileList[i], osImageFilename))
    4891             :         {
    4892           7 :             continue;
    4893             :         }
    4894         273 :         if (VSIUnlink(papszFileList[i]) != 0)
    4895             :         {
    4896           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Deleting %s failed:\n%s",
    4897           0 :                      papszFileList[i], VSIStrerror(errno));
    4898           0 :             eErr = CE_Failure;
    4899             :         }
    4900             :     }
    4901             : 
    4902         108 :     CSLDestroy(papszFileList);
    4903             : 
    4904         108 :     return eErr;
    4905             : }
    4906             : 
    4907             : /************************************************************************/
    4908             : /*                         GDALRegister_PDS4()                          */
    4909             : /************************************************************************/
    4910             : 
    4911        1935 : void GDALRegister_PDS4()
    4912             : 
    4913             : {
    4914        1935 :     if (GDALGetDriverByName(PDS4_DRIVER_NAME) != nullptr)
    4915         282 :         return;
    4916             : 
    4917        1653 :     GDALDriver *poDriver = new GDALDriver();
    4918        1653 :     PDS4DriverSetCommonMetadata(poDriver);
    4919             : 
    4920        1653 :     poDriver->pfnOpen = PDS4Dataset::Open;
    4921        1653 :     poDriver->pfnCreate = PDS4Dataset::Create;
    4922        1653 :     poDriver->pfnCreateCopy = PDS4Dataset::CreateCopy;
    4923        1653 :     poDriver->pfnDelete = PDS4Dataset::Delete;
    4924             : 
    4925        1653 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    4926             : }

Generated by: LCOV version 1.14