/******************************************************************************

 *

 * Project:  Earth Engine Data API Images driver

 * Purpose:  Earth Engine Data API Images driver

 * Author:   Even Rouault, even dot rouault at spatialys.com

 *

 ******************************************************************************

 * Copyright (c) 2017-2018, Planet Labs

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "gdal_priv.h"

#include "cpl_http.h"

#include "cpl_conv.h"

#include "ogrlibjsonutils.h"

#include "eeda.h"

#include <algorithm>

#include <vector>

#include <map>

#include <limits>

extern "C" void GDALRegister_EEDAI();

static const int DEFAULT_BLOCK_SIZE = 256;

const GUInt32 RETRY_PER_BAND = 1;

const GUInt32 RETRY_SPATIAL_SPLIT = 2;

// Eart engine server only allows up to 16 MB per request

const int SERVER_BYTE_LIMIT = 16 * 1024 * 1024;

const int SERVER_SIMUTANEOUS_BAND_LIMIT = 100;

const int SERVER_DIMENSION_LIMIT = 10000;

/************************************************************************/

/*                          GDALEEDAIDataset                            */

/************************************************************************/

class GDALEEDAIDataset final : public GDALEEDABaseDataset

{

    CPL_DISALLOW_COPY_ASSIGN(GDALEEDAIDataset)

    friend class GDALEEDAIRasterBand;

    int m_nBlockSize;

    CPLString m_osAsset{};

    CPLString m_osAssetName{};

    GDALEEDAIDataset *m_poParentDS;

#ifdef DEBUG_VERBOSE

    int m_iOvrLevel;

#endif

    CPLString m_osPixelEncoding{};

    bool m_bQueryMultipleBands;

    OGRSpatialReference m_oSRS{};

    double m_adfGeoTransform[6];

    std::vector<GDALEEDAIDataset *> m_apoOverviewDS{};

    GDALEEDAIDataset(GDALEEDAIDataset *poParentDS, int iOvrLevel);

    void

    SetMetadataFromProperties(json_object *poProperties,

                              const std::map<CPLString, int> &aoMapBandNames);

  public:

    GDALEEDAIDataset();

    virtual ~GDALEEDAIDataset();

    const OGRSpatialReference *GetSpatialRef() const override;

    virtual CPLErr GetGeoTransform(double *) override;

    virtual CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,

                             int nXSize, int nYSize, void *pData, int nBufXSize,

                             int nBufYSize, GDALDataType eBufType,

                             int nBandCount, BANDMAP_TYPE panBandMap,

                             GSpacing nPixelSpace, GSpacing nLineSpace,

                             GSpacing nBandSpace,

                             GDALRasterIOExtraArg *psExtraArg) override;

    bool ComputeQueryStrategy();

    bool Open(GDALOpenInfo *poOpenInfo);

};

/************************************************************************/

/*                        GDALEEDAIRasterBand                           */

/************************************************************************/

class GDALEEDAIRasterBand final : public GDALRasterBand

{

    CPL_DISALLOW_COPY_ASSIGN(GDALEEDAIRasterBand)

    friend class GDALEEDAIDataset;

    GDALColorInterp m_eInterp;

    bool DecodeNPYArray(const GByte *pabyData, int nDataLen,

                        bool bQueryAllBands, void *pDstBuffer, int nBlockXOff,

                        int nBlockYOff, int nXBlocks, int nYBlocks,

                        int nReqXSize, int nReqYSize) const;

    bool DecodeGDALDataset(const GByte *pabyData, int nDataLen,

                           bool bQueryAllBands, void *pDstBuffer,

                           int nBlockXOff, int nBlockYOff, int nXBlocks,

                           int nYBlocks, int nReqXSize, int nReqYSize);

    CPLErr GetBlocks(int nBlockXOff, int nBlockYOff, int nXBlocks, int nYBlocks,

                     bool bQueryAllBands, void *pBuffer);

    GUInt32 PrefetchBlocks(int nXOff, int nYOff, int nXSize, int nYSize,

                           int nBufXSize, int nBufYSize, bool bQueryAllBands);

  public:

    GDALEEDAIRasterBand(GDALEEDAIDataset *poDSIn, GDALDataType eDT);

    virtual ~GDALEEDAIRasterBand();

    virtual CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,

                             int nXSize, int nYSize, void *pData, int nBufXSize,

                             int nBufYSize, GDALDataType eBufType,

