text: /******************************************************************************
   *
   * Project:  GDAL TileDB Driver
   * Purpose:  Implement GDAL TileDB Support based on https://www.tiledb.io
   * Author:   TileDB, Inc
   *
   ******************************************************************************
   * Copyright (c) 2019, TileDB, Inc
   *
   * Permission is hereby granted, free of charge, to any person obtaining a
   * copy of this software and associated documentation files (the "Software"),
   * to deal in the Software without restriction, including without limitation
   * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   * and/or sell copies of the Software, and to permit persons to whom the
   * Software is furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included
   * in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
   * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
   * DEALINGS IN THE SOFTWARE.
   ****************************************************************************/
  
  #include <cassert>
  #include <list>
  
  #include "cpl_string.h"
  #include "gdal_frmts.h"
  #include "gdal_pam.h"
  
  #include "tiledb_headers.h"
  
  CPL_CVSID("$Id$")
  
  #ifdef _MSC_VER
  #pragma warning( push )
  // 'tiledb::Query::set_buffer': was declared deprecated
  // 'tiledb::Array::Array': was declared deprecated
  #pragma warning( disable : 4996 ) /* XXXX was deprecated */
  #endif
  
  const CPLString TILEDB_VALUES( "TDB_VALUES" );
  
  /************************************************************************/
  /* ==================================================================== */
  /*                               TileDBDataset                          */
  /* ==================================================================== */
  /************************************************************************/
  
  class TileDBRasterBand;
  
  class TileDBDataset final: public GDALPamDataset
  {
      friend class TileDBRasterBand;
  
      protected:
          uint64_t      nTimestamp = 0;
          int           nBitsPerSample = 8;
          GDALDataType  eDataType = GDT_Unknown;
          int           nBlockXSize = -1;
          int           nBlockYSize = -1;
          int           nBlocksX = 0;
          int           nBlocksY = 0;
          int           nBandStart = 1;
          bool          bGlobalOrder = false;
          bool          bHasSubDatasets = false;
          int           nSubDataCount = 0;
          char          **papszSubDatasets = nullptr;
          CPLStringList m_osSubdatasetMD{};
          CPLXMLNode*   psSubDatasetsTree = nullptr;
          CPLString     osMetaDoc;
          TILEDB_INTERLEAVE_MODE eIndexMode = BAND;
  
          std::unique_ptr<tiledb::Context> m_ctx;
          std::unique_ptr<tiledb::Context> m_roCtx;
          std::unique_ptr<tiledb::Array> m_array;
          std::unique_ptr<tiledb::Array> m_roArray;
          std::unique_ptr<tiledb::ArraySchema> m_schema;
          std::unique_ptr<tiledb::FilterList> m_filterList;
  
          char          **papszAttributes = nullptr;
          std::list<std::unique_ptr<GDALDataset>> lpoAttributeDS = {};
  
          bool bStats = FALSE;
          CPLErr AddFilter( const char* pszFilterName, const int level );
          CPLErr CreateAttribute( GDALDataType eType, const CPLString& osAttrName,
                                  const int nSubRasterCount=1 );
          CPLErr AddDimensions( tiledb::Domain& domain,
                               tiledb::Dimension& y,
                               tiledb::Dimension& x,
                               tiledb::Dimension* poBands=nullptr );
          CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
                          void *, int, int, GDALDataType,
                          int, int *,
                          GSpacing nPixelSpace, GSpacing nLineSpace,
                          GSpacing nBandSpace,
                          GDALRasterIOExtraArg* psExtraArg ) override;
      public:
          virtual ~TileDBDataset();
  
          CPLErr TryLoadCachedXML(char **papszSiblingFiles = nullptr, bool bReload=true);
          CPLErr TryLoadXML(char **papszSiblingFiles = nullptr) override;
          CPLErr TrySaveXML() override;
          char** GetMetadata(const char *pszDomain) override;
          static GDALDataset      *Open( GDALOpenInfo * );
          static int              Identify( GDALOpenInfo * );
          static CPLErr           Delete( const char * pszFilename );
          static CPLErr           CopySubDatasets( GDALDataset* poSrcDS,
                                                  TileDBDataset* poDstDS,
                                                  GDALProgressFunc pfnProgress,
                                                  void *pProgressData );
          static TileDBDataset    *CreateLL( const char * pszFilename,
                                      int nXSize, int nYSize, int nBands,
                                      GDALDataType eType, char ** papszOptions );
          static GDALDataset      *Create( const char * pszFilename,
                                      int nXSize, int nYSize, int nBands,
                                      GDALDataType eType, char ** papszOptions );
          static GDALDataset      *CreateCopy( const char * pszFilename,
                                      GDALDataset *poSrcDS,
                                      int bStrict, char ** papszOptions,
                                      GDALProgressFunc pfnProgress,
                                      void *pProgressData);
          static void             ErrorHandler( const std::string& msg );
          static void             SetBlockSize( GDALRasterBand* poBand,
                                                  char ** &papszOptions );
  
          virtual void            FlushCache( bool bAtClosing ) override;
  
  
  };
  
  /************************************************************************/
  /* ==================================================================== */
  /*                            TileDBRasterBand                          */
  /* ==================================================================== */
  /************************************************************************/
  
  class TileDBRasterBand final: public GDALPamRasterBand
  {
      friend class TileDBDataset;
      protected:
          TileDBDataset  *poGDS;
          bool bStats;
          CPLString osAttrName;
          int nCurrBlockX = 0;
          int nCurrBlockY = 0;
          std::unique_ptr<tiledb::Query> m_query;
          std::unique_ptr<tiledb::Query> m_roQuery;
          void   Finalize( );
      public:
          TileDBRasterBand( TileDBDataset *, int, CPLString = TILEDB_VALUES );
          virtual CPLErr IReadBlock( int, int, void * ) override;
          virtual CPLErr IWriteBlock( int, int, void * ) override;
          virtual GDALColorInterp GetColorInterpretation() override;
  
  };
  
  static CPLString vsi_to_tiledb_uri( const char* pszUri )
  {
      CPLString osUri;
  
      if ( STARTS_WITH_CI( pszUri, "/VSIS3/") )
          osUri.Printf("s3://%s", pszUri + 7);
      else if ( STARTS_WITH_CI( pszUri, "/VSIGS/") )
          osUri.Printf("gcs://%s", pszUri + 7);
      else
      {
          osUri = pszUri;
          // tiledb (at least at 2.4.2 on Conda) wrongly interprets relative
          // directories on Windows as absolute ones.
          if( CPLIsFilenameRelative(pszUri) )
          {
              char* pszCurDir = CPLGetCurrentDir();
              if( pszCurDir )
                  osUri = CPLFormFilename(pszCurDir, pszUri, nullptr);
              CPLFree(pszCurDir);
          }
      }
  
      return osUri;
  }
  
  static CPLErr option_to_index_type( const char* pszIndexingType, TILEDB_INTERLEAVE_MODE& eMode )
  {
      CPLErr eErr = CE_None;
  
      if ( pszIndexingType )
      {
          if EQUAL( pszIndexingType, "BAND" )
              eMode = BAND;
          else if EQUAL( pszIndexingType, "ATTRIBUTES" )
              eMode = ATTRIBUTES;
          else if EQUAL( pszIndexingType, "PIXEL" )
              eMode = PIXEL;
          else
          {
              eErr = CE_Failure;
              CPLError( eErr, CPLE_AppDefined,
                  "Unable to identify TileDB index mode %s.",
                  pszIndexingType );
          }
      }
      else
      {
          eMode = BAND;
      }
  
      return eErr;
  }
  
  static const char* index_type_name( TILEDB_INTERLEAVE_MODE eMode )
  {
      switch ( eMode )
      {
          case PIXEL:
              return "PIXEL";
          case ATTRIBUTES:
              return "ATTRIBUTES";
          case BAND:
              return "BAND";
          default:
              return nullptr;
      }
  }
  
  /************************************************************************/
  /*                             SetBuffer()                              */
  /************************************************************************/
  
  static CPLErr SetBuffer( tiledb::Query* poQuery, GDALDataType eType,
                          const CPLString& osAttrName, void * pImage, int nSize )
  {
      switch (eType)
      {
          case GDT_Byte:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<unsigned char*>( pImage ), nSize );
              break;
          case GDT_UInt16:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<unsigned short*>( pImage ), nSize );
              break;
          case GDT_UInt32:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<unsigned int*>( pImage ), nSize );
              break;
          case GDT_Int16:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<short*>( pImage ), nSize );
              break;
          case GDT_Int32:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<int*>( pImage ), nSize );
              break;
          case GDT_Float32:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<float*>( pImage ), nSize );
              break;
          case GDT_Float64:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<double*>( pImage ), nSize );
              break;
          case GDT_CInt16:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<short*>( pImage ), nSize * 2 );
              break;
          case GDT_CInt32:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<int*>( pImage ), nSize * 2 );
              break;
          case GDT_CFloat32:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<float*>( pImage ), nSize * 2 );
              break;
          case GDT_CFloat64:
              poQuery->set_buffer(
                  osAttrName, reinterpret_cast<double*>( pImage ), nSize * 2 );
              break;
          default:
              return CE_Failure;
      }
      return CE_None;
  }
  
  /************************************************************************/
  /*                          TileDBRasterBand()                          */
  /************************************************************************/
  
  TileDBRasterBand::TileDBRasterBand(
          TileDBDataset *poDSIn, int nBandIn, CPLString osAttr ) :
      poGDS( poDSIn ),
      bStats( poDSIn->bStats ),
      osAttrName(osAttr)
  {
      poDS = poDSIn;
      nBand = nBandIn;
      eDataType = poGDS->eDataType;
      eAccess = poGDS->eAccess;
      nRasterXSize = poGDS->nRasterXSize;
      nRasterYSize = poGDS->nRasterYSize;
      nBlockXSize = poGDS->nBlockXSize;
      nBlockYSize = poGDS->nBlockYSize;
  
      m_query.reset(new tiledb::Query( *poGDS->m_ctx, *poGDS->m_array ) );
  
      if ( ( eAccess == GA_Update ) && ( poGDS->m_roArray ) )
      {
          m_roQuery.reset( new tiledb::Query( *poGDS->m_roCtx, *poGDS->m_roArray ) );
          m_roQuery->set_layout( TILEDB_ROW_MAJOR );
      }
  
      if ( poGDS->bGlobalOrder )
      {
          m_query->set_layout( TILEDB_GLOBAL_ORDER );
          int nBandIdx = poGDS->nBandStart + nBand - 1;
          // initialize to complete image block layout
          std::vector<uint64_t> oaSubarray = { uint64_t( nBandIdx ),
                                          uint64_t( nBandIdx ),
                                          0,
                                          uint64_t(poDSIn->nBlocksY) * nBlockYSize - 1,
                                          0,
                                          uint64_t(poDSIn->nBlocksX) * nBlockXSize - 1,
                                          };
  
          if ( poGDS->m_array->schema().domain().ndim() == 3 )
          {
              if ( poGDS->eIndexMode == PIXEL )
                  std::rotate( oaSubarray.begin(),
                              oaSubarray.begin() + 2,
                              oaSubarray.end() );
              m_query->set_subarray( oaSubarray );
          }
          else
          {
              m_query->set_subarray( std::vector<uint64_t> (
                                          oaSubarray.cbegin() + 2,
                                          oaSubarray.cend() ) );
          }
      }
      else
          m_query->set_layout( TILEDB_ROW_MAJOR );
  }
  
  /************************************************************************/
  /*                          Finalize()                                  */
  /************************************************************************/
  
  void TileDBRasterBand::Finalize()
  
  {
      if ( poGDS->bGlobalOrder )
      {
          m_query->finalize();
      }
  }
  
  /************************************************************************/
  /*                             IReadBlock()                             */
  /************************************************************************/
  
  CPLErr TileDBRasterBand::IReadBlock( int nBlockXOff,
                                      int nBlockYOff,
                                      void * pImage )
  {
      int nStartX = nBlockXSize * nBlockXOff;
      int nStartY = nBlockYSize * nBlockYOff;
      uint64_t nEndX =  nStartX + nBlockXSize;
      uint64_t nEndY =  nStartY + nBlockYSize;
      int nBandIdx = poGDS->nBandStart + nBand - 1;
  
      std::vector<uint64_t> oaSubarray = {
                                      uint64_t( nBandIdx ),
                                      uint64_t( nBandIdx ),
                                      (uint64_t) nStartY,
                                      (uint64_t) nEndY - 1,
                                      (uint64_t) nStartX,
                                      (uint64_t) nEndX - 1 };
  
      if ( poGDS->eIndexMode == PIXEL )
          std::rotate( oaSubarray.begin(),
                      oaSubarray.begin() + 2,
                      oaSubarray.end() );
  
      tiledb::Query* q;
      if ( ( eAccess == GA_Update ) && ( poGDS->m_roArray ) )
          q = m_roQuery.get();
      else
          q = m_query.get();
  
      if ( poGDS->m_array->schema().domain().ndim() == 3 )
      {
          q->set_subarray( oaSubarray );
      }
      else
      {
          q->set_subarray( std::vector<uint64_t> (
                                      oaSubarray.cbegin() + 2,
                                      oaSubarray.cend() ) );
      }
  
      SetBuffer(q, eDataType, osAttrName,
                       pImage, nBlockXSize * nBlockYSize );
  
      if ( bStats )
          tiledb::Stats::enable();
  
      auto status = q->submit();
  
      if ( bStats )
      {
          tiledb::Stats::dump(stdout);
          tiledb::Stats::disable();
      }
  
      if ( status == tiledb::Query::Status::FAILED )
          return CE_Failure;
      else
          return CE_None;
  }
  
  /************************************************************************/
  /*                             IWriteBlock()                            */
  /************************************************************************/
  CPLErr TileDBRasterBand::IWriteBlock( int nBlockXOff,
                                        int nBlockYOff, void * pImage )
  
  {
      if( eAccess == GA_ReadOnly )
      {
          CPLError( CE_Failure, CPLE_NoWriteAccess,
                    "Unable to write block, dataset is opened read only.\n" );
          return CE_Failure;
      }
  
      CPLAssert( poGDS != nullptr
                 && nBlockXOff >= 0
                 && nBlockYOff >= 0
                 && pImage != nullptr );
  
      int nBandIdx = poGDS->nBandStart + nBand - 1;
      if ( !poGDS->bGlobalOrder )
      {
          int nStartX = nBlockXSize * nBlockXOff;
          int nStartY = nBlockYSize * nBlockYOff;
          uint64_t nEndX =  nStartX + nBlockXSize;
          uint64_t nEndY =  nStartY + nBlockYSize;
  
          std::vector<uint64_t> oaSubarray = {
                                          uint64_t( nBandIdx ),
                                          uint64_t( nBandIdx ),
                                          (uint64_t) nStartY,
                                          (uint64_t) nEndY - 1,
                                          (uint64_t) nStartX,
                                          (uint64_t) nEndX - 1 };
  
          if ( poGDS->eIndexMode == PIXEL )
              std::rotate( oaSubarray.begin(),
                          oaSubarray.begin() + 2,
                          oaSubarray.end() );
  
          if ( poGDS->eIndexMode == ATTRIBUTES )
          {
              m_query->set_subarray( std::vector<uint64_t> (
                                          oaSubarray.cbegin() + 2,
                                          oaSubarray.cend() ) );
          }
          else
          {
              m_query->set_subarray( oaSubarray );
          }
      }
      else
      {
          // global order requires ordered blocks (see FlushCache(bool bAtClosing))
          if ( ( nCurrBlockX != nBlockXOff ) || ( nCurrBlockY != nBlockYOff ) )
          {   CPLError( CE_Failure, CPLE_AppDefined,
                      "Non-sequential global write to TileDB.\n");
              return CE_Failure;
          }
          else
          {
              if ( ++nCurrBlockX == poGDS->nBlocksX )
              {
                  nCurrBlockX = 0;
                  nCurrBlockY++;
              }
          }
      }
  
      std::vector<std::unique_ptr<void, decltype(&VSIFree)>> aBlocks;
  
      if ( poGDS->lpoAttributeDS.size() > 0 )
      {
          for ( auto const& poAttrDS: poGDS->lpoAttributeDS )
          {
              GDALRasterBand* poAttrBand = poAttrDS->GetRasterBand( nBandIdx );
              GDALDataType eAttrType = poAttrBand->GetRasterDataType();
              int nBytes = GDALGetDataTypeSizeBytes( eAttrType );
              int nValues = nBlockXSize * nBlockYSize;
              void* pAttrBlock = VSIMalloc( nBytes * nValues );
              aBlocks.emplace_back(pAttrBlock, &VSIFree);
  
              if ( pAttrBlock == nullptr )
              {
                  CPLError(CE_Failure, CPLE_OutOfMemory,
                              "Cannot allocate attribute buffer");
                  return CE_Failure;
              }
  
              int nXSize = nBlockXSize;
              int nYSize = nBlockYSize;
              if( nBlockXOff + nXSize > nRasterXSize )
                  nXSize = nRasterXSize - nBlockXOff;
              if( nBlockYOff + nYSize > nRasterYSize )
                  nYSize = nRasterYSize - nBlockYOff;
  
              poAttrBand->AdviseRead(
                  nBlockXOff, nBlockYOff, nXSize, nYSize,
                  nBlockXSize, nBlockYSize, eAttrType, nullptr
              );
  
              CPLErr eErr = poAttrBand->RasterIO(
                  GF_Read,
                  nBlockXOff, nBlockYOff, nXSize, nYSize,
                  pAttrBlock, nBlockXSize, nBlockYSize,
                  eAttrType, 0, 0, nullptr );
  
              if ( eErr == CE_None )
              {
                  CPLString osName = CPLString().Printf( "%s",
                      CPLGetBasename( poAttrDS->GetDescription() ) );
  
                  SetBuffer(m_query.get(), eAttrType,
                      osName,
                      pAttrBlock, nBlockXSize * nBlockYSize );
              }
              else
              {
                  return eErr;
              }
          }
      }
  
      if ( poGDS->eIndexMode == ATTRIBUTES )
      {
          if ( nBand != poGDS->GetRasterCount() )
          {
              tiledb::Query* q = ( ( TileDBRasterBand* ) poGDS->GetRasterBand( poGDS->nBands ) )->m_query.get();
              SetBuffer( q, eDataType, osAttrName,
                         pImage, nBlockXSize * nBlockYSize );
          }
          else
          {
              SetBuffer( m_query.get(), eDataType, osAttrName, pImage, nBlockXSize * nBlockYSize );
          }
      }
      else
      {
          SetBuffer( m_query.get(), eDataType, osAttrName,
                      pImage, nBlockXSize * nBlockYSize );
      }
  
      if ( ( poGDS->eIndexMode != ATTRIBUTES) ||
           ( nBand == poGDS->GetRasterCount() ) )
      {
          if ( bStats )
              tiledb::Stats::enable();
  
          auto status = m_query->submit();
  
          if ( bStats )
          {
              tiledb::Stats::dump(stdout);
              tiledb::Stats::disable();
          }
  
          if (status == tiledb::Query::Status::FAILED)
              return CE_Failure;
      }
  
      return CE_None;
  }
  
  /************************************************************************/
  /*                       GetColorInterpretation()                       */
  /************************************************************************/
  
  GDALColorInterp TileDBRasterBand::GetColorInterpretation()
  
  {
      if (poGDS->nBands == 1)
          return GCI_GrayIndex;
  
      if ( nBand == 1 )
          return GCI_RedBand;
  
      else if( nBand == 2 )
          return GCI_GreenBand;
  
      else if ( nBand == 3 )
          return GCI_BlueBand;
  
      return GCI_AlphaBand;
  }
  
  /************************************************************************/
  /* ==================================================================== */
  /*                             TileDBDataset                            */
  /* ==================================================================== */
  /************************************************************************/
  
  /************************************************************************/
  /*                           ~TileDBDataset()                           */
  /************************************************************************/
  
  TileDBDataset::~TileDBDataset()
  
  {
      TileDBDataset::FlushCache(true);
  
      // important to finalize arrays before closing array when updating
      if ( bGlobalOrder )
      {
          for( auto&& poBand: GetBands() )
          {
              static_cast<TileDBRasterBand*>( poBand )->Finalize();
          }
      }
  
      try
      {
          if( m_array )
          {
              m_array->close();
          }
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
      }
      CPLDestroyXMLNode( psSubDatasetsTree );
      CSLDestroy( papszSubDatasets );
      CSLDestroy( papszAttributes );
  }
  
  /************************************************************************/
  /*                             AddDimensions()                          */
  /************************************************************************/
  
  CPLErr TileDBDataset::AddDimensions( tiledb::Domain& domain,
                                       tiledb::Dimension& y,
                                       tiledb::Dimension& x,
                                       tiledb::Dimension* poBands )
  
  {
      CPLErr eErr = CE_None;
      switch ( eIndexMode )
      {
      case ATTRIBUTES:
          domain.add_dimensions( y, x );
          CreateAttribute( eDataType, TILEDB_VALUES, nBands );
          break;
      case PIXEL:
          assert( poBands );
          domain.add_dimensions( y, x, *poBands );
          CreateAttribute( eDataType, TILEDB_VALUES, 1 );
          break;
      default: // BAND
          assert( poBands );
          domain.add_dimensions( *poBands, y, x );
          CreateAttribute( eDataType, TILEDB_VALUES, 1 );
          break;
      }
  
      return eErr;
  }
  
  /************************************************************************/
  /*                             FlushCache()                             */
  /************************************************************************/
  
  void TileDBDataset::FlushCache(bool bAtClosing)
  
  {
      BlockBasedFlushCache(bAtClosing);
  
      if( nPamFlags & GPF_DIRTY )
          TrySaveXML();
  }
  
  /************************************************************************/
  /*                           TrySaveXML()                               */
  /************************************************************************/
  
  CPLErr TileDBDataset::TrySaveXML()
  
  {
      if( m_array == nullptr )
          return CE_None;
  
      CPLXMLNode *psTree = nullptr;
      try
      {
          tiledb::VFS vfs( *m_ctx, m_ctx->config() );
  
          nPamFlags &= ~GPF_DIRTY;
  
          if( psPam == nullptr || (nPamFlags & GPF_NOSAVE) )
              return CE_None;
  
  /* -------------------------------------------------------------------- */
  /*      Make sure we know the filename we want to store in.             */
  /* -------------------------------------------------------------------- */
          if( !BuildPamFilename() )
              return CE_None;
  
  /* -------------------------------------------------------------------- */
  /*      Build the XML representation of the auxiliary metadata.          */
  /* -------------------------------------------------------------------- */
          psTree = SerializeToXML( nullptr );
  
          if( psTree == nullptr )
          {
              /* If we have unset all metadata, we have to delete the PAM file */
  #if TILEDB_VERSION_MAJOR == 1 && TILEDB_VERSION_MINOR < 7
              vfs.remove_file(psPam->pszPamFilename);
  #else
              m_array->delete_metadata("_gdal");
  #endif
              return CE_None;
          }
  
          if ( psSubDatasetsTree != nullptr )
          {
              CPLAddXMLChild( psTree, CPLCloneXMLTree( psSubDatasetsTree->psChild ) );
          }
  
  /* -------------------------------------------------------------------- */
  /*      If we are working with a subdataset, we need to integrate       */
  /*      the subdataset tree within the whole existing pam tree,         */
  /*      after removing any old version of the same subdataset.          */
  /* -------------------------------------------------------------------- */
          if( !psPam->osSubdatasetName.empty() )
          {
              CPLXMLNode *psOldTree, *psSubTree;
  
              CPLErrorReset();
              {
                  CPLErrorHandlerPusher oQuietError( CPLQuietErrorHandler );
                  psOldTree = CPLParseXMLFile( psPam->pszPamFilename );
              }
  
              if( psOldTree == nullptr )
                  psOldTree = CPLCreateXMLNode( nullptr, CXT_Element, "PAMDataset" );
  
              for( psSubTree = psOldTree->psChild;
                  psSubTree != nullptr;
                  psSubTree = psSubTree->psNext )
              {
                  if( psSubTree->eType != CXT_Element
                      || !EQUAL(psSubTree->pszValue,"Subdataset") )
                      continue;
  
                  if( !EQUAL(CPLGetXMLValue( psSubTree, "name", "" ),
                          psPam->osSubdatasetName) )
                      continue;
  
                  break;
              }
  
              if( psSubTree == nullptr )
              {
                  psSubTree = CPLCreateXMLNode( psOldTree, CXT_Element,
                                              "Subdataset" );
                  CPLCreateXMLNode(
                      CPLCreateXMLNode( psSubTree, CXT_Attribute, "name" ),
                      CXT_Text, psPam->osSubdatasetName );
              }
  
              CPLXMLNode *psOldPamDataset = CPLGetXMLNode( psSubTree, "PAMDataset");
              if( psOldPamDataset != nullptr )
              {
                  CPLRemoveXMLChild( psSubTree, psOldPamDataset );
                  CPLDestroyXMLNode( psOldPamDataset );
              }
  
              CPLAddXMLChild( psSubTree, psTree );
              psTree = psOldTree;
          }
  
  /* -------------------------------------------------------------------- */
  /*      Try saving the auxiliary metadata.                               */
  /* -------------------------------------------------------------------- */
          bool bSaved = false;
          CPLErrorHandlerPusher oQuietError( CPLQuietErrorHandler );
          char* pszTree = CPLSerializeXMLTree( psTree );
  
  #if TILEDB_VERSION_MAJOR > 1 || TILEDB_VERSION_MINOR >= 7
          if ( eAccess == GA_ReadOnly )
          {
              if ( nTimestamp )
              {
                  auto oMeta = std::unique_ptr<tiledb::Array>(
                      new tiledb::Array( *m_ctx, m_array->uri(), TILEDB_WRITE, nTimestamp )
                  );
                  oMeta->put_metadata("_gdal", TILEDB_UINT8, static_cast<int>(strlen( pszTree )), pszTree);
                  oMeta->close();
              }
              else
              {
                  auto oMeta = std::unique_ptr<tiledb::Array>(
                          new tiledb::Array( *m_ctx, m_array->uri(), TILEDB_WRITE )
                      );
                  oMeta->put_metadata("_gdal", TILEDB_UINT8, static_cast<int>(strlen( pszTree )), pszTree);
                  oMeta->close();
              }
          }
          else
          {
              m_array->put_metadata("_gdal", TILEDB_UINT8, static_cast<int>(strlen( pszTree )), pszTree);
          }
  
          bSaved = true;
  #endif
  
          // cppcheck-suppress knownConditionTrueFalse
          if (!bSaved)
          {
              vfs.touch( psPam->pszPamFilename );
              tiledb::VFS::filebuf fbuf( vfs );
              fbuf.open( psPam->pszPamFilename, std::ios::out );
              std::ostream os(&fbuf);
  
              if (os.good())
              {
                  os.write( pszTree, strlen(pszTree));
                  bSaved = true;
              }
  
              fbuf.close();
          }
  
          CPLFree( pszTree );
  
  /* -------------------------------------------------------------------- */
  /*      If it fails, check if we have a proxy directory for auxiliary    */
  /*      metadata to be stored in, and try to save there.                */
  /* -------------------------------------------------------------------- */
          CPLErr eErr = CE_None;
  
          if( bSaved )
              eErr = CE_None;
          else
          {
              const char *pszBasename = GetDescription();
  
              if( psPam->osPhysicalFilename.length() > 0 )
                  pszBasename = psPam->osPhysicalFilename;
  
              const char *pszNewPam = nullptr;
              if( PamGetProxy(pszBasename) == nullptr
                  && ((pszNewPam = PamAllocateProxy(pszBasename)) != nullptr))
              {
                  CPLErrorReset();
                  CPLFree( psPam->pszPamFilename );
                  psPam->pszPamFilename = CPLStrdup(pszNewPam);
                  eErr = TrySaveXML();
              }
              /* No way we can save into a /vsicurl resource */
              else if( !STARTS_WITH(psPam->pszPamFilename, "/vsicurl") )
              {
                  CPLError( CE_Warning, CPLE_AppDefined,
                          "Unable to save auxiliary information in %s.",
                          psPam->pszPamFilename );
                  eErr = CE_Warning;
              }
          }
  
  /* -------------------------------------------------------------------- */
  /*      Cleanup                                                         */
  /* -------------------------------------------------------------------- */
          if ( psTree )
              CPLDestroyXMLNode( psTree );
  
          return eErr;
      }
      catch(const std::exception& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          if ( psTree )
              CPLDestroyXMLNode( psTree );
          return CE_Failure;
      }
  }
  
  /************************************************************************/
  /*                           TryLoadXML()                               */
  /************************************************************************/
  
  CPLErr TileDBDataset::TryLoadXML( char ** papszSiblingFiles )
  
  {
      return TryLoadCachedXML( papszSiblingFiles, true );
  }
  
  /************************************************************************/
  /*                           TryLoadCachedXML()                               */
  /************************************************************************/
  
  CPLErr TileDBDataset::TryLoadCachedXML( char ** /*papszSiblingFiles*/, bool bReload )
  
  {
      CPLXMLNode *psTree = nullptr;
      try
      {
          PamInitialize();
          tiledb::VFS vfs( *m_ctx, m_ctx->config() );
  
  /* -------------------------------------------------------------------- */
  /*      Clear dirty flag.  Generally when we get to this point is       */
  /*      from a call at the end of the Open() method, and some calls     */
  /*      may have already marked the PAM info as dirty (for instance     */
  /*      setting metadata), but really everything to this point is       */
  /*      reproducible, and so the PAM info should not really be          */
  /*      thought of as dirty.                                            */
  /* -------------------------------------------------------------------- */
          nPamFlags &= ~GPF_DIRTY;
  
  /* -------------------------------------------------------------------- */
  /*      Try reading the file.                                           */
  /* -------------------------------------------------------------------- */
          if( !BuildPamFilename() )
              return CE_None;
  
  /* -------------------------------------------------------------------- */
  /*      In case the PAM filename is a .aux.xml file next to the         */
  /*      physical file and we have a siblings list, then we can skip     */
  /*      stat'ing the filesystem.                                        */
  /* -------------------------------------------------------------------- */
          CPLErr eLastErr = CPLGetLastErrorType();
          int nLastErrNo = CPLGetLastErrorNo();
          CPLString osLastErrorMsg = CPLGetLastErrorMsg();
  
          CPLErrorReset();
          {
              CPLErrorHandlerPusher oQuietError( CPLQuietErrorHandler );
  
  #if TILEDB_VERSION_MAJOR > 1 || TILEDB_VERSION_MINOR >= 7
              if ( bReload )
              {
                  tiledb_datatype_t v_type = TILEDB_UINT8; // CPLSerializeXMLTree returns char*
                  const void* v_r = nullptr;
                  uint32_t v_num = 0;
                  if ( ( eAccess == GA_Update ) && ( m_roArray ) )
                  {
                      m_roArray->get_metadata("_gdal", &v_type, &v_num, &v_r);
                      if ( v_r )
                      {
                          osMetaDoc = CPLString( static_cast<const char*>( v_r ), v_num);
                      }
                  }
                  else
                  {
                      m_array->get_metadata("_gdal", &v_type, &v_num, &v_r);
  
                      if ( v_r )
                      {
                          osMetaDoc = CPLString( static_cast<const char*>( v_r ), v_num);
                      }
                  }
                  psTree = CPLParseXMLString( osMetaDoc );
              }
  #endif
              if ( bReload &&
                   psTree == nullptr &&
                   vfs.is_file( psPam->pszPamFilename ) )
              {
                  auto nBytes = vfs.file_size( psPam->pszPamFilename );
                  tiledb::VFS::filebuf fbuf( vfs );
                  fbuf.open( psPam->pszPamFilename, std::ios::in );
                  std::istream is ( &fbuf );
                  osMetaDoc.resize(nBytes);
                  is.read( ( char* ) osMetaDoc.data(), nBytes );
                  fbuf.close();
                  psTree = CPLParseXMLString( osMetaDoc );
              }
  
              if ( !bReload )
              {
                  psTree = CPLParseXMLString( osMetaDoc );
              }
          }
          CPLErrorReset();
  
          if( eLastErr != CE_None )
              CPLErrorSetState( eLastErr, nLastErrNo, osLastErrorMsg.c_str() );
  
  /* -------------------------------------------------------------------- */
  /*      If we are looking for a subdataset, search for its subtree not. */
  /* -------------------------------------------------------------------- */
          if( psTree && !psPam->osSubdatasetName.empty() )
          {
              CPLXMLNode *psSubTree = psTree->psChild;
  
              for( ;
                  psSubTree != nullptr;
                  psSubTree = psSubTree->psNext )
              {
                  if( psSubTree->eType != CXT_Element
                      || !EQUAL(psSubTree->pszValue,"Subdataset") )
                      continue;
  
                  if( !EQUAL(CPLGetXMLValue( psSubTree, "name", "" ),
                          psPam->osSubdatasetName) )
                      continue;
  
                  psSubTree = CPLGetXMLNode( psSubTree, "PAMDataset" );
                  break;
              }
  
              if( psSubTree != nullptr )
                  psSubTree = CPLCloneXMLTree( psSubTree );
  
              CPLDestroyXMLNode( psTree );
              psTree = psSubTree;
          }
          if( psTree == nullptr )
              return CE_Failure;
  
  /* -------------------------------------------------------------------- */
  /*      Initialize ourselves from this XML tree.                        */
  /* -------------------------------------------------------------------- */
  
          CPLString osPath(CPLGetPath(psPam->pszPamFilename));
          const CPLErr eErr = XMLInit( psTree, osPath );
  
          CPLDestroyXMLNode( psTree );
  
          if( eErr != CE_None )
              PamClear();
  
          return eErr;
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          if ( psTree )
              CPLDestroyXMLNode( psTree );
          return CE_Failure;
      }
  }
  
  /************************************************************************/
  /*                             GetMetadata()                            */
  /************************************************************************/
  
  char **TileDBDataset::GetMetadata(const char *pszDomain)
  
  {
      if( pszDomain != nullptr && EQUAL(pszDomain, "SUBDATASETS") )
      {
          char** papszMeta = CSLDuplicate(GDALPamDataset::GetMetadata( pszDomain ));
          for ( int i = 0; papszMeta && papszMeta[i];i++ )
          {
              if( STARTS_WITH(papszMeta[i], "SUBDATASET_") &&
                  strstr(papszMeta[i], "_NAME=") )
              {
                  char* pszKey = nullptr;
                  const char* pszAttr = CPLParseNameValue(papszMeta[i], &pszKey);
                  if ( pszAttr && !STARTS_WITH( pszAttr, "TILEDB:" ) )
                  {
                      CPLString osAttr(pszAttr);
                      CPLFree( papszMeta[i] );
                      papszMeta[i] =  CPLStrdup(
                          CPLString().Printf(
                              "%s=TILEDB:\"%s\":%s",
                              pszKey,
                              GetDescription(),
                              osAttr.c_str()
                          )
                      );
                  }
                  CPLFree(pszKey);
              }
          }
          m_osSubdatasetMD.Assign(papszMeta);
          return m_osSubdatasetMD.List();
      }
      else
      {
          return GDALPamDataset::GetMetadata( pszDomain );
      }
  }
  
  /************************************************************************/
  /*                           AddFilter()                                */
  /************************************************************************/
  
  CPLErr TileDBDataset::AddFilter( const char* pszFilterName, const int level )
  
  {
      try
      {
          if (pszFilterName == nullptr)
              m_filterList->add_filter( tiledb::Filter( *m_ctx, TILEDB_FILTER_NONE )
                  .set_option( TILEDB_COMPRESSION_LEVEL, level ) );
          else if EQUAL(pszFilterName, "GZIP")
              m_filterList->add_filter( tiledb::Filter( *m_ctx, TILEDB_FILTER_GZIP )
                  .set_option( TILEDB_COMPRESSION_LEVEL, level ) );
          else if EQUAL(pszFilterName, "ZSTD")
              m_filterList->add_filter( tiledb::Filter( *m_ctx, TILEDB_FILTER_ZSTD )
                  .set_option( TILEDB_COMPRESSION_LEVEL, level ) );
          else if EQUAL(pszFilterName, "LZ4")
              m_filterList->add_filter( tiledb::Filter( *m_ctx, TILEDB_FILTER_LZ4 )
                  .set_option( TILEDB_COMPRESSION_LEVEL, level ) );
          else if EQUAL(pszFilterName, "RLE")
              m_filterList->add_filter( tiledb::Filter( *m_ctx, TILEDB_FILTER_RLE )
                  .set_option( TILEDB_COMPRESSION_LEVEL, level ) );
          else if EQUAL(pszFilterName, "BZIP2")
              m_filterList->add_filter( tiledb::Filter( *m_ctx, TILEDB_FILTER_BZIP2 )
                  .set_option( TILEDB_COMPRESSION_LEVEL, level ) );
          else if EQUAL(pszFilterName, "DOUBLE-DELTA")
              m_filterList->add_filter( tiledb::Filter( *m_ctx,
                  TILEDB_FILTER_DOUBLE_DELTA ) );
          else if EQUAL(pszFilterName, "POSITIVE-DELTA")
              m_filterList->add_filter( tiledb::Filter( *m_ctx,
                  TILEDB_FILTER_POSITIVE_DELTA ) );
          else
              return CE_Failure;
  
          return CE_None;
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          return CE_Failure;
      }
  }
  
  /************************************************************************/
  /*                              Delete()                                */
  /************************************************************************/
  
  CPLErr TileDBDataset::Delete( const char * pszFilename )
  
  {
      try
      {
          tiledb::Context ctx;
          tiledb::VFS vfs( ctx );
          CPLString osArrayPath = vsi_to_tiledb_uri( pszFilename );
  
          if ( vfs.is_dir( osArrayPath ) )
          {
              vfs.remove_dir( osArrayPath );
              return CE_None;
          }
          else
              return CE_Failure;
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          return CE_Failure;
      }
  }
  
  /************************************************************************/
  /*                              Identify()                              */
  /************************************************************************/
  
  int TileDBDataset::Identify( GDALOpenInfo * poOpenInfo )
  
  {
      if( STARTS_WITH_CI(poOpenInfo->pszFilename, "TILEDB:") )
      {
          return TRUE;
      }
  
      try
      {
          const char* pszConfig = CSLFetchNameValue(
              poOpenInfo->papszOpenOptions, "TILEDB_CONFIG" );
  
          if ( pszConfig != nullptr )
          {
              return TRUE;
          }
  
          if( poOpenInfo->bIsDirectory ||
              ( (STARTS_WITH_CI( poOpenInfo->pszFilename, "/VSIS3/" ) ||
                 STARTS_WITH_CI( poOpenInfo->pszFilename, "/VSIGS/" )) &&
                  !EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "tif")) )
          {
              tiledb::Context ctx;
              CPLString osArrayPath = vsi_to_tiledb_uri( poOpenInfo->pszFilename );
              return tiledb::Object::object( ctx, osArrayPath ).type() == tiledb::Object::Type::Array;
          }
  
          return FALSE;
      }
      catch( ... )
      {
          return FALSE;
      }
  }
  
  /************************************************************************/
  /*                                Open()                                */
  /************************************************************************/
  
  GDALDataset *TileDBDataset::Open( GDALOpenInfo * poOpenInfo )
  
  {
      try
      {
          if( !Identify( poOpenInfo ) )
              return nullptr;
  
          auto poDS = cpl::make_unique<TileDBDataset>();
  
          const char* pszConfig = CSLFetchNameValue(
                                      poOpenInfo->papszOpenOptions,
                                      "TILEDB_CONFIG" );
  
          const char* pszTimestamp = CSLFetchNameValue(
                                          poOpenInfo->papszOpenOptions,
                                          "TILEDB_TIMESTAMP" );
  
          if( pszConfig != nullptr )
          {
              tiledb::Config cfg( pszConfig );
              poDS->m_ctx.reset( new tiledb::Context( cfg ) );
          }
          else
          {
              poDS->m_ctx.reset( new tiledb::Context() );
          }
  
          if ( pszTimestamp )
              poDS->nTimestamp = atoi ( pszTimestamp );
  
          CPLString osArrayPath;
          CPLString osAux;
          CPLString osSubdataset;
  
          const char* pszAttr = CSLFetchNameValue(
              poOpenInfo->papszOpenOptions,
              "TILEDB_ATTRIBUTE" );
  
          if( STARTS_WITH_CI(poOpenInfo->pszFilename, "TILEDB:") &&
              !STARTS_WITH_CI(poOpenInfo->pszFilename, "TILEDB://") )
          {
              // form required read attributes and open file
              // Create a corresponding GDALDataset.
              CPLStringList apszName(
                  CSLTokenizeString2(poOpenInfo->pszFilename, ":",
                                  CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
  
              if( apszName.size() != 3 )
              {
                  return nullptr;
              }
  
              osArrayPath = vsi_to_tiledb_uri(apszName[1]);
              osSubdataset = apszName[2];
              poDS->SetSubdatasetName( osSubdataset.c_str() );
          }
          else
          {
              if ( pszAttr != nullptr )
              {
                  poDS->SetSubdatasetName( pszAttr );
              }
  
              osArrayPath = vsi_to_tiledb_uri( poOpenInfo->pszFilename );
          }
  
          const char* pszArrayName = CPLGetBasename( osArrayPath );
          osAux.Printf( "%s.tdb", pszArrayName );
  
          // aux file is in array folder
          poDS->SetPhysicalFilename( CPLFormFilename( osArrayPath, osAux, nullptr ) );
          // Initialize any PAM information.
          poDS->SetDescription( osArrayPath );
  
          tiledb_query_type_t eMode = TILEDB_READ;
          if ( poOpenInfo->eAccess == GA_Update )
          {
              eMode = TILEDB_WRITE;
              poDS->m_roCtx.reset(new tiledb::Context( poDS->m_ctx->config() ) );
              if ( poDS->nTimestamp )
                  poDS->m_roArray.reset(
                      new tiledb::Array( *poDS->m_roCtx, osArrayPath,
                                          TILEDB_READ, poDS->nTimestamp ) );
              else
                  poDS->m_roArray.reset(
                      new tiledb::Array( *poDS->m_roCtx, osArrayPath, TILEDB_READ )
                  );
          }
  
          if ( poDS->nTimestamp )
              poDS->m_array.reset(
                  new tiledb::Array( *poDS->m_ctx, osArrayPath,
                                      eMode, poDS->nTimestamp ) );
          else
              poDS->m_array.reset(
                  new tiledb::Array( *poDS->m_ctx, osArrayPath, eMode ) );
  
          poDS->eAccess = poOpenInfo->eAccess;
  
          // dependent on PAM metadata for information about array
          poDS->TryLoadXML();
  
          tiledb::ArraySchema schema = poDS->m_array->schema();
  
          char ** papszStructMeta = poDS->GetMetadata( "IMAGE_STRUCTURE" );
          const char* pszXSize = CSLFetchNameValue( papszStructMeta, "X_SIZE");
          if ( pszXSize )
          {
              poDS->nRasterXSize = atoi( pszXSize );
              if ( poDS->nRasterXSize <= 0 )
              {
                  CPLError( CE_Failure, CPLE_AppDefined,
                      "Width %i should be greater than zero.", poDS->nRasterXSize );
                  return nullptr;
              }
          }
  
          const char* pszYSize = CSLFetchNameValue( papszStructMeta, "Y_SIZE");
          if ( pszYSize )
          {
              poDS->nRasterYSize = atoi( pszYSize );
              if ( poDS->nRasterYSize <= 0 )
              {
                  CPLError( CE_Failure, CPLE_AppDefined,
                      "Height %i should be greater than zero.", poDS->nRasterYSize );
                  return nullptr;
              }
          }
  
          const char* pszNBits = CSLFetchNameValue( papszStructMeta, "NBITS");
          if ( pszNBits )
          {
              poDS->nBitsPerSample = atoi( pszNBits );
          }
  
          const char* pszDataType = CSLFetchNameValue( papszStructMeta, "DATA_TYPE");
          if ( pszDataType )
          {
              // handle the case where arrays have been written with int type (2.5.0)
              GDALDataType eDT = GDALGetDataTypeByName( pszDataType );
              if ( eDT == GDT_Unknown )
              {
                  int t = atoi( pszDataType );
                  if ( ( t > 0 ) && ( t < GDT_TypeCount ) )
                      poDS->eDataType = static_cast<GDALDataType>( atoi( pszDataType ) );
                  else
                  {
                      CPLError( CE_Failure, CPLE_AppDefined,
                          "Unknown data type %s.", pszDataType );
                      return nullptr;
                  }
              }
              else
              {
                  poDS->eDataType = eDT;
              }
          }
  
          const char* pszIndexMode = CSLFetchNameValue( papszStructMeta, "INTERLEAVE");
  
          if ( pszIndexMode )
              option_to_index_type( pszIndexMode, poDS->eIndexMode );
  
          std::vector<tiledb::Dimension> dims = schema.domain().dimensions();
  
          if( ( dims.size() == 2 ) || ( dims.size() == 3) )
          {
              if ( dims.size() == 3 )
              {
                  if ( ( pszAttr != nullptr ) && ( schema.attributes().count( pszAttr ) == 0 ) )
                  {
                      CPLError( CE_Failure, CPLE_NotSupported,
                          "%s attribute is not found in TileDB schema.",
                          pszAttr );
                      return nullptr;
                  }
  
                  if ( poDS->eIndexMode == PIXEL )
                      std::rotate( dims.begin(), dims.begin() + 2, dims.end() );
  
                  poDS->nBandStart = static_cast<int>(dims[0].domain<uint64_t>().first);
                  poDS->nBands = static_cast<int>(dims[0].domain<uint64_t>().second
                                  - dims[0].domain<uint64_t>().first + 1);
                  poDS->nBlockYSize = static_cast<int>(dims[1].tile_extent<uint64_t>());
                  poDS->nBlockXSize = static_cast<int>(dims[2].tile_extent<uint64_t>());
              }
              else
              {
                  const char* pszBands = poDS->GetMetadataItem("NUM_BANDS",
                                                              "IMAGE_STRUCTURE");
                  if ( pszBands )
                  {
                      poDS->nBands = atoi( pszBands );
                      poDS->nBlockYSize = static_cast<int>(dims[0].tile_extent<uint64_t>());
                      poDS->nBlockXSize = static_cast<int>(dims[1].tile_extent<uint64_t>());
                  }
  
                  poDS->eIndexMode = ATTRIBUTES;
              }
          }
          else
          {
              CPLError( CE_Failure, CPLE_AppDefined,
                  "Wrong number of dimensions %li expected 2 or 3.", dims.size() );
              return nullptr;
          }
  
          poDS->nBlocksX = DIV_ROUND_UP( poDS->nRasterXSize, poDS->nBlockXSize );
          poDS->nBlocksY = DIV_ROUND_UP( poDS->nRasterYSize, poDS->nBlockYSize );
  
          if ( dims.size() == 3 )
          {
              // Create band information objects.
              for ( int i = 1;i <= poDS->nBands;++i )
              {
                  if ( pszAttr == nullptr )
                      poDS->SetBand( i, new TileDBRasterBand( poDS.get(), i ) );
                  else
                      poDS->SetBand( i,
                          new TileDBRasterBand(
                              poDS.get(), i, CPLString( pszAttr ) ) );
              }
          }
          else // subdatasets or only attributes
          {
              if ( ( poOpenInfo->eAccess == GA_Update ) &&
                   ( poDS->GetMetadata( "SUBDATASETS" ) != nullptr ) )
              {
                  CPLError( CE_Failure, CPLE_NotSupported,
                          "The TileDB driver does not support update access "
                          "to subdatasets." );
                  return nullptr;
              }
  
              if( !osSubdataset.empty() )
              {
                  // do we have the attribute in the schema
                  if ( schema.attributes().count( osSubdataset ) )
                  {
                      poDS->SetBand(1,
                          new TileDBRasterBand(
                              poDS.get(), 1, osSubdataset ) );
                  }
                  else
                  {
                      if ( schema.attributes().count( osSubdataset + "_1" ) )
                      {
                          // Create band information objects.
                          for ( int i = 1;i <= poDS->nBands;++i )
                          {
                              CPLString osAttr = CPLString().Printf("%s_%d",
                                                  osSubdataset.c_str(), i);
                              poDS->SetBand( i,
                                  new TileDBRasterBand(
                                      poDS.get(), i,
                                      osAttr ) );
                          }
                      }
                      else
                      {
                          CPLError( CE_Failure, CPLE_NotSupported,
                              "%s attribute is not found in TileDB schema.",
                              osSubdataset.c_str() );
                          return nullptr;
                      }
                  }
              }
              else
              {
                  char** papszMeta = poDS->GetMetadata( "SUBDATASETS" );
                  if ( papszMeta != nullptr )
                  {
                      if ( ( CSLCount( papszMeta ) / 2 ) == 1 )
                      {
                          CPLString osDSName =
                              CSLFetchNameValueDef(poDS->papszSubDatasets,
                                                  "SUBDATASET_1_NAME", "");
                          return (GDALDataset *)GDALOpen( osDSName, poOpenInfo->eAccess );
                      }
                  }
                  else if ( poDS->eIndexMode == ATTRIBUTES )
                  {
                      poDS->nBands = schema.attribute_num();
                      poDS->nBlockYSize = static_cast<int>(dims[0].tile_extent<uint64_t>());
                      poDS->nBlockXSize = static_cast<int>(dims[1].tile_extent<uint64_t>());
  
                      // Create band information objects.
                      for ( int i = 1;i <= poDS->nBands;++i )
                      {
                          CPLString osAttr = TILEDB_VALUES + CPLString().Printf( "_%i", i+1 );
                          poDS->SetBand( i,
                              new TileDBRasterBand(
                                  poDS.get(), i,
                                  osAttr ) );
                      }
                  }
                  else
                  {
                      CPLError( CE_Failure, CPLE_NotSupported,
                          "%s is missing required TileDB subdataset metadata.",
                          osArrayPath.c_str() );
                      return nullptr;
                  }
              }
          }
  
          // reload metadata now that bands are created to populate band metadata
          poDS->TryLoadCachedXML( nullptr, false );
  
          tiledb::VFS vfs( *poDS->m_ctx, poDS->m_ctx->config() );
  
          if ( !STARTS_WITH_CI(osArrayPath, "TILEDB:")
                  &&  vfs.is_dir( osArrayPath ) )
              poDS->oOvManager.Initialize( poDS.get(), ":::VIRTUAL:::" );
          else
              CPLError( CE_Warning, CPLE_AppDefined,
                  "Overviews not supported for network writes." );
  
        return poDS.release();
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          return nullptr;
      }
  }
  
  /************************************************************************/
  /*                              CreateAttribute()                       */
  /************************************************************************/
  
  CPLErr TileDBDataset::CreateAttribute( GDALDataType eType,
              const CPLString& osAttrName, const int nSubRasterCount )
  {
      try
      {
          for ( int i = 0; i < nSubRasterCount;++i )
          {
              CPLString osName( osAttrName );
              // a few special cases
              // remove any leading slashes or
              // additional slashes as in the case of hdf5
              if STARTS_WITH( osName, "//" )
              {
                  osName = osName.substr(2);
              }
  
              osName.replaceAll( "/", "_" );
              CPLString osPrettyName = osName;
  
              if ( ( eIndexMode == ATTRIBUTES ) ||
                 ( ( bHasSubDatasets ) && ( nSubRasterCount > 1 ) ) )
              {
                  osName = CPLString().Printf( "%s_%d", osName.c_str(), i + 1);
              }
  
              switch (eType)
              {
                  case GDT_Byte:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<unsigned char>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 8;
                      break;
                  }
                  case GDT_UInt16:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<unsigned short>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 16;
                      break;
                  }
                  case GDT_UInt32:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<unsigned int>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 32;
                      break;
                  }
                  case GDT_Int16:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<short>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 16;
                      break;
                  }
                  case GDT_Int32:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<int>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 32;
                      break;
                  }
                  case GDT_Float32:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<float>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 32;
                      break;
                  }
                  case GDT_Float64:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<double>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 64;
                      break;
                  }
                  case GDT_CInt16:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<short[2]>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 16;
                      break;
                  }
                  case GDT_CInt32:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<int[2]>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 32;
                      break;
                  }
                  case GDT_CFloat32:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<float[2]>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 32;
                      break;
                  }
                  case GDT_CFloat64:
                  {
                      m_schema->add_attribute( tiledb::Attribute::create<double[2]>( *m_ctx, osName, *m_filterList ) );
                      nBitsPerSample = 64;
                      break;
                  }
                  default:
                      return CE_Failure;
              }
  
              if ( ( bHasSubDatasets ) && ( i == 0 ) )
              {
                  ++nSubDataCount;
  
                  CPLString osDim;
                  switch( nSubRasterCount )
                  {
                  case 2:
                      osDim.Printf( "%dx%d", nRasterXSize, nRasterYSize );
                      break;
                  default:
                      osDim.Printf( "%dx%dx%d", nSubRasterCount,
                          nRasterXSize, nRasterYSize );
                      break;
                  }
  
                  papszSubDatasets = CSLSetNameValue( papszSubDatasets,
                                  CPLString().Printf( "SUBDATASET_%d_NAME",
                                      nSubDataCount ),
                                  CPLString().Printf( "%s", osPrettyName.c_str() ) );
  
                  papszSubDatasets = CSLSetNameValue( papszSubDatasets,
                      CPLString().Printf( "SUBDATASET_%d_DESC", nSubDataCount ),
                      CPLString().Printf( "[%s] %s (%s)",
                          osDim.c_str(),
                          osPrettyName.c_str(),
                          GDALGetDataTypeName( eType )
                      ) );
  
  
                  // add to PAM metadata
                  if ( psSubDatasetsTree == nullptr )
                  {
                      psSubDatasetsTree = CPLCreateXMLNode( nullptr,
                                              CXT_Element, "PAMDataset" );
                  }
  
                  CPLXMLNode* psSubNode = CPLCreateXMLNode(
                          psSubDatasetsTree, CXT_Element, "Subdataset" );
                  CPLAddXMLAttributeAndValue( psSubNode, "name", osPrettyName.c_str() );
  
                  CPLXMLNode* psMetaNode = CPLCreateXMLNode(
                              CPLCreateXMLNode(
                              psSubNode, CXT_Element, "PAMDataset" ),
                          CXT_Element, "Metadata");
                  CPLAddXMLAttributeAndValue( psMetaNode, "domain", "IMAGE_STRUCTURE" );
  
                  CPLAddXMLAttributeAndValue(
                      CPLCreateXMLElementAndValue( psMetaNode, "MDI",
                          CPLString().Printf( "%d", nRasterXSize ) ),
                      "KEY",
                      "X_SIZE"
                  );
  
                  CPLAddXMLAttributeAndValue(
                      CPLCreateXMLElementAndValue( psMetaNode, "MDI",
                          CPLString().Printf( "%d", nRasterYSize ) ),
                      "KEY",
                      "Y_SIZE"
                  );
  
                  CPLAddXMLAttributeAndValue(
                      CPLCreateXMLElementAndValue( psMetaNode, "MDI",
                          CPLString().Printf( "%s", GDALGetDataTypeName( eType ) ) ),
                      "KEY",
                      "DATA_TYPE"
                  );
  
                  if ( lpoAttributeDS.size() > 0 )
                  {
                      CPLAddXMLAttributeAndValue(
                          CPLCreateXMLElementAndValue( psMetaNode, "MDI",
                              CPLString().Printf( "%d", nBands ) ),
                          "KEY",
                          "NUM_BANDS"
                      );
                  }
                  else
                  {
                      CPLAddXMLAttributeAndValue(
                          CPLCreateXMLElementAndValue( psMetaNode, "MDI",
                              CPLString().Printf( "%d", nSubRasterCount ) ),
                          "KEY",
                          "NUM_BANDS"
                      );
                  }
  
                  CPLAddXMLAttributeAndValue(
                      CPLCreateXMLElementAndValue( psMetaNode, "MDI",
                          CPLString().Printf( "%d", nBitsPerSample ) ),
                      "KEY",
                      "NBITS"
                  );
              }
          }
          return CE_None;
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          return CE_Failure;
      }
  }
  
  /************************************************************************/
  /*                              SetBlockSize()                          */
  /************************************************************************/
  
  void TileDBDataset::SetBlockSize( GDALRasterBand* poBand, char** &papszOptions)
  
  {
      int nX = 0;
      int nY = 0;
      poBand->GetBlockSize( &nX, &nY );
  
      if ( CSLFetchNameValue(papszOptions, "BLOCKXSIZE") == nullptr )
      {
          papszOptions = CSLSetNameValue(
                              papszOptions, "BLOCKXSIZE",
                              CPLString().Printf( "%d", nX ) );
      }
  
      if ( CSLFetchNameValue(papszOptions, "BLOCKYSIZE") == nullptr )
      {
          papszOptions = CSLSetNameValue(
                              papszOptions, "BLOCKYSIZE",
                              CPLString().Printf( "%d", nY ) );
      }
  }
  
  /************************************************************************/
  /*                              CreateLL()                              */
  /*                                                                      */
  /*      Shared functionality between TileDBDataset::Create() and        */
  /*      TileDBDataset::CreateCopy() for creating TileDB array based on  */
  /*      a set of options and a configuration.                           */
  /************************************************************************/
  
  TileDBDataset* TileDBDataset::CreateLL( const char *pszFilename,
                           int nXSize, int nYSize, int nBands,
                           GDALDataType eType,
                           char **papszOptions )
  {
      try
      {
          if( ( nXSize <= 0 && nYSize <= 0 ) )
          {
              return nullptr;
          }
  
          auto poDS = cpl::make_unique<TileDBDataset>();
          poDS->nRasterXSize = nXSize;
          poDS->nRasterYSize = nYSize;
          poDS->eDataType = eType;
          poDS->nBands = nBands;
          poDS->eAccess = GA_Update;
          poDS->bGlobalOrder = true;
  
          if ( nBands == 0 ) // subdatasets
          {
              poDS->eIndexMode = ATTRIBUTES;
          }
          else
          {
              const char* pszIndexMode = CSLFetchNameValue(
                                      papszOptions, "INTERLEAVE" );
  
              if ( option_to_index_type( pszIndexMode, poDS->eIndexMode ) )
                  return nullptr;
          }
  
          const char* pszConfig = CSLFetchNameValue( papszOptions, "TILEDB_CONFIG" );
          if( pszConfig != nullptr )
          {
              tiledb::Config cfg( pszConfig );
              poDS->m_ctx.reset( new tiledb::Context( cfg ) );
          }
          else
          {
              poDS->m_ctx.reset( new tiledb::Context() );
          }
  
          const char* pszCompression = CSLFetchNameValue(
                                          papszOptions, "COMPRESSION" );
          const char* pszCompressionLevel = CSLFetchNameValue(
                                          papszOptions, "COMPRESSION_LEVEL" );
  
          const char* pszBlockXSize = CSLFetchNameValue( papszOptions, "BLOCKXSIZE" );
          poDS->nBlockXSize = ( pszBlockXSize ) ? atoi( pszBlockXSize ) : 256;
          const char* pszBlockYSize = CSLFetchNameValue( papszOptions, "BLOCKYSIZE" );
          poDS->nBlockYSize = ( pszBlockYSize ) ? atoi( pszBlockYSize ) : 256;
          poDS->bStats = CSLFetchBoolean( papszOptions, "STATS", FALSE );
  
          const char* pszTimestamp = CSLFetchNameValue(
                                          papszOptions,
                                          "TILEDB_TIMESTAMP" );
          if ( pszTimestamp != nullptr )
              poDS->nTimestamp = atoi( pszTimestamp );
  
          // set dimensions and attribute type for schema
          poDS->m_schema.reset( new tiledb::ArraySchema( *poDS->m_ctx, TILEDB_DENSE ) );
          poDS->m_schema->set_tile_order( TILEDB_ROW_MAJOR );
          poDS->m_schema->set_cell_order( TILEDB_ROW_MAJOR );
  
          poDS->m_filterList.reset(new tiledb::FilterList(*poDS->m_ctx));
  
          if (pszCompression != nullptr)
          {
              int nLevel = ( pszCompressionLevel ) ? atoi( pszCompressionLevel ) : -1;
              if ( poDS->AddFilter( pszCompression, nLevel ) == CE_None )
              {
                  poDS->SetMetadataItem( "COMPRESSION", pszCompression,
                                      "IMAGE_STRUCTURE" );
                  poDS->m_schema->set_coords_filter_list( *poDS->m_filterList );
              }
          }
  
          CPLString osAux;
          const char* pszArrayName = CPLGetBasename( pszFilename );
          osAux.Printf( "%s.tdb", pszArrayName );
  
          poDS->SetPhysicalFilename( CPLFormFilename(
                                      pszFilename, osAux.c_str(), nullptr ) );
  
          // Initialize PAM information.
          poDS->SetDescription( pszFilename );
  
          // this driver enforces that all subdatasets are the same size
          tiledb::Domain domain( *poDS->m_ctx );
  
          // Note the dimension bounds are inclusive and are expanded to the match the block size
          poDS->nBlocksX = DIV_ROUND_UP( nXSize, poDS->nBlockXSize );
          poDS->nBlocksY = DIV_ROUND_UP( nYSize, poDS->nBlockYSize );
  
          uint64_t w = (uint64_t)poDS->nBlocksX * poDS->nBlockXSize - 1;
          uint64_t h = (uint64_t)poDS->nBlocksY * poDS->nBlockYSize - 1;
  
          auto d1 = tiledb::Dimension::create<uint64_t>(
                      *poDS->m_ctx, "X", {0, w},
                      uint64_t( poDS->nBlockXSize ) );
          auto d2 = tiledb::Dimension::create<uint64_t>(
                      *poDS->m_ctx, "Y",{0, h},
                      uint64_t( poDS->nBlockYSize ) );
  
          if ( ( nBands == 0 ) || ( poDS->eIndexMode == ATTRIBUTES ) )
          {
              poDS->AddDimensions( domain, d2, d1, nullptr );
          }
          else
          {
              auto d3 = tiledb::Dimension::create<uint64_t>( *poDS->m_ctx, "BANDS", {1, uint64_t( nBands )}, 1);
              poDS->AddDimensions( domain, d2, d1, &d3 );
          }
  
          poDS->m_schema->set_domain( domain ).set_order({{TILEDB_ROW_MAJOR, TILEDB_ROW_MAJOR}});;
  
          // register additional attributes to the pixel value, these will be
          // be reported as subdatasets on future reads
          poDS->papszAttributes = CSLFetchNameValueMultiple(
                                                      papszOptions, "TILEDB_ATTRIBUTE" );
  
          for( int i = 0;
              poDS->papszAttributes != nullptr &&
              poDS->papszAttributes[i] != nullptr;
              i++ )
          {
              // modeling additional attributes as subdatasets
              poDS->bHasSubDatasets = true;
              // check each attribute is a GDAL source
              std::unique_ptr<GDALDataset> poAttrDS(
                  GDALDataset::Open( poDS->papszAttributes[i],
                  GA_ReadOnly ));
  
              if ( poAttrDS != nullptr )
              {
                  const char* pszAttrName = CPLGetBasename(
                                              poAttrDS->GetDescription());
                  // check each is co-registered
                  // candidate band
                  int nAttrBands = poAttrDS->GetRasterCount();
                  if ( nAttrBands > 0 )
                  {
                      GDALRasterBand* poAttrBand = poAttrDS->GetRasterBand( 1 );
  
                      if ( ( poAttrBand->GetXSize() == poDS->nRasterXSize )
                          && ( poAttrBand->GetYSize() == poDS->nRasterYSize )
                          && ( poDS->nBands == nAttrBands ) )
                      {
                          // could check geotransform, but it is sufficient
                          // that cartesian dimensions are equal
                          poDS->lpoAttributeDS.push_back( std::move( poAttrDS ) );
                          poDS->CreateAttribute( poAttrBand->GetRasterDataType(),
                              pszAttrName,
                              1
                          );
                      }
                      else
                      {
                          CPLError( CE_Warning, CPLE_AppDefined,
                                  "Skipping %s as it has a different dimension\n",
                                  poDS->papszAttributes[i]);
                      }
                  }
                  else
                  {
                      CPLError( CE_Warning, CPLE_AppDefined,
                              "Skipping %s as it doesn't have any bands\n",
                              poDS->papszAttributes[i]);
                  }
              }
              else
              {
                  CPLError( CE_Warning, CPLE_AppDefined,
                          "Skipping %s, not recognized as a GDAL dataset\n",
                          poDS->papszAttributes[i]);
              }
          }
  
          return poDS.release();
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "TileDB: %s", e.what() );
          return nullptr;
      }
  }
  
  /************************************************************************/
  /*                              CopySubDatasets()                       */
  /*                                                                      */
  /*      Copy SubDatasets from src to a TileDBDataset                    */
  /*                                                                      */
  /************************************************************************/
  
  CPLErr TileDBDataset::CopySubDatasets( GDALDataset* poSrcDS,
                                      TileDBDataset* poDstDS,
                                      GDALProgressFunc pfnProgress,
                                      void *pProgressData )
  
  {
      try
      {
          std::vector<std::unique_ptr<GDALDataset>> apoDatasets;
          poDstDS->bHasSubDatasets = true;
          char ** papszSrcSubDatasets = poSrcDS->GetMetadata( "SUBDATASETS" );
          if( !papszSrcSubDatasets )
              return CE_Failure;
          const char* pszSubDSName = CSLFetchNameValue(papszSrcSubDatasets,
                                                          "SUBDATASET_1_NAME");
          if( !pszSubDSName )
              return CE_Failure;
  
          CPLStringList apszTokens(
                  CSLTokenizeString2(pszSubDSName, ":",
                                  CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
          // FIXME? this is tailored for HDF5-like subdataset names
          // HDF5:foo.hdf5:attrname.
          if( apszTokens.size() != 3 )
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                      "Cannot guess attribute name in %s", pszSubDSName);
              return CE_Failure;
          }
  
          std::unique_ptr<GDALDataset> poSubDataset(GDALDataset::Open(pszSubDSName));
          if( poSubDataset.get() == nullptr || poSubDataset->GetRasterCount() == 0 )
          {
              return CE_Failure;
          }
  
          uint64_t nSubXSize = poSubDataset->GetRasterXSize();
          uint64_t nSubYSize = poSubDataset->GetRasterYSize();
  
          const char* pszAttrName = apszTokens[2];
  
          poDstDS->CreateAttribute( poSubDataset->GetRasterBand( 1 )
                                                      ->GetRasterDataType(),
                                                  pszAttrName,
                                                  poSubDataset->GetRasterCount() );
          apoDatasets.push_back( std::move(poSubDataset) );
  
          for( int i = 0; papszSrcSubDatasets[i] != nullptr;i ++ )
          {
              if( STARTS_WITH_CI(papszSrcSubDatasets[i], "SUBDATASET_1_NAME=") ||
                  strstr(papszSrcSubDatasets[i], "_DESC=") )
              {
                  continue;
              }
              pszSubDSName = CPLParseNameValue(papszSrcSubDatasets[i], nullptr);
              apszTokens =
                  CSLTokenizeString2(pszSubDSName, ":",
                                  CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES);
              if( apszTokens.size() != 3 )
              {
                  CPLError(CE_Failure, CPLE_AppDefined,
                          "Cannot guess attribute name in %s", pszSubDSName);
                  continue;
              }
  
              std::unique_ptr<GDALDataset> poSubDS(GDALDataset::Open(pszSubDSName));
              if ( ( poSubDS != nullptr ) && poSubDS->GetRasterCount() > 0 )
              {
                  GDALRasterBand* poBand = poSubDS->GetRasterBand( 1 );
                  int nBlockXSize, nBlockYSize;
                  poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
  
                  if ( ( poSubDS->GetRasterXSize() != (int) nSubXSize ) ||
                      ( poSubDS->GetRasterYSize() != (int) nSubYSize ) ||
                      ( nBlockXSize != poDstDS->nBlockXSize) ||
                      ( nBlockYSize != poDstDS->nBlockYSize ))
                  {
                      CPLError( CE_Warning, CPLE_AppDefined,
                          "Sub-datasets must have the same dimension,"
                          " and block sizes, skipping %s\n",
                          pszSubDSName );
                  }
                  else
                  {
                      pszAttrName = apszTokens[2];
                      poDstDS->CreateAttribute(
                          poSubDS->GetRasterBand(1)->GetRasterDataType(),
                          pszAttrName, poSubDS->GetRasterCount() );
                      apoDatasets.push_back( std::move(poSubDS) );
                  }
              }
              else
              {
                  CPLError( CE_Warning, CPLE_AppDefined,
                      "Sub-datasets must be not null and contain data in bands,"
                      "skipping %s\n",
                      pszSubDSName );
              }
          }
  
          poDstDS->SetMetadata( poDstDS->papszSubDatasets, "SUBDATASETS" );
          tiledb::Array::create( poDstDS->GetDescription(), *poDstDS->m_schema );
  
          if ( poDstDS->nTimestamp )
              poDstDS->m_array.reset( new tiledb::Array(
                  *poDstDS->m_ctx, poDstDS->GetDescription(),
                  TILEDB_WRITE, poDstDS->nTimestamp ) );
          else
              poDstDS->m_array.reset( new tiledb::Array(
                      *poDstDS->m_ctx, poDstDS->GetDescription(), TILEDB_WRITE ) );
  
          /* --------------------------------------------------------  */
          /*      Report preliminary (0) progress.                     */
          /* --------------------------------------------------------- */
          if( !pfnProgress( 0.0, nullptr, pProgressData ) )
          {
              CPLError( CE_Failure, CPLE_UserInterrupt,
                      "User terminated CreateCopy()" );
              return CE_Failure;
          }
  
          // copy over subdatasets by block
          tiledb::Query query( *poDstDS->m_ctx, *poDstDS->m_array );
          query.set_layout( TILEDB_GLOBAL_ORDER );
          int nTotalBlocks = poDstDS->nBlocksX * poDstDS->nBlocksY;
  
          // row-major
          for ( int j = 0; j < poDstDS->nBlocksY;++j )
          {
              for ( int i = 0; i < poDstDS->nBlocksX;++i)
              {
                  std::vector<std::unique_ptr<void, decltype(&VSIFree)>> aBlocks;
                  // have to write set all tiledb attributes on write
                  int iAttr = 0;
                  for ( auto& poSubDS : apoDatasets )
                  {
                      GDALDataType eDT = poSubDS->GetRasterBand( 1 )->
                                                  GetRasterDataType();
  
                      for ( int b = 1;b <= poSubDS->GetRasterCount();++b )
                      {
                          int nBytes = GDALGetDataTypeSizeBytes( eDT );
                          int nValues = nBytes * poDstDS->nBlockXSize * poDstDS->nBlockYSize;
                          void* pBlock = VSIMalloc( nBytes * nValues );
                          aBlocks.emplace_back( pBlock, &VSIFree );
                          GDALRasterBand* poBand = poSubDS->GetRasterBand( b );
                          if ( poBand->ReadBlock( i, j, pBlock ) == CE_None )
                          {
                              SetBuffer( &query, eDT,
                                  poDstDS->m_schema->attribute( iAttr++ ).name(),
                                  pBlock, poDstDS->nBlockXSize * poDstDS->nBlockYSize );
                          }
                      }
                  }
  
                  if ( poDstDS->bStats )
                      tiledb::Stats::enable();
  
                  auto status = query.submit();
  
                  if ( poDstDS->bStats )
                  {
                      tiledb::Stats::dump(stdout);
                      tiledb::Stats::disable();
                  }
  
                  if ( status == tiledb::Query::Status::FAILED )
                  {
                      return CE_Failure;
                  }
  
                  int nBlocks = ( (j + 1) * poDstDS->nBlocksX );
  
                  if( !pfnProgress( nBlocks /  static_cast<double>( nTotalBlocks),
                                      nullptr, pProgressData ) )
                  {
                      CPLError( CE_Failure, CPLE_UserInterrupt,
                              "User terminated CreateCopy()" );
                      return CE_Failure;
                  }
  
              }
          }
  
          query.finalize();
  
          return CE_None;
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          return CE_Failure;
      }
  }
  
  
  /************************************************************************/
  /*                              Create()                                */
  /************************************************************************/
  
  GDALDataset *
  TileDBDataset::Create( const char * pszFilename, int nXSize, int nYSize, int nBands,
          GDALDataType eType, char ** papszOptions )
  
  {
      try
      {
          CPLString osArrayPath = vsi_to_tiledb_uri( pszFilename );
  
          std::unique_ptr<TileDBDataset> poDS(
              TileDBDataset::CreateLL( osArrayPath, nXSize, nYSize,
                                      nBands, eType, papszOptions ));
  
          if( !poDS )
              return nullptr;
  
          tiledb::Array::create( osArrayPath, *poDS->m_schema );
  
          if ( poDS->nTimestamp )
              poDS->m_array.reset( new tiledb::Array( *poDS->m_ctx, osArrayPath,
                                      TILEDB_WRITE, poDS->nTimestamp ) );
          else
              poDS->m_array.reset( new tiledb::Array( *poDS->m_ctx,
                                   osArrayPath, TILEDB_WRITE ) );
  
          for( int i = 0; i < poDS->nBands;i++ )
          {
              if ( poDS->eIndexMode == ATTRIBUTES )
              {
                  poDS->SetBand( i+1, new TileDBRasterBand( poDS.get(), i+1,
                                 TILEDB_VALUES + CPLString().Printf( "_%i", i+1 ) ) );
              }
              else
                  poDS->SetBand( i+1, new TileDBRasterBand( poDS.get(), i+1 ) );
          }
  
          char** papszImageStruct = nullptr;
          papszImageStruct = CSLAddNameValue(papszImageStruct, "NBITS",
                          CPLString().Printf( "%d", poDS->nBitsPerSample ) );
          papszImageStruct = CSLAddNameValue(papszImageStruct, "DATA_TYPE",
                          CPLString().Printf( "%s", GDALGetDataTypeName( poDS->eDataType ) ) );
          papszImageStruct = CSLAddNameValue(papszImageStruct, "X_SIZE",
                          CPLString().Printf( "%d", poDS->nRasterXSize ) );
          papszImageStruct = CSLAddNameValue(papszImageStruct, "Y_SIZE",
                          CPLString().Printf( "%d", poDS->nRasterYSize ) );
          papszImageStruct = CSLAddNameValue(papszImageStruct, "INTERLEAVE",
                                              index_type_name( poDS->eIndexMode ) );
  
          if ( poDS->lpoAttributeDS.size() > 0 )
          {
              int i = 0;
              for ( auto const& poAttrDS: poDS->lpoAttributeDS )
              {
                  papszImageStruct = CSLAddNameValue( papszImageStruct,
                                  CPLString().Printf( "TILEDB_ATTRIBUTE_%i", ++i ),
                                  CPLGetBasename( poAttrDS->GetDescription() ) );
              }
          }
          poDS->SetMetadata(papszImageStruct, "IMAGE_STRUCTURE");
  
          CSLDestroy( papszImageStruct );
  
          return poDS.release();
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          return nullptr;
      }
  }
  
  /************************************************************************/
  /*                              CreateCopy()                            */
  /************************************************************************/
  
  GDALDataset *
  TileDBDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
                  int bStrict, char ** papszOptions,
                  GDALProgressFunc pfnProgress,
                  void *pProgressData )
  
  {
      char** papszCopyOptions = CSLDuplicate( papszOptions );
      CPLString osArrayPath = vsi_to_tiledb_uri( pszFilename );
  
      try
      {
          std::unique_ptr<TileDBDataset> poDstDS;
  
          if ( CSLFetchNameValue(papszOptions, "APPEND_SUBDATASET" ) )
          {
              // TileDB schemas are fixed
              CPLError(CE_Failure, CPLE_NotSupported,
                  "TileDB driver does not support "
                  "appending to an existing schema.");
              CSLDestroy( papszCopyOptions );
              return nullptr;
          }
  
          char** papszSrcSubDatasets = poSrcDS->GetMetadata( "SUBDATASETS" );
  
          if ( papszSrcSubDatasets == nullptr )
          {
              const int nBands = poSrcDS->GetRasterCount();
  
              if ( nBands > 0 )
              {
                  GDALRasterBand* poBand = poSrcDS->GetRasterBand( 1 );
                  GDALDataType eType = poBand->GetRasterDataType();
  
                  for (int i = 2; i <= nBands; ++i)
                  {
                      if (eType != poSrcDS->GetRasterBand( i )->GetRasterDataType( ) )
                      {
                          CPLError(CE_Failure, CPLE_NotSupported,
                              "TileDB driver does not support "
                              "source dataset with different band data types.");
                          CSLDestroy( papszCopyOptions );
                          return nullptr;
                      }
                  }
  
                  poDstDS.reset(
                      static_cast<TileDBDataset*>(TileDBDataset::Create( osArrayPath,
                              poSrcDS->GetRasterXSize(),
                              poSrcDS->GetRasterYSize(),
                              nBands, eType, papszOptions )));
  
                  if( !poDstDS )
                  {
                      CSLDestroy( papszCopyOptions );
                      return nullptr;
                  }
  
                  CPLErr eErr = GDALDatasetCopyWholeRaster( poSrcDS, poDstDS.get(),
                                                  papszOptions, pfnProgress,
                                                  pProgressData );
  
                  if ( eErr != CE_None )
                  {
                      CPLError(eErr, CPLE_AppDefined,
                              "Error copying raster to TileDB.");
                  }
              }
              else
              {
                  CPLError(CE_Failure, CPLE_NotSupported,
                      "TileDB driver does not support "
                      "source dataset with zero bands.");
              }
          }
          else
          {
              if ( bStrict )
              {
                  CPLError(CE_Failure, CPLE_NotSupported,
                      "TileDB driver does not support copying "
                      "subdatasets in strict mode.");
                  CSLDestroy( papszCopyOptions );
                  return nullptr;
              }
  
              if ( CSLFetchNameValue(papszOptions, "BLOCKXSIZE" ) ||
                  CSLFetchNameValue(papszOptions, "BLOCKYSIZE" ) )
              {
                  CPLError(CE_Failure, CPLE_NotSupported,
                      "Changing block size is not supported when copying subdatasets.");
                  CSLDestroy( papszCopyOptions );
                  return nullptr;
              }
  
              const int nSubDatasetCount = CSLCount( papszSrcSubDatasets ) / 2;
              const int nMaxFiles =
                  atoi(CPLGetConfigOption( "GDAL_READDIR_LIMIT_ON_OPEN", "1000" ) );
  
              if ( nSubDatasetCount <= nMaxFiles )
              {
                  const char* pszSource = CSLFetchNameValue(
                      papszSrcSubDatasets, "SUBDATASET_1_NAME");
                  if( pszSource )
                  {
                      std::unique_ptr<GDALDataset> poSubDataset(GDALDataset::Open( pszSource ));
                      if ( poSubDataset && poSubDataset->GetRasterCount() > 0 )
                      {
                          GDALRasterBand* poBand = poSubDataset->GetRasterBand( 1 );
  
                          TileDBDataset::SetBlockSize( poBand, papszCopyOptions );
                          poDstDS.reset(TileDBDataset::CreateLL(
                                      osArrayPath, poBand->GetXSize(),
                                      poBand->GetYSize(), 0, poBand->GetRasterDataType(),
                                      papszCopyOptions ));
  
                          if ( poDstDS &&
                              TileDBDataset::CopySubDatasets( poSrcDS, poDstDS.get(),
                                                  pfnProgress, pProgressData ) != CE_None )
                          {
                              poDstDS.reset();
                              CPLError( CE_Failure, CPLE_AppDefined,
                                  "Unable to copy subdatasets.");
                          }
                      }
                  }
              }
              else
              {
              CPLError( CE_Failure, CPLE_AppDefined,
                  "Please increase GDAL_READDIR_LIMIT_ON_OPEN variable.");
              }
          }
  
          // TODO Supporting mask bands is a possible future task
          if ( poDstDS != nullptr )
          {
              int nCloneFlags = GCIF_PAM_DEFAULT & ~GCIF_MASK;
              poDstDS->CloneInfo( poSrcDS, nCloneFlags );
          }
  
          CSLDestroy( papszCopyOptions );
  
          return poDstDS.release();
      }
      catch(const tiledb::TileDBError& e)
      {
          CPLError( CE_Failure, CPLE_AppDefined, "%s", e.what() );
          CSLDestroy( papszCopyOptions );
          return nullptr;
      }
  }
  
  /************************************************************************/
  /*                             IRasterIO()                              */
  /*                                                                      */
  /*      Multi-band raster io handler.  We will use block based          */
  /*      loading for attribute only arrays and for                       */
  /*      arrays that are indexed band last as they                       */
  /*      are effectively pixel interleaved, so processing all bands      */
  /*      for a given block together avoid extra seeks.                   */
  /************************************************************************/
  
  CPLErr TileDBDataset::IRasterIO( GDALRWFlag eRWFlag,
                                int nXOff, int nYOff, int nXSize, int nYSize,
                                void *pData, int nBufXSize, int nBufYSize,
                                GDALDataType eBufType,
                                int nBandCount, int *panBandMap,
                                GSpacing nPixelSpace, GSpacing nLineSpace,
                                GSpacing nBandSpace,
                                GDALRasterIOExtraArg* psExtraArg )
  
  {
      if ( ( eIndexMode == PIXEL ) ||
           ( eIndexMode == ATTRIBUTES ))
      {
          return GDALDataset::BlockBasedRasterIO(
              eRWFlag, nXOff, nYOff, nXSize, nYSize,
              pData, nBufXSize, nBufYSize, eBufType,
              nBandCount, panBandMap, nPixelSpace, nLineSpace, nBandSpace, psExtraArg );
      }
      else
          return
              GDALDataset::IRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize,
                                      pData, nBufXSize, nBufYSize, eBufType,
                                      nBandCount, panBandMap,
                                      nPixelSpace, nLineSpace, nBandSpace, psExtraArg );
  }
  
  /************************************************************************/
  /*                         GDALRegister_TILEDB()                        */
  /************************************************************************/
  
  void GDALRegister_TileDB()
  
  {
      if( GDALGetDriverByName( "TileDB" ) != nullptr )
          return;
  
      GDALDriver *poDriver = new GDALDriver();
  
      poDriver->SetDescription( "TileDB" );
      poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
      poDriver->SetMetadataItem( GDAL_DCAP_SUBCREATECOPY, "YES" );
      poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "TileDB" );
      poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/tiledb.html" );
      poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
                                 "Byte UInt16 Int16 UInt32 Int32 Float32 "
                                 "Float64 CInt16 CInt32 CFloat32 CFloat64" );
      poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
      poDriver->SetMetadataItem( GDAL_DMD_SUBDATASETS, "YES" );
      poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
  "<CreationOptionList>\n"
  "   <Option name='COMPRESSION' type='string-select' description='image compression to use' default='NONE'>\n"
  "       <Value>NONE</Value>\n"
  "       <Value>GZIP</Value>\n"
  "       <Value>ZSTD</Value>\n"
  "       <Value>LZ4</Value>\n"
  "       <Value>RLE</Value>\n"
  "       <Value>BZIP2</Value>\n"
  "       <Value>DOUBLE-DELTA</Value>\n"
  "       <Value>POSITIVE-DELTA</Value>\n"
  "   </Option>\n"
  "   <Option name='COMPRESSION_LEVEL' type='int' description='Compression level'/>\n"
  "   <Option name='BLOCKXSIZE' type='int' description='Tile Width'/>"
  "   <Option name='BLOCKYSIZE' type='int' description='Tile Height'/>"
  "   <Option name='STATS' type='boolean' description='Dump TileDB stats'/>"
  "   <Option name='TILEDB_CONFIG' type='string' description='location of configuration file for TileDB'/>"
  "   <Option name='TILEDB_ATTRIBUTE' type='string' description='co-registered file to add as TileDB attributes, only applicable for interleave types of band or pixel'/>"
  "   <Option name='INTERLEAVE' type='string-select' description='Indexing order' default='BAND'>\n"
  "        <Value>BAND</Value>\n"
  "        <Value>PIXEL</Value>\n"
  "        <Value>ATTRIBUTES</Value>\n"
  "   </Option>\n"
  "   <Option name='TILEDB_TIMESTAMP' type='int' description='Create array at this timestamp, the timestamp should be > 0'/>\n"
  "</CreationOptionList>\n" );
  
      poDriver->SetMetadataItem( GDAL_DMD_OPENOPTIONLIST,
  "<OpenOptionList>"
  "   <Option name='STATS' type='boolean' description='Dump TileDB stats'/>"
  "   <Option name='TILEDB_ATTRIBUTE' type='string' description='Attribute to read from each band'/>"
  "   <Option name='TILEDB_CONFIG' type='string' description='location of configuration file for TileDB'/>"
  "   <Option name='TILEDB_TIMESTAMP' type='int' description='Open array at this timestamp, the timestamp should be > 0'/>"
  "</OpenOptionList>" );
  
      poDriver->pfnIdentify = TileDBDataset::Identify;
      poDriver->pfnOpen = TileDBDataset::Open;
      poDriver->pfnCreate = TileDBDataset::Create;
      poDriver->pfnCreateCopy = TileDBDataset::CreateCopy;
      poDriver->pfnDelete = TileDBDataset::Delete;
  
      GetGDALDriverManager()->RegisterDriver( poDriver );
  }
  
  
  #ifdef _MSC_VER
  #pragma warning( pop )
  #endif
  