1 /******************************************************************************
2  *
3  * Project:  Earth Engine Data API Images driver
4  * Purpose:  Earth Engine Data API Images driver
5  * Author:   Even Rouault, even dot rouault at spatialys.com
6  *
7  ******************************************************************************
8  * Copyright (c) 2017-2018, Planet Labs
9  *
10  * Permission is hereby granted, free of charge, to any person obtaining a
11  * copy of this software and associated documentation files (the "Software"),
12  * to deal in the Software without restriction, including without limitation
13  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  * and/or sell copies of the Software, and to permit persons to whom the
15  * Software is furnished to do so, subject to the following conditions:
16  *
17  * The above copyright notice and this permission notice shall be included
18  * in all copies or substantial portions of the Software.
19  *
20  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  * DEALINGS IN THE SOFTWARE.
27  ****************************************************************************/
28 
29 #include "gdal_priv.h"
30 #include "cpl_http.h"
31 #include "cpl_conv.h"
32 #include "ogrgeojsonreader.h"
33 #include "ogrgeojsonwriter.h"
34 #include "eeda.h"
35 
36 #include <algorithm>
37 #include <vector>
38 #include <map>
39 #include <limits>
40 
41 extern "C" void GDALRegister_EEDAI();
42 
43 static const int DEFAULT_BLOCK_SIZE = 256;
44 
45 const GUInt32 RETRY_PER_BAND = 1;
46 const GUInt32 RETRY_SPATIAL_SPLIT = 2;
47 
48 // Eart engine server only allows up to 16 MB per request
49 const int SERVER_BYTE_LIMIT = 16 * 1024 * 1024;
50 const int SERVER_SIMUTANEOUS_BAND_LIMIT = 100;
51 const int SERVER_DIMENSION_LIMIT = 10000;
52 
53 /************************************************************************/
54 /*                          GDALEEDAIDataset                            */
55 /************************************************************************/
56 
57 class GDALEEDAIDataset final: public GDALEEDABaseDataset
58 {
59             CPL_DISALLOW_COPY_ASSIGN(GDALEEDAIDataset)
60 
61             friend class GDALEEDAIRasterBand;
62 
63             int         m_nBlockSize;
64             CPLString   m_osAsset{};
65             CPLString   m_osAssetName{};
66             GDALEEDAIDataset* m_poParentDS;
67 #ifdef DEBUG_VERBOSE
68             int         m_iOvrLevel;
69 #endif
70             CPLString   m_osPixelEncoding{};
71             bool        m_bQueryMultipleBands;
72             CPLString   m_osWKT{};
73             double      m_adfGeoTransform[6];
74             std::vector<GDALEEDAIDataset*> m_apoOverviewDS{};
75 
76                         GDALEEDAIDataset(GDALEEDAIDataset* poParentDS,
77                                          int iOvrLevel);
78 
79             void        SetMetadataFromProperties(
80                             json_object* poProperties,
81                             const std::map<CPLString, int>& aoMapBandNames );
82     public:
83                 GDALEEDAIDataset();
84                 virtual ~GDALEEDAIDataset();
85 
86                 virtual const char* _GetProjectionRef() override;
GetSpatialRef() const87                 const OGRSpatialReference* GetSpatialRef() const override {
88                     return GetSpatialRefFromOldGetProjectionRef();
89                 }
90                 virtual CPLErr GetGeoTransform( double* ) override;
91 
92                 virtual CPLErr IRasterIO( GDALRWFlag eRWFlag,
93                                  int nXOff, int nYOff, int nXSize, int nYSize,
94                                  void * pData, int nBufXSize, int nBufYSize,
95                                  GDALDataType eBufType,
96                                  int nBandCount, int *panBandMap,
97                                  GSpacing nPixelSpace, GSpacing nLineSpace,
98                                  GSpacing nBandSpace,
99                                  GDALRasterIOExtraArg* psExtraArg ) override;
100 
101                 bool ComputeQueryStrategy();
102 
103                 bool Open(GDALOpenInfo* poOpenInfo);
104 };
105 
106 /************************************************************************/
107 /*                        GDALEEDAIRasterBand                           */
108 /************************************************************************/
109 
110 class GDALEEDAIRasterBand final: public GDALRasterBand
111 {
112                 CPL_DISALLOW_COPY_ASSIGN(GDALEEDAIRasterBand)
113 
114                 friend class GDALEEDAIDataset;
115 
116                 GDALColorInterp m_eInterp;
117 
118                 bool    DecodeNPYArray( const GByte* pabyData,
119                                         int nDataLen,
120                                         bool bQueryAllBands,
121                                         void* pDstBuffer,
122                                         int nBlockXOff, int nBlockYOff,
123                                         int nXBlocks, int nYBlocks,
124                                         int nReqXSize, int nReqYSize ) const;
125                 bool    DecodeGDALDataset( const GByte* pabyData,
126                                          int nDataLen,
127                                          bool bQueryAllBands,
128                                          void* pDstBuffer,
129                                          int nBlockXOff, int nBlockYOff,
130                                          int nXBlocks, int nYBlocks,
131                                          int nReqXSize, int nReqYSize );
132 
133                 CPLErr  GetBlocks(    int nBlockXOff, int nBlockYOff,
134                                       int nXBlocks, int nYBlocks,
135                                       bool bQueryAllBands,
136                                       void* pBuffer);
137                 GUInt32 PrefetchBlocks(int nXOff, int nYOff,
138                                        int nXSize, int nYSize,
139                                        int nBufXSize, int nBufYSize,
140                                        bool bQueryAllBands);
141 
142     public:
143                 GDALEEDAIRasterBand(GDALEEDAIDataset *poDSIn,
144                                     GDALDataType eDT,
145                                     bool bSignedByte);
146                 virtual ~GDALEEDAIRasterBand();
147 
148                 virtual CPLErr IRasterIO( GDALRWFlag eRWFlag,
149                                   int nXOff, int nYOff, int nXSize, int nYSize,
150                                   void * pData, int nBufXSize, int nBufYSize,
151                                   GDALDataType eBufType,
152                                   GSpacing nPixelSpace, GSpacing nLineSpace,
153                                   GDALRasterIOExtraArg* psExtraArg ) CPL_OVERRIDE;
154 
155                 virtual CPLErr IReadBlock( int, int, void * ) CPL_OVERRIDE;
156                 virtual int GetOverviewCount() CPL_OVERRIDE;
157                 virtual GDALRasterBand* GetOverview(int) CPL_OVERRIDE;
SetColorInterpretation(GDALColorInterp eInterp)158                 virtual CPLErr SetColorInterpretation(GDALColorInterp eInterp) CPL_OVERRIDE
159                                         { m_eInterp = eInterp; return CE_None;}
GetColorInterpretation()160                 virtual GDALColorInterp GetColorInterpretation() CPL_OVERRIDE
161                                         { return m_eInterp; }
162 };
163 
164 /************************************************************************/
165 /*                         GDALEEDAIDataset()                           */
166 /************************************************************************/
167 
GDALEEDAIDataset()168 GDALEEDAIDataset::GDALEEDAIDataset() :
169     m_nBlockSize(DEFAULT_BLOCK_SIZE),
170     m_poParentDS(nullptr),
171 #ifdef DEBUG_VERBOSE
172     m_iOvrLevel(0),
173 #endif
174     m_bQueryMultipleBands(false)
175 {
176     m_adfGeoTransform[0] = 0.0;
177     m_adfGeoTransform[1] = 1.0;
178     m_adfGeoTransform[2] = 0.0;
179     m_adfGeoTransform[3] = 0.0;
180     m_adfGeoTransform[4] = 0.0;
181     m_adfGeoTransform[5] = 1.0;
182 }
183 
184 /************************************************************************/
185 /*                         GDALEEDAIDataset()                           */
186 /************************************************************************/
187 
GDALEEDAIDataset(GDALEEDAIDataset * poParentDS,int iOvrLevel)188 GDALEEDAIDataset::GDALEEDAIDataset(GDALEEDAIDataset* poParentDS,
189                                    int iOvrLevel ) :
190     m_nBlockSize(poParentDS->m_nBlockSize),
191     m_osAsset(poParentDS->m_osAsset),
192     m_osAssetName(poParentDS->m_osAssetName),
193     m_poParentDS(poParentDS),
194 #ifdef DEBUG_VERBOSE
195     m_iOvrLevel(iOvrLevel),
196 #endif
197     m_osPixelEncoding(poParentDS->m_osPixelEncoding),
198     m_bQueryMultipleBands(poParentDS->m_bQueryMultipleBands),
199     m_osWKT(poParentDS->m_osWKT)
200 {
201     m_osBaseURL = poParentDS->m_osBaseURL;
202     nRasterXSize = m_poParentDS->nRasterXSize >> iOvrLevel;
203     nRasterYSize = m_poParentDS->nRasterYSize >> iOvrLevel;
204     m_adfGeoTransform[0] = m_poParentDS->m_adfGeoTransform[0];
205     m_adfGeoTransform[1] = m_poParentDS->m_adfGeoTransform[1] *
206                                     m_poParentDS->nRasterXSize / nRasterXSize;
207     m_adfGeoTransform[2] = m_poParentDS->m_adfGeoTransform[2];
208     m_adfGeoTransform[3] = m_poParentDS->m_adfGeoTransform[3];
209     m_adfGeoTransform[4] = m_poParentDS->m_adfGeoTransform[4];
210     m_adfGeoTransform[5] = m_poParentDS->m_adfGeoTransform[5] *
211                                     m_poParentDS->nRasterYSize / nRasterYSize;
212 }
213 
214 /************************************************************************/
215 /*                        ~GDALEEDAIDataset()                           */
216 /************************************************************************/
217 
~GDALEEDAIDataset()218 GDALEEDAIDataset::~GDALEEDAIDataset()
219 {
220     for(size_t i = 0; i < m_apoOverviewDS.size(); i++ )
221     {
222         delete m_apoOverviewDS[i];
223     }
224 }
225 
226 /************************************************************************/
227 /*                        GDALEEDAIRasterBand()                         */
228 /************************************************************************/
229 
GDALEEDAIRasterBand(GDALEEDAIDataset * poDSIn,GDALDataType eDT,bool bSignedByte)230 GDALEEDAIRasterBand::GDALEEDAIRasterBand(GDALEEDAIDataset* poDSIn,
231                                          GDALDataType eDT,
232                                          bool bSignedByte) :
233     m_eInterp(GCI_Undefined)
234 {
235     eDataType = eDT;
236     nBlockXSize = poDSIn->m_nBlockSize;
237     nBlockYSize = poDSIn->m_nBlockSize;
238     if( bSignedByte )
239     {
240         SetMetadataItem("PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE");
241     }
242 }
243 
244 /************************************************************************/
245 /*                       ~GDALEEDAIRasterBand()                         */
246 /************************************************************************/
247 
~GDALEEDAIRasterBand()248 GDALEEDAIRasterBand::~GDALEEDAIRasterBand()
249 {
250 }
251 
252 /************************************************************************/
253 /*                           GetOverviewCount()                         */
254 /************************************************************************/
255 
GetOverviewCount()256 int GDALEEDAIRasterBand::GetOverviewCount()
257 {
258     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
259     return static_cast<int>(poGDS->m_apoOverviewDS.size());
260 }
261 
262 /************************************************************************/
263 /*                              GetOverview()                           */
264 /************************************************************************/
265 
GetOverview(int iIndex)266 GDALRasterBand* GDALEEDAIRasterBand::GetOverview(int iIndex)
267 {
268     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
269     if( iIndex >= 0 &&
270         iIndex < static_cast<int>(poGDS->m_apoOverviewDS.size()) )
271     {
272         return poGDS->m_apoOverviewDS[iIndex]->GetRasterBand(nBand);
273     }
274     return nullptr;
275 }
276 
277 
278 /************************************************************************/
279 /*                            DecodeNPYArray()                          */
280 /************************************************************************/
281 
DecodeNPYArray(const GByte * pabyData,int nDataLen,bool bQueryAllBands,void * pDstBuffer,int nBlockXOff,int nBlockYOff,int nXBlocks,int nYBlocks,int nReqXSize,int nReqYSize) const282 bool GDALEEDAIRasterBand::DecodeNPYArray( const GByte* pabyData,
283                                           int nDataLen,
284                                           bool bQueryAllBands,
285                                           void* pDstBuffer,
286                                           int nBlockXOff, int nBlockYOff,
287                                           int nXBlocks, int nYBlocks,
288                                           int nReqXSize, int nReqYSize ) const
289 {
290     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
291 
292     // See https://docs.scipy.org/doc/numpy-1.13.0/neps/npy-format.html
293     // for description of NPY array serialization format
294     if( nDataLen < 10 )
295     {
296         CPLError(CE_Failure, CPLE_AppDefined, "Non NPY array returned");
297         return false;
298     }
299 
300     if( memcmp(pabyData, "\x93NUMPY", 6) != 0 )
301     {
302         CPLError(CE_Failure, CPLE_AppDefined, "Non NPY array returned");
303         return false;
304     }
305     const int nVersionMajor = pabyData[6];
306     if( nVersionMajor != 1 )
307     {
308         CPLError(CE_Failure, CPLE_AppDefined,
309                  "Only version 1 of NPY array supported. Here found %d",
310                  nVersionMajor);
311         return false;
312     }
313     // Ignore version minor
314     const int nHeaderLen = pabyData[8] | (pabyData[9] << 8);
315     if( nDataLen < 10 + nHeaderLen )
316     {
317         CPLError(CE_Failure, CPLE_AppDefined,
318                  "Corrupted NPY array returned: not enough bytes for header");
319         return false;
320     }
321 
322 #ifdef DEBUG
323     CPLString osDescr;
324     osDescr.assign(reinterpret_cast<const char *>(pabyData) + 10, nHeaderLen);
325     // Should be something like
326     // {'descr': [('B2', '<u2'), ('B3', '<u2'), ('B4', '<u2'), ('B8', '<u2'),
327     // ('QA10', '<u2')], 'fortran_order': False, 'shape': (256, 256), }
328     CPLDebug("EEDAI", "NPY descr: %s", osDescr.c_str());
329     // TODO: validate that the descr is the one expected
330 #endif
331 
332     int nTotalDataTypeSize = 0;
333     for( int i = 1; i <= poGDS->GetRasterCount(); i++ )
334     {
335         if( bQueryAllBands || i == nBand )
336         {
337             nTotalDataTypeSize += GDALGetDataTypeSizeBytes(
338                         poGDS->GetRasterBand(i)->GetRasterDataType());
339         }
340     }
341     int nDataSize = nTotalDataTypeSize * nReqXSize * nReqYSize;
342     if( nDataLen < 10 + nHeaderLen + nDataSize )
343     {
344         CPLError(CE_Failure, CPLE_AppDefined,
345                  "Corrupted NPY array returned: not enough bytes for payload. "
346                  "%d needed, only %d found",
347                  10 + nHeaderLen + nDataSize, nDataLen);
348         return false;
349     }
350     else if( nDataLen > 10 + nHeaderLen + nDataSize )
351     {
352         CPLError(CE_Warning, CPLE_AppDefined,
353                  "Possibly corrupted NPY array returned: "
354                  "expected bytes for payload. "
355                  "%d needed, got %d found",
356                  10 + nHeaderLen + nDataSize, nDataLen);
357     }
358 
359     for( int iYBlock = 0; iYBlock < nYBlocks; iYBlock++ )
360     {
361         int nBlockActualYSize = nBlockYSize;
362         if( (iYBlock + nBlockYOff + 1) * nBlockYSize > nRasterYSize )
363         {
364             nBlockActualYSize = nRasterYSize -
365                                     (iYBlock + nBlockYOff) * nBlockYSize;
366         }
367 
368         for( int iXBlock = 0; iXBlock < nXBlocks; iXBlock++ )
369         {
370             int nBlockActualXSize = nBlockXSize;
371             if( (iXBlock + nBlockXOff + 1) * nBlockXSize > nRasterXSize )
372             {
373                 nBlockActualXSize = nRasterXSize -
374                                         (iXBlock + nBlockXOff) * nBlockXSize;
375             }
376 
377             int nOffsetBand = 10 + nHeaderLen +
378                 (iYBlock * nBlockYSize * nReqXSize + iXBlock * nBlockXSize) *
379                         nTotalDataTypeSize;
380 
381             for( int i = 1; i <= poGDS->GetRasterCount(); i++ )
382             {
383                 GDALRasterBlock* poBlock = nullptr;
384                 GByte* pabyDstBuffer;
385                 if( i == nBand && pDstBuffer != nullptr )
386                     pabyDstBuffer = reinterpret_cast<GByte*>(pDstBuffer);
387                 else if( bQueryAllBands || (i == nBand && pDstBuffer == nullptr) )
388                 {
389                     GDALEEDAIRasterBand* poOtherBand =
390                         reinterpret_cast<GDALEEDAIRasterBand*>(
391                                                     poGDS->GetRasterBand(i) );
392                     poBlock = poOtherBand->TryGetLockedBlockRef(
393                         nBlockXOff + iXBlock, nBlockYOff + iYBlock);
394                     if (poBlock != nullptr)
395                     {
396                         poBlock->DropLock();
397                         continue;
398                     }
399                     poBlock = poOtherBand->GetLockedBlockRef(
400                         nBlockXOff + iXBlock, nBlockYOff + iYBlock, TRUE);
401                     if (poBlock == nullptr)
402                     {
403                         continue;
404                     }
405                     pabyDstBuffer =
406                             reinterpret_cast<GByte*>(poBlock->GetDataRef());
407                 }
408                 else
409                 {
410                     continue;
411                 }
412 
413                 GDALDataType eDT = poGDS->GetRasterBand(i)->GetRasterDataType();
414                 const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
415 
416                 for( int iLine = 0; iLine < nBlockActualYSize; iLine++ )
417                 {
418                     GByte* pabyLineDest =
419                         pabyDstBuffer + iLine * nDTSize * nBlockXSize;
420                     GDALCopyWords(
421                         const_cast<GByte*>(pabyData) + nOffsetBand +
422                                     iLine * nTotalDataTypeSize * nReqXSize,
423                         eDT, nTotalDataTypeSize,
424                         pabyLineDest,
425                         eDT, nDTSize,
426                         nBlockActualXSize);
427 #ifdef CPL_MSB
428                     if( nDTSize > 1 )
429                     {
430                         GDALSwapWords(pabyLineDest, nDTSize,
431                                       nBlockActualXSize, nDTSize);
432                     }
433 #endif
434                 }
435 
436                 nOffsetBand += nDTSize;
437 
438                 if( poBlock )
439                     poBlock->DropLock();
440             }
441         }
442     }
443     return true;
444 }
445 /************************************************************************/
446 /*                            DecodeGDALDataset()                         */
447 /************************************************************************/
448 
DecodeGDALDataset(const GByte * pabyData,int nDataLen,bool bQueryAllBands,void * pDstBuffer,int nBlockXOff,int nBlockYOff,int nXBlocks,int nYBlocks,int nReqXSize,int nReqYSize)449 bool GDALEEDAIRasterBand::DecodeGDALDataset( const GByte* pabyData,
450                                            int nDataLen,
451                                            bool bQueryAllBands,
452                                            void* pDstBuffer,
453                                            int nBlockXOff, int nBlockYOff,
454                                            int nXBlocks, int nYBlocks,
455                                            int nReqXSize, int nReqYSize )
456 {
457     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
458 
459     CPLString osTmpFilename(CPLSPrintf("/vsimem/eeai/%p", this));
460     VSIFCloseL(VSIFileFromMemBuffer(osTmpFilename,
461                                     const_cast<GByte*>(pabyData),
462                                     nDataLen,
463                                     false));
464     const char* const apszDrivers[] = { "PNG", "JPEG", "GTIFF", nullptr };
465     GDALDataset* poTileDS = reinterpret_cast<GDALDataset*>(
466         GDALOpenEx(osTmpFilename, GDAL_OF_RASTER, apszDrivers, nullptr, nullptr));
467     if( poTileDS == nullptr )
468     {
469         CPLError(CE_Failure, CPLE_AppDefined,
470                     "Cannot decode buffer returned by the "
471                     "server as a PNG, JPEG or GeoTIFF image");
472         VSIUnlink(osTmpFilename);
473         return false;
474     }
475     if( poTileDS->GetRasterXSize() != nReqXSize ||
476         poTileDS->GetRasterYSize() != nReqYSize ||
477         // The server might return a RGBA image even if only 3 bands are requested
478         poTileDS->GetRasterCount() <
479                 (bQueryAllBands ? poGDS->GetRasterCount() : 1) )
480     {
481         CPLError(CE_Failure, CPLE_AppDefined,
482                  "Bad dimensions/band count for image returned "
483                  "by server: %dx%dx%d",
484                  poTileDS->GetRasterXSize(),
485                  poTileDS->GetRasterYSize(),
486                  poTileDS->GetRasterCount());
487         delete poTileDS;
488         VSIUnlink(osTmpFilename);
489         return false;
490     }
491 
492     for( int iYBlock = 0; iYBlock < nYBlocks; iYBlock++ )
493     {
494         int nBlockActualYSize = nBlockYSize;
495         if( (iYBlock + nBlockYOff + 1) * nBlockYSize > nRasterYSize )
496         {
497             nBlockActualYSize = nRasterYSize -
498                                         (iYBlock + nBlockYOff) * nBlockYSize;
499         }
500 
501         for( int iXBlock = 0; iXBlock < nXBlocks; iXBlock++ )
502         {
503             int nBlockActualXSize = nBlockXSize;
504             if( (iXBlock + nBlockXOff + 1) * nBlockXSize > nRasterXSize )
505             {
506                 nBlockActualXSize = nRasterXSize -
507                                         (iXBlock + nBlockXOff) * nBlockXSize;
508             }
509 
510             for(int i=1; i <= poGDS->GetRasterCount(); i++)
511             {
512                 GDALRasterBlock* poBlock = nullptr;
513                 GByte* pabyDstBuffer;
514                 if( i == nBand && pDstBuffer != nullptr )
515                     pabyDstBuffer = reinterpret_cast<GByte*>(pDstBuffer);
516                 else if( bQueryAllBands || (i == nBand && pDstBuffer == nullptr) )
517                 {
518                     GDALEEDAIRasterBand* poOtherBand =
519                         reinterpret_cast<GDALEEDAIRasterBand*>(
520                                                     poGDS->GetRasterBand(i) );
521                     poBlock = poOtherBand->TryGetLockedBlockRef(
522                         nBlockXOff + iXBlock, nBlockYOff + iYBlock);
523                     if (poBlock != nullptr)
524                     {
525                         poBlock->DropLock();
526                         continue;
527                     }
528                     poBlock = poOtherBand->GetLockedBlockRef(
529                         nBlockXOff + iXBlock, nBlockYOff + iYBlock, TRUE);
530                     if (poBlock == nullptr)
531                     {
532                         continue;
533                     }
534                     pabyDstBuffer =
535                             reinterpret_cast<GByte*>(poBlock->GetDataRef());
536                 }
537                 else
538                 {
539                     continue;
540                 }
541 
542                 GDALDataType eDT = poGDS->GetRasterBand(i)->GetRasterDataType();
543                 const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
544                 const int nTileBand = bQueryAllBands ? i : 1;
545                 CPLErr eErr = poTileDS->GetRasterBand(nTileBand)->
546                     RasterIO(GF_Read,
547                              iXBlock * nBlockXSize,
548                              iYBlock * nBlockYSize,
549                              nBlockActualXSize, nBlockActualYSize,
550                              pabyDstBuffer,
551                              nBlockActualXSize, nBlockActualYSize,
552                              eDT,
553                              nDTSize, nDTSize * nBlockXSize, nullptr);
554 
555                 if( poBlock )
556                     poBlock->DropLock();
557                 if( eErr != CE_None )
558                 {
559                     delete poTileDS;
560                     VSIUnlink(osTmpFilename);
561                     return false;
562                 }
563             }
564         }
565     }
566 
567     delete poTileDS;
568     VSIUnlink(osTmpFilename);
569     return true;
570 }
571 
GetBlocks(int nBlockXOff,int nBlockYOff,int nXBlocks,int nYBlocks,bool bQueryAllBands,void * pBuffer)572 CPLErr GDALEEDAIRasterBand::GetBlocks(int nBlockXOff, int nBlockYOff,
573                                       int nXBlocks, int nYBlocks,
574                                       bool bQueryAllBands,
575                                       void* pBuffer)
576 {
577     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
578 
579     // Build request content
580     json_object* poReq = json_object_new_object();
581     json_object_object_add(poReq, "fileFormat",
582                            json_object_new_string(poGDS->m_osPixelEncoding));
583     json_object* poBands = json_object_new_array();
584     for( int i = 1; i <= poGDS->GetRasterCount(); i++ )
585     {
586         if( bQueryAllBands || i == nBand )
587         {
588             json_object_array_add(poBands,
589                 json_object_new_string(
590                     poGDS->GetRasterBand(i)->GetDescription()));
591         }
592     }
593     json_object_object_add(poReq, "bandIds",
594                            poBands);
595 
596     int nReqXSize = nBlockXSize * nXBlocks;
597     if( (nBlockXOff + nXBlocks) * nBlockXSize > nRasterXSize )
598         nReqXSize = nRasterXSize - nBlockXOff * nBlockXSize;
599     int nReqYSize = nBlockYSize * nYBlocks;
600     if( (nBlockYOff + nYBlocks) * nBlockYSize > nRasterYSize )
601         nReqYSize = nRasterYSize - nBlockYOff * nBlockYSize;
602     const double dfX0 = poGDS->m_adfGeoTransform[0] +
603         nBlockXOff * nBlockXSize * poGDS->m_adfGeoTransform[1];
604     const double dfY0 = poGDS->m_adfGeoTransform[3] +
605         nBlockYOff * nBlockYSize * poGDS->m_adfGeoTransform[5];
606 #ifdef DEBUG_VERBOSE
607     CPLDebug("EEDAI", "nBlockYOff=%d nBlockYOff=%d "
608              "nXBlocks=%d nYBlocks=%d nReqXSize=%d nReqYSize=%d",
609              nBlockYOff, nBlockYOff, nXBlocks, nYBlocks, nReqXSize, nReqYSize);
610 #endif
611 
612     json_object* poPixelGrid = json_object_new_object();
613 
614     json_object* poAffineTransform = json_object_new_object();
615     json_object_object_add(poAffineTransform, "translateX",
616         json_object_new_double_with_significant_figures(dfX0, 18));
617     json_object_object_add(poAffineTransform, "translateY",
618         json_object_new_double_with_significant_figures(dfY0, 18));
619     json_object_object_add(poAffineTransform, "scaleX",
620         json_object_new_double_with_significant_figures(
621             poGDS->m_adfGeoTransform[1], 18));
622     json_object_object_add(poAffineTransform, "scaleY",
623         json_object_new_double_with_significant_figures(
624             poGDS->m_adfGeoTransform[5], 18));
625     json_object_object_add(poAffineTransform, "shearX",
626         json_object_new_double_with_significant_figures(0.0, 18));
627     json_object_object_add(poAffineTransform, "shearY",
628         json_object_new_double_with_significant_figures(0.0, 18));
629     json_object_object_add(poPixelGrid, "affineTransform", poAffineTransform);
630 
631     json_object* poDimensions = json_object_new_object();
632     json_object_object_add(poDimensions, "width",
633                            json_object_new_int(nReqXSize));
634     json_object_object_add(poDimensions, "height",
635                            json_object_new_int(nReqYSize));
636     json_object_object_add(poPixelGrid, "dimensions", poDimensions);
637     json_object_object_add(poReq, "grid", poPixelGrid);
638 
639     CPLString osPostContent = json_object_get_string(poReq);
640     json_object_put(poReq);
641 
642     // Issue request
643     char** papszOptions = (poGDS->m_poParentDS) ?
644         poGDS->m_poParentDS->GetBaseHTTPOptions() :
645         poGDS->GetBaseHTTPOptions();
646     papszOptions = CSLSetNameValue(papszOptions, "CUSTOMREQUEST", "POST");
647     CPLString osHeaders = CSLFetchNameValueDef(papszOptions, "HEADERS", "");
648     if( !osHeaders.empty() )
649         osHeaders += "\r\n";
650     osHeaders += "Content-Type: application/json";
651     papszOptions = CSLSetNameValue(papszOptions, "HEADERS", osHeaders);
652     papszOptions = CSLSetNameValue(papszOptions, "POSTFIELDS", osPostContent);
653     CPLHTTPResult* psResult = EEDAHTTPFetch(
654         (poGDS->m_osBaseURL + poGDS->m_osAssetName + ":getPixels").c_str(),
655         papszOptions);
656     CSLDestroy(papszOptions);
657     if( psResult == nullptr )
658         return CE_Failure;
659 
660     if( psResult->pszErrBuf != nullptr )
661     {
662         if( psResult->pabyData )
663         {
664             CPLError(CE_Failure, CPLE_AppDefined, "%s: %s",
665                      psResult->pszErrBuf,
666                      reinterpret_cast<const char*>(psResult->pabyData));
667         }
668         else
669         {
670             CPLError(CE_Failure, CPLE_AppDefined, "%s",
671                      psResult->pabyData ? reinterpret_cast<const char*>(psResult->pabyData) :
672                      psResult->pszErrBuf);
673         }
674         CPLHTTPDestroyResult(psResult);
675         return CE_Failure;
676     }
677 
678     if( psResult->pabyData == nullptr )
679     {
680         CPLError(CE_Failure, CPLE_AppDefined,
681                  "Empty content returned by server");
682         CPLHTTPDestroyResult(psResult);
683         return CE_Failure;
684     }
685 #ifdef DEBUG_VERBOSE
686     CPLDebug("EEADI", "Result: %s (%d bytes)",
687              reinterpret_cast<const char*>(psResult->pabyData),
688              psResult->nDataLen);
689 #endif
690 
691     if( EQUAL(poGDS->m_osPixelEncoding, "NPY") )
692     {
693         if( !DecodeNPYArray( psResult->pabyData, psResult->nDataLen,
694                              bQueryAllBands,
695                              pBuffer,
696                              nBlockXOff, nBlockYOff,
697                              nXBlocks, nYBlocks,
698                              nReqXSize, nReqYSize ) )
699         {
700             CPLHTTPDestroyResult(psResult);
701             return CE_Failure;
702         }
703     }
704     else
705     {
706         if( !DecodeGDALDataset( psResult->pabyData, psResult->nDataLen,
707                               bQueryAllBands,
708                               pBuffer,
709                               nBlockXOff, nBlockYOff,
710                               nXBlocks, nYBlocks,
711                               nReqXSize, nReqYSize ) )
712         {
713             CPLHTTPDestroyResult(psResult);
714             return CE_Failure;
715         }
716     }
717 
718     CPLHTTPDestroyResult(psResult);
719 
720     return CE_None;
721 }
722 
723 
724 /************************************************************************/
725 /*                             IReadBlock()                             */
726 /************************************************************************/
727 
IReadBlock(int nBlockXOff,int nBlockYOff,void * pBuffer)728 CPLErr GDALEEDAIRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
729                                         void *pBuffer )
730 {
731     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
732 #ifdef DEBUG_VERBOSE
733     CPLDebug("EEDAI", "ReadBlock x=%d y=%d band=%d level=%d",
734              nBlockXOff, nBlockYOff, nBand, poGDS->m_iOvrLevel);
735 #endif
736 
737     return GetBlocks(nBlockXOff, nBlockYOff, 1, 1,
738                      poGDS->m_bQueryMultipleBands,
739                      pBuffer);
740 }
741 
742 /************************************************************************/
743 /*                          PrefetchBlocks()                            */
744 /************************************************************************/
745 
746 // Return or'ed flags among 0, RETRY_PER_BAND, RETRY_SPATIAL_SPLIT if the user
747 // should try to split the request in smaller chunks
748 
PrefetchBlocks(int nXOff,int nYOff,int nXSize,int nYSize,int nBufXSize,int nBufYSize,bool bQueryAllBands)749 GUInt32 GDALEEDAIRasterBand::PrefetchBlocks(int nXOff, int nYOff,
750                                          int nXSize, int nYSize,
751                                          int nBufXSize, int nBufYSize,
752                                          bool bQueryAllBands)
753 {
754     CPL_IGNORE_RET_VAL(nBufXSize);
755     CPL_IGNORE_RET_VAL(nBufYSize);
756 
757     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
758     int nBlockXOff = nXOff / nBlockXSize;
759     int nBlockYOff = nYOff / nBlockYSize;
760     int nXBlocks = (nXOff + nXSize - 1) / nBlockXSize - nBlockXOff + 1;
761     int nYBlocks = (nYOff + nYSize - 1) / nBlockYSize - nBlockYOff + 1;
762 
763     const int nThisDTSize = GDALGetDataTypeSizeBytes(GetRasterDataType());
764     int nTotalDataTypeSize = 0;
765     int nQueriedBands = 0;
766     for( int i = 1; i <= poGDS->GetRasterCount(); i++ )
767     {
768         if( bQueryAllBands || i == nBand )
769         {
770             nQueriedBands ++;
771             nTotalDataTypeSize += GDALGetDataTypeSizeBytes(
772                         poGDS->GetRasterBand(i)->GetRasterDataType());
773         }
774     }
775 
776     // Check the number of already cached blocks, and remove fully
777     // cached lines at the top of the area of interest from the queried
778     // blocks
779     int nBlocksCached = 0;
780     int nBlocksCachedForThisBand = 0;
781     bool bAllLineCached = true;
782     for( int iYBlock = 0; iYBlock < nYBlocks; )
783     {
784         for( int iXBlock = 0; iXBlock < nXBlocks; iXBlock++ )
785         {
786             for(int i=1; i <= poGDS->GetRasterCount(); i++)
787             {
788                 GDALRasterBlock* poBlock = nullptr;
789                 if( bQueryAllBands || i == nBand )
790                 {
791                     GDALEEDAIRasterBand* poOtherBand =
792                         reinterpret_cast<GDALEEDAIRasterBand*>(
793                                                 poGDS->GetRasterBand(i) );
794                     poBlock = poOtherBand->TryGetLockedBlockRef(
795                         nBlockXOff + iXBlock, nBlockYOff + iYBlock);
796                     if (poBlock != nullptr)
797                     {
798                         nBlocksCached ++;
799                         if( i == nBand )
800                             nBlocksCachedForThisBand ++;
801                         poBlock->DropLock();
802                         continue;
803                     }
804                     else
805                     {
806                         bAllLineCached = false;
807                     }
808                 }
809             }
810         }
811 
812         if( bAllLineCached )
813         {
814             nBlocksCached -= nXBlocks * nQueriedBands;
815             nBlocksCachedForThisBand -= nXBlocks;
816             nBlockYOff ++;
817             nYBlocks --;
818         }
819         else
820         {
821             iYBlock ++;
822         }
823     }
824 
825     if( nXBlocks > 0 && nYBlocks > 0 )
826     {
827         bool bMustReturn = false;
828         GUInt32 nRetryFlags = 0;
829 
830         // Get the blocks if the number of already cached blocks is lesser
831         // than 25% of the to be queried blocks
832         if( nBlocksCached > (nQueriedBands * nXBlocks * nYBlocks) / 4 )
833         {
834             if( nBlocksCachedForThisBand <= (nXBlocks * nYBlocks) / 4 )
835             {
836                 nRetryFlags |= RETRY_PER_BAND;
837             }
838             else
839             {
840                 bMustReturn = true;
841             }
842         }
843 
844         // Don't request too many pixels in one dimension
845         if( nXBlocks * nBlockXSize > SERVER_DIMENSION_LIMIT ||
846             nYBlocks * nBlockYSize > SERVER_DIMENSION_LIMIT )
847         {
848             bMustReturn = true;
849             nRetryFlags |= RETRY_SPATIAL_SPLIT;
850         }
851 
852         // Make sure that we have enough cache (with a margin of 50%)
853         // and the number of queried pixels isn't too big w.r.t server
854         // limit
855         const GIntBig nUncompressedSize =
856             static_cast<GIntBig>(nXBlocks) * nYBlocks *
857                         nBlockXSize * nBlockYSize * nTotalDataTypeSize;
858         const GIntBig nCacheMax = GDALGetCacheMax64()/2;
859         if( nUncompressedSize > nCacheMax ||
860             nUncompressedSize > SERVER_BYTE_LIMIT )
861         {
862             if( bQueryAllBands && poGDS->GetRasterCount() > 1 )
863             {
864                 const GIntBig nUncompressedSizeThisBand =
865                     static_cast<GIntBig>(nXBlocks) * nYBlocks *
866                             nBlockXSize * nBlockYSize * nThisDTSize;
867                 if( nUncompressedSizeThisBand <= SERVER_BYTE_LIMIT &&
868                     nUncompressedSizeThisBand <= nCacheMax )
869                 {
870                     nRetryFlags |= RETRY_PER_BAND;
871                 }
872             }
873             if( nXBlocks > 1 || nYBlocks > 1 )
874             {
875                 nRetryFlags |= RETRY_SPATIAL_SPLIT;
876             }
877             return nRetryFlags;
878         }
879         if( bMustReturn )
880             return nRetryFlags;
881 
882         GetBlocks(nBlockXOff, nBlockYOff, nXBlocks, nYBlocks,
883                    bQueryAllBands, nullptr);
884     }
885 
886     return 0;
887 }
888 
889 /************************************************************************/
890 /*                              IRasterIO()                             */
891 /************************************************************************/
892 
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)893 CPLErr GDALEEDAIRasterBand::IRasterIO( GDALRWFlag eRWFlag,
894                                   int nXOff, int nYOff, int nXSize, int nYSize,
895                                   void * pData, int nBufXSize, int nBufYSize,
896                                   GDALDataType eBufType,
897                                   GSpacing nPixelSpace, GSpacing nLineSpace,
898                                   GDALRasterIOExtraArg* psExtraArg )
899 
900 {
901 
902 /* ==================================================================== */
903 /*      Do we have overviews that would be appropriate to satisfy       */
904 /*      this request?                                                   */
905 /* ==================================================================== */
906     if( (nBufXSize < nXSize || nBufYSize < nYSize)
907         && GetOverviewCount() > 0 && eRWFlag == GF_Read )
908     {
909         GDALRasterIOExtraArg sExtraArg;
910         GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
911 
912         const int nOverview =
913             GDALBandGetBestOverviewLevel2( this, nXOff, nYOff, nXSize, nYSize,
914                                            nBufXSize, nBufYSize, &sExtraArg );
915         if (nOverview >= 0)
916         {
917             GDALRasterBand* poOverviewBand = GetOverview(nOverview);
918             if (poOverviewBand == nullptr)
919                 return CE_Failure;
920 
921             return poOverviewBand->RasterIO(
922                 eRWFlag, nXOff, nYOff, nXSize, nYSize,
923                 pData, nBufXSize, nBufYSize, eBufType,
924                 nPixelSpace, nLineSpace, &sExtraArg );
925         }
926     }
927 
928     GDALEEDAIDataset* poGDS = reinterpret_cast<GDALEEDAIDataset*>(poDS);
929     GUInt32 nRetryFlags = PrefetchBlocks(
930         nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
931         poGDS->m_bQueryMultipleBands);
932     if( (nRetryFlags & RETRY_SPATIAL_SPLIT) &&
933         nXSize == nBufXSize && nYSize == nBufYSize && nYSize > nBlockYSize )
934     {
935         GDALRasterIOExtraArg sExtraArg;
936         INIT_RASTERIO_EXTRA_ARG(sExtraArg);
937 
938         int nHalf = std::max(nBlockYSize,
939                              ((nYSize / 2 ) / nBlockYSize) * nBlockYSize);
940         CPLErr eErr = IRasterIO(eRWFlag, nXOff, nYOff,
941                                 nXSize, nHalf,
942                                 pData,
943                                 nXSize, nHalf,
944                                 eBufType,
945                                 nPixelSpace, nLineSpace,
946                                 &sExtraArg);
947         if( eErr == CE_None )
948         {
949             eErr = IRasterIO(eRWFlag,
950                                 nXOff, nYOff + nHalf,
951                                 nXSize, nYSize - nHalf,
952                                 static_cast<GByte*>(pData) +
953                                         nHalf * nLineSpace,
954                                 nXSize, nYSize - nHalf,
955                                 eBufType,
956                                 nPixelSpace, nLineSpace,
957                                 &sExtraArg);
958         }
959         return eErr;
960     }
961     else if( (nRetryFlags & RETRY_SPATIAL_SPLIT) &&
962         nXSize == nBufXSize && nYSize == nBufYSize && nXSize > nBlockXSize )
963     {
964         GDALRasterIOExtraArg sExtraArg;
965         INIT_RASTERIO_EXTRA_ARG(sExtraArg);
966 
967         int nHalf = std::max(nBlockXSize,
968                              ((nXSize / 2 ) / nBlockXSize) * nBlockXSize);
969         CPLErr eErr = IRasterIO(eRWFlag, nXOff, nYOff,
970                                 nHalf, nYSize,
971                                 pData,
972                                 nHalf, nYSize,
973                                 eBufType,
974                                 nPixelSpace, nLineSpace,
975                                 &sExtraArg);
976         if( eErr == CE_None )
977         {
978             eErr = IRasterIO(eRWFlag,
979                                 nXOff + nHalf, nYOff,
980                                 nXSize - nHalf, nYSize,
981                                 static_cast<GByte*>(pData) +
982                                         nHalf * nPixelSpace,
983                                 nXSize - nHalf, nYSize,
984                                 eBufType,
985                                 nPixelSpace, nLineSpace,
986                                 &sExtraArg);
987         }
988         return eErr;
989     }
990     else if( (nRetryFlags & RETRY_PER_BAND) &&
991              poGDS->m_bQueryMultipleBands && poGDS->nBands > 1 )
992     {
993         CPL_IGNORE_RET_VAL(PrefetchBlocks(
994             nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, false));
995     }
996 
997     return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
998                                      pData, nBufXSize, nBufYSize,
999                                      eBufType,
1000                                      nPixelSpace, nLineSpace,
1001                                      psExtraArg);
1002 }
1003 
1004 
1005 /************************************************************************/
1006 /*                              IRasterIO()                             */
1007 /************************************************************************/
1008 
1009 
1010 CPLErr
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)1011 GDALEEDAIDataset::IRasterIO( GDALRWFlag eRWFlag,
1012                                  int nXOff, int nYOff, int nXSize, int nYSize,
1013                                  void * pData, int nBufXSize, int nBufYSize,
1014                                  GDALDataType eBufType,
1015                                  int nBandCount, int *panBandMap,
1016                                  GSpacing nPixelSpace, GSpacing nLineSpace,
1017                                  GSpacing nBandSpace,
1018                                  GDALRasterIOExtraArg* psExtraArg )
1019 {
1020 
1021 /* ==================================================================== */
1022 /*      Do we have overviews that would be appropriate to satisfy       */
1023 /*      this request?                                                   */
1024 /* ==================================================================== */
1025     if( (nBufXSize < nXSize || nBufYSize < nYSize)
1026         && GetRasterBand(1)->GetOverviewCount() > 0 && eRWFlag == GF_Read )
1027     {
1028         GDALRasterIOExtraArg sExtraArg;
1029         GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
1030 
1031         const int nOverview =
1032             GDALBandGetBestOverviewLevel2( GetRasterBand(1),
1033                                             nXOff, nYOff, nXSize, nYSize,
1034                                            nBufXSize, nBufYSize, &sExtraArg );
1035         if (nOverview >= 0)
1036         {
1037             GDALRasterBand* poOverviewBand =
1038                         GetRasterBand(1)->GetOverview(nOverview);
1039             if (poOverviewBand == nullptr ||
1040                 poOverviewBand->GetDataset() == nullptr)
1041             {
1042                 return CE_Failure;
1043             }
1044 
1045             return poOverviewBand->GetDataset()->RasterIO(
1046                 eRWFlag, nXOff, nYOff, nXSize, nYSize,
1047                 pData, nBufXSize, nBufYSize, eBufType,
1048                 nBandCount, panBandMap,
1049                 nPixelSpace, nLineSpace, nBandSpace, &sExtraArg );
1050         }
1051     }
1052 
1053     GDALEEDAIRasterBand* poBand =
1054         cpl::down_cast<GDALEEDAIRasterBand*>(GetRasterBand(1));
1055 
1056     GUInt32 nRetryFlags = poBand->PrefetchBlocks(
1057                                 nXOff, nYOff, nXSize, nYSize,
1058                                 nBufXSize, nBufYSize,
1059                                 m_bQueryMultipleBands);
1060     int nBlockXSize, nBlockYSize;
1061     poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
1062     if( (nRetryFlags & RETRY_SPATIAL_SPLIT) &&
1063         nXSize == nBufXSize && nYSize == nBufYSize && nYSize > nBlockYSize )
1064     {
1065         GDALRasterIOExtraArg sExtraArg;
1066         INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1067 
1068         int nHalf = std::max(nBlockYSize,
1069                          ((nYSize / 2 ) / nBlockYSize) * nBlockYSize);
1070         CPLErr eErr = IRasterIO(eRWFlag, nXOff, nYOff,
1071                                 nXSize, nHalf,
1072                                 pData,
1073                                 nXSize, nHalf,
1074                                 eBufType,
1075                                 nBandCount, panBandMap,
1076                                 nPixelSpace, nLineSpace, nBandSpace,
1077                                 &sExtraArg);
1078         if( eErr == CE_None )
1079         {
1080             eErr = IRasterIO(eRWFlag,
1081                                 nXOff, nYOff + nHalf,
1082                                 nXSize, nYSize - nHalf,
1083                                 static_cast<GByte*>(pData) +
1084                                     nHalf * nLineSpace,
1085                                 nXSize, nYSize - nHalf,
1086                                 eBufType,
1087                                 nBandCount, panBandMap,
1088                                 nPixelSpace, nLineSpace, nBandSpace,
1089                                 &sExtraArg);
1090         }
1091         return eErr;
1092     }
1093     else if( (nRetryFlags & RETRY_SPATIAL_SPLIT) &&
1094         nXSize == nBufXSize && nYSize == nBufYSize && nXSize > nBlockXSize )
1095     {
1096         GDALRasterIOExtraArg sExtraArg;
1097         INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1098 
1099         int nHalf = std::max(nBlockXSize,
1100                          ((nXSize / 2 ) / nBlockXSize) * nBlockXSize);
1101         CPLErr eErr = IRasterIO(eRWFlag, nXOff, nYOff,
1102                                 nHalf, nYSize,
1103                                 pData,
1104                                 nHalf, nYSize,
1105                                 eBufType,
1106                                 nBandCount, panBandMap,
1107                                 nPixelSpace, nLineSpace, nBandSpace,
1108                                 &sExtraArg);
1109         if( eErr == CE_None )
1110         {
1111             eErr = IRasterIO(eRWFlag,
1112                                 nXOff + nHalf, nYOff,
1113                                 nXSize - nHalf, nYSize,
1114                                 static_cast<GByte*>(pData) +
1115                                         nHalf * nPixelSpace,
1116                                 nXSize - nHalf, nYSize,
1117                                 eBufType,
1118                                 nBandCount, panBandMap,
1119                                 nPixelSpace, nLineSpace, nBandSpace,
1120                                 &sExtraArg);
1121         }
1122         return eErr;
1123     }
1124     else if( (nRetryFlags & RETRY_PER_BAND) &&
1125              m_bQueryMultipleBands && nBands > 1 )
1126     {
1127         for( int iBand = 1; iBand <= nBands; iBand++ )
1128         {
1129             poBand =
1130                 cpl::down_cast<GDALEEDAIRasterBand*>(GetRasterBand(iBand));
1131             CPL_IGNORE_RET_VAL(poBand->PrefetchBlocks(
1132                                     nXOff, nYOff, nXSize, nYSize,
1133                                     nBufXSize, nBufYSize, false));
1134         }
1135     }
1136 
1137     return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
1138                                   pData, nBufXSize, nBufYSize,
1139                                   eBufType,
1140                                   nBandCount, panBandMap,
1141                                   nPixelSpace, nLineSpace, nBandSpace,
1142                                   psExtraArg);
1143 }
1144 
1145 /************************************************************************/
1146 /*                        ComputeQueryStrategy()                        */
1147 /************************************************************************/
1148 
ComputeQueryStrategy()1149 bool GDALEEDAIDataset::ComputeQueryStrategy()
1150 {
1151     m_bQueryMultipleBands = true;
1152     m_osPixelEncoding.toupper();
1153 
1154     bool bHeterogeneousDataTypes = false;
1155     if( nBands >= 2 )
1156     {
1157         GDALDataType eDTFirstBand = GetRasterBand(1)->GetRasterDataType();
1158         for( int i = 2; i <= nBands; i++ )
1159         {
1160             if( GetRasterBand(i)->GetRasterDataType() != eDTFirstBand )
1161             {
1162                 bHeterogeneousDataTypes = true;
1163                 break;
1164             }
1165         }
1166     }
1167 
1168     if( EQUAL(m_osPixelEncoding, "AUTO") )
1169     {
1170         if( bHeterogeneousDataTypes )
1171         {
1172             m_osPixelEncoding = "NPY";
1173         }
1174         else
1175         {
1176             m_osPixelEncoding = "PNG";
1177             for( int i = 1; i <= nBands; i++ )
1178             {
1179                 if( GetRasterBand(i)->GetRasterDataType() != GDT_Byte )
1180                 {
1181                     m_osPixelEncoding = "GEO_TIFF";
1182                 }
1183             }
1184         }
1185     }
1186 
1187     if( EQUAL(m_osPixelEncoding, "PNG") ||
1188         EQUAL(m_osPixelEncoding, "JPEG") ||
1189         EQUAL(m_osPixelEncoding, "AUTO_JPEG_PNG") )
1190     {
1191         if( nBands != 1 && nBands != 3 )
1192         {
1193             m_bQueryMultipleBands = false;
1194         }
1195         for( int i = 1; i <= nBands; i++ )
1196         {
1197             if( GetRasterBand(i)->GetRasterDataType() != GDT_Byte )
1198             {
1199                 CPLError(CE_Failure, CPLE_NotSupported,
1200                  "This dataset has non-Byte bands, which is incompatible "
1201                  "with PIXEL_ENCODING=%s", m_osPixelEncoding.c_str());
1202                 return false;
1203             }
1204         }
1205     }
1206 
1207     if(nBands > SERVER_SIMUTANEOUS_BAND_LIMIT )
1208     {
1209         m_bQueryMultipleBands = false;
1210     }
1211 
1212     if( m_bQueryMultipleBands && m_osPixelEncoding != "NPY" &&
1213         bHeterogeneousDataTypes )
1214     {
1215         CPLDebug("EEDAI",
1216                  "%s PIXEL_ENCODING does not support heterogeneous data types. "
1217                  "Falling back to querying band per band",
1218                  m_osPixelEncoding.c_str());
1219         m_bQueryMultipleBands = false;
1220     }
1221 
1222     return true;
1223 }
1224 
1225 /************************************************************************/
1226 /*                          GetProjectionRef()                          */
1227 /************************************************************************/
1228 
_GetProjectionRef()1229 const char* GDALEEDAIDataset::_GetProjectionRef()
1230 {
1231     return m_osWKT.c_str();
1232 }
1233 
1234 /************************************************************************/
1235 /*                          GetGeoTransform()                           */
1236 /************************************************************************/
1237 
GetGeoTransform(double * adfGeoTransform)1238 CPLErr GDALEEDAIDataset::GetGeoTransform( double* adfGeoTransform )
1239 {
1240     memcpy( adfGeoTransform, m_adfGeoTransform, 6 * sizeof(double) );
1241     return CE_None;
1242 }
1243 
1244 /************************************************************************/
1245 /*                               Open()                                 */
1246 /************************************************************************/
1247 
Open(GDALOpenInfo * poOpenInfo)1248 bool GDALEEDAIDataset::Open(GDALOpenInfo* poOpenInfo)
1249 {
1250     m_osBaseURL = CPLGetConfigOption("EEDA_URL",
1251                             "https://earthengine-highvolume.googleapis.com/v1alpha/");
1252 
1253     m_osAsset =
1254             CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "ASSET", "");
1255     CPLString osBandList(
1256             CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "BANDS", "") );
1257     if( m_osAsset.empty() )
1258     {
1259         char** papszTokens =
1260                 CSLTokenizeString2(poOpenInfo->pszFilename, ":", 0);
1261         if( CSLCount(papszTokens) < 2 )
1262         {
1263             CPLError(CE_Failure, CPLE_AppDefined,
1264                      "No asset specified in connection string or "
1265                      "ASSET open option");
1266             CSLDestroy(papszTokens);
1267             return false;
1268         }
1269         if( CSLCount(papszTokens) == 3 )
1270         {
1271             osBandList = papszTokens[2];
1272         }
1273 
1274         m_osAsset = papszTokens[1];
1275         CSLDestroy(papszTokens);
1276     }
1277     m_osAssetName = ConvertPathToName(m_osAsset);
1278 
1279     m_osPixelEncoding =
1280         CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
1281                              "PIXEL_ENCODING", "AUTO");
1282     m_nBlockSize = atoi(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
1283                         "BLOCK_SIZE", CPLSPrintf("%d", DEFAULT_BLOCK_SIZE)));
1284     if( m_nBlockSize < 128 &&
1285         !CPLTestBool(CPLGetConfigOption("EEDA_FORCE_BLOCK_SIZE", "FALSE")) )
1286     {
1287         CPLError(CE_Failure, CPLE_NotSupported, "Invalid BLOCK_SIZE");
1288         return false;
1289     }
1290 
1291     std::set<CPLString> oSetUserBandNames;
1292     {
1293         char** papszTokens = CSLTokenizeString2(osBandList, ",", 0);
1294         for( int i = 0; papszTokens && papszTokens[i]; i++ )
1295             oSetUserBandNames.insert(papszTokens[i]);
1296         CSLDestroy(papszTokens);
1297     }
1298 
1299     // Issue request to get image metadata
1300     char** papszOptions = GetBaseHTTPOptions();
1301     if( papszOptions == nullptr )
1302         return false;
1303     CPLHTTPResult* psResult = EEDAHTTPFetch(
1304                 (m_osBaseURL + m_osAssetName).c_str(), papszOptions);
1305     CSLDestroy(papszOptions);
1306     if( psResult == nullptr )
1307         return false;
1308     if( psResult->pszErrBuf != nullptr )
1309     {
1310         if( psResult->pabyData )
1311         {
1312             CPLError(CE_Failure, CPLE_AppDefined, "%s: %s",
1313                      psResult->pszErrBuf,
1314                      reinterpret_cast<const char*>(psResult->pabyData));
1315         }
1316         else
1317         {
1318             CPLError(CE_Failure, CPLE_AppDefined, "%s",
1319                      psResult->pabyData ? reinterpret_cast<const char*>(psResult->pabyData) :
1320                      psResult->pszErrBuf);
1321         }
1322         CPLHTTPDestroyResult(psResult);
1323         return false;
1324     }
1325 
1326     if( psResult->pabyData == nullptr )
1327     {
1328         CPLError(CE_Failure, CPLE_AppDefined,
1329                  "Empty content returned by server");
1330         CPLHTTPDestroyResult(psResult);
1331         return false;
1332     }
1333 
1334     const char* pszText = reinterpret_cast<const char*>(psResult->pabyData);
1335 #ifdef DEBUG_VERBOSE
1336     CPLDebug("EEDAI", "%s", pszText);
1337 #endif
1338 
1339     json_object* poObj = nullptr;
1340     if( !OGRJSonParse(pszText, &poObj, true) )
1341     {
1342         CPLHTTPDestroyResult(psResult);
1343         return false;
1344     }
1345 
1346     CPLHTTPDestroyResult(psResult);
1347 
1348     if( json_object_get_type(poObj) != json_type_object )
1349     {
1350         CPLError( CE_Failure, CPLE_AppDefined,
1351                   "Return is not a JSON dictionary");
1352         json_object_put(poObj);
1353         return false;
1354     }
1355 
1356     json_object* poType = CPL_json_object_object_get(poObj, "type");
1357     const char* pszType = json_object_get_string(poType);
1358     if( pszType == nullptr || !EQUAL(pszType, "IMAGE") )
1359     {
1360         CPLError( CE_Failure, CPLE_AppDefined,
1361                   "Asset is not an image, but %s",
1362                   pszType ? pszType : "(null)" );
1363         json_object_put(poObj);
1364         return false;
1365     }
1366 
1367     json_object* poBands = CPL_json_object_object_get(poObj, "bands");
1368     if( poBands == nullptr ||  json_object_get_type(poBands) != json_type_array )
1369     {
1370         CPLError(CE_Failure, CPLE_AppDefined, "No band found");
1371         json_object_put(poObj);
1372         return false;
1373     }
1374 
1375     std::map<CPLString, CPLString> oMapCodeToWKT;
1376     std::vector<EEDAIBandDesc> aoBandDesc = BuildBandDescArray(poBands,
1377                                                                oMapCodeToWKT);
1378     std::map<CPLString, int> aoMapBandNames;
1379 
1380     if( aoBandDesc.empty() )
1381     {
1382         CPLError(CE_Failure, CPLE_AppDefined, "No band found");
1383         json_object_put(poObj);
1384         return false;
1385     }
1386 
1387     // Indices are aoBandDesc indices
1388     std::map< int, std::vector<int> > oMapSimilarBands;
1389 
1390     size_t iIdxFirstBand = 0;
1391     for( size_t i = 0; i < aoBandDesc.size(); i++ )
1392     {
1393         // Instantiate bands if they are compatible between them, and
1394         // if they are requested by the user (when user explicitly requested
1395         // them)
1396         if( (oSetUserBandNames.empty() ||
1397              oSetUserBandNames.find(aoBandDesc[i].osName) !=
1398                                         oSetUserBandNames.end() ) &&
1399             (nBands == 0 ||
1400              aoBandDesc[i].IsSimilar(aoBandDesc[iIdxFirstBand])) )
1401         {
1402             if( nBands == 0 )
1403             {
1404                 iIdxFirstBand = i;
1405                 nRasterXSize = aoBandDesc[i].nWidth;
1406                 nRasterYSize = aoBandDesc[i].nHeight;
1407                 memcpy(m_adfGeoTransform, aoBandDesc[i].adfGeoTransform.data(),
1408                        6 * sizeof(double));
1409                 m_osWKT = aoBandDesc[i].osWKT;
1410                 int iOvr = 0;
1411                 while( (nRasterXSize >> iOvr) > 256 ||
1412                        (nRasterYSize >> iOvr) > 256 )
1413                 {
1414                     iOvr ++;
1415                     m_apoOverviewDS.push_back(
1416                         new GDALEEDAIDataset(this, iOvr));
1417                 }
1418             }
1419 
1420             GDALRasterBand* poBand =
1421                 new GDALEEDAIRasterBand(this, aoBandDesc[i].eDT,
1422                                         aoBandDesc[i].bSignedByte);
1423             const int iBand = nBands + 1;
1424             SetBand( iBand, poBand );
1425             poBand->SetDescription( aoBandDesc[i].osName );
1426 
1427             // as images in USDA/NAIP/DOQQ catalog
1428             if( EQUAL(aoBandDesc[i].osName, "R") )
1429                 poBand->SetColorInterpretation(GCI_RedBand);
1430             else if( EQUAL(aoBandDesc[i].osName, "G") )
1431                 poBand->SetColorInterpretation(GCI_GreenBand);
1432             else if( EQUAL(aoBandDesc[i].osName, "B") )
1433                 poBand->SetColorInterpretation(GCI_BlueBand);
1434 
1435             for(size_t iOvr = 0; iOvr < m_apoOverviewDS.size(); iOvr++ )
1436             {
1437                 GDALRasterBand* poOvrBand =
1438                     new GDALEEDAIRasterBand(m_apoOverviewDS[iOvr],
1439                                             aoBandDesc[i].eDT,
1440                                             aoBandDesc[i].bSignedByte);
1441                 m_apoOverviewDS[iOvr]->SetBand( iBand, poOvrBand );
1442                 poOvrBand->SetDescription( aoBandDesc[i].osName );
1443             }
1444 
1445             aoMapBandNames[ aoBandDesc[i].osName ] = iBand;
1446         }
1447         else
1448         {
1449             if( oSetUserBandNames.find(aoBandDesc[i].osName) !=
1450                                                 oSetUserBandNames.end() )
1451             {
1452                 CPLError(CE_Warning, CPLE_AppDefined,
1453                          "Band %s is not compatible of other bands",
1454                          aoBandDesc[i].osName.c_str());
1455             }
1456             aoMapBandNames[ aoBandDesc[i].osName ] = -1;
1457         }
1458 
1459         // Group similar bands to be able to build subdataset list
1460         std::map< int, std::vector<int> >::iterator oIter =
1461                                                     oMapSimilarBands.begin();
1462         for( ; oIter != oMapSimilarBands.end(); ++oIter )
1463         {
1464             if( aoBandDesc[i].IsSimilar(aoBandDesc[oIter->first]) )
1465             {
1466                 oIter->second.push_back(static_cast<int>(i));
1467                 break;
1468             }
1469         }
1470         if( oIter == oMapSimilarBands.end() )
1471         {
1472             oMapSimilarBands[static_cast<int>(i)].push_back(static_cast<int>(i));
1473         }
1474     }
1475 
1476     if( !ComputeQueryStrategy() )
1477     {
1478         json_object_put(poObj);
1479         return false;
1480     }
1481     for( size_t i = 0; i < m_apoOverviewDS.size(); i++ )
1482     {
1483         m_apoOverviewDS[i]->ComputeQueryStrategy();
1484     }
1485 
1486     if( nBands > 1 )
1487     {
1488         SetMetadataItem("INTERLEAVE",
1489                         m_bQueryMultipleBands ? "PIXEL" : "BAND",
1490                         "IMAGE_STRUCTURE");
1491     }
1492 
1493     // Build subdataset list
1494     if( oSetUserBandNames.empty() && oMapSimilarBands.size() > 1 )
1495     {
1496         CPLStringList aoSubDSList;
1497         std::map< int, std::vector<int> >::iterator oIter =
1498                                                     oMapSimilarBands.begin();
1499         for( ; oIter != oMapSimilarBands.end(); ++oIter )
1500         {
1501             CPLString osSubDSSuffix;
1502             for( size_t i = 0; i < oIter->second.size(); ++i )
1503             {
1504                 if( !osSubDSSuffix.empty() )
1505                     osSubDSSuffix += ",";
1506                 osSubDSSuffix += aoBandDesc[oIter->second[i]].osName;
1507             }
1508             aoSubDSList.AddNameValue(
1509                 CPLSPrintf("SUBDATASET_%d_NAME", aoSubDSList.size() / 2 + 1),
1510                 CPLSPrintf("EEDAI:%s:%s",
1511                            m_osAsset.c_str(), osSubDSSuffix.c_str()) );
1512             aoSubDSList.AddNameValue(
1513                 CPLSPrintf("SUBDATASET_%d_DESC", aoSubDSList.size() / 2 + 1),
1514                 CPLSPrintf("Band%s %s of %s",
1515                            oIter->second.size() > 1 ? "s" : "",
1516                            osSubDSSuffix.c_str(), m_osAsset.c_str()) );
1517         }
1518         SetMetadata( aoSubDSList.List(), "SUBDATASETS" );
1519     }
1520 
1521     // Attach metadata to dataset or bands
1522     json_object* poProperties = CPL_json_object_object_get(poObj,
1523                                                            "properties");
1524     if( poProperties &&
1525         json_object_get_type(poProperties) == json_type_object )
1526     {
1527         SetMetadataFromProperties(poProperties, aoMapBandNames);
1528     }
1529     json_object_put(poObj);
1530 
1531     SetDescription( poOpenInfo->pszFilename );
1532 
1533     return true;
1534 }
1535 
1536 /************************************************************************/
1537 /*                       SetMetadataFromProperties()                    */
1538 /************************************************************************/
1539 
SetMetadataFromProperties(json_object * poProperties,const std::map<CPLString,int> & aoMapBandNames)1540 void GDALEEDAIDataset::SetMetadataFromProperties(
1541                 json_object* poProperties,
1542                 const std::map<CPLString, int>& aoMapBandNames )
1543 {
1544     json_object_iter it;
1545     it.key = nullptr;
1546     it.val = nullptr;
1547     it.entry = nullptr;
1548     json_object_object_foreachC(poProperties, it)
1549     {
1550         if( it.val )
1551         {
1552             CPLString osKey(it.key);
1553             int nBandForMD = 0;
1554             std::map<CPLString, int>::const_iterator oIter =
1555                                                 aoMapBandNames.begin();
1556             for( ; oIter != aoMapBandNames.end(); ++oIter )
1557             {
1558                 CPLString osBandName(oIter->first);
1559                 CPLString osNeedle("_" + osBandName);
1560                 size_t nPos = osKey.find(osNeedle);
1561                 if( nPos != std::string::npos &&
1562                     nPos + osNeedle.size() == osKey.size() )
1563                 {
1564                     nBandForMD = oIter->second;
1565                     osKey.resize(nPos);
1566                     break;
1567                 }
1568 
1569                 // Landsat bands are named Bxxx, must their metadata
1570                 // are _BAND_xxxx ...
1571                 if( osBandName.size() > 1 && osBandName[0] == 'B' &&
1572                     atoi(osBandName.c_str() + 1) > 0 )
1573                 {
1574                     osNeedle = "_BAND_" + osBandName.substr(1);
1575                     nPos = osKey.find(osNeedle);
1576                     if( nPos != std::string::npos &&
1577                         nPos + osNeedle.size() == osKey.size() )
1578                     {
1579                         nBandForMD = oIter->second;
1580                         osKey.resize(nPos);
1581                         break;
1582                     }
1583                 }
1584             }
1585 
1586             if( nBandForMD > 0 )
1587             {
1588                 GetRasterBand(nBandForMD)->SetMetadataItem(
1589                     osKey, json_object_get_string(it.val));
1590             }
1591             else if( nBandForMD == 0 )
1592             {
1593                 SetMetadataItem(osKey, json_object_get_string(it.val));
1594             }
1595         }
1596     }
1597 }
1598 
1599 /************************************************************************/
1600 /*                          GDALEEDAIIdentify()                         */
1601 /************************************************************************/
1602 
GDALEEDAIIdentify(GDALOpenInfo * poOpenInfo)1603 static int GDALEEDAIIdentify(GDALOpenInfo* poOpenInfo)
1604 {
1605     return STARTS_WITH_CI(poOpenInfo->pszFilename, "EEDAI:");
1606 }
1607 
1608 
1609 /************************************************************************/
1610 /*                            GDALEEDAIOpen()                           */
1611 /************************************************************************/
1612 
GDALEEDAIOpen(GDALOpenInfo * poOpenInfo)1613 static GDALDataset* GDALEEDAIOpen(GDALOpenInfo* poOpenInfo)
1614 {
1615     if(! GDALEEDAIIdentify(poOpenInfo) )
1616         return nullptr;
1617 
1618     GDALEEDAIDataset* poDS = new GDALEEDAIDataset();
1619     if( !poDS->Open(poOpenInfo) )
1620     {
1621         delete poDS;
1622         return nullptr;
1623     }
1624     return poDS;
1625 }
1626 
1627 /************************************************************************/
1628 /*                         GDALRegister_EEDAI()                         */
1629 /************************************************************************/
1630 
GDALRegister_EEDAI()1631 void GDALRegister_EEDAI()
1632 
1633 {
1634     if( GDALGetDriverByName( "EEDAI" ) != nullptr )
1635         return;
1636 
1637     GDALDriver *poDriver = new GDALDriver();
1638 
1639     poDriver->SetDescription( "EEDAI" );
1640     poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
1641     poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1642                                "Earth Engine Data API Image" );
1643     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/eedai.html" );
1644     poDriver->SetMetadataItem( GDAL_DMD_CONNECTION_PREFIX, "EEDAI:" );
1645     poDriver->SetMetadataItem( GDAL_DMD_OPENOPTIONLIST,
1646 "<OpenOptionList>"
1647 "  <Option name='ASSET' type='string' description='Asset name'/>"
1648 "  <Option name='BANDS' type='string' "
1649                         "description='Comma separated list of band names'/>"
1650 "  <Option name='PIXEL_ENCODING' type='string-select' "
1651                         "description='Format in which pixls are queried'>"
1652 "       <Value>AUTO</Value>"
1653 "       <Value>PNG</Value>"
1654 "       <Value>JPEG</Value>"
1655 "       <Value>GEO_TIFF</Value>"
1656 "       <Value>AUTO_JPEG_PNG</Value>"
1657 "       <Value>NPY</Value>"
1658 "   </Option>"
1659 "  <Option name='BLOCK_SIZE' type='integer' "
1660                                 "description='Size of a block' default='256'/>"
1661 "</OpenOptionList>");
1662     poDriver->SetMetadataItem( GDAL_DMD_SUBDATASETS, "YES" );
1663 
1664     poDriver->pfnOpen = GDALEEDAIOpen;
1665     poDriver->pfnIdentify = GDALEEDAIIdentify;
1666 
1667     GetGDALDriverManager()->RegisterDriver( poDriver );
1668 }
1669