                             GSpacing nPixelSpace, GSpacing nLineSpace,

                             GDALRasterIOExtraArg *psExtraArg) CPL_OVERRIDE;

    virtual CPLErr IReadBlock(int, int, void *) CPL_OVERRIDE;

    virtual int GetOverviewCount() CPL_OVERRIDE;

    virtual GDALRasterBand *GetOverview(int) CPL_OVERRIDE;

    virtual CPLErr SetColorInterpretation(GDALColorInterp eInterp) CPL_OVERRIDE

    {

        m_eInterp = eInterp;

        return CE_None;

    }

    virtual GDALColorInterp GetColorInterpretation() CPL_OVERRIDE

    {

        return m_eInterp;

    }

};

/************************************************************************/

/*                         GDALEEDAIDataset()                           */

/************************************************************************/

GDALEEDAIDataset::GDALEEDAIDataset()

    : m_nBlockSize(DEFAULT_BLOCK_SIZE), m_poParentDS(nullptr),

#ifdef DEBUG_VERBOSE

      m_iOvrLevel(0),

#endif

      m_bQueryMultipleBands(false)

{

    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

    m_adfGeoTransform[0] = 0.0;

    m_adfGeoTransform[1] = 1.0;

    m_adfGeoTransform[2] = 0.0;

    m_adfGeoTransform[3] = 0.0;

    m_adfGeoTransform[4] = 0.0;

    m_adfGeoTransform[5] = 1.0;

}

/************************************************************************/

/*                         GDALEEDAIDataset()                           */

/************************************************************************/

GDALEEDAIDataset::GDALEEDAIDataset(GDALEEDAIDataset *poParentDS, int iOvrLevel)

    : m_nBlockSize(poParentDS->m_nBlockSize), m_osAsset(poParentDS->m_osAsset),

      m_osAssetName(poParentDS->m_osAssetName), m_poParentDS(poParentDS),

#ifdef DEBUG_VERBOSE

      m_iOvrLevel(iOvrLevel),

#endif

      m_osPixelEncoding(poParentDS->m_osPixelEncoding),

      m_bQueryMultipleBands(poParentDS->m_bQueryMultipleBands),

      m_oSRS(poParentDS->m_oSRS)

{

    m_osBaseURL = poParentDS->m_osBaseURL;

    nRasterXSize = m_poParentDS->nRasterXSize >> iOvrLevel;

    nRasterYSize = m_poParentDS->nRasterYSize >> iOvrLevel;

    m_adfGeoTransform[0] = m_poParentDS->m_adfGeoTransform[0];

    m_adfGeoTransform[1] = m_poParentDS->m_adfGeoTransform[1] *

                           m_poParentDS->nRasterXSize / nRasterXSize;

    m_adfGeoTransform[2] = m_poParentDS->m_adfGeoTransform[2];

    m_adfGeoTransform[3] = m_poParentDS->m_adfGeoTransform[3];

    m_adfGeoTransform[4] = m_poParentDS->m_adfGeoTransform[4];

    m_adfGeoTransform[5] = m_poParentDS->m_adfGeoTransform[5] *

                           m_poParentDS->nRasterYSize / nRasterYSize;

}

/************************************************************************/

/*                        ~GDALEEDAIDataset()                           */

/************************************************************************/

GDALEEDAIDataset::~GDALEEDAIDataset()

{

    for (size_t i = 0; i < m_apoOverviewDS.size(); i++)

    {

        delete m_apoOverviewDS[i];

    }

}

/************************************************************************/

/*                        GDALEEDAIRasterBand()                         */

/************************************************************************/

GDALEEDAIRasterBand::GDALEEDAIRasterBand(GDALEEDAIDataset *poDSIn,

                                         GDALDataType eDT)

    : m_eInterp(GCI_Undefined)

{

    eDataType = eDT;

    nBlockXSize = poDSIn->m_nBlockSize;

    nBlockYSize = poDSIn->m_nBlockSize;

}

/************************************************************************/

/*                       ~GDALEEDAIRasterBand()                         */

/************************************************************************/

GDALEEDAIRasterBand::~GDALEEDAIRasterBand()

{

}

/************************************************************************/

/*                           GetOverviewCount()                         */

/************************************************************************/

int GDALEEDAIRasterBand::GetOverviewCount()

{

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    return static_cast<int>(poGDS->m_apoOverviewDS.size());

}

/************************************************************************/

/*                              GetOverview()                           */

/************************************************************************/

