1 /******************************************************************************
2  *
3  * Project:  ECRG TOC read Translator
4  * Purpose:  Implementation of ECRGTOCDataset and ECRGTOCSubDataset.
5  * Author:   Even Rouault, even.rouault at spatialys.com
6  *
7  ******************************************************************************
8  * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
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 // g++ -g -Wall -fPIC frmts/nitf/ecrgtocdataset.cpp -shared -o gdal_ECRGTOC.so -Iport -Igcore -Iogr -Ifrmts/vrt -L. -lgdal
30 
31 #include "cpl_port.h"
32 
33 #include <cassert>
34 #include <cmath>
35 #include <cstddef>
36 #include <cstdio>
37 #include <cstdlib>
38 #include <cstring>
39 #include <memory>
40 #include <string>
41 #include <vector>
42 
43 #include "cpl_conv.h"
44 #include "cpl_error.h"
45 #include "cpl_minixml.h"
46 #include "cpl_string.h"
47 #include "gdal.h"
48 #include "gdal_frmts.h"
49 #include "gdal_pam.h"
50 #include "gdal_priv.h"
51 #include "gdal_proxy.h"
52 #include "ogr_srs_api.h"
53 #include "vrtdataset.h"
54 
55 CPL_CVSID("$Id: ecrgtocdataset.cpp f6099e5ed704166bf5cc113a053dd1b2725cb391 2020-03-22 11:20:10 +0100 Kai Pastor $")
56 
57 /** Overview of used classes :
58    - ECRGTOCDataset : lists the different subdatasets, listed in the .xml,
59                       as subdatasets
60    - ECRGTOCSubDataset : one of these subdatasets, implemented as a VRT, of
61                          the relevant NITF tiles
62    - ECRGTOCProxyRasterDataSet : a "proxy" dataset that maps to a NITF tile
63 */
64 
65 typedef struct
66 {
67     const char* pszName;
68     const char* pszPath;
69     int         nScale;
70     int         nZone;
71 } FrameDesc;
72 
73 /************************************************************************/
74 /* ==================================================================== */
75 /*                            ECRGTOCDataset                            */
76 /* ==================================================================== */
77 /************************************************************************/
78 
79 class ECRGTOCDataset final: public GDALPamDataset
80 {
81   char      **papszSubDatasets;
82   double      adfGeoTransform[6];
83 
84   char      **papszFileList;
85 
86   public:
87     ECRGTOCDataset()
88     {
89         papszSubDatasets = nullptr;
90         papszFileList = nullptr;
91         memset( adfGeoTransform, 0, sizeof(adfGeoTransform) );
92     }
93 
94     virtual ~ECRGTOCDataset()
95     {
96         CSLDestroy( papszSubDatasets );
97         CSLDestroy(papszFileList);
98     }
99 
100     virtual char      **GetMetadata( const char * pszDomain = "" ) override;
101 
102     virtual char      **GetFileList() override { return CSLDuplicate(papszFileList); }
103 
104     void                AddSubDataset(const char* pszFilename,
105                                       const char* pszProductTitle,
106                                       const char* pszDiscId,
107                                       const char* pszScale);
108 
109     virtual CPLErr GetGeoTransform( double * padfGeoTransform) override
110     {
111         memcpy(padfGeoTransform, adfGeoTransform, 6 * sizeof(double));
112         return CE_None;
113     }
114 
115     virtual const char *_GetProjectionRef(void) override
116     {
117         return SRS_WKT_WGS84_LAT_LONG;
118     }
119     const OGRSpatialReference* GetSpatialRef() const override {
120         return GetSpatialRefFromOldGetProjectionRef();
121     }
122 
123     static GDALDataset* Build(  const char* pszTOCFilename,
124                                 CPLXMLNode* psXML,
125                                 CPLString osProduct,
126                                 CPLString osDiscId,
127                                 CPLString osScale,
128                                 const char* pszFilename);
129 
130     static int Identify( GDALOpenInfo * poOpenInfo );
131     static GDALDataset* Open( GDALOpenInfo * poOpenInfo );
132 };
133 
134 /************************************************************************/
135 /* ==================================================================== */
136 /*                            ECRGTOCSubDataset                          */
137 /* ==================================================================== */
138 /************************************************************************/
139 
140 class ECRGTOCSubDataset final: public VRTDataset
141 {
142   char**       papszFileList;
143 
144   public:
145     ECRGTOCSubDataset(int nXSize, int nYSize) : VRTDataset(nXSize, nYSize)
146     {
147         /* Don't try to write a VRT file */
148         SetWritable(FALSE);
149 
150         /* The driver is set to VRT in VRTDataset constructor. */
151         /* We have to set it to the expected value ! */
152         poDriver = reinterpret_cast<GDALDriver *>(
153             GDALGetDriverByName( "ECRGTOC" ) );
154 
155         papszFileList = nullptr;
156     }
157 
158     ~ECRGTOCSubDataset()
159     {
160         CSLDestroy(papszFileList);
161     }
162 
163     virtual char      **GetFileList() override { return CSLDuplicate(papszFileList); }
164 
165     static GDALDataset* Build(  const char* pszProductTitle,
166                                 const char* pszDiscId,
167                                 int nScale,
168                                 int nCountSubDataset,
169                                 const char* pszTOCFilename,
170                                 const std::vector<FrameDesc>& aosFrameDesc,
171                                 double dfGlobalMinX,
172                                 double dfGlobalMinY,
173                                 double dfGlobalMaxX,
174                                 double dfGlobalMaxY,
175                                 double dfGlobalPixelXSize,
176                                 double dfGlobalPixelYSize);
177 };
178 
179 /************************************************************************/
180 /*                           LaunderString()                            */
181 /************************************************************************/
182 
183 static CPLString LaunderString(const char* pszStr)
184 {
185     CPLString osRet(pszStr);
186     for(size_t i=0;i<osRet.size();i++)
187     {
188         if( osRet[i] == ':' || osRet[i] == ' ' )
189             osRet[i] = '_';
190     }
191     return osRet;
192 }
193 
194 /************************************************************************/
195 /*                           AddSubDataset()                            */
196 /************************************************************************/
197 
198 void ECRGTOCDataset::AddSubDataset( const char* pszFilename,
199                                     const char* pszProductTitle,
200                                     const char* pszDiscId,
201                                     const char* pszScale)
202 
203 {
204     char szName[80];
205     const int nCount = CSLCount(papszSubDatasets ) / 2;
206 
207     snprintf( szName, sizeof(szName), "SUBDATASET_%d_NAME", nCount+1 );
208     papszSubDatasets =
209         CSLSetNameValue( papszSubDatasets, szName,
210               CPLSPrintf( "ECRG_TOC_ENTRY:%s:%s:%s:%s",
211                           LaunderString(pszProductTitle).c_str(),
212                           LaunderString(pszDiscId).c_str(),
213                           LaunderString(pszScale).c_str(), pszFilename ) );
214 
215     snprintf( szName, sizeof(szName), "SUBDATASET_%d_DESC", nCount+1 );
216     papszSubDatasets =
217         CSLSetNameValue( papszSubDatasets, szName,
218             CPLSPrintf( "Product %s, disc %s, scale %s", pszProductTitle, pszDiscId, pszScale));
219 }
220 
221 /************************************************************************/
222 /*                            GetMetadata()                             */
223 /************************************************************************/
224 
225 char **ECRGTOCDataset::GetMetadata( const char *pszDomain )
226 
227 {
228     if( pszDomain != nullptr && EQUAL(pszDomain,"SUBDATASETS") )
229         return papszSubDatasets;
230 
231     return GDALPamDataset::GetMetadata( pszDomain );
232 }
233 
234 /************************************************************************/
235 /*                         GetScaleFromString()                         */
236 /************************************************************************/
237 
238 static int GetScaleFromString(const char* pszScale)
239 {
240     const char* pszPtr = strstr(pszScale, "1:");
241     if (pszPtr)
242         pszPtr = pszPtr + 2;
243     else
244         pszPtr = pszScale;
245 
246     int nScale = 0;
247     char ch;
248     while((ch = *pszPtr) != '\0')
249     {
250         if (ch >= '0' && ch <= '9')
251             nScale = nScale * 10 + ch - '0';
252         else if (ch == ' ')
253             ;
254         else if (ch == 'k' || ch == 'K')
255             return nScale * 1000;
256         else if (ch == 'm' || ch == 'M')
257             return nScale * 1000000;
258         else
259             return 0;
260         pszPtr ++;
261     }
262     return nScale;
263 }
264 
265 /************************************************************************/
266 /*                            GetFromBase34()                           */
267 /************************************************************************/
268 
269 static GIntBig GetFromBase34(const char* pszVal, int nMaxSize)
270 {
271     GIntBig nFrameNumber = 0;
272     for(int i=0;i<nMaxSize;i++)
273     {
274         char ch = pszVal[i];
275         if (ch == '\0')
276             break;
277         int chVal;
278         if (ch >= 'A' && ch <= 'Z')
279             ch += 'a' - 'A';
280         /* i and o letters are excluded, */
281         if (ch >= '0' && ch <= '9')
282             chVal = ch - '0';
283         else if (ch >= 'a' && ch <= 'h')
284             chVal = ch - 'a' + 10;
285         else if (ch >= 'j' && ch <= 'n')
286             chVal = ch - 'a' + 10 - 1;
287         else if (ch >= 'p' && ch <= 'z')
288             chVal = ch - 'a' + 10 - 2;
289         else
290         {
291             CPLDebug("ECRG", "Invalid base34 value : %s", pszVal);
292             break;
293         }
294         nFrameNumber = nFrameNumber * 34 + chVal;
295     }
296 
297     return nFrameNumber;
298 }
299 
300 /************************************************************************/
301 /*                             GetExtent()                              */
302 /************************************************************************/
303 
304 /* MIL-PRF-32283 - Table II. ECRG zone limits. */
305 /* starting with a fake zone 0 for convenience. */
306 constexpr int anZoneUpperLat[] = { 0, 32, 48, 56, 64, 68, 72, 76, 80 };
307 
308 /* APPENDIX 70, TABLE III of MIL-A-89007 */
309 constexpr int anACst_ADRG[] =
310     { 369664, 302592, 245760, 199168, 163328, 137216, 110080, 82432 };
311 constexpr int nBCst_ADRG = 400384;
312 
313 // TODO: Why are these two functions done this way?
314 static int CEIL_ROUND(double a, double b)
315 {
316     return static_cast<int>( ceil( a / b ) * b );
317 }
318 
319 static int NEAR_ROUND(double a, double b)
320 {
321     return static_cast<int>( floor( ( a / b ) + 0.5 ) * b );
322 }
323 
324 constexpr int ECRG_PIXELS = 2304;
325 
326 static
327 int GetExtent(const char* pszFrameName, int nScale, int nZone,
328               double& dfMinX, double& dfMaxX, double& dfMinY, double& dfMaxY,
329               double& dfPixelXSize, double& dfPixelYSize)
330 {
331     const int nAbsZone = abs(nZone);
332 #ifdef DEBUG
333     assert( nAbsZone > 0 && nAbsZone <= 8 );
334 #endif
335 
336 /************************************************************************/
337 /*  Compute east-west constant                                          */
338 /************************************************************************/
339     /* MIL-PRF-89038 - 60.1.2 - East-west pixel constant. */
340     const int nEW_ADRG = CEIL_ROUND(anACst_ADRG[nAbsZone-1] * (1e6 / nScale), 512);
341     const int nEW_CADRG = NEAR_ROUND(nEW_ADRG / (150. / 100.), 256);
342     /* MIL-PRF-32283 - D.2.1.2 - East-west pixel constant. */
343     const int nEW = nEW_CADRG / 256 * 384;
344 
345 /************************************************************************/
346 /*  Compute number of longitudinal frames                               */
347 /************************************************************************/
348     /* MIL-PRF-32283 - D.2.1.7 - Longitudinal frames and subframes */
349     const int nCols = static_cast<int>(
350         ceil(static_cast<double>(nEW) / ECRG_PIXELS) );
351 
352 /************************************************************************/
353 /*  Compute north-south constant                                        */
354 /************************************************************************/
355     /* MIL-PRF-89038 - 60.1.1 -  North-south. pixel constant */
356     const int nNS_ADRG = CEIL_ROUND(nBCst_ADRG * (1e6 / nScale), 512) / 4;
357     const int nNS_CADRG = NEAR_ROUND(nNS_ADRG / (150. / 100.), 256);
358     /* MIL-PRF-32283 - D.2.1.1 - North-south. pixel constant and Frame Width/Height */
359     const int nNS = nNS_CADRG / 256 * 384;
360 
361 /************************************************************************/
362 /*  Compute number of latitudinal frames and latitude of top of zone    */
363 /************************************************************************/
364     dfPixelYSize = 90.0 / nNS;
365 
366     const double dfFrameLatHeight = dfPixelYSize * ECRG_PIXELS;
367 
368     /* MIL-PRF-32283 - D.2.1.5 - Equatorward and poleward zone extents. */
369     int nUpperZoneFrames = static_cast<int>(
370         ceil(anZoneUpperLat[nAbsZone] / dfFrameLatHeight) );
371     int nBottomZoneFrames = static_cast<int>(
372         floor(anZoneUpperLat[nAbsZone-1] / dfFrameLatHeight) );
373     const int nRows = nUpperZoneFrames - nBottomZoneFrames;
374 
375     /* Not sure to really understand D.2.1.5.a. Testing needed */
376     if (nZone < 0)
377     {
378         nUpperZoneFrames = -nBottomZoneFrames;
379         /*nBottomZoneFrames = nUpperZoneFrames - nRows;*/
380     }
381 
382     const double dfUpperZoneTopLat = dfFrameLatHeight * nUpperZoneFrames;
383 
384 /************************************************************************/
385 /*  Compute coordinates of the frame in the zone                        */
386 /************************************************************************/
387 
388     /* Converts the first 10 characters into a number from base 34 */
389     const GIntBig nFrameNumber = GetFromBase34(pszFrameName, 10);
390 
391     /*  MIL-PRF-32283 - A.2.6.1 */
392     const GIntBig nY = nFrameNumber / nCols;
393     const GIntBig nX = nFrameNumber % nCols;
394 
395 /************************************************************************/
396 /*  Compute extent of the frame                                         */
397 /************************************************************************/
398 
399     /* The nY is counted from the bottom of the zone... Pfff */
400     dfMaxY = dfUpperZoneTopLat - (nRows - 1 - nY) * dfFrameLatHeight;
401     dfMinY = dfMaxY - dfFrameLatHeight;
402 
403     dfPixelXSize = 360.0 / nEW;
404 
405     const double dfFrameLongWidth = dfPixelXSize * ECRG_PIXELS;
406     dfMinX = -180.0 + nX * dfFrameLongWidth;
407     dfMaxX = dfMinX + dfFrameLongWidth;
408 
409 #ifdef DEBUG_VERBOSE
410     CPLDebug("ECRG", "Frame %s : minx=%.16g, maxy=%.16g, maxx=%.16g, miny=%.16g",
411              pszFrameName, dfMinX, dfMaxY, dfMaxX, dfMinY);
412 #endif
413 
414     return TRUE;
415 }
416 
417 /************************************************************************/
418 /* ==================================================================== */
419 /*                        ECRGTOCProxyRasterDataSet                       */
420 /* ==================================================================== */
421 /************************************************************************/
422 
423 class ECRGTOCProxyRasterDataSet final: public GDALProxyPoolDataset
424 {
425     /* The following parameters are only for sanity checking */
426     mutable int checkDone;
427     mutable int checkOK;
428     double dfMinX;
429     double dfMaxY;
430     double dfPixelXSize;
431     double dfPixelYSize;
432 
433     public:
434         ECRGTOCProxyRasterDataSet( ECRGTOCSubDataset* /* poSubDataset */,
435                                    const char* fileName,
436                                    int nXSize, int nYSize,
437                                    double dfMinX, double dfMaxY,
438                                    double dfPixelXSize, double dfPixelYSize );
439 
440         GDALDataset* RefUnderlyingDataset() const override
441         {
442             GDALDataset* poSourceDS = GDALProxyPoolDataset::RefUnderlyingDataset();
443             if (poSourceDS)
444             {
445                 if (!checkDone)
446                     SanityCheckOK(poSourceDS);
447                 if (!checkOK)
448                 {
449                     GDALProxyPoolDataset::UnrefUnderlyingDataset(poSourceDS);
450                     poSourceDS = nullptr;
451                 }
452             }
453             return poSourceDS;
454         }
455 
456         void UnrefUnderlyingDataset(GDALDataset* poUnderlyingDataset) const override
457         {
458             GDALProxyPoolDataset::UnrefUnderlyingDataset(poUnderlyingDataset);
459         }
460 
461         int SanityCheckOK(GDALDataset* poSourceDS) const;
462 };
463 
464 /************************************************************************/
465 /*                    ECRGTOCProxyRasterDataSet()                       */
466 /************************************************************************/
467 
468 ECRGTOCProxyRasterDataSet::ECRGTOCProxyRasterDataSet(
469     ECRGTOCSubDataset* /* poSubDatasetIn */,
470     const char* fileNameIn,
471     int nXSizeIn, int nYSizeIn,
472     double dfMinXIn, double dfMaxYIn,
473     double dfPixelXSizeIn, double dfPixelYSizeIn ) :
474     // Mark as shared since the VRT will take several references if we are in
475     // RGBA mode (4 bands for this dataset).
476     GDALProxyPoolDataset(fileNameIn, nXSizeIn, nYSizeIn, GA_ReadOnly,
477                          TRUE, SRS_WKT_WGS84_LAT_LONG),
478     checkDone(FALSE),
479     checkOK(FALSE),
480     dfMinX(dfMinXIn),
481     dfMaxY(dfMaxYIn),
482     dfPixelXSize(dfPixelXSizeIn),
483     dfPixelYSize(dfPixelYSizeIn)
484 {
485 
486     for( int i = 0; i < 3; i++ )
487     {
488         SetBand(i + 1,
489                 new GDALProxyPoolRasterBand(this, i+1, GDT_Byte, nXSizeIn, 1));
490     }
491 }
492 
493 /************************************************************************/
494 /*                    SanityCheckOK()                                   */
495 /************************************************************************/
496 
497 #define WARN_CHECK_DS(x) do { if (!(x)) { \
498     CPLError(CE_Warning, CPLE_AppDefined,                             \
499              "For %s, assert '" #x "' failed",                        \
500              GetDescription()); checkOK = FALSE; } } while( false )
501 
502 int ECRGTOCProxyRasterDataSet::SanityCheckOK( GDALDataset* poSourceDS ) const
503 {
504     // int nSrcBlockXSize;
505     // int nSrcBlockYSize;
506     // int nBlockXSize;
507     // int nBlockYSize;
508     double l_adfGeoTransform[6] = {};
509     if( checkDone )
510         return checkOK;
511 
512     checkOK = TRUE;
513     checkDone = TRUE;
514 
515     poSourceDS->GetGeoTransform(l_adfGeoTransform);
516     WARN_CHECK_DS(fabs(l_adfGeoTransform[0] - dfMinX) < 1e-10);
517     WARN_CHECK_DS(fabs(l_adfGeoTransform[3] - dfMaxY) < 1e-10);
518     WARN_CHECK_DS(fabs(l_adfGeoTransform[1] - dfPixelXSize) < 1e-10);
519     WARN_CHECK_DS(fabs(l_adfGeoTransform[5] - (-dfPixelYSize)) < 1e-10);
520     WARN_CHECK_DS(l_adfGeoTransform[2] == 0 &&
521                   l_adfGeoTransform[4] == 0);  // No rotation.
522     WARN_CHECK_DS(poSourceDS->GetRasterCount() == 3);
523     WARN_CHECK_DS(poSourceDS->GetRasterXSize() == nRasterXSize);
524     WARN_CHECK_DS(poSourceDS->GetRasterYSize() == nRasterYSize);
525     WARN_CHECK_DS(EQUAL(poSourceDS->GetProjectionRef(), SRS_WKT_WGS84_LAT_LONG));
526     // poSourceDS->GetRasterBand(1)->GetBlockSize(&nSrcBlockXSize,
527     //                                            &nSrcBlockYSize);
528     // GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
529     // WARN_CHECK_DS(nSrcBlockXSize == nBlockXSize);
530     // WARN_CHECK_DS(nSrcBlockYSize == nBlockYSize);
531     WARN_CHECK_DS(
532         poSourceDS->GetRasterBand(1)->GetRasterDataType() == GDT_Byte);
533 
534     return checkOK;
535 }
536 
537 /************************************************************************/
538 /*                           BuildFullName()                            */
539 /************************************************************************/
540 
541 static const char* BuildFullName(const char* pszTOCFilename,
542                                  const char* pszFramePath,
543                                  const char* pszFrameName)
544 {
545     char* pszPath = nullptr;
546     if (pszFramePath[0] == '.' &&
547         (pszFramePath[1] == '/' ||pszFramePath[1] == '\\'))
548         pszPath = CPLStrdup(pszFramePath + 2);
549     else
550         pszPath = CPLStrdup(pszFramePath);
551     for(int i=0;pszPath[i] != '\0';i++)
552     {
553         if (pszPath[i] == '\\')
554             pszPath[i] = '/';
555     }
556     const char* pszName = CPLFormFilename(pszPath, pszFrameName, nullptr);
557     CPLFree(pszPath);
558     pszPath = nullptr;
559     const char* pszTOCPath = CPLGetDirname(pszTOCFilename);
560     const char* pszFirstSlashInName = strchr(pszName, '/');
561     if (pszFirstSlashInName != nullptr)
562     {
563         int nFirstDirLen = static_cast<int>(pszFirstSlashInName - pszName);
564         if (static_cast<int>( strlen(pszTOCPath) ) >= nFirstDirLen + 1 &&
565             (pszTOCPath[strlen(pszTOCPath) - (nFirstDirLen + 1)] == '/' ||
566                 pszTOCPath[strlen(pszTOCPath) - (nFirstDirLen + 1)] == '\\') &&
567             strncmp(pszTOCPath + strlen(pszTOCPath) - nFirstDirLen, pszName, nFirstDirLen) == 0)
568         {
569             pszTOCPath = CPLGetDirname(pszTOCPath);
570         }
571     }
572     return CPLProjectRelativeFilename(pszTOCPath, pszName);
573 }
574 
575 /************************************************************************/
576 /*                              Build()                                 */
577 /************************************************************************/
578 
579 /* Builds a ECRGTOCSubDataset from the set of files of the toc entry */
580 GDALDataset* ECRGTOCSubDataset::Build(  const char* pszProductTitle,
581                                         const char* pszDiscId,
582                                         int nScale,
583                                         int nCountSubDataset,
584                                         const char* pszTOCFilename,
585                                         const std::vector<FrameDesc>& aosFrameDesc,
586                                         double dfGlobalMinX,
587                                         double dfGlobalMinY,
588                                         double dfGlobalMaxX,
589                                         double dfGlobalMaxY,
590                                         double dfGlobalPixelXSize,
591                                         double dfGlobalPixelYSize)
592 {
593     GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("VRT");
594     if( poDriver == nullptr )
595         return nullptr;
596 
597     const int nSizeX = static_cast<int>(
598         (dfGlobalMaxX - dfGlobalMinX) / dfGlobalPixelXSize + 0.5);
599     const int nSizeY = static_cast<int>(
600         (dfGlobalMaxY - dfGlobalMinY) / dfGlobalPixelYSize + 0.5);
601 
602     /* ------------------------------------ */
603     /* Create the VRT with the overall size */
604     /* ------------------------------------ */
605     ECRGTOCSubDataset *poVirtualDS = new ECRGTOCSubDataset( nSizeX, nSizeY );
606 
607     poVirtualDS->SetProjection(SRS_WKT_WGS84_LAT_LONG);
608 
609     double adfGeoTransform[6] = {
610       dfGlobalMinX,
611       dfGlobalPixelXSize,
612       0,
613       dfGlobalMaxY,
614       0,
615       -dfGlobalPixelYSize
616     };
617     poVirtualDS->SetGeoTransform(adfGeoTransform);
618 
619     for (int i=0;i<3;i++)
620     {
621         poVirtualDS->AddBand(GDT_Byte, nullptr);
622         GDALRasterBand *poBand = poVirtualDS->GetRasterBand( i + 1 );
623         poBand->SetColorInterpretation((GDALColorInterp)(GCI_RedBand+i));
624     }
625 
626     poVirtualDS->SetDescription(pszTOCFilename);
627 
628     poVirtualDS->SetMetadataItem("PRODUCT_TITLE", pszProductTitle);
629     poVirtualDS->SetMetadataItem("DISC_ID", pszDiscId);
630     if (nScale != -1)
631         poVirtualDS->SetMetadataItem("SCALE", CPLString().Printf("%d", nScale));
632 
633     /* -------------------------------------------------------------------- */
634     /*      Check for overviews.                                            */
635     /* -------------------------------------------------------------------- */
636 
637     poVirtualDS->oOvManager.Initialize( poVirtualDS,
638                                         CPLString().Printf("%s.%d", pszTOCFilename, nCountSubDataset));
639 
640     poVirtualDS->papszFileList = poVirtualDS->GDALDataset::GetFileList();
641 
642     for( int i=0; i < static_cast<int>( aosFrameDesc.size() ); i++)
643     {
644         const char* pszName = BuildFullName(pszTOCFilename,
645                                             aosFrameDesc[i].pszPath,
646                                             aosFrameDesc[i].pszName);
647 
648         double dfMinX = 0.0;
649         double dfMaxX = 0.0;
650         double dfMinY = 0.0;
651         double dfMaxY = 0.0;
652         double dfPixelXSize = 0.0;
653         double dfPixelYSize = 0.0;
654         GetExtent(aosFrameDesc[i].pszName,
655                   aosFrameDesc[i].nScale, aosFrameDesc[i].nZone,
656                   dfMinX, dfMaxX, dfMinY, dfMaxY, dfPixelXSize, dfPixelYSize);
657 
658         const int nFrameXSize = static_cast<int>(
659             (dfMaxX - dfMinX) / dfPixelXSize + 0.5);
660         const int nFrameYSize = static_cast<int>(
661             (dfMaxY - dfMinY) / dfPixelYSize + 0.5);
662 
663         poVirtualDS->papszFileList = CSLAddString(poVirtualDS->papszFileList, pszName);
664 
665         /* We create proxy datasets and raster bands */
666         /* Using real datasets and raster bands is possible in theory */
667         /* However for large datasets, a TOC entry can include several hundreds of files */
668         /* and we finally reach the limit of maximum file descriptors open at the same time ! */
669         /* So the idea is to warp the datasets into a proxy and open the underlying dataset only when it is */
670         /* needed (IRasterIO operation). To improve a bit efficiency, we have a cache of opened */
671         /* underlying datasets */
672         ECRGTOCProxyRasterDataSet* poDS = new ECRGTOCProxyRasterDataSet(
673             reinterpret_cast<ECRGTOCSubDataset *>( poVirtualDS), pszName,
674             nFrameXSize, nFrameYSize,
675             dfMinX, dfMaxY, dfPixelXSize, dfPixelYSize);
676 
677         for( int j=0; j<3; j++)
678         {
679             VRTSourcedRasterBand *poBand = reinterpret_cast<VRTSourcedRasterBand *>(
680                 poVirtualDS->GetRasterBand( j + 1 ) );
681             /* Place the raster band at the right position in the VRT */
682             poBand->AddSimpleSource(
683                 poDS->GetRasterBand(j + 1),
684                 0, 0, nFrameXSize, nFrameYSize,
685                 static_cast<int>((dfMinX - dfGlobalMinX) / dfGlobalPixelXSize + 0.5),
686                 static_cast<int>((dfGlobalMaxY - dfMaxY) / dfGlobalPixelYSize + 0.5),
687                 static_cast<int>((dfMaxX - dfMinX) / dfGlobalPixelXSize + 0.5),
688                 static_cast<int>((dfMaxY - dfMinY) / dfGlobalPixelYSize + 0.5));
689         }
690 
691         /* The ECRGTOCProxyRasterDataSet will be destroyed when its last raster band will be */
692         /* destroyed */
693         poDS->Dereference();
694     }
695 
696     poVirtualDS->SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
697 
698     return poVirtualDS;
699 }
700 
701 /************************************************************************/
702 /*                             Build()                                  */
703 /************************************************************************/
704 
705 GDALDataset* ECRGTOCDataset::Build(const char* pszTOCFilename,
706                                    CPLXMLNode* psXML,
707                                    CPLString osProduct,
708                                    CPLString osDiscId,
709                                    CPLString osScale,
710                                    const char* pszOpenInfoFilename)
711 {
712     CPLXMLNode* psTOC = CPLGetXMLNode(psXML, "=Table_of_Contents");
713     if (psTOC == nullptr)
714     {
715         CPLError(CE_Warning, CPLE_AppDefined,
716                  "Cannot find Table_of_Contents element");
717         return nullptr;
718     }
719 
720     double dfGlobalMinX = 0.0;
721     double dfGlobalMinY = 0.0;
722     double dfGlobalMaxX = 0.0;
723     double dfGlobalMaxY = 0.0;
724     double dfGlobalPixelXSize = 0.0;
725     double dfGlobalPixelYSize = 0.0;
726     bool bGlobalExtentValid = false;
727 
728     ECRGTOCDataset* poDS = new ECRGTOCDataset();
729     int nSubDatasets = 0;
730 
731     int bLookForSubDataset = !osProduct.empty() && !osDiscId.empty();
732 
733     int nCountSubDataset = 0;
734 
735     poDS->SetDescription( pszOpenInfoFilename );
736     poDS->papszFileList = poDS->GDALDataset::GetFileList();
737 
738     for(CPLXMLNode* psIter1 = psTOC->psChild;
739                     psIter1 != nullptr;
740                     psIter1 = psIter1->psNext)
741     {
742         if (!(psIter1->eType == CXT_Element && psIter1->pszValue != nullptr &&
743               strcmp(psIter1->pszValue, "product") == 0))
744             continue;
745 
746         const char* pszProductTitle =
747             CPLGetXMLValue(psIter1, "product_title", nullptr);
748         if (pszProductTitle == nullptr)
749         {
750             CPLError(CE_Warning, CPLE_AppDefined,
751                         "Cannot find product_title attribute");
752             continue;
753         }
754 
755         if (bLookForSubDataset && strcmp(LaunderString(pszProductTitle), osProduct.c_str()) != 0)
756             continue;
757 
758         for(CPLXMLNode* psIter2 = psIter1->psChild;
759                         psIter2 != nullptr;
760                         psIter2 = psIter2->psNext)
761         {
762             if (!(psIter2->eType == CXT_Element && psIter2->pszValue != nullptr &&
763                   strcmp(psIter2->pszValue, "disc") == 0))
764                 continue;
765 
766             const char* pszDiscId = CPLGetXMLValue(psIter2, "id", nullptr);
767             if (pszDiscId == nullptr)
768             {
769                 CPLError(CE_Warning, CPLE_AppDefined,
770                             "Cannot find id attribute");
771                 continue;
772             }
773 
774             if (bLookForSubDataset && strcmp(LaunderString(pszDiscId), osDiscId.c_str()) != 0)
775                 continue;
776 
777             CPLXMLNode* psFrameList = CPLGetXMLNode(psIter2, "frame_list");
778             if (psFrameList == nullptr)
779             {
780                 CPLError(CE_Warning, CPLE_AppDefined,
781                             "Cannot find frame_list element");
782                 continue;
783             }
784 
785             for(CPLXMLNode* psIter3 = psFrameList->psChild;
786                             psIter3 != nullptr;
787                             psIter3 = psIter3->psNext)
788             {
789                 if (!(psIter3->eType == CXT_Element &&
790                       psIter3->pszValue != nullptr &&
791                       strcmp(psIter3->pszValue, "scale") == 0))
792                     continue;
793 
794                 const char* pszSize = CPLGetXMLValue(psIter3, "size", nullptr);
795                 if (pszSize == nullptr)
796                 {
797                     CPLError(CE_Warning, CPLE_AppDefined,
798                                 "Cannot find size attribute");
799                     continue;
800                 }
801 
802                 int nScale = GetScaleFromString(pszSize);
803                 if (nScale <= 0)
804                 {
805                     CPLError(CE_Warning, CPLE_AppDefined,
806                              "Invalid scale %s", pszSize);
807                     continue;
808                 }
809 
810                 if( bLookForSubDataset )
811                 {
812                     if( !osScale.empty() )
813                     {
814                         if( strcmp(LaunderString(pszSize), osScale.c_str()) != 0 )
815                         {
816                             continue;
817                         }
818                     }
819                     else
820                     {
821                         int nCountScales = 0;
822                         for(CPLXMLNode* psIter4 = psFrameList->psChild;
823                                 psIter4 != nullptr;
824                                 psIter4 = psIter4->psNext)
825                         {
826                             if (!(psIter4->eType == CXT_Element &&
827                                 psIter4->pszValue != nullptr &&
828                                 strcmp(psIter4->pszValue, "scale") == 0))
829                                 continue;
830                             nCountScales ++;
831                         }
832                         if( nCountScales > 1 )
833                         {
834                             CPLError( CE_Failure, CPLE_AppDefined,
835                                       "Scale should be mentioned in "
836                                       "subdatasets syntax since this disk "
837                                       "contains several scales" );
838                             delete poDS;
839                             return nullptr;
840                         }
841                     }
842                 }
843 
844                 nCountSubDataset ++;
845 
846                 std::vector<FrameDesc> aosFrameDesc;
847                 int nValidFrames = 0;
848 
849                 for(CPLXMLNode* psIter4 = psIter3->psChild;
850                                 psIter4 != nullptr;
851                                 psIter4 = psIter4->psNext)
852                 {
853                     if (!(psIter4->eType == CXT_Element &&
854                           psIter4->pszValue != nullptr &&
855                           strcmp(psIter4->pszValue, "frame") == 0))
856                         continue;
857 
858                     const char* pszFrameName =
859                         CPLGetXMLValue(psIter4, "name", nullptr);
860                     if (pszFrameName == nullptr)
861                     {
862                         CPLError(CE_Warning, CPLE_AppDefined,
863                                  "Cannot find name element");
864                         continue;
865                     }
866 
867                     if (strlen(pszFrameName) != 18)
868                     {
869                         CPLError(CE_Warning, CPLE_AppDefined,
870                                  "Invalid value for name element : %s",
871                                  pszFrameName);
872                         continue;
873                     }
874 
875                     const char* pszFramePath =
876                         CPLGetXMLValue(psIter4, "frame_path", nullptr);
877                     if (pszFramePath == nullptr)
878                     {
879                         CPLError(CE_Warning, CPLE_AppDefined,
880                                  "Cannot find frame_path element");
881                         continue;
882                     }
883 
884                     const char* pszFrameZone =
885                         CPLGetXMLValue(psIter4, "frame_zone", nullptr);
886                     if (pszFrameZone == nullptr)
887                     {
888                         CPLError(CE_Warning, CPLE_AppDefined,
889                                  "Cannot find frame_zone element");
890                         continue;
891                     }
892                     if (strlen(pszFrameZone) != 1)
893                     {
894                         CPLError(CE_Warning, CPLE_AppDefined,
895                                  "Invalid value for frame_zone element : %s",
896                                  pszFrameZone);
897                         continue;
898                     }
899                     char chZone = pszFrameZone[0];
900                     int nZone = 0;
901                     if (chZone >= '1' && chZone <= '9')
902                         nZone = chZone - '0';
903                     else if (chZone >= 'a' && chZone <= 'h')
904                         nZone = -(chZone - 'a' + 1);
905                     else if (chZone >= 'A' && chZone <= 'H')
906                         nZone = -(chZone - 'A' + 1);
907                     else if (chZone == 'j' || chZone == 'J')
908                         nZone = -9;
909                     else
910                     {
911                         CPLError(CE_Warning, CPLE_AppDefined,
912                                  "Invalid value for frame_zone element : %s",
913                                   pszFrameZone);
914                         continue;
915                     }
916                     if (nZone == 9 || nZone == -9)
917                     {
918                         CPLError(CE_Warning, CPLE_AppDefined,
919                                  "Polar zones unhandled by current implementation");
920                         continue;
921                     }
922 
923                     double dfMinX = 0.0;
924                     double dfMaxX = 0.0;
925                     double dfMinY = 0.0;
926                     double dfMaxY = 0.0;
927                     double dfPixelXSize = 0.0;
928                     double dfPixelYSize = 0.0;
929                     if (!GetExtent(pszFrameName, nScale, nZone,
930                                    dfMinX, dfMaxX, dfMinY, dfMaxY,
931                                    dfPixelXSize, dfPixelYSize))
932                     {
933                         continue;
934                     }
935 
936                     nValidFrames ++;
937 
938                     const char* pszFullName = BuildFullName(pszTOCFilename,
939                                                         pszFramePath,
940                                                         pszFrameName);
941                     poDS->papszFileList = CSLAddString(poDS->papszFileList, pszFullName);
942 
943                     if (!bGlobalExtentValid)
944                     {
945                         dfGlobalMinX = dfMinX;
946                         dfGlobalMinY = dfMinY;
947                         dfGlobalMaxX = dfMaxX;
948                         dfGlobalMaxY = dfMaxY;
949                         dfGlobalPixelXSize = dfPixelXSize;
950                         dfGlobalPixelYSize = dfPixelYSize;
951                         bGlobalExtentValid = true;
952                     }
953                     else
954                     {
955                         if (dfMinX < dfGlobalMinX) dfGlobalMinX = dfMinX;
956                         if (dfMinY < dfGlobalMinY) dfGlobalMinY = dfMinY;
957                         if (dfMaxX > dfGlobalMaxX) dfGlobalMaxX = dfMaxX;
958                         if (dfMaxY > dfGlobalMaxY) dfGlobalMaxY = dfMaxY;
959                         if (dfPixelXSize < dfGlobalPixelXSize)
960                             dfGlobalPixelXSize = dfPixelXSize;
961                         if (dfPixelYSize < dfGlobalPixelYSize)
962                             dfGlobalPixelYSize = dfPixelYSize;
963                     }
964 
965                     nValidFrames ++;
966 
967                     if (bLookForSubDataset)
968                     {
969                         FrameDesc frameDesc;
970                         frameDesc.pszName = pszFrameName;
971                         frameDesc.pszPath = pszFramePath;
972                         frameDesc.nScale = nScale;
973                         frameDesc.nZone = nZone;
974                         aosFrameDesc.push_back(frameDesc);
975                     }
976                 }
977 
978                 if (bLookForSubDataset)
979                 {
980                     delete poDS;
981                     if (nValidFrames == 0)
982                         return nullptr;
983                     return ECRGTOCSubDataset::Build(pszProductTitle,
984                                                     pszDiscId,
985                                                     nScale,
986                                                     nCountSubDataset,
987                                                     pszTOCFilename,
988                                                     aosFrameDesc,
989                                                     dfGlobalMinX,
990                                                     dfGlobalMinY,
991                                                     dfGlobalMaxX,
992                                                     dfGlobalMaxY,
993                                                     dfGlobalPixelXSize,
994                                                     dfGlobalPixelYSize);
995                 }
996 
997                 if (nValidFrames)
998                 {
999                     poDS->AddSubDataset(pszOpenInfoFilename,
1000                                         pszProductTitle, pszDiscId, pszSize);
1001                     nSubDatasets ++;
1002                 }
1003             }
1004         }
1005     }
1006 
1007     if (!bGlobalExtentValid)
1008     {
1009         delete poDS;
1010         return nullptr;
1011     }
1012 
1013     if (nSubDatasets == 1)
1014     {
1015         const char* pszSubDatasetName = CSLFetchNameValue(
1016             poDS->GetMetadata("SUBDATASETS"), "SUBDATASET_1_NAME");
1017         GDALOpenInfo oOpenInfo(pszSubDatasetName, GA_ReadOnly);
1018         delete poDS;
1019         GDALDataset* poRetDS = Open(&oOpenInfo);
1020         if (poRetDS)
1021             poRetDS->SetDescription(pszOpenInfoFilename);
1022         return poRetDS;
1023     }
1024 
1025     poDS->adfGeoTransform[0] = dfGlobalMinX;
1026     poDS->adfGeoTransform[1] = dfGlobalPixelXSize;
1027     poDS->adfGeoTransform[2] = 0.0;
1028     poDS->adfGeoTransform[3] = dfGlobalMaxY;
1029     poDS->adfGeoTransform[4] = 0.0;
1030     poDS->adfGeoTransform[5] = - dfGlobalPixelYSize;
1031 
1032     poDS->nRasterXSize = static_cast<int>(
1033         0.5 + (dfGlobalMaxX - dfGlobalMinX) / dfGlobalPixelXSize);
1034     poDS->nRasterYSize = static_cast<int>(
1035         0.5 + (dfGlobalMaxY - dfGlobalMinY) / dfGlobalPixelYSize);
1036 
1037 /* -------------------------------------------------------------------- */
1038 /*      Initialize any PAM information.                                 */
1039 /* -------------------------------------------------------------------- */
1040     poDS->TryLoadXML();
1041 
1042     return poDS;
1043 }
1044 
1045 /************************************************************************/
1046 /*                              Identify()                              */
1047 /************************************************************************/
1048 
1049 int ECRGTOCDataset::Identify( GDALOpenInfo * poOpenInfo )
1050 
1051 {
1052     const char *pszFilename = poOpenInfo->pszFilename;
1053 
1054 /* -------------------------------------------------------------------- */
1055 /*      Is this a sub-dataset selector? If so, it is obviously ECRGTOC. */
1056 /* -------------------------------------------------------------------- */
1057     if( STARTS_WITH_CI(pszFilename, "ECRG_TOC_ENTRY:"))
1058         return TRUE;
1059 
1060 /* -------------------------------------------------------------------- */
1061 /*  First we check to see if the file has the expected header           */
1062 /*  bytes.                                                              */
1063 /* -------------------------------------------------------------------- */
1064     const char *pabyHeader
1065         = reinterpret_cast<const char *>( poOpenInfo->pabyHeader );
1066     if( pabyHeader == nullptr )
1067         return FALSE;
1068 
1069     if ( strstr(pabyHeader, "<Table_of_Contents") != nullptr &&
1070          strstr(pabyHeader, "<file_header ") != nullptr)
1071         return TRUE;
1072 
1073     if ( strstr(pabyHeader, "<!DOCTYPE Table_of_Contents [") != nullptr)
1074         return TRUE;
1075 
1076     return FALSE;
1077 }
1078 
1079 /************************************************************************/
1080 /*                                Open()                                */
1081 /************************************************************************/
1082 
1083 GDALDataset *ECRGTOCDataset::Open( GDALOpenInfo * poOpenInfo )
1084 
1085 {
1086     if( !Identify( poOpenInfo ) )
1087         return nullptr;
1088 
1089     const char *pszFilename = poOpenInfo->pszFilename;
1090     CPLString osFilename;
1091     CPLString osProduct, osDiscId, osScale;
1092 
1093     if( STARTS_WITH_CI(pszFilename, "ECRG_TOC_ENTRY:") )
1094     {
1095         pszFilename += strlen("ECRG_TOC_ENTRY:");
1096 
1097         /* PRODUCT:DISK:SCALE:FILENAME (or PRODUCT:DISK:FILENAME historically) */
1098         /* with FILENAME potentially C:\BLA... */
1099         char** papszTokens = CSLTokenizeString2(pszFilename, ":", 0);
1100         int nTokens = CSLCount(papszTokens);
1101         if( nTokens != 3 && nTokens != 4 && nTokens != 5 )
1102         {
1103             CSLDestroy(papszTokens);
1104             return nullptr;
1105         }
1106 
1107         osProduct = papszTokens[0];
1108         osDiscId = papszTokens[1];
1109 
1110         if( nTokens == 3 )
1111             osFilename = papszTokens[2];
1112         else if( nTokens == 4 )
1113         {
1114             if( strlen(papszTokens[2]) == 1 &&
1115                 (papszTokens[3][0] == '\\' ||
1116                  papszTokens[3][0] == '/') )
1117             {
1118                 osFilename = papszTokens[2];
1119                 osFilename += ":";
1120                 osFilename += papszTokens[3];
1121             }
1122             else
1123             {
1124                 osScale = papszTokens[2];
1125                 osFilename = papszTokens[3];
1126             }
1127         }
1128         else if( nTokens == 5 &&
1129                 strlen(papszTokens[3]) == 1 &&
1130                 (papszTokens[4][0] == '\\' ||
1131                  papszTokens[4][0] == '/') )
1132         {
1133             osScale = papszTokens[2];
1134             osFilename = papszTokens[3];
1135             osFilename += ":";
1136             osFilename += papszTokens[4];
1137         }
1138         else
1139         {
1140             CSLDestroy(papszTokens);
1141             return nullptr;
1142         }
1143 
1144         CSLDestroy(papszTokens);
1145         pszFilename = osFilename.c_str();
1146     }
1147 
1148 /* -------------------------------------------------------------------- */
1149 /*      Parse the XML file                                              */
1150 /* -------------------------------------------------------------------- */
1151     CPLXMLNode* psXML = CPLParseXMLFile(pszFilename);
1152     if (psXML == nullptr)
1153     {
1154         return nullptr;
1155     }
1156 
1157     GDALDataset* poDS = Build( pszFilename, psXML, osProduct, osDiscId,
1158                                osScale,
1159                                poOpenInfo->pszFilename);
1160     CPLDestroyXMLNode(psXML);
1161 
1162     if (poDS && poOpenInfo->eAccess == GA_Update)
1163     {
1164         CPLError(CE_Failure, CPLE_NotSupported,
1165                  "ECRGTOC driver does not support update mode");
1166         delete poDS;
1167         return nullptr;
1168     }
1169 
1170     return poDS;
1171 }
1172 
1173 /************************************************************************/
1174 /*                         GDALRegister_ECRGTOC()                       */
1175 /************************************************************************/
1176 
1177 void GDALRegister_ECRGTOC()
1178 
1179 {
1180     if( GDALGetDriverByName( "ECRGTOC" ) != nullptr )
1181         return;
1182 
1183     GDALDriver *poDriver = new GDALDriver();
1184 
1185     poDriver->SetDescription( "ECRGTOC" );
1186     poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
1187     poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "ECRG TOC format" );
1188 
1189     poDriver->pfnIdentify = ECRGTOCDataset::Identify;
1190     poDriver->pfnOpen = ECRGTOCDataset::Open;
1191 
1192     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1193                                "drivers/raster/ecrgtoc.html" );
1194     poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "xml" );
1195     poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1196     poDriver->SetMetadataItem( GDAL_DMD_SUBDATASETS, "YES" );
1197 
1198     GetGDALDriverManager()->RegisterDriver( poDriver );
1199 }
1200