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