GDALRasterBand *GDALEEDAIRasterBand::GetOverview(int iIndex)

{

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    if (iIndex >= 0 && iIndex < static_cast<int>(poGDS->m_apoOverviewDS.size()))

    {

        return poGDS->m_apoOverviewDS[iIndex]->GetRasterBand(nBand);

    }

    return nullptr;

}

/************************************************************************/

/*                            DecodeNPYArray()                          */

/************************************************************************/

bool GDALEEDAIRasterBand::DecodeNPYArray(const GByte *pabyData, int nDataLen,

                                         bool bQueryAllBands, void *pDstBuffer,

                                         int nBlockXOff, int nBlockYOff,

                                         int nXBlocks, int nYBlocks,

                                         int nReqXSize, int nReqYSize) const

{

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    // See https://docs.scipy.org/doc/numpy-1.13.0/neps/npy-format.html

    // for description of NPY array serialization format

    if (nDataLen < 10)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "Non NPY array returned");

        return false;

    }

    if (memcmp(pabyData, "\x93NUMPY", 6) != 0)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "Non NPY array returned");

        return false;

    }

    const int nVersionMajor = pabyData[6];

    if (nVersionMajor != 1)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Only version 1 of NPY array supported. Here found %d",

                 nVersionMajor);

        return false;

    }

    // Ignore version minor

    const int nHeaderLen = pabyData[8] | (pabyData[9] << 8);

    if (nDataLen < 10 + nHeaderLen)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Corrupted NPY array returned: not enough bytes for header");

        return false;

    }

#ifdef DEBUG

    CPLString osDescr;

    osDescr.assign(reinterpret_cast<const char *>(pabyData) + 10, nHeaderLen);

    // Should be something like

    // {'descr': [('B2', '<u2'), ('B3', '<u2'), ('B4', '<u2'), ('B8', '<u2'),

    // ('QA10', '<u2')], 'fortran_order': False, 'shape': (256, 256), }

    CPLDebug("EEDAI", "NPY descr: %s", osDescr.c_str());

    // TODO: validate that the descr is the one expected

#endif

    int nTotalDataTypeSize = 0;

    for (int i = 1; i <= poGDS->GetRasterCount(); i++)

    {

        if (bQueryAllBands || i == nBand)

        {

            nTotalDataTypeSize += GDALGetDataTypeSizeBytes(

                poGDS->GetRasterBand(i)->GetRasterDataType());

        }

    }

    int nDataSize = nTotalDataTypeSize * nReqXSize * nReqYSize;

    if (nDataLen < 10 + nHeaderLen + nDataSize)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Corrupted NPY array returned: not enough bytes for payload. "

                 "%d needed, only %d found",

                 10 + nHeaderLen + nDataSize, nDataLen);

        return false;

    }

    else if (nDataLen > 10 + nHeaderLen + nDataSize)

    {

        CPLError(CE_Warning, CPLE_AppDefined,

                 "Possibly corrupted NPY array returned: "

                 "expected bytes for payload. "

                 "%d needed, got %d found",

                 10 + nHeaderLen + nDataSize, nDataLen);

    }

    for (int iYBlock = 0; iYBlock < nYBlocks; iYBlock++)

    {

        int nBlockActualYSize = nBlockYSize;

        if ((iYBlock + nBlockYOff + 1) * nBlockYSize > nRasterYSize)

        {

            nBlockActualYSize =

                nRasterYSize - (iYBlock + nBlockYOff) * nBlockYSize;

        }

        for (int iXBlock = 0; iXBlock < nXBlocks; iXBlock++)

        {

            int nBlockActualXSize = nBlockXSize;

            if ((iXBlock + nBlockXOff + 1) * nBlockXSize > nRasterXSize)

            {

                nBlockActualXSize =

                    nRasterXSize - (iXBlock + nBlockXOff) * nBlockXSize;

            }

            int nOffsetBand =

                10 + nHeaderLen +

                (iYBlock * nBlockYSize * nReqXSize + iXBlock * nBlockXSize) *

                    nTotalDataTypeSize;

            for (int i = 1; i <= poGDS->GetRasterCount(); i++)

            {

                GDALRasterBlock *poBlock = nullptr;

                GByte *pabyDstBuffer;

                if (i == nBand && pDstBuffer != nullptr)

                    pabyDstBuffer = reinterpret_cast<GByte *>(pDstBuffer);

                else if (bQueryAllBands ||

                         (i == nBand && pDstBuffer == nullptr))

                {

                    GDALEEDAIRasterBand *poOtherBand =

                        reinterpret_cast<GDALEEDAIRasterBand *>(

                            poGDS->GetRasterBand(i));

                    poBlock = poOtherBand->TryGetLockedBlockRef(

                        nBlockXOff + iXBlock, nBlockYOff + iYBlock);

                    if (poBlock != nullptr)

                    {

                        poBlock->DropLock();

                        continue;

                    }

                    poBlock = poOtherBand->GetLockedBlockRef(

                        nBlockXOff + iXBlock, nBlockYOff + iYBlock, TRUE);

                    if (poBlock == nullptr)

                    {

                        continue;

                    }

                    pabyDstBuffer =

                        reinterpret_cast<GByte *>(poBlock->GetDataRef());

                }

                else

                {

                    continue;

                }

                GDALDataType eDT = poGDS->GetRasterBand(i)->GetRasterDataType();

                const int nDTSize = GDALGetDataTypeSizeBytes(eDT);

                for (int iLine = 0; iLine < nBlockActualYSize; iLine++)

                {

                    GByte *pabyLineDest =

                        pabyDstBuffer + iLine * nDTSize * nBlockXSize;

                    GDALCopyWords(const_cast<GByte *>(pabyData) + nOffsetBand +

                                      iLine * nTotalDataTypeSize * nReqXSize,

                                  eDT, nTotalDataTypeSize, pabyLineDest, eDT,

                                  nDTSize, nBlockActualXSize);

#ifdef CPL_MSB

                    if (nDTSize > 1)

                    {

                        GDALSwapWords(pabyLineDest, nDTSize, nBlockActualXSize,

                                      nDTSize);

                    }

#endif

                }

                nOffsetBand += nDTSize;

                if (poBlock)

                    poBlock->DropLock();

            }

        }

    }

    return true;

}

/************************************************************************/

/*                            DecodeGDALDataset()                         */

/************************************************************************/

bool GDALEEDAIRasterBand::DecodeGDALDataset(const GByte *pabyData, int nDataLen,

                                            bool bQueryAllBands,

                                            void *pDstBuffer, int nBlockXOff,

                                            int nBlockYOff, int nXBlocks,

                                            int nYBlocks, int nReqXSize,

                                            int nReqYSize)

{

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    const CPLString osTmpFilename(VSIMemGenerateHiddenFilename("eedai"));

    VSIFCloseL(VSIFileFromMemBuffer(

        osTmpFilename, const_cast<GByte *>(pabyData), nDataLen, false));

    const char *const apszDrivers[] = {"PNG", "JPEG", "GTIFF", nullptr};

    GDALDataset *poTileDS = GDALDataset::FromHandle(GDALOpenEx(

        osTmpFilename, GDAL_OF_RASTER, apszDrivers, nullptr, nullptr));

    if (poTileDS == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Cannot decode buffer returned by the "

                 "server as a PNG, JPEG or GeoTIFF image");

        VSIUnlink(osTmpFilename);

        return false;

    }

    if (poTileDS->GetRasterXSize() != nReqXSize ||

        poTileDS->GetRasterYSize() != nReqYSize ||

        // The server might return a RGBA image even if only 3 bands are

        // requested

        poTileDS->GetRasterCount() <

            (bQueryAllBands ? poGDS->GetRasterCount() : 1))

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Bad dimensions/band count for image returned "

                 "by server: %dx%dx%d",

                 poTileDS->GetRasterXSize(), poTileDS->GetRasterYSize(),

                 poTileDS->GetRasterCount());

        delete poTileDS;

        VSIUnlink(osTmpFilename);

        return false;

    }

    for (int iYBlock = 0; iYBlock < nYBlocks; iYBlock++)

    {

        int nBlockActualYSize = nBlockYSize;

        if ((iYBlock + nBlockYOff + 1) * nBlockYSize > nRasterYSize)

        {

            nBlockActualYSize =

                nRasterYSize - (iYBlock + nBlockYOff) * nBlockYSize;

        }

        for (int iXBlock = 0; iXBlock < nXBlocks; iXBlock++)

        {

            int nBlockActualXSize = nBlockXSize;

            if ((iXBlock + nBlockXOff + 1) * nBlockXSize > nRasterXSize)

            {

                nBlockActualXSize =

                    nRasterXSize - (iXBlock + nBlockXOff) * nBlockXSize;

            }

            for (int i = 1; i <= poGDS->GetRasterCount(); i++)

            {

                GDALRasterBlock *poBlock = nullptr;

                GByte *pabyDstBuffer;

                if (i == nBand && pDstBuffer != nullptr)

                    pabyDstBuffer = reinterpret_cast<GByte *>(pDstBuffer);

                else if (bQueryAllBands ||

                         (i == nBand && pDstBuffer == nullptr))

                {

                    GDALEEDAIRasterBand *poOtherBand =

                        reinterpret_cast<GDALEEDAIRasterBand *>(

                            poGDS->GetRasterBand(i));

                    poBlock = poOtherBand->TryGetLockedBlockRef(

                        nBlockXOff + iXBlock, nBlockYOff + iYBlock);

                    if (poBlock != nullptr)

                    {

                        poBlock->DropLock();

                        continue;

                    }

                    poBlock = poOtherBand->GetLockedBlockRef(

                        nBlockXOff + iXBlock, nBlockYOff + iYBlock, TRUE);

                    if (poBlock == nullptr)

                    {

                        continue;

                    }

                    pabyDstBuffer =

                        reinterpret_cast<GByte *>(poBlock->GetDataRef());

                }

                else

                {

                    continue;

                }

                GDALDataType eDT = poGDS->GetRasterBand(i)->GetRasterDataType();

                const int nDTSize = GDALGetDataTypeSizeBytes(eDT);

                const int nTileBand = bQueryAllBands ? i : 1;

                CPLErr eErr = poTileDS->GetRasterBand(nTileBand)->RasterIO(

                    GF_Read, iXBlock * nBlockXSize, iYBlock * nBlockYSize,

                    nBlockActualXSize, nBlockActualYSize, pabyDstBuffer,

                    nBlockActualXSize, nBlockActualYSize, eDT, nDTSize,

                    static_cast<GSpacing>(nDTSize) * nBlockXSize, nullptr);

                if (poBlock)

                    poBlock->DropLock();

                if (eErr != CE_None)

                {

                    delete poTileDS;

                    VSIUnlink(osTmpFilename);

                    return false;

                }

            }

        }

    }

    delete poTileDS;

    VSIUnlink(osTmpFilename);

    return true;

}

CPLErr GDALEEDAIRasterBand::GetBlocks(int nBlockXOff, int nBlockYOff,

                                      int nXBlocks, int nYBlocks,

                                      bool bQueryAllBands, void *pBuffer)

{

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    // Build request content

    json_object *poReq = json_object_new_object();

    json_object_object_add(poReq, "fileFormat",

                           json_object_new_string(poGDS->m_osPixelEncoding));

    json_object *poBands = json_object_new_array();

    for (int i = 1; i <= poGDS->GetRasterCount(); i++)

    {

        if (bQueryAllBands || i == nBand)

        {

            json_object_array_add(

                poBands, json_object_new_string(

                             poGDS->GetRasterBand(i)->GetDescription()));

        }

    }

    json_object_object_add(poReq, "bandIds", poBands);

    int nReqXSize = nBlockXSize * nXBlocks;

    if ((nBlockXOff + nXBlocks) * nBlockXSize > nRasterXSize)

        nReqXSize = nRasterXSize - nBlockXOff * nBlockXSize;

    int nReqYSize = nBlockYSize * nYBlocks;

    if ((nBlockYOff + nYBlocks) * nBlockYSize > nRasterYSize)

        nReqYSize = nRasterYSize - nBlockYOff * nBlockYSize;

    const double dfX0 = poGDS->m_adfGeoTransform[0] +

                        nBlockXOff * nBlockXSize * poGDS->m_adfGeoTransform[1];

    const double dfY0 = poGDS->m_adfGeoTransform[3] +

                        nBlockYOff * nBlockYSize * poGDS->m_adfGeoTransform[5];

#ifdef DEBUG_VERBOSE

    CPLDebug("EEDAI",

             "nBlockYOff=%d nBlockYOff=%d "

             "nXBlocks=%d nYBlocks=%d nReqXSize=%d nReqYSize=%d",

             nBlockYOff, nBlockYOff, nXBlocks, nYBlocks, nReqXSize, nReqYSize);

#endif

    json_object *poPixelGrid = json_object_new_object();

    json_object *poAffineTransform = json_object_new_object();

    json_object_object_add(

        poAffineTransform, "translateX",

        json_object_new_double_with_significant_figures(dfX0, 18));

    json_object_object_add(

        poAffineTransform, "translateY",

        json_object_new_double_with_significant_figures(dfY0, 18));

    json_object_object_add(poAffineTransform, "scaleX",

                           json_object_new_double_with_significant_figures(

                               poGDS->m_adfGeoTransform[1], 18));

    json_object_object_add(poAffineTransform, "scaleY",

                           json_object_new_double_with_significant_figures(

                               poGDS->m_adfGeoTransform[5], 18));

    json_object_object_add(

        poAffineTransform, "shearX",

        json_object_new_double_with_significant_figures(0.0, 18));

    json_object_object_add(

        poAffineTransform, "shearY",

        json_object_new_double_with_significant_figures(0.0, 18));

    json_object_object_add(poPixelGrid, "affineTransform", poAffineTransform);

    json_object *poDimensions = json_object_new_object();

    json_object_object_add(poDimensions, "width",

                           json_object_new_int(nReqXSize));

    json_object_object_add(poDimensions, "height",

                           json_object_new_int(nReqYSize));

    json_object_object_add(poPixelGrid, "dimensions", poDimensions);

    json_object_object_add(poReq, "grid", poPixelGrid);

    CPLString osPostContent = json_object_get_string(poReq);

    json_object_put(poReq);

    // Issue request

    char **papszOptions = (poGDS->m_poParentDS)

                              ? poGDS->m_poParentDS->GetBaseHTTPOptions()

                              : poGDS->GetBaseHTTPOptions();

    papszOptions = CSLSetNameValue(papszOptions, "CUSTOMREQUEST", "POST");

    CPLString osHeaders = CSLFetchNameValueDef(papszOptions, "HEADERS", "");

    if (!osHeaders.empty())

        osHeaders += "\r\n";

    osHeaders += "Content-Type: application/json";

    papszOptions = CSLSetNameValue(papszOptions, "HEADERS", osHeaders);

    papszOptions = CSLSetNameValue(papszOptions, "POSTFIELDS", osPostContent);

    CPLHTTPResult *psResult = EEDAHTTPFetch(

        (poGDS->m_osBaseURL + poGDS->m_osAssetName + ":getPixels").c_str(),

        papszOptions);

    CSLDestroy(papszOptions);

    if (psResult == nullptr)

        return CE_Failure;

    if (psResult->pszErrBuf != nullptr)

    {

        if (psResult->pabyData)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "%s: %s", psResult->pszErrBuf,

                     reinterpret_cast<const char *>(psResult->pabyData));

        }

        else

        {

            CPLError(CE_Failure, CPLE_AppDefined, "%s", psResult->pszErrBuf);

        }

        CPLHTTPDestroyResult(psResult);

        return CE_Failure;

    }

    if (psResult->pabyData == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Empty content returned by server");

        CPLHTTPDestroyResult(psResult);

        return CE_Failure;

    }

#ifdef DEBUG_VERBOSE

    CPLDebug("EEADI", "Result: %s (%d bytes)",

             reinterpret_cast<const char *>(psResult->pabyData),

             psResult->nDataLen);

#endif

    if (EQUAL(poGDS->m_osPixelEncoding, "NPY"))

    {

        if (!DecodeNPYArray(psResult->pabyData, psResult->nDataLen,

                            bQueryAllBands, pBuffer, nBlockXOff, nBlockYOff,

                            nXBlocks, nYBlocks, nReqXSize, nReqYSize))

        {

            CPLHTTPDestroyResult(psResult);

            return CE_Failure;

        }

    }

    else

    {

        if (!DecodeGDALDataset(psResult->pabyData, psResult->nDataLen,

                               bQueryAllBands, pBuffer, nBlockXOff, nBlockYOff,

                               nXBlocks, nYBlocks, nReqXSize, nReqYSize))

        {

            CPLHTTPDestroyResult(psResult);

            return CE_Failure;

        }

    }

    CPLHTTPDestroyResult(psResult);

    return CE_None;

}

/************************************************************************/

/*                             IReadBlock()                             */

/************************************************************************/

CPLErr GDALEEDAIRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,

                                       void *pBuffer)

{

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

#ifdef DEBUG_VERBOSE

    CPLDebug("EEDAI", "ReadBlock x=%d y=%d band=%d level=%d", nBlockXOff,

             nBlockYOff, nBand, poGDS->m_iOvrLevel);

#endif

    return GetBlocks(nBlockXOff, nBlockYOff, 1, 1, poGDS->m_bQueryMultipleBands,

                     pBuffer);

}

/************************************************************************/

/*                          PrefetchBlocks()                            */

/************************************************************************/

// Return or'ed flags among 0, RETRY_PER_BAND, RETRY_SPATIAL_SPLIT if the user

// should try to split the request in smaller chunks

GUInt32 GDALEEDAIRasterBand::PrefetchBlocks(int nXOff, int nYOff, int nXSize,

                                            int nYSize, int nBufXSize,

                                            int nBufYSize, bool bQueryAllBands)

{

    CPL_IGNORE_RET_VAL(nBufXSize);

    CPL_IGNORE_RET_VAL(nBufYSize);

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    int nBlockXOff = nXOff / nBlockXSize;

    int nBlockYOff = nYOff / nBlockYSize;

    int nXBlocks = (nXOff + nXSize - 1) / nBlockXSize - nBlockXOff + 1;

    int nYBlocks = (nYOff + nYSize - 1) / nBlockYSize - nBlockYOff + 1;

    const int nThisDTSize = GDALGetDataTypeSizeBytes(GetRasterDataType());

    int nTotalDataTypeSize = 0;

    int nQueriedBands = 0;

    for (int i = 1; i <= poGDS->GetRasterCount(); i++)

    {

        if (bQueryAllBands || i == nBand)

        {

            nQueriedBands++;

            nTotalDataTypeSize += GDALGetDataTypeSizeBytes(

                poGDS->GetRasterBand(i)->GetRasterDataType());

        }

    }

    // Check the number of already cached blocks, and remove fully

    // cached lines at the top of the area of interest from the queried

    // blocks

    int nBlocksCached = 0;

    int nBlocksCachedForThisBand = 0;

    bool bAllLineCached = true;

    for (int iYBlock = 0; iYBlock < nYBlocks;)

    {

        for (int iXBlock = 0; iXBlock < nXBlocks; iXBlock++)

        {

            for (int i = 1; i <= poGDS->GetRasterCount(); i++)

            {

                GDALRasterBlock *poBlock = nullptr;

                if (bQueryAllBands || i == nBand)

                {

                    GDALEEDAIRasterBand *poOtherBand =

                        reinterpret_cast<GDALEEDAIRasterBand *>(

                            poGDS->GetRasterBand(i));

                    poBlock = poOtherBand->TryGetLockedBlockRef(

                        nBlockXOff + iXBlock, nBlockYOff + iYBlock);

                    if (poBlock != nullptr)

                    {

                        nBlocksCached++;

                        if (i == nBand)

                            nBlocksCachedForThisBand++;

                        poBlock->DropLock();

                        continue;

                    }

                    else

                    {

                        bAllLineCached = false;

                    }

                }

            }

        }

        if (bAllLineCached)

        {

            nBlocksCached -= nXBlocks * nQueriedBands;

            nBlocksCachedForThisBand -= nXBlocks;

            nBlockYOff++;

            nYBlocks--;

        }

        else

        {

            iYBlock++;

        }

    }

    if (nXBlocks > 0 && nYBlocks > 0)

    {

        bool bMustReturn = false;

        GUInt32 nRetryFlags = 0;

        // Get the blocks if the number of already cached blocks is lesser

        // than 25% of the to be queried blocks

        if (nBlocksCached > (nQueriedBands * nXBlocks * nYBlocks) / 4)

        {

            if (nBlocksCachedForThisBand <= (nXBlocks * nYBlocks) / 4)

            {

                nRetryFlags |= RETRY_PER_BAND;

            }

            else

            {

                bMustReturn = true;

            }

        }

        // Don't request too many pixels in one dimension

        if (nXBlocks * nBlockXSize > SERVER_DIMENSION_LIMIT ||

            nYBlocks * nBlockYSize > SERVER_DIMENSION_LIMIT)

        {

            bMustReturn = true;

            nRetryFlags |= RETRY_SPATIAL_SPLIT;

        }

        // Make sure that we have enough cache (with a margin of 50%)

        // and the number of queried pixels isn't too big w.r.t server

        // limit

        const GIntBig nUncompressedSize = static_cast<GIntBig>(nXBlocks) *

                                          nYBlocks * nBlockXSize * nBlockYSize *

                                          nTotalDataTypeSize;

        const GIntBig nCacheMax = GDALGetCacheMax64() / 2;

        if (nUncompressedSize > nCacheMax ||

            nUncompressedSize > SERVER_BYTE_LIMIT)

        {

            if (bQueryAllBands && poGDS->GetRasterCount() > 1)

            {

                const GIntBig nUncompressedSizeThisBand =

                    static_cast<GIntBig>(nXBlocks) * nYBlocks * nBlockXSize *

                    nBlockYSize * nThisDTSize;

                if (nUncompressedSizeThisBand <= SERVER_BYTE_LIMIT &&

                    nUncompressedSizeThisBand <= nCacheMax)

                {

                    nRetryFlags |= RETRY_PER_BAND;

                }

            }

            if (nXBlocks > 1 || nYBlocks > 1)

            {

                nRetryFlags |= RETRY_SPATIAL_SPLIT;

            }

            return nRetryFlags;

        }

        if (bMustReturn)

            return nRetryFlags;

        GetBlocks(nBlockXOff, nBlockYOff, nXBlocks, nYBlocks, bQueryAllBands,

                  nullptr);

    }

    return 0;

}

/************************************************************************/

/*                              IRasterIO()                             */

/************************************************************************/

CPLErr GDALEEDAIRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,

                                      int nXSize, int nYSize, void *pData,

                                      int nBufXSize, int nBufYSize,

                                      GDALDataType eBufType,

                                      GSpacing nPixelSpace, GSpacing nLineSpace,

                                      GDALRasterIOExtraArg *psExtraArg)

{

    /* ==================================================================== */

    /*      Do we have overviews that would be appropriate to satisfy       */

    /*      this request?                                                   */

    /* ==================================================================== */

    if ((nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0 &&

        eRWFlag == GF_Read)

    {

        GDALRasterIOExtraArg sExtraArg;

        GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);

        const int nOverview =

            GDALBandGetBestOverviewLevel2(this, nXOff, nYOff, nXSize, nYSize,

                                          nBufXSize, nBufYSize, &sExtraArg);

        if (nOverview >= 0)

        {

            GDALRasterBand *poOverviewBand = GetOverview(nOverview);

            if (poOverviewBand == nullptr)

                return CE_Failure;

            return poOverviewBand->RasterIO(

                eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,

                nBufYSize, eBufType, nPixelSpace, nLineSpace, &sExtraArg);

        }

    }

    GDALEEDAIDataset *poGDS = reinterpret_cast<GDALEEDAIDataset *>(poDS);

    GUInt32 nRetryFlags =

        PrefetchBlocks(nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,

                       poGDS->m_bQueryMultipleBands);

    if ((nRetryFlags & RETRY_SPATIAL_SPLIT) && nXSize == nBufXSize &&

        nYSize == nBufYSize && nYSize > nBlockYSize)

    {

        GDALRasterIOExtraArg sExtraArg;

        INIT_RASTERIO_EXTRA_ARG(sExtraArg);

        int nHalf =

            std::max(nBlockYSize, ((nYSize / 2) / nBlockYSize) * nBlockYSize);

        CPLErr eErr =

            IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nHalf, pData, nXSize,

                      nHalf, eBufType, nPixelSpace, nLineSpace, &sExtraArg);

        if (eErr == CE_None)

        {

            eErr = IRasterIO(

                eRWFlag, nXOff, nYOff + nHalf, nXSize, nYSize - nHalf,

                static_cast<GByte *>(pData) + nHalf * nLineSpace, nXSize,

                nYSize - nHalf, eBufType, nPixelSpace, nLineSpace, &sExtraArg);

        }

        return eErr;

    }

    else if ((nRetryFlags & RETRY_SPATIAL_SPLIT) && nXSize == nBufXSize &&

             nYSize == nBufYSize && nXSize > nBlockXSize)

    {

        GDALRasterIOExtraArg sExtraArg;

        INIT_RASTERIO_EXTRA_ARG(sExtraArg);

        int nHalf =

            std::max(nBlockXSize, ((nXSize / 2) / nBlockXSize) * nBlockXSize);

        CPLErr eErr =

            IRasterIO(eRWFlag, nXOff, nYOff, nHalf, nYSize, pData, nHalf,

                      nYSize, eBufType, nPixelSpace, nLineSpace, &sExtraArg);

        if (eErr == CE_None)

        {

            eErr =

                IRasterIO(eRWFlag, nXOff + nHalf, nYOff, nXSize - nHalf, nYSize,

                          static_cast<GByte *>(pData) + nHalf * nPixelSpace,