1 /******************************************************************************
2 *
3 * Project: Erdas Imagine (.img) Translator
4 * Purpose: Supporting functions for HFA (.img) ... main (C callable) API
5 * that is not dependent on GDAL (just CPL).
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 *
8 ******************************************************************************
9 * Copyright (c) 1999, Intergraph Corporation
10 * Copyright (c) 2007-2011, Even Rouault <even dot rouault at spatialys.com>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ******************************************************************************
30 *
31 * hfaopen.cpp
32 *
33 * Supporting routines for reading Erdas Imagine (.imf) Hierarchical
34 * File Architecture files. This is intended to be a library independent
35 * of the GDAL core, but dependent on the Common Portability Library.
36 *
37 */
38
39 #include "cpl_port.h"
40 #include "hfa_p.h"
41
42 #include <cerrno>
43 #include <climits>
44 #include <cmath>
45 #include <cstddef>
46 #include <cstdio>
47 #include <cstdlib>
48 #include <cstring>
49 #if HAVE_FCNTL_H
50 # include <fcntl.h>
51 #endif
52 #include <algorithm>
53 #include <memory>
54 #include <string>
55 #include <vector>
56
57 #include "cpl_conv.h"
58 #include "cpl_error.h"
59 #include "cpl_string.h"
60 #include "cpl_vsi.h"
61 #include "gdal_priv.h"
62 #include "hfa.h"
63
64 CPL_CVSID("$Id: hfaopen.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
65
66 static const char * const apszAuxMetadataItems[] = {
67 // node/entry field_name metadata_key type
68 "Statistics", "dminimum", "STATISTICS_MINIMUM", "Esta_Statistics",
69 "Statistics", "dmaximum", "STATISTICS_MAXIMUM", "Esta_Statistics",
70 "Statistics", "dmean", "STATISTICS_MEAN", "Esta_Statistics",
71 "Statistics", "dmedian", "STATISTICS_MEDIAN", "Esta_Statistics",
72 "Statistics", "dmode", "STATISTICS_MODE", "Esta_Statistics",
73 "Statistics", "dstddev", "STATISTICS_STDDEV", "Esta_Statistics",
74 "HistogramParameters", "lBinFunction.numBins", "STATISTICS_HISTONUMBINS","Eimg_StatisticsParameters830",
75 "HistogramParameters", "dBinFunction.minLimit", "STATISTICS_HISTOMIN", "Eimg_StatisticsParameters830",
76 "HistogramParameters", "dBinFunction.maxLimit", "STATISTICS_HISTOMAX", "Eimg_StatisticsParameters830",
77 "StatisticsParameters", "lSkipFactorX", "STATISTICS_SKIPFACTORX", "",
78 "StatisticsParameters", "lSkipFactorY", "STATISTICS_SKIPFACTORY", "",
79 "StatisticsParameters", "dExcludedValues", "STATISTICS_EXCLUDEDVALUES","",
80 "", "elayerType", "LAYER_TYPE", "",
81 "RRDInfoList", "salgorithm.string", "OVERVIEWS_ALGORITHM", "Emif_String",
82 nullptr
83 };
84
GetHFAAuxMetaDataList()85 const char * const * GetHFAAuxMetaDataList()
86 {
87 return apszAuxMetadataItems;
88 }
89
90 /************************************************************************/
91 /* HFAGetDictionary() */
92 /************************************************************************/
93
HFAGetDictionary(HFAHandle hHFA)94 static char * HFAGetDictionary( HFAHandle hHFA )
95
96 {
97 int nDictMax = 100;
98 char *pszDictionary = static_cast<char *>(CPLMalloc(nDictMax));
99 int nDictSize = 0;
100
101 if( VSIFSeekL(hHFA->fp, hHFA->nDictionaryPos, SEEK_SET) < 0 )
102 {
103 pszDictionary[nDictSize] = '\0';
104 return pszDictionary;
105 }
106
107 while( true )
108 {
109 if( nDictSize >= nDictMax - 1 )
110 {
111 nDictMax = nDictSize * 2 + 100;
112 pszDictionary =
113 static_cast<char *>(CPLRealloc(pszDictionary, nDictMax));
114 }
115
116 if( VSIFReadL(pszDictionary + nDictSize, 1, 1, hHFA->fp) < 1 ||
117 pszDictionary[nDictSize] == '\0' ||
118 (nDictSize > 2 && pszDictionary[nDictSize - 2] == ',' &&
119 pszDictionary[nDictSize-1] == '.') )
120 break;
121
122 nDictSize++;
123 }
124
125 pszDictionary[nDictSize] = '\0';
126
127 return pszDictionary;
128 }
129
130 /************************************************************************/
131 /* HFAOpen() */
132 /************************************************************************/
133
HFAOpen(const char * pszFilename,const char * pszAccess)134 HFAHandle HFAOpen( const char *pszFilename, const char *pszAccess )
135
136 {
137 VSILFILE *fp = VSIFOpenL(
138 pszFilename,
139 (EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb")) ? "rb" : "r+b");
140
141 // Should this be changed to use some sort of CPLFOpen() which will
142 // set the error?
143 if( fp == nullptr )
144 {
145 CPLError(CE_Failure, CPLE_OpenFailed, "File open of %s failed.",
146 pszFilename);
147
148 return nullptr;
149 }
150
151 // Read and verify the header.
152 char szHeader[16] = {};
153 if( VSIFReadL(szHeader, 16, 1, fp) < 1 )
154 {
155 CPLError(CE_Failure, CPLE_AppDefined,
156 "Attempt to read 16 byte header failed for\n%s.", pszFilename);
157 CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
158 return nullptr;
159 }
160
161 if( !STARTS_WITH_CI(szHeader, "EHFA_HEADER_TAG") )
162 {
163 CPLError(CE_Failure, CPLE_AppDefined,
164 "File %s is not an Imagine HFA file ... header wrong.",
165 pszFilename);
166 CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
167 return nullptr;
168 }
169
170 // Create the HFAInfo_t.
171 HFAInfo_t *psInfo =
172 static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
173
174 psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
175 psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
176 psInfo->fp = fp;
177 if( EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb") )
178 psInfo->eAccess = HFA_ReadOnly;
179 else
180 psInfo->eAccess = HFA_Update;
181 psInfo->bTreeDirty = false;
182
183 // Where is the header?
184 GUInt32 nHeaderPos = 0;
185 bool bRet = VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, fp) > 0;
186 HFAStandard(4, &nHeaderPos);
187
188 // Read the header.
189 bRet &= VSIFSeekL(fp, nHeaderPos, SEEK_SET) >= 0;
190
191 bRet &= VSIFReadL(&(psInfo->nVersion), sizeof(GInt32), 1, fp) > 0;
192 HFAStandard(4, &(psInfo->nVersion));
193
194 bRet &= VSIFReadL(szHeader, 4, 1, fp) > 0; // Skip freeList.
195
196 bRet &= VSIFReadL(&(psInfo->nRootPos), sizeof(GInt32), 1, fp) > 0;
197 HFAStandard(4, &(psInfo->nRootPos));
198
199 bRet &= VSIFReadL(&(psInfo->nEntryHeaderLength), sizeof(GInt16), 1, fp) > 0;
200 HFAStandard(2, &(psInfo->nEntryHeaderLength));
201
202 bRet &= VSIFReadL(&(psInfo->nDictionaryPos), sizeof(GInt32), 1, fp) > 0;
203 HFAStandard(4, &(psInfo->nDictionaryPos));
204
205 // Collect file size.
206 bRet &= VSIFSeekL(fp, 0, SEEK_END) >= 0;
207 if( !bRet )
208 {
209 CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
210 CPLFree(psInfo->pszFilename);
211 CPLFree(psInfo->pszPath);
212 CPLFree(psInfo);
213 return nullptr;
214 }
215 psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
216
217 // Instantiate the root entry.
218 psInfo->poRoot = HFAEntry::New(psInfo, psInfo->nRootPos, nullptr, nullptr);
219 if( psInfo->poRoot == nullptr )
220 {
221 CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
222 CPLFree(psInfo->pszFilename);
223 CPLFree(psInfo->pszPath);
224 CPLFree(psInfo);
225 return nullptr;
226 }
227
228 // Read the dictionary.
229 psInfo->pszDictionary = HFAGetDictionary(psInfo);
230 psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
231
232 // Collect band definitions.
233 HFAParseBandInfo(psInfo);
234
235 return psInfo;
236 }
237
238 /************************************************************************/
239 /* HFACreateDependent() */
240 /* */
241 /* Create a .rrd file for the named file if it does not exist, */
242 /* or return the existing dependent if it already exists. */
243 /************************************************************************/
244
HFACreateDependent(HFAInfo_t * psBase)245 HFAInfo_t *HFACreateDependent( HFAInfo_t *psBase )
246
247 {
248 if( psBase->psDependent != nullptr )
249 return psBase->psDependent;
250
251 // Create desired RRD filename.
252 const CPLString oBasename = CPLGetBasename(psBase->pszFilename);
253 const CPLString oRRDFilename =
254 CPLFormFilename(psBase->pszPath, oBasename, "rrd");
255
256 // Does this file already exist? If so, re-use it.
257 VSILFILE *fp = VSIFOpenL(oRRDFilename, "rb");
258 if( fp != nullptr )
259 {
260 CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
261 psBase->psDependent = HFAOpen(oRRDFilename, "rb");
262 // FIXME? this is not going to be reused but recreated...
263 }
264
265 // Otherwise create it now.
266 HFAInfo_t *psDep = HFACreateLL(oRRDFilename);
267 psBase->psDependent = psDep;
268 if( psDep == nullptr )
269 return nullptr;
270
271 /* -------------------------------------------------------------------- */
272 /* Add the DependentFile node with the pointer back to the */
273 /* parent. When working from an .aux file we really want the */
274 /* .rrd to point back to the original file, not the .aux file. */
275 /* -------------------------------------------------------------------- */
276 HFAEntry *poEntry = psBase->poRoot->GetNamedChild("DependentFile");
277 const char *pszDependentFile = nullptr;
278 if( poEntry != nullptr )
279 pszDependentFile = poEntry->GetStringField("dependent.string");
280 if( pszDependentFile == nullptr )
281 pszDependentFile = psBase->pszFilename;
282
283 HFAEntry *poDF = HFAEntry::New(psDep, "DependentFile",
284 "Eimg_DependentFile", psDep->poRoot);
285
286 poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
287 poDF->SetPosition();
288 poDF->SetStringField("dependent.string", pszDependentFile);
289
290 return psDep;
291 }
292
293 /************************************************************************/
294 /* HFAGetDependent() */
295 /************************************************************************/
296
HFAGetDependent(HFAInfo_t * psBase,const char * pszFilename)297 HFAInfo_t *HFAGetDependent(HFAInfo_t *psBase, const char *pszFilename)
298
299 {
300 if( EQUAL(pszFilename, psBase->pszFilename) )
301 return psBase;
302
303 if( psBase->psDependent != nullptr )
304 {
305 if( EQUAL(pszFilename, psBase->psDependent->pszFilename) )
306 return psBase->psDependent;
307 else
308 return nullptr;
309 }
310
311 // Try to open the dependent file.
312 const char *pszMode = psBase->eAccess == HFA_Update ? "r+b" : "rb";
313
314 char *pszDependent =
315 CPLStrdup(CPLFormFilename(psBase->pszPath, pszFilename, nullptr));
316
317 VSILFILE *fp = VSIFOpenL(pszDependent, pszMode);
318 if( fp != nullptr )
319 {
320 CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
321 psBase->psDependent = HFAOpen(pszDependent, pszMode);
322 }
323
324 CPLFree(pszDependent);
325
326 return psBase->psDependent;
327 }
328
329 /************************************************************************/
330 /* HFAParseBandInfo() */
331 /* */
332 /* This is used by HFAOpen() and HFACreate() to initialize the */
333 /* band structures. */
334 /************************************************************************/
335
HFAParseBandInfo(HFAInfo_t * psInfo)336 CPLErr HFAParseBandInfo( HFAInfo_t *psInfo )
337
338 {
339 // Find the first band node.
340 psInfo->nBands = 0;
341 HFAEntry *poNode = psInfo->poRoot->GetChild();
342 while( poNode != nullptr )
343 {
344 if( EQUAL(poNode->GetType(), "Eimg_Layer") &&
345 poNode->GetIntField("width") > 0 &&
346 poNode->GetIntField("height") > 0 )
347 {
348 if( psInfo->nBands == 0 )
349 {
350 psInfo->nXSize = poNode->GetIntField("width");
351 psInfo->nYSize = poNode->GetIntField("height");
352 }
353 else if( poNode->GetIntField("width") != psInfo->nXSize ||
354 poNode->GetIntField("height") != psInfo->nYSize )
355 {
356 return CE_Failure;
357 }
358
359 psInfo->papoBand = static_cast<HFABand **>(CPLRealloc(
360 psInfo->papoBand, sizeof(HFABand *) * (psInfo->nBands + 1)));
361 psInfo->papoBand[psInfo->nBands] = new HFABand(psInfo, poNode);
362 if( psInfo->papoBand[psInfo->nBands]->nWidth == 0 )
363 {
364 delete psInfo->papoBand[psInfo->nBands];
365 return CE_Failure;
366 }
367 psInfo->nBands++;
368 }
369
370 poNode = poNode->GetNext();
371 }
372
373 return CE_None;
374 }
375
376 /************************************************************************/
377 /* HFAClose() */
378 /************************************************************************/
379
HFAClose(HFAHandle hHFA)380 int HFAClose( HFAHandle hHFA )
381
382 {
383 if( hHFA->eAccess == HFA_Update && (hHFA->bTreeDirty ||
384 (hHFA->poDictionary != nullptr &&
385 hHFA->poDictionary->bDictionaryTextDirty)) )
386 HFAFlush(hHFA);
387
388 int nRet = 0;
389 if( hHFA->psDependent != nullptr )
390 {
391 if( HFAClose(hHFA->psDependent) != 0 )
392 nRet = -1;
393 }
394
395 delete hHFA->poRoot;
396
397 if( VSIFCloseL(hHFA->fp) != 0 )
398 nRet = -1;
399
400 if( hHFA->poDictionary != nullptr )
401 delete hHFA->poDictionary;
402
403 CPLFree(hHFA->pszDictionary);
404 CPLFree(hHFA->pszFilename);
405 CPLFree(hHFA->pszIGEFilename);
406 CPLFree(hHFA->pszPath);
407
408 for( int i = 0; i < hHFA->nBands; i++ )
409 {
410 delete hHFA->papoBand[i];
411 }
412
413 CPLFree(hHFA->papoBand);
414
415 if( hHFA->pProParameters != nullptr )
416 {
417 Eprj_ProParameters *psProParams = (Eprj_ProParameters *)
418 hHFA->pProParameters;
419
420 CPLFree(psProParams->proExeName);
421 CPLFree(psProParams->proName);
422 CPLFree(psProParams->proSpheroid.sphereName);
423
424 CPLFree(psProParams);
425 }
426
427 if( hHFA->pDatum != nullptr )
428 {
429 CPLFree(((Eprj_Datum *)hHFA->pDatum)->datumname);
430 CPLFree(((Eprj_Datum *)hHFA->pDatum)->gridname);
431 CPLFree(hHFA->pDatum);
432 }
433
434 if( hHFA->pMapInfo != nullptr )
435 {
436 CPLFree(((Eprj_MapInfo *)hHFA->pMapInfo)->proName);
437 CPLFree(((Eprj_MapInfo *)hHFA->pMapInfo)->units);
438 CPLFree(hHFA->pMapInfo);
439 }
440
441 CPLFree(hHFA);
442 return nRet;
443 }
444
445 /************************************************************************/
446 /* HFARemove() */
447 /* Used from HFADelete() function. */
448 /************************************************************************/
449
HFARemove(const char * pszFilename)450 static CPLErr HFARemove( const char *pszFilename )
451
452 {
453 VSIStatBufL sStat;
454
455 if( VSIStatL(pszFilename, &sStat) == 0 && VSI_ISREG(sStat.st_mode) )
456 {
457 if( VSIUnlink(pszFilename) == 0 )
458 return CE_None;
459 else
460 {
461 CPLError(CE_Failure, CPLE_AppDefined,
462 "Attempt to unlink %s failed.", pszFilename);
463 return CE_Failure;
464 }
465 }
466
467 CPLError(CE_Failure, CPLE_AppDefined, "Unable to delete %s, not a file.",
468 pszFilename);
469 return CE_Failure;
470 }
471
472 /************************************************************************/
473 /* HFADelete() */
474 /************************************************************************/
475
HFADelete(const char * pszFilename)476 CPLErr HFADelete( const char *pszFilename )
477
478 {
479 HFAInfo_t *psInfo = HFAOpen(pszFilename, "rb");
480 HFAEntry *poDMS = nullptr;
481 HFAEntry *poLayer = nullptr;
482 HFAEntry *poNode = nullptr;
483
484 if( psInfo != nullptr )
485 {
486 poNode = psInfo->poRoot->GetChild();
487 while( (poNode != nullptr) && (poLayer == nullptr) )
488 {
489 if( EQUAL(poNode->GetType(), "Eimg_Layer") )
490 {
491 poLayer = poNode;
492 }
493 poNode = poNode->GetNext();
494 }
495
496 if( poLayer != nullptr )
497 poDMS = poLayer->GetNamedChild("ExternalRasterDMS");
498
499 if( poDMS )
500 {
501 const char *pszRawFilename =
502 poDMS->GetStringField("fileName.string");
503
504 if( pszRawFilename != nullptr )
505 HFARemove(
506 CPLFormFilename(psInfo->pszPath, pszRawFilename, nullptr));
507 }
508
509 CPL_IGNORE_RET_VAL(HFAClose(psInfo));
510 }
511 return HFARemove(pszFilename);
512 }
513
514 /************************************************************************/
515 /* HFAGetRasterInfo() */
516 /************************************************************************/
517
HFAGetRasterInfo(HFAHandle hHFA,int * pnXSize,int * pnYSize,int * pnBands)518 CPLErr HFAGetRasterInfo( HFAHandle hHFA, int *pnXSize, int *pnYSize,
519 int *pnBands )
520
521 {
522 if( pnXSize != nullptr )
523 *pnXSize = hHFA->nXSize;
524 if( pnYSize != nullptr )
525 *pnYSize = hHFA->nYSize;
526 if( pnBands != nullptr )
527 *pnBands = hHFA->nBands;
528 return CE_None;
529 }
530
531 /************************************************************************/
532 /* HFAGetBandInfo() */
533 /************************************************************************/
534
HFAGetBandInfo(HFAHandle hHFA,int nBand,EPTType * peDataType,int * pnBlockXSize,int * pnBlockYSize,int * pnCompressionType)535 CPLErr HFAGetBandInfo( HFAHandle hHFA, int nBand, EPTType *peDataType,
536 int *pnBlockXSize, int *pnBlockYSize,
537 int *pnCompressionType )
538
539 {
540 if( nBand < 0 || nBand > hHFA->nBands )
541 {
542 CPLAssert(false);
543 return CE_Failure;
544 }
545
546 HFABand *poBand = hHFA->papoBand[nBand - 1];
547
548 if( peDataType != nullptr )
549 *peDataType = poBand->eDataType;
550
551 if( pnBlockXSize != nullptr )
552 *pnBlockXSize = poBand->nBlockXSize;
553
554 if( pnBlockYSize != nullptr )
555 *pnBlockYSize = poBand->nBlockYSize;
556
557 // Get compression code from RasterDMS.
558 if( pnCompressionType != nullptr )
559 {
560 *pnCompressionType = 0;
561
562 HFAEntry *poDMS = poBand->poNode->GetNamedChild("RasterDMS");
563
564 if( poDMS != nullptr )
565 *pnCompressionType = poDMS->GetIntField("compressionType");
566 }
567
568 return CE_None;
569 }
570
571 /************************************************************************/
572 /* HFAGetBandNoData() */
573 /* */
574 /* returns TRUE if value is set, otherwise FALSE. */
575 /************************************************************************/
576
HFAGetBandNoData(HFAHandle hHFA,int nBand,double * pdfNoData)577 int HFAGetBandNoData( HFAHandle hHFA, int nBand, double *pdfNoData )
578
579 {
580 if( nBand < 0 || nBand > hHFA->nBands )
581 {
582 CPLAssert(false);
583 return CE_Failure;
584 }
585
586 HFABand *poBand = hHFA->papoBand[nBand - 1];
587
588 if( !poBand->bNoDataSet && poBand->nOverviews > 0 )
589 {
590 poBand = poBand->papoOverviews[0];
591 if( poBand == nullptr )
592 return FALSE;
593 }
594
595 *pdfNoData = poBand->dfNoData;
596 return poBand->bNoDataSet;
597 }
598
599 /************************************************************************/
600 /* HFASetBandNoData() */
601 /* */
602 /* attempts to set a no-data value on the given band */
603 /************************************************************************/
604
HFASetBandNoData(HFAHandle hHFA,int nBand,double dfValue)605 CPLErr HFASetBandNoData( HFAHandle hHFA, int nBand, double dfValue )
606
607 {
608 if( nBand < 0 || nBand > hHFA->nBands )
609 {
610 CPLAssert(false);
611 return CE_Failure;
612 }
613
614 HFABand *poBand = hHFA->papoBand[nBand - 1];
615
616 return poBand->SetNoDataValue( dfValue );
617 }
618
619 /************************************************************************/
620 /* HFAGetOverviewCount() */
621 /************************************************************************/
622
HFAGetOverviewCount(HFAHandle hHFA,int nBand)623 int HFAGetOverviewCount( HFAHandle hHFA, int nBand )
624
625 {
626 if( nBand < 0 || nBand > hHFA->nBands )
627 {
628 CPLAssert(false);
629 return CE_Failure;
630 }
631
632 HFABand *poBand = hHFA->papoBand[nBand - 1];
633 poBand->LoadOverviews();
634
635 return poBand->nOverviews;
636 }
637
638 /************************************************************************/
639 /* HFAGetOverviewInfo() */
640 /************************************************************************/
641
HFAGetOverviewInfo(HFAHandle hHFA,int nBand,int iOverview,int * pnXSize,int * pnYSize,int * pnBlockXSize,int * pnBlockYSize,EPTType * peHFADataType)642 CPLErr HFAGetOverviewInfo( HFAHandle hHFA, int nBand, int iOverview,
643 int *pnXSize, int *pnYSize,
644 int *pnBlockXSize, int *pnBlockYSize,
645 EPTType *peHFADataType )
646
647 {
648 if( nBand < 0 || nBand > hHFA->nBands )
649 {
650 CPLAssert(false);
651 return CE_Failure;
652 }
653
654 HFABand *poBand = hHFA->papoBand[nBand - 1];
655 poBand->LoadOverviews();
656
657 if( iOverview < 0 || iOverview >= poBand->nOverviews )
658 {
659 CPLAssert(false);
660 return CE_Failure;
661 }
662 poBand = poBand->papoOverviews[iOverview];
663 if( poBand == nullptr )
664 {
665 return CE_Failure;
666 }
667
668 if( pnXSize != nullptr )
669 *pnXSize = poBand->nWidth;
670
671 if( pnYSize != nullptr )
672 *pnYSize = poBand->nHeight;
673
674 if( pnBlockXSize != nullptr )
675 *pnBlockXSize = poBand->nBlockXSize;
676
677 if( pnBlockYSize != nullptr )
678 *pnBlockYSize = poBand->nBlockYSize;
679
680 if( peHFADataType != nullptr )
681 *peHFADataType = poBand->eDataType;
682
683 return CE_None;
684 }
685
686 /************************************************************************/
687 /* HFAGetRasterBlock() */
688 /************************************************************************/
689
HFAGetRasterBlock(HFAHandle hHFA,int nBand,int nXBlock,int nYBlock,void * pData)690 CPLErr HFAGetRasterBlock( HFAHandle hHFA, int nBand,
691 int nXBlock, int nYBlock, void *pData )
692
693 {
694 return HFAGetRasterBlockEx(hHFA, nBand, nXBlock, nYBlock, pData, -1);
695 }
696
697 /************************************************************************/
698 /* HFAGetRasterBlockEx() */
699 /************************************************************************/
700
HFAGetRasterBlockEx(HFAHandle hHFA,int nBand,int nXBlock,int nYBlock,void * pData,int nDataSize)701 CPLErr HFAGetRasterBlockEx( HFAHandle hHFA, int nBand,
702 int nXBlock, int nYBlock, void *pData,
703 int nDataSize )
704
705 {
706 if( nBand < 1 || nBand > hHFA->nBands )
707 return CE_Failure;
708
709 return hHFA->papoBand[nBand - 1]->
710 GetRasterBlock(nXBlock, nYBlock, pData,nDataSize);
711 }
712
713 /************************************************************************/
714 /* HFAGetOverviewRasterBlock() */
715 /************************************************************************/
716
HFAGetOverviewRasterBlock(HFAHandle hHFA,int nBand,int iOverview,int nXBlock,int nYBlock,void * pData)717 CPLErr HFAGetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
718 int nXBlock, int nYBlock, void *pData )
719
720 {
721 return HFAGetOverviewRasterBlockEx(hHFA, nBand, iOverview, nXBlock, nYBlock,
722 pData, -1);
723 }
724
725 /************************************************************************/
726 /* HFAGetOverviewRasterBlockEx() */
727 /************************************************************************/
728
HFAGetOverviewRasterBlockEx(HFAHandle hHFA,int nBand,int iOverview,int nXBlock,int nYBlock,void * pData,int nDataSize)729 CPLErr HFAGetOverviewRasterBlockEx( HFAHandle hHFA, int nBand, int iOverview,
730 int nXBlock, int nYBlock, void *pData,
731 int nDataSize )
732
733 {
734 if( nBand < 1 || nBand > hHFA->nBands )
735 return CE_Failure;
736
737 if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand - 1]->nOverviews )
738 return CE_Failure;
739
740 return hHFA->papoBand[nBand - 1]->papoOverviews[iOverview]->
741 GetRasterBlock(nXBlock, nYBlock, pData, nDataSize);
742 }
743
744 /************************************************************************/
745 /* HFASetRasterBlock() */
746 /************************************************************************/
747
HFASetRasterBlock(HFAHandle hHFA,int nBand,int nXBlock,int nYBlock,void * pData)748 CPLErr HFASetRasterBlock( HFAHandle hHFA, int nBand,
749 int nXBlock, int nYBlock, void *pData )
750
751 {
752 if( nBand < 1 || nBand > hHFA->nBands )
753 return CE_Failure;
754
755 return hHFA->papoBand[nBand - 1]->SetRasterBlock(nXBlock, nYBlock, pData);
756 }
757
758 /************************************************************************/
759 /* HFASetRasterBlock() */
760 /************************************************************************/
761
HFASetOverviewRasterBlock(HFAHandle hHFA,int nBand,int iOverview,int nXBlock,int nYBlock,void * pData)762 CPLErr HFASetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
763 int nXBlock, int nYBlock, void *pData )
764
765 {
766 if( nBand < 1 || nBand > hHFA->nBands )
767 return CE_Failure;
768
769 if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand - 1]->nOverviews )
770 return CE_Failure;
771
772 return hHFA->papoBand[nBand - 1]->papoOverviews[iOverview]->
773 SetRasterBlock(nXBlock, nYBlock, pData);
774 }
775
776 /************************************************************************/
777 /* HFAGetBandName() */
778 /************************************************************************/
779
HFAGetBandName(HFAHandle hHFA,int nBand)780 const char *HFAGetBandName( HFAHandle hHFA, int nBand )
781 {
782 if( nBand < 1 || nBand > hHFA->nBands )
783 return "";
784
785 return hHFA->papoBand[nBand - 1]->GetBandName();
786 }
787
788 /************************************************************************/
789 /* HFASetBandName() */
790 /************************************************************************/
791
HFASetBandName(HFAHandle hHFA,int nBand,const char * pszName)792 void HFASetBandName( HFAHandle hHFA, int nBand, const char *pszName )
793 {
794 if( nBand < 1 || nBand > hHFA->nBands )
795 return;
796
797 hHFA->papoBand[nBand - 1]->SetBandName(pszName);
798 }
799
800 /************************************************************************/
801 /* HFAGetDataTypeBits() */
802 /************************************************************************/
803
HFAGetDataTypeBits(EPTType eDataType)804 int HFAGetDataTypeBits( EPTType eDataType )
805
806 {
807 switch( eDataType )
808 {
809 case EPT_u1:
810 return 1;
811
812 case EPT_u2:
813 return 2;
814
815 case EPT_u4:
816 return 4;
817
818 case EPT_u8:
819 case EPT_s8:
820 return 8;
821
822 case EPT_u16:
823 case EPT_s16:
824 return 16;
825
826 case EPT_u32:
827 case EPT_s32:
828 case EPT_f32:
829 return 32;
830
831 case EPT_f64:
832 case EPT_c64:
833 return 64;
834
835 case EPT_c128:
836 return 128;
837 }
838
839 CPLAssert(false);
840 return 1;
841 }
842
843 /************************************************************************/
844 /* HFAGetDataTypeName() */
845 /************************************************************************/
846
HFAGetDataTypeName(EPTType eDataType)847 const char *HFAGetDataTypeName( EPTType eDataType )
848
849 {
850 switch( eDataType )
851 {
852 case EPT_u1:
853 return "u1";
854
855 case EPT_u2:
856 return "u2";
857
858 case EPT_u4:
859 return "u4";
860
861 case EPT_u8:
862 return "u8";
863
864 case EPT_s8:
865 return "s8";
866
867 case EPT_u16:
868 return "u16";
869
870 case EPT_s16:
871 return "s16";
872
873 case EPT_u32:
874 return "u32";
875
876 case EPT_s32:
877 return "s32";
878
879 case EPT_f32:
880 return "f32";
881
882 case EPT_f64:
883 return "f64";
884
885 case EPT_c64:
886 return "c64";
887
888 case EPT_c128:
889 return "c128";
890
891 default:
892 CPLAssert(false);
893 return "unknown";
894 }
895 }
896
897 /************************************************************************/
898 /* HFAGetMapInfo() */
899 /************************************************************************/
900
HFAGetMapInfo(HFAHandle hHFA)901 const Eprj_MapInfo *HFAGetMapInfo( HFAHandle hHFA )
902
903 {
904 if( hHFA->nBands < 1 )
905 return nullptr;
906
907 // Do we already have it?
908 if( hHFA->pMapInfo != nullptr )
909 return (Eprj_MapInfo *)hHFA->pMapInfo;
910
911 // Get the HFA node. If we don't find it under the usual name
912 // we search for any node of the right type (#3338).
913 HFAEntry *poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild("Map_Info");
914 if( poMIEntry == nullptr )
915 {
916 for( HFAEntry *poChild = hHFA->papoBand[0]->poNode->GetChild();
917 poChild != nullptr && poMIEntry == nullptr;
918 poChild = poChild->GetNext() )
919 {
920 if( EQUAL(poChild->GetType(), "Eprj_MapInfo") )
921 poMIEntry = poChild;
922 }
923 }
924
925 if( poMIEntry == nullptr )
926 {
927 return nullptr;
928 }
929
930 // Allocate the structure.
931 Eprj_MapInfo *psMapInfo =
932 static_cast<Eprj_MapInfo *>(CPLCalloc(sizeof(Eprj_MapInfo), 1));
933
934 // Fetch the fields.
935 psMapInfo->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
936
937 psMapInfo->upperLeftCenter.x =
938 poMIEntry->GetDoubleField("upperLeftCenter.x");
939 psMapInfo->upperLeftCenter.y =
940 poMIEntry->GetDoubleField("upperLeftCenter.y");
941
942 psMapInfo->lowerRightCenter.x =
943 poMIEntry->GetDoubleField("lowerRightCenter.x");
944 psMapInfo->lowerRightCenter.y =
945 poMIEntry->GetDoubleField("lowerRightCenter.y");
946
947 CPLErr eErr = CE_None;
948 psMapInfo->pixelSize.width =
949 poMIEntry->GetDoubleField("pixelSize.width", &eErr);
950 psMapInfo->pixelSize.height =
951 poMIEntry->GetDoubleField("pixelSize.height", &eErr);
952
953 // The following is basically a hack to get files with
954 // non-standard MapInfo's that misname the pixelSize fields. (#3338)
955 if( eErr != CE_None )
956 {
957 psMapInfo->pixelSize.width = poMIEntry->GetDoubleField("pixelSize.x");
958 psMapInfo->pixelSize.height = poMIEntry->GetDoubleField("pixelSize.y");
959 }
960
961 psMapInfo->units = CPLStrdup(poMIEntry->GetStringField("units"));
962
963 hHFA->pMapInfo = (void *)psMapInfo;
964
965 return psMapInfo;
966 }
967
968 /************************************************************************/
969 /* HFAInvGeoTransform() */
970 /************************************************************************/
971
HFAInvGeoTransform(const double * gt_in,double * gt_out)972 static bool HFAInvGeoTransform( const double *gt_in, double *gt_out )
973
974 {
975 // Assume a 3rd row that is [1 0 0].
976 // Compute determinate.
977 const double det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
978
979 if( fabs(det) < 1.0e-15 )
980 return false;
981
982 const double inv_det = 1.0 / det;
983
984 // Compute adjoint, and divide by determinate.
985 gt_out[1] = gt_in[5] * inv_det;
986 gt_out[4] = -gt_in[4] * inv_det;
987
988 gt_out[2] = -gt_in[2] * inv_det;
989 gt_out[5] = gt_in[1] * inv_det;
990
991 gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
992 gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
993
994 return true;
995 }
996
997 /************************************************************************/
998 /* HFAGetGeoTransform() */
999 /************************************************************************/
1000
HFAGetGeoTransform(HFAHandle hHFA,double * padfGeoTransform)1001 int HFAGetGeoTransform( HFAHandle hHFA, double *padfGeoTransform )
1002
1003 {
1004 const Eprj_MapInfo *psMapInfo = HFAGetMapInfo(hHFA);
1005
1006 padfGeoTransform[0] = 0.0;
1007 padfGeoTransform[1] = 1.0;
1008 padfGeoTransform[2] = 0.0;
1009 padfGeoTransform[3] = 0.0;
1010 padfGeoTransform[4] = 0.0;
1011 padfGeoTransform[5] = 1.0;
1012
1013 // Simple (north up) MapInfo approach.
1014 if( psMapInfo != nullptr )
1015 {
1016 padfGeoTransform[0] =
1017 psMapInfo->upperLeftCenter.x - psMapInfo->pixelSize.width * 0.5;
1018 padfGeoTransform[1] = psMapInfo->pixelSize.width;
1019 if( padfGeoTransform[1] == 0.0 )
1020 padfGeoTransform[1] = 1.0;
1021 padfGeoTransform[2] = 0.0;
1022 if( psMapInfo->upperLeftCenter.y >= psMapInfo->lowerRightCenter.y )
1023 padfGeoTransform[5] = -psMapInfo->pixelSize.height;
1024 else
1025 padfGeoTransform[5] = psMapInfo->pixelSize.height;
1026 if( padfGeoTransform[5] == 0.0 )
1027 padfGeoTransform[5] = 1.0;
1028
1029 padfGeoTransform[3] =
1030 psMapInfo->upperLeftCenter.y - padfGeoTransform[5] * 0.5;
1031 padfGeoTransform[4] = 0.0;
1032
1033 // Special logic to fixup odd angular units.
1034 if( EQUAL(psMapInfo->units, "ds") )
1035 {
1036 padfGeoTransform[0] /= 3600.0;
1037 padfGeoTransform[1] /= 3600.0;
1038 padfGeoTransform[2] /= 3600.0;
1039 padfGeoTransform[3] /= 3600.0;
1040 padfGeoTransform[4] /= 3600.0;
1041 padfGeoTransform[5] /= 3600.0;
1042 }
1043
1044 return TRUE;
1045 }
1046
1047 // Try for a MapToPixelXForm affine polynomial supporting
1048 // rotated and sheared affine transformations.
1049 if( hHFA->nBands == 0 )
1050 return FALSE;
1051
1052 HFAEntry *poXForm0 =
1053 hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
1054
1055 if( poXForm0 == nullptr )
1056 return FALSE;
1057
1058 if( poXForm0->GetIntField("order") != 1 ||
1059 poXForm0->GetIntField("numdimtransform") != 2 ||
1060 poXForm0->GetIntField("numdimpolynomial") != 2 ||
1061 poXForm0->GetIntField("termcount") != 3 )
1062 return FALSE;
1063
1064 // Verify that there aren't any further xform steps.
1065 if( hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm1")
1066 != nullptr )
1067 return FALSE;
1068
1069 // We should check that the exponent list is 0 0 1 0 0 1, but
1070 // we don't because we are lazy.
1071
1072 // Fetch geotransform values.
1073 double adfXForm[6] = {
1074 poXForm0->GetDoubleField("polycoefvector[0]"),
1075 poXForm0->GetDoubleField("polycoefmtx[0]"),
1076 poXForm0->GetDoubleField("polycoefmtx[2]"),
1077 poXForm0->GetDoubleField("polycoefvector[1]"),
1078 poXForm0->GetDoubleField("polycoefmtx[1]"),
1079 poXForm0->GetDoubleField("polycoefmtx[3]")
1080 };
1081
1082 // Invert.
1083
1084 if( !HFAInvGeoTransform(adfXForm, padfGeoTransform) )
1085 memset(padfGeoTransform, 0, 6 * sizeof(double));
1086
1087 // Adjust origin from center of top left pixel to top left corner
1088 // of top left pixel.
1089 padfGeoTransform[0] -= padfGeoTransform[1] * 0.5;
1090 padfGeoTransform[0] -= padfGeoTransform[2] * 0.5;
1091 padfGeoTransform[3] -= padfGeoTransform[4] * 0.5;
1092 padfGeoTransform[3] -= padfGeoTransform[5] * 0.5;
1093
1094 return TRUE;
1095 }
1096
1097 /************************************************************************/
1098 /* HFASetMapInfo() */
1099 /************************************************************************/
1100
HFASetMapInfo(HFAHandle hHFA,const Eprj_MapInfo * poMapInfo)1101 CPLErr HFASetMapInfo( HFAHandle hHFA, const Eprj_MapInfo *poMapInfo )
1102
1103 {
1104 // Loop over bands, setting information on each one.
1105 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1106 {
1107 // Create a new Map_Info if there isn't one present already.
1108 HFAEntry *poMIEntry =
1109 hHFA->papoBand[iBand]->poNode->GetNamedChild("Map_Info");
1110 if( poMIEntry == nullptr )
1111 {
1112 poMIEntry = HFAEntry::New(hHFA, "Map_Info", "Eprj_MapInfo",
1113 hHFA->papoBand[iBand]->poNode);
1114 }
1115
1116 poMIEntry->MarkDirty();
1117
1118 // Ensure we have enough space for all the data.
1119 // TODO(schwehr): Explain 48 and 40 constants.
1120 const int nSize = static_cast<int>(
1121 48 + 40
1122 + strlen(poMapInfo->proName) + 1
1123 + strlen(poMapInfo->units) + 1);
1124
1125 GByte *pabyData = poMIEntry->MakeData(nSize);
1126 memset(pabyData, 0, nSize);
1127
1128 poMIEntry->SetPosition();
1129
1130 // Write the various fields.
1131 poMIEntry->SetStringField("proName", poMapInfo->proName);
1132
1133 poMIEntry->SetDoubleField("upperLeftCenter.x",
1134 poMapInfo->upperLeftCenter.x);
1135 poMIEntry->SetDoubleField("upperLeftCenter.y",
1136 poMapInfo->upperLeftCenter.y);
1137
1138 poMIEntry->SetDoubleField("lowerRightCenter.x",
1139 poMapInfo->lowerRightCenter.x);
1140 poMIEntry->SetDoubleField("lowerRightCenter.y",
1141 poMapInfo->lowerRightCenter.y);
1142
1143 poMIEntry->SetDoubleField("pixelSize.width",
1144 poMapInfo->pixelSize.width);
1145 poMIEntry->SetDoubleField("pixelSize.height",
1146 poMapInfo->pixelSize.height);
1147
1148 poMIEntry->SetStringField("units", poMapInfo->units);
1149 }
1150
1151 return CE_None;
1152 }
1153
1154 /************************************************************************/
1155 /* HFAGetPEString() */
1156 /* */
1157 /* Some files have a ProjectionX node containing the ESRI style */
1158 /* PE_STRING. This function allows fetching from it. */
1159 /************************************************************************/
1160
HFAGetPEString(HFAHandle hHFA)1161 char *HFAGetPEString( HFAHandle hHFA )
1162
1163 {
1164 if( hHFA->nBands == 0 )
1165 return nullptr;
1166
1167 // Get the HFA node.
1168 HFAEntry *poProX = hHFA->papoBand[0]->poNode->GetNamedChild("ProjectionX");
1169 if( poProX == nullptr )
1170 return nullptr;
1171
1172 const char *pszType = poProX->GetStringField("projection.type.string");
1173 if( pszType == nullptr || !EQUAL(pszType, "PE_COORDSYS") )
1174 return nullptr;
1175
1176 // Use a gross hack to scan ahead to the actual projection
1177 // string. We do it this way because we don't have general
1178 // handling for MIFObjects.
1179 GByte *pabyData = poProX->GetData();
1180 int nDataSize = poProX->GetDataSize();
1181
1182 while( nDataSize > 10 &&
1183 !STARTS_WITH_CI((const char *)pabyData, "PE_COORDSYS,.") ) {
1184 pabyData++;
1185 nDataSize--;
1186 }
1187
1188 if( nDataSize < 31 )
1189 return nullptr;
1190
1191 // Skip ahead to the actual string.
1192 pabyData += 30;
1193 // nDataSize -= 30;
1194
1195 return CPLStrdup((const char *)pabyData);
1196 }
1197
1198 /************************************************************************/
1199 /* HFASetPEString() */
1200 /************************************************************************/
1201
HFASetPEString(HFAHandle hHFA,const char * pszPEString)1202 CPLErr HFASetPEString( HFAHandle hHFA, const char *pszPEString )
1203
1204 {
1205 if( !CPLTestBool(CPLGetConfigOption("HFA_WRITE_PE_STRING", "YES")) )
1206 return CE_None;
1207
1208 // Loop over bands, setting information on each one.
1209 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1210 {
1211 // Verify we don't already have the node, since update-in-place
1212 // is likely to be more complicated.
1213 HFAEntry *poProX =
1214 hHFA->papoBand[iBand]->poNode->GetNamedChild("ProjectionX");
1215
1216 // If we are setting an empty string then a missing entry is equivalent.
1217 if( strlen(pszPEString) == 0 && poProX == nullptr )
1218 continue;
1219
1220 // Create the node.
1221 if( poProX == nullptr )
1222 {
1223 poProX = HFAEntry::New(hHFA, "ProjectionX", "Eprj_MapProjection842",
1224 hHFA->papoBand[iBand]->poNode);
1225 if( poProX->GetTypeObject() == nullptr )
1226 return CE_Failure;
1227 }
1228
1229 // Prepare the data area with some extra space just in case.
1230 GByte *pabyData =
1231 poProX->MakeData(static_cast<int>(700 + strlen(pszPEString)));
1232 if( !pabyData )
1233 return CE_Failure;
1234
1235 memset(pabyData, 0, 250 + strlen(pszPEString));
1236
1237 poProX->SetPosition();
1238
1239 poProX->SetStringField("projection.type.string", "PE_COORDSYS");
1240 poProX->SetStringField("projection.MIFDictionary.string",
1241 "{0:pcstring,}Emif_String,{1:x{0:pcstring,}"
1242 "Emif_String,coordSys,}PE_COORDSYS,.");
1243
1244 // Use a gross hack to scan ahead to the actual projection
1245 // string. We do it this way because we don't have general
1246 // handling for MIFObjects
1247 pabyData = poProX->GetData();
1248 int nDataSize = poProX->GetDataSize();
1249 GUInt32 iOffset = poProX->GetDataPos();
1250
1251 while( nDataSize > 10 &&
1252 !STARTS_WITH_CI((const char *) pabyData, "PE_COORDSYS,.") )
1253 {
1254 pabyData++;
1255 nDataSize--;
1256 iOffset++;
1257 }
1258
1259 CPLAssert(nDataSize > static_cast<int>(strlen(pszPEString)) + 10);
1260
1261 pabyData += 14;
1262 iOffset += 14;
1263
1264 // Set the size and offset of the mifobject.
1265 iOffset += 8;
1266
1267 GUInt32 nSize = static_cast<GUInt32>(strlen(pszPEString) + 9);
1268
1269 HFAStandard(4, &nSize);
1270 memcpy(pabyData, &nSize, 4);
1271 pabyData += 4;
1272
1273 HFAStandard(4, &iOffset);
1274 memcpy(pabyData, &iOffset, 4);
1275 pabyData += 4;
1276
1277 // Set the size and offset of the string value.
1278 nSize = static_cast<GUInt32>(strlen(pszPEString) + 1);
1279
1280 HFAStandard(4, &nSize);
1281 memcpy(pabyData, &nSize, 4);
1282 pabyData += 4;
1283
1284 iOffset = 8;
1285 HFAStandard(4, &iOffset);
1286 memcpy(pabyData, &iOffset, 4);
1287 pabyData += 4;
1288
1289 // Place the string itself.
1290 memcpy(pabyData, pszPEString, strlen(pszPEString) + 1);
1291
1292 poProX->SetStringField("title.string", "PE");
1293 }
1294
1295 return CE_None;
1296 }
1297
1298 /************************************************************************/
1299 /* HFAGetProParameters() */
1300 /************************************************************************/
1301
HFAGetProParameters(HFAHandle hHFA)1302 const Eprj_ProParameters *HFAGetProParameters( HFAHandle hHFA )
1303
1304 {
1305 if( hHFA->nBands < 1 )
1306 return nullptr;
1307
1308 // Do we already have it?
1309 if( hHFA->pProParameters != nullptr )
1310 return (Eprj_ProParameters *)hHFA->pProParameters;
1311
1312 // Get the HFA node.
1313 HFAEntry *poMIEntry =
1314 hHFA->papoBand[0]->poNode->GetNamedChild("Projection");
1315 if( poMIEntry == nullptr )
1316 return nullptr;
1317
1318 // Allocate the structure.
1319 Eprj_ProParameters *psProParams = static_cast<Eprj_ProParameters *>(
1320 CPLCalloc(sizeof(Eprj_ProParameters), 1));
1321
1322 // Fetch the fields.
1323 const int proType = poMIEntry->GetIntField("proType");
1324 if( proType != EPRJ_INTERNAL && proType != EPRJ_EXTERNAL )
1325 {
1326 CPLError(CE_Failure, CPLE_AppDefined, "Wrong value for proType");
1327 CPLFree(psProParams);
1328 return nullptr;
1329 }
1330 psProParams->proType = static_cast<Eprj_ProType>(proType);
1331 psProParams->proNumber = poMIEntry->GetIntField("proNumber");
1332 psProParams->proExeName = CPLStrdup(poMIEntry->GetStringField("proExeName"));
1333 psProParams->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
1334 psProParams->proZone = poMIEntry->GetIntField("proZone");
1335
1336 for( int i = 0; i < 15; i++ )
1337 {
1338 char szFieldName[40] = {};
1339
1340 snprintf(szFieldName, sizeof(szFieldName), "proParams[%d]", i);
1341 psProParams->proParams[i] = poMIEntry->GetDoubleField(szFieldName);
1342 }
1343
1344 psProParams->proSpheroid.sphereName =
1345 CPLStrdup(poMIEntry->GetStringField("proSpheroid.sphereName"));
1346 psProParams->proSpheroid.a = poMIEntry->GetDoubleField("proSpheroid.a");
1347 psProParams->proSpheroid.b = poMIEntry->GetDoubleField("proSpheroid.b");
1348 psProParams->proSpheroid.eSquared =
1349 poMIEntry->GetDoubleField("proSpheroid.eSquared");
1350 psProParams->proSpheroid.radius =
1351 poMIEntry->GetDoubleField("proSpheroid.radius");
1352
1353 hHFA->pProParameters = (void *)psProParams;
1354
1355 return psProParams;
1356 }
1357
1358 /************************************************************************/
1359 /* HFASetProParameters() */
1360 /************************************************************************/
1361
HFASetProParameters(HFAHandle hHFA,const Eprj_ProParameters * poPro)1362 CPLErr HFASetProParameters( HFAHandle hHFA, const Eprj_ProParameters *poPro )
1363
1364 {
1365 // Loop over bands, setting information on each one.
1366 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1367 {
1368 // Create a new Projection if there isn't one present already.
1369 HFAEntry *poMIEntry =
1370 hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1371 if( poMIEntry == nullptr )
1372 {
1373 poMIEntry = HFAEntry::New(hHFA, "Projection", "Eprj_ProParameters",
1374 hHFA->papoBand[iBand]->poNode);
1375 }
1376
1377 poMIEntry->MarkDirty();
1378
1379 // Ensure we have enough space for all the data.
1380 // TODO(schwehr): Explain all these constants.
1381 int nSize =
1382 static_cast<int>(34 + 15 * 8 + 8 + strlen(poPro->proName) + 1 + 32 +
1383 8 + strlen(poPro->proSpheroid.sphereName) + 1);
1384
1385 if( poPro->proExeName != nullptr )
1386 nSize += static_cast<int>(strlen(poPro->proExeName) + 1);
1387
1388 GByte *pabyData = poMIEntry->MakeData(nSize);
1389 if( !pabyData )
1390 return CE_Failure;
1391
1392 poMIEntry->SetPosition();
1393
1394 // Initialize the whole thing to zeros for a clean start.
1395 memset(poMIEntry->GetData(), 0, poMIEntry->GetDataSize());
1396
1397 // Write the various fields.
1398 poMIEntry->SetIntField("proType", poPro->proType);
1399
1400 poMIEntry->SetIntField("proNumber", poPro->proNumber);
1401
1402 poMIEntry->SetStringField("proExeName", poPro->proExeName);
1403 poMIEntry->SetStringField("proName", poPro->proName);
1404 poMIEntry->SetIntField("proZone", poPro->proZone);
1405 poMIEntry->SetDoubleField("proParams[0]", poPro->proParams[0]);
1406 poMIEntry->SetDoubleField("proParams[1]", poPro->proParams[1]);
1407 poMIEntry->SetDoubleField("proParams[2]", poPro->proParams[2]);
1408 poMIEntry->SetDoubleField("proParams[3]", poPro->proParams[3]);
1409 poMIEntry->SetDoubleField("proParams[4]", poPro->proParams[4]);
1410 poMIEntry->SetDoubleField("proParams[5]", poPro->proParams[5]);
1411 poMIEntry->SetDoubleField("proParams[6]", poPro->proParams[6]);
1412 poMIEntry->SetDoubleField("proParams[7]", poPro->proParams[7]);
1413 poMIEntry->SetDoubleField("proParams[8]", poPro->proParams[8]);
1414 poMIEntry->SetDoubleField("proParams[9]", poPro->proParams[9]);
1415 poMIEntry->SetDoubleField("proParams[10]", poPro->proParams[10]);
1416 poMIEntry->SetDoubleField("proParams[11]", poPro->proParams[11]);
1417 poMIEntry->SetDoubleField("proParams[12]", poPro->proParams[12]);
1418 poMIEntry->SetDoubleField("proParams[13]", poPro->proParams[13]);
1419 poMIEntry->SetDoubleField("proParams[14]", poPro->proParams[14]);
1420 poMIEntry->SetStringField("proSpheroid.sphereName",
1421 poPro->proSpheroid.sphereName);
1422 poMIEntry->SetDoubleField("proSpheroid.a", poPro->proSpheroid.a);
1423 poMIEntry->SetDoubleField("proSpheroid.b", poPro->proSpheroid.b);
1424 poMIEntry->SetDoubleField("proSpheroid.eSquared",
1425 poPro->proSpheroid.eSquared);
1426 poMIEntry->SetDoubleField("proSpheroid.radius",
1427 poPro->proSpheroid.radius);
1428 }
1429
1430 return CE_None;
1431 }
1432
1433 /************************************************************************/
1434 /* HFAGetDatum() */
1435 /************************************************************************/
1436
HFAGetDatum(HFAHandle hHFA)1437 const Eprj_Datum *HFAGetDatum( HFAHandle hHFA )
1438
1439 {
1440 if( hHFA->nBands < 1 )
1441 return nullptr;
1442
1443 // Do we already have it?
1444 if( hHFA->pDatum != nullptr )
1445 return (Eprj_Datum *)hHFA->pDatum;
1446
1447 // Get the HFA node.
1448 HFAEntry *poMIEntry =
1449 hHFA->papoBand[0]->poNode->GetNamedChild("Projection.Datum");
1450 if( poMIEntry == nullptr )
1451 return nullptr;
1452
1453 // Allocate the structure.
1454 Eprj_Datum *psDatum =
1455 static_cast<Eprj_Datum *>(CPLCalloc(sizeof(Eprj_Datum), 1));
1456
1457 // Fetch the fields.
1458 psDatum->datumname = CPLStrdup(poMIEntry->GetStringField("datumname"));
1459 const int nDatumType = poMIEntry->GetIntField("type");
1460 if( nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE )
1461 {
1462 CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
1463 psDatum->type = EPRJ_DATUM_NONE;
1464 }
1465 else
1466 psDatum->type = static_cast<Eprj_DatumType>(nDatumType);
1467
1468 for( int i = 0; i < 7; i++ )
1469 {
1470 char szFieldName[30] = {};
1471 snprintf(szFieldName, sizeof(szFieldName), "params[%d]", i);
1472 psDatum->params[i] = poMIEntry->GetDoubleField(szFieldName);
1473 }
1474
1475 psDatum->gridname = CPLStrdup(poMIEntry->GetStringField("gridname"));
1476
1477 hHFA->pDatum = (void *)psDatum;
1478
1479 return psDatum;
1480 }
1481
1482 /************************************************************************/
1483 /* HFASetDatum() */
1484 /************************************************************************/
1485
HFASetDatum(HFAHandle hHFA,const Eprj_Datum * poDatum)1486 CPLErr HFASetDatum( HFAHandle hHFA, const Eprj_Datum *poDatum )
1487
1488 {
1489 // Loop over bands, setting information on each one.
1490 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1491 {
1492 // Create a new Projection if there isn't one present already.
1493 HFAEntry *poProParams =
1494 hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1495 if( poProParams == nullptr )
1496 {
1497 CPLError(CE_Failure, CPLE_AppDefined,
1498 "Can't add Eprj_Datum with no Eprj_ProjParameters.");
1499 return CE_Failure;
1500 }
1501
1502 HFAEntry *poDatumEntry = poProParams->GetNamedChild("Datum");
1503 if( poDatumEntry == nullptr )
1504 {
1505 poDatumEntry =
1506 HFAEntry::New(hHFA, "Datum", "Eprj_Datum", poProParams);
1507 }
1508
1509 poDatumEntry->MarkDirty();
1510
1511 // Ensure we have enough space for all the data.
1512 // TODO(schwehr): Explain constants.
1513 int nSize =
1514 static_cast<int>(26 + strlen(poDatum->datumname) + 1 + 7 * 8);
1515
1516 if( poDatum->gridname != nullptr )
1517 nSize += static_cast<int>(strlen(poDatum->gridname) + 1);
1518
1519 GByte *pabyData = poDatumEntry->MakeData(nSize);
1520 if( !pabyData )
1521 return CE_Failure;
1522
1523 poDatumEntry->SetPosition();
1524
1525 // Initialize the whole thing to zeros for a clean start.
1526 memset(poDatumEntry->GetData(), 0, poDatumEntry->GetDataSize());
1527
1528 // Write the various fields.
1529 poDatumEntry->SetStringField("datumname", poDatum->datumname);
1530 poDatumEntry->SetIntField("type", poDatum->type);
1531
1532 poDatumEntry->SetDoubleField("params[0]", poDatum->params[0]);
1533 poDatumEntry->SetDoubleField("params[1]", poDatum->params[1]);
1534 poDatumEntry->SetDoubleField("params[2]", poDatum->params[2]);
1535 poDatumEntry->SetDoubleField("params[3]", poDatum->params[3]);
1536 poDatumEntry->SetDoubleField("params[4]", poDatum->params[4]);
1537 poDatumEntry->SetDoubleField("params[5]", poDatum->params[5]);
1538 poDatumEntry->SetDoubleField("params[6]", poDatum->params[6]);
1539
1540 poDatumEntry->SetStringField("gridname", poDatum->gridname);
1541 }
1542
1543 return CE_None;
1544 }
1545
1546 /************************************************************************/
1547 /* HFAGetPCT() */
1548 /* */
1549 /* Read the PCT from a band, if it has one. */
1550 /************************************************************************/
1551
HFAGetPCT(HFAHandle hHFA,int nBand,int * pnColors,double ** ppadfRed,double ** ppadfGreen,double ** ppadfBlue,double ** ppadfAlpha,double ** ppadfBins)1552 CPLErr HFAGetPCT( HFAHandle hHFA, int nBand, int *pnColors,
1553 double **ppadfRed, double **ppadfGreen,
1554 double **ppadfBlue, double **ppadfAlpha,
1555 double **ppadfBins )
1556
1557 {
1558 if( nBand < 1 || nBand > hHFA->nBands )
1559 return CE_Failure;
1560
1561 return hHFA->papoBand[nBand - 1]->GetPCT(pnColors, ppadfRed, ppadfGreen,
1562 ppadfBlue, ppadfAlpha, ppadfBins);
1563 }
1564
1565 /************************************************************************/
1566 /* HFASetPCT() */
1567 /* */
1568 /* Set the PCT on a band. */
1569 /************************************************************************/
1570
HFASetPCT(HFAHandle hHFA,int nBand,int nColors,double * padfRed,double * padfGreen,double * padfBlue,double * padfAlpha)1571 CPLErr HFASetPCT( HFAHandle hHFA, int nBand, int nColors,
1572 double *padfRed, double *padfGreen, double *padfBlue,
1573 double *padfAlpha )
1574
1575 {
1576 if( nBand < 1 || nBand > hHFA->nBands )
1577 return CE_Failure;
1578
1579 return hHFA->papoBand[nBand - 1]->SetPCT(nColors, padfRed, padfGreen,
1580 padfBlue, padfAlpha);
1581 }
1582
1583 /************************************************************************/
1584 /* HFAGetDataRange() */
1585 /************************************************************************/
1586
HFAGetDataRange(HFAHandle hHFA,int nBand,double * pdfMin,double * pdfMax)1587 CPLErr HFAGetDataRange( HFAHandle hHFA, int nBand,
1588 double *pdfMin, double *pdfMax )
1589
1590 {
1591 if( nBand < 1 || nBand > hHFA->nBands )
1592 return CE_Failure;
1593
1594 HFAEntry *poBinInfo =
1595 hHFA->papoBand[nBand-1]->poNode->GetNamedChild("Statistics" );
1596
1597 if( poBinInfo == nullptr )
1598 return CE_Failure;
1599
1600 *pdfMin = poBinInfo->GetDoubleField("minimum");
1601 *pdfMax = poBinInfo->GetDoubleField("maximum");
1602
1603 if( *pdfMax > *pdfMin )
1604 return CE_None;
1605 else
1606 return CE_Failure;
1607 }
1608
1609 /************************************************************************/
1610 /* HFADumpNode() */
1611 /************************************************************************/
1612
HFADumpNode(HFAEntry * poEntry,int nIndent,bool bVerbose,FILE * fp)1613 static void HFADumpNode( HFAEntry *poEntry, int nIndent, bool bVerbose,
1614 FILE *fp )
1615
1616 {
1617 std::string osSpaces(nIndent * 2, ' ');
1618
1619 fprintf(fp, "%s%s(%s) @ %u + %u @ %u\n", osSpaces.c_str(),
1620 poEntry->GetName(), poEntry->GetType(),
1621 poEntry->GetFilePos(),
1622 poEntry->GetDataSize(), poEntry->GetDataPos());
1623
1624 if( bVerbose )
1625 {
1626 osSpaces += "+ ";
1627 poEntry->DumpFieldValues(fp, osSpaces.c_str());
1628 fprintf(fp, "\n");
1629 }
1630
1631 if( poEntry->GetChild() != nullptr )
1632 HFADumpNode(poEntry->GetChild(), nIndent+1, bVerbose, fp);
1633
1634 if( poEntry->GetNext() != nullptr )
1635 HFADumpNode(poEntry->GetNext(), nIndent, bVerbose, fp);
1636 }
1637
1638 /************************************************************************/
1639 /* HFADumpTree() */
1640 /* */
1641 /* Dump the tree of information in a HFA file. */
1642 /************************************************************************/
1643
HFADumpTree(HFAHandle hHFA,FILE * fpOut)1644 void HFADumpTree( HFAHandle hHFA, FILE *fpOut )
1645
1646 {
1647 HFADumpNode(hHFA->poRoot, 0, true, fpOut);
1648 }
1649
1650 /************************************************************************/
1651 /* HFADumpDictionary() */
1652 /* */
1653 /* Dump the dictionary (in raw, and parsed form) to the named */
1654 /* device. */
1655 /************************************************************************/
1656
HFADumpDictionary(HFAHandle hHFA,FILE * fpOut)1657 void HFADumpDictionary( HFAHandle hHFA, FILE *fpOut )
1658
1659 {
1660 fprintf(fpOut, "%s\n", hHFA->pszDictionary);
1661
1662 hHFA->poDictionary->Dump(fpOut);
1663 }
1664
1665 /************************************************************************/
1666 /* HFAStandard() */
1667 /* */
1668 /* Swap byte order on MSB systems. */
1669 /************************************************************************/
1670
1671 #ifdef CPL_MSB
HFAStandard(int nBytes,void * pData)1672 void HFAStandard( int nBytes, void *pData )
1673
1674 {
1675 GByte *pabyData = static_cast<GByte *>(pData);
1676
1677 for( int i = nBytes / 2 - 1; i >= 0; i-- )
1678 {
1679 GByte byTemp = pabyData[i];
1680 pabyData[i] = pabyData[nBytes - i - 1];
1681 pabyData[nBytes - i - 1] = byTemp;
1682 }
1683 }
1684 #endif
1685
1686 /* ==================================================================== */
1687 /* Default data dictionary. Emitted verbatim into the imagine */
1688 /* file. */
1689 /* ==================================================================== */
1690
1691 static const char * const aszDefaultDD[] = {
1692 "{1:lversion,1:LfreeList,1:LrootEntryPtr,1:sentryHeaderLength,1:LdictionaryPtr,}Ehfa_File,{1:Lnext,1:Lprev,1:Lparent,1:Lchild,1:Ldata,1:ldataSize,64:cname,32:ctype,1:tmodTime,}Ehfa_Entry,{16:clabel,1:LheaderPtr,}Ehfa_HeaderTag,{1:LfreeList,1:lfreeSize,}Ehfa_FreeListNode,{1:lsize,1:Lptr,}Ehfa_Data,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,",
1693 "1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer_SubSample,{1:e2:raster,vector,type,1:LdictionaryPtr,}Ehfa_Layer,{1:LspaceUsedForRasterData,}ImgFormatInfo831,{1:sfileCode,1:Loffset,1:lsize,1:e2:false,true,logvalid,",
1694 "1:e2:no compression,ESRI GRID compression,compressionType,}Edms_VirtualBlockInfo,{1:lmin,1:lmax,}Edms_FreeIDList,{1:lnumvirtualblocks,1:lnumobjectsperblock,1:lnextobjectnum,1:e2:no compression,RLC compression,compressionType,0:poEdms_VirtualBlockInfo,blockinfo,0:poEdms_FreeIDList,freelist,1:tmodTime,}Edms_State,{0:pcstring,}Emif_String,{1:oEmif_String,fileName,2:LlayerStackValidFlagsOffset,2:LlayerStackDataOffset,1:LlayerStackCount,1:LlayerStackIndex,}ImgExternalRaster,{1:oEmif_String,algorithm,0:poEmif_String,nameList,}Eimg_RRDNamesList,{1:oEmif_String,projection,1:oEmif_String,units,}Eimg_MapInformation,",
1695 "{1:oEmif_String,dependent,}Eimg_DependentFile,{1:oEmif_String,ImageLayerName,}Eimg_DependentLayerName,{1:lnumrows,1:lnumcolumns,1:e13:EGDA_TYPE_U1,EGDA_TYPE_U2,EGDA_TYPE_U4,EGDA_TYPE_U8,EGDA_TYPE_S8,EGDA_TYPE_U16,EGDA_TYPE_S16,EGDA_TYPE_U32,EGDA_TYPE_S32,EGDA_TYPE_F32,EGDA_TYPE_F64,EGDA_TYPE_C64,EGDA_TYPE_C128,datatype,1:e4:EGDA_SCALAR_OBJECT,EGDA_TABLE_OBJECT,EGDA_MATRIX_OBJECT,EGDA_RASTER_OBJECT,objecttype,}Egda_BaseData,{1:*bvalueBD,}Eimg_NonInitializedValue,{1:dx,1:dy,}Eprj_Coordinate,{1:dwidth,1:dheight,}Eprj_Size,{0:pcproName,1:*oEprj_Coordinate,upperLeftCenter,",
1696 "1:*oEprj_Coordinate,lowerRightCenter,1:*oEprj_Size,pixelSize,0:pcunits,}Eprj_MapInfo,{0:pcdatumname,1:e3:EPRJ_DATUM_PARAMETRIC,EPRJ_DATUM_GRID,EPRJ_DATUM_REGRESSION,type,0:pdparams,0:pcgridname,}Eprj_Datum,{0:pcsphereName,1:da,1:db,1:deSquared,1:dradius,}Eprj_Spheroid,{1:e2:EPRJ_INTERNAL,EPRJ_EXTERNAL,proType,1:lproNumber,0:pcproExeName,0:pcproName,1:lproZone,0:pdproParams,1:*oEprj_Spheroid,proSpheroid,}Eprj_ProParameters,{1:dminimum,1:dmaximum,1:dmean,1:dmedian,1:dmode,1:dstddev,}Esta_Statistics,{1:lnumBins,1:e4:direct,linear,logarithmic,explicit,binFunctionType,1:dminLimit,1:dmaxLimit,1:*bbinLimits,}Edsc_BinFunction,{0:poEmif_String,LayerNames,1:*bExcludedValues,1:oEmif_String,AOIname,",
1697 "1:lSkipFactorX,1:lSkipFactorY,1:*oEdsc_BinFunction,BinFunction,}Eimg_StatisticsParameters830,{1:lnumrows,}Edsc_Table,{1:lnumRows,1:LcolumnDataPtr,1:e4:integer,real,complex,string,dataType,1:lmaxNumChars,}Edsc_Column,{1:lposition,0:pcname,1:e2:EMSC_FALSE,EMSC_TRUE,editable,1:e3:LEFT,CENTER,RIGHT,alignment,0:pcformat,1:e3:DEFAULT,APPLY,AUTO-APPLY,formulamode,0:pcformula,1:dcolumnwidth,0:pcunits,1:e5:NO_COLOR,RED,GREEN,BLUE,COLOR,colorflag,0:pcgreenname,0:pcbluename,}Eded_ColumnAttributes_1,{1:lversion,1:lnumobjects,1:e2:EAOI_UNION,EAOI_INTERSECTION,operation,}Eaoi_AreaOfInterest,",
1698 "{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,",
1699 "{1:x{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,projection,1:x{0:pcstring,}Emif_String,title,}Eprj_MapProjection842,",
1700 "{0:poEmif_String,titleList,}Exfr_GenericXFormHeader,{1:lorder,1:lnumdimtransform,1:lnumdimpolynomial,1:ltermcount,0:plexponentlist,1:*bpolycoefmtx,1:*bpolycoefvector,}Efga_Polynomial,",
1701 ".",
1702 nullptr
1703 };
1704
1705 /************************************************************************/
1706 /* HFACreateLL() */
1707 /* */
1708 /* Low level creation of an Imagine file. Writes out the */
1709 /* Ehfa_HeaderTag, dictionary and Ehfa_File. */
1710 /************************************************************************/
1711
HFACreateLL(const char * pszFilename)1712 HFAHandle HFACreateLL( const char *pszFilename )
1713
1714 {
1715 // Create the file in the file system.
1716 VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
1717 if( fp == nullptr )
1718 {
1719 CPLError(CE_Failure, CPLE_OpenFailed, "Creation of file %s failed.",
1720 pszFilename);
1721 return nullptr;
1722 }
1723
1724 // Create the HFAInfo_t.
1725 HFAInfo_t *psInfo =
1726 static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
1727
1728 psInfo->fp = fp;
1729 psInfo->eAccess = HFA_Update;
1730 psInfo->nXSize = 0;
1731 psInfo->nYSize = 0;
1732 psInfo->nBands = 0;
1733 psInfo->papoBand = nullptr;
1734 psInfo->pMapInfo = nullptr;
1735 psInfo->pDatum = nullptr;
1736 psInfo->pProParameters = nullptr;
1737 psInfo->bTreeDirty = false;
1738 psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
1739 psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
1740
1741 // Write out the Ehfa_HeaderTag.
1742 bool bRet = VSIFWriteL((void *)"EHFA_HEADER_TAG", 1, 16, fp) > 0;
1743
1744 GInt32 nHeaderPos = 20;
1745 HFAStandard(4, &nHeaderPos);
1746 bRet &= VSIFWriteL(&nHeaderPos, 4, 1, fp) > 0;
1747
1748 // Write the Ehfa_File node, locked in at offset 20.
1749 GInt32 nVersion = 1;
1750 GInt32 nFreeList = 0;
1751 GInt32 nRootEntry = 0;
1752 GInt16 nEntryHeaderLength = 128;
1753 GInt32 nDictionaryPtr = 38;
1754
1755 psInfo->nEntryHeaderLength = nEntryHeaderLength;
1756 psInfo->nRootPos = 0;
1757 psInfo->nDictionaryPos = nDictionaryPtr;
1758 psInfo->nVersion = nVersion;
1759
1760 HFAStandard(4, &nVersion);
1761 HFAStandard(4, &nFreeList);
1762 HFAStandard(4, &nRootEntry);
1763 HFAStandard(2, &nEntryHeaderLength);
1764 HFAStandard(4, &nDictionaryPtr);
1765
1766 bRet &= VSIFWriteL(&nVersion, 4, 1, fp) > 0;
1767 bRet &= VSIFWriteL(&nFreeList, 4, 1, fp) > 0;
1768 bRet &= VSIFWriteL(&nRootEntry, 4, 1, fp) > 0;
1769 bRet &= VSIFWriteL(&nEntryHeaderLength, 2, 1, fp) > 0;
1770 bRet &= VSIFWriteL(&nDictionaryPtr, 4, 1, fp) > 0;
1771
1772 // Write the dictionary, locked in at location 38. Note that
1773 // we jump through a bunch of hoops to operate on the
1774 // dictionary in chunks because some compiles (such as VC++)
1775 // don't allow particularly large static strings.
1776 int nDictLen = 0;
1777
1778 for( int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++ )
1779 nDictLen += static_cast<int>(strlen(aszDefaultDD[iChunk]));
1780
1781 psInfo->pszDictionary = static_cast<char *>(CPLMalloc(nDictLen + 1));
1782 psInfo->pszDictionary[0] = '\0';
1783
1784 for( int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++ )
1785 strcat(psInfo->pszDictionary, aszDefaultDD[iChunk]);
1786
1787 bRet &= VSIFWriteL((void *)psInfo->pszDictionary,
1788 strlen(psInfo->pszDictionary) + 1, 1, fp) > 0;
1789 if( !bRet )
1790 {
1791 CPL_IGNORE_RET_VAL(HFAClose(psInfo));
1792 return nullptr;
1793 }
1794
1795 psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
1796
1797 psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
1798
1799 // Create a root entry.
1800 psInfo->poRoot = new HFAEntry(psInfo, "root", "root", nullptr);
1801
1802 // If an .ige or .rrd file exists with the same base name,
1803 // delete them. (#1784)
1804 CPLString osExtension = CPLGetExtension(pszFilename);
1805 if( !EQUAL(osExtension, "rrd") && !EQUAL(osExtension, "aux") )
1806 {
1807 CPLString osPath = CPLGetPath(pszFilename);
1808 CPLString osBasename = CPLGetBasename(pszFilename);
1809 VSIStatBufL sStatBuf;
1810 CPLString osSupFile = CPLFormCIFilename(osPath, osBasename, "rrd");
1811
1812 if( VSIStatL(osSupFile, &sStatBuf) == 0 )
1813 VSIUnlink(osSupFile);
1814
1815 osSupFile = CPLFormCIFilename(osPath, osBasename, "ige");
1816
1817 if( VSIStatL(osSupFile, &sStatBuf) == 0 )
1818 VSIUnlink(osSupFile);
1819 }
1820
1821 return psInfo;
1822 }
1823
1824 /************************************************************************/
1825 /* HFAAllocateSpace() */
1826 /* */
1827 /* Return an area in the file to the caller to write the */
1828 /* requested number of bytes. Currently this is always at the */
1829 /* end of the file, but eventually we might actually keep track */
1830 /* of free space. The HFAInfo_t's concept of file size is */
1831 /* updated, even if nothing ever gets written to this region. */
1832 /* */
1833 /* Returns the offset to the requested space, or zero one */
1834 /* failure. */
1835 /************************************************************************/
1836
HFAAllocateSpace(HFAInfo_t * psInfo,GUInt32 nBytes)1837 GUInt32 HFAAllocateSpace( HFAInfo_t *psInfo, GUInt32 nBytes )
1838
1839 {
1840 // TODO(schwehr): Check if this will wrap over 2GB limit.
1841
1842 psInfo->nEndOfFile += nBytes;
1843 return psInfo->nEndOfFile - nBytes;
1844 }
1845
1846 /************************************************************************/
1847 /* HFAFlush() */
1848 /* */
1849 /* Write out any dirty tree information to disk, putting the */
1850 /* disk file in a consistent state. */
1851 /************************************************************************/
1852
HFAFlush(HFAHandle hHFA)1853 CPLErr HFAFlush( HFAHandle hHFA )
1854
1855 {
1856 if( !hHFA->bTreeDirty && !hHFA->poDictionary->bDictionaryTextDirty )
1857 return CE_None;
1858
1859 CPLAssert(hHFA->poRoot != nullptr);
1860
1861 // Flush HFAEntry tree to disk.
1862 if( hHFA->bTreeDirty )
1863 {
1864 const CPLErr eErr = hHFA->poRoot->FlushToDisk();
1865 if( eErr != CE_None )
1866 return eErr;
1867
1868 hHFA->bTreeDirty = false;
1869 }
1870
1871 // Flush Dictionary to disk.
1872 GUInt32 nNewDictionaryPos = hHFA->nDictionaryPos;
1873 bool bRet = true;
1874 if( hHFA->poDictionary->bDictionaryTextDirty )
1875 {
1876 bRet &= VSIFSeekL(hHFA->fp, 0, SEEK_END) >= 0;
1877 nNewDictionaryPos = static_cast<GUInt32>(VSIFTellL(hHFA->fp));
1878 bRet &=
1879 VSIFWriteL(hHFA->poDictionary->osDictionaryText.c_str(),
1880 strlen(hHFA->poDictionary->osDictionaryText.c_str()) + 1,
1881 1, hHFA->fp) > 0;
1882 hHFA->poDictionary->bDictionaryTextDirty = false;
1883 }
1884
1885 // Do we need to update the Ehfa_File pointer to the root node?
1886 if( hHFA->nRootPos != hHFA->poRoot->GetFilePos() ||
1887 nNewDictionaryPos != hHFA->nDictionaryPos )
1888 {
1889 GUInt32 nHeaderPos = 0;
1890
1891 bRet &= VSIFSeekL(hHFA->fp, 16, SEEK_SET) >= 0;
1892 bRet &= VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, hHFA->fp) > 0;
1893 HFAStandard(4, &nHeaderPos);
1894
1895 GUInt32 nOffset = hHFA->poRoot->GetFilePos();
1896 hHFA->nRootPos = nOffset;
1897 HFAStandard(4, &nOffset);
1898 bRet &= VSIFSeekL(hHFA->fp, nHeaderPos + 8, SEEK_SET) >= 0;
1899 bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
1900
1901 nOffset = nNewDictionaryPos;
1902 hHFA->nDictionaryPos = nNewDictionaryPos;
1903 HFAStandard(4, &nOffset);
1904 bRet &= VSIFSeekL(hHFA->fp, nHeaderPos + 14, SEEK_SET) >= 0;
1905 bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
1906 }
1907
1908 return bRet ? CE_None : CE_Failure;
1909 }
1910
1911 /************************************************************************/
1912 /* HFACreateLayer() */
1913 /* */
1914 /* Create a layer object, and corresponding RasterDMS. */
1915 /* Suitable for use with primary layers, and overviews. */
1916 /************************************************************************/
1917
1918 int
HFACreateLayer(HFAHandle psInfo,HFAEntry * poParent,const char * pszLayerName,int bOverview,int nBlockSize,int bCreateCompressed,int bCreateLargeRaster,int bDependentLayer,int nXSize,int nYSize,EPTType eDataType,char **,GIntBig nStackValidFlagsOffset,GIntBig nStackDataOffset,int nStackCount,int nStackIndex)1919 HFACreateLayer( HFAHandle psInfo, HFAEntry *poParent,
1920 const char *pszLayerName,
1921 int bOverview, int nBlockSize,
1922 int bCreateCompressed, int bCreateLargeRaster,
1923 int bDependentLayer,
1924 int nXSize, int nYSize, EPTType eDataType,
1925 char ** /* papszOptions */,
1926 // These are only related to external (large) files.
1927 GIntBig nStackValidFlagsOffset,
1928 GIntBig nStackDataOffset,
1929 int nStackCount, int nStackIndex )
1930
1931 {
1932 const char *pszLayerType =
1933 bOverview ? "Eimg_Layer_SubSample" : "Eimg_Layer";
1934
1935 if( nBlockSize <= 0 )
1936 {
1937 CPLError(CE_Failure, CPLE_IllegalArg,
1938 "HFACreateLayer: nBlockXSize < 0");
1939 return FALSE;
1940 }
1941
1942 // Work out some details about the tiling scheme.
1943 const int nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
1944 const int nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
1945 const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
1946 const int nBytesPerBlock =
1947 (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
1948
1949 // Create the Eimg_Layer for the band.
1950 HFAEntry *poEimg_Layer =
1951 HFAEntry::New(psInfo, pszLayerName, pszLayerType, poParent);
1952
1953 poEimg_Layer->SetIntField("width", nXSize);
1954 poEimg_Layer->SetIntField("height", nYSize);
1955 poEimg_Layer->SetStringField("layerType", "athematic");
1956 poEimg_Layer->SetIntField("pixelType", eDataType);
1957 poEimg_Layer->SetIntField("blockWidth", nBlockSize);
1958 poEimg_Layer->SetIntField("blockHeight", nBlockSize);
1959
1960 // Create the RasterDMS (block list). This is a complex type
1961 // with pointers, and variable size. We set the superstructure
1962 // ourselves rather than trying to have the HFA type management
1963 // system do it for us (since this would be hard to implement).
1964 if( !bCreateLargeRaster && !bDependentLayer )
1965 {
1966 HFAEntry *poEdms_State =
1967 HFAEntry::New(psInfo, "RasterDMS", "Edms_State", poEimg_Layer);
1968
1969 // TODO(schwehr): Explain constants.
1970 const int nDmsSize = 14 * nBlocks + 38;
1971 GByte *pabyData = poEdms_State->MakeData(nDmsSize);
1972
1973 // Set some simple values.
1974 poEdms_State->SetIntField("numvirtualblocks", nBlocks);
1975 poEdms_State->SetIntField("numobjectsperblock",
1976 nBlockSize * nBlockSize);
1977 poEdms_State->SetIntField("nextobjectnum",
1978 nBlockSize * nBlockSize * nBlocks);
1979
1980 // Is file compressed or not?
1981 if( bCreateCompressed )
1982 {
1983 poEdms_State->SetStringField("compressionType", "RLC compression");
1984 }
1985 else
1986 {
1987 poEdms_State->SetStringField("compressionType", "no compression");
1988 }
1989
1990 // We need to hardcode file offset into the data, so locate it now.
1991 poEdms_State->SetPosition();
1992
1993 // Set block info headers.
1994
1995 // Blockinfo count.
1996 GUInt32 nValue = nBlocks;
1997 HFAStandard(4, &nValue);
1998 memcpy(pabyData + 14, &nValue, 4);
1999
2000 // Blockinfo position.
2001 nValue = poEdms_State->GetDataPos() + 22;
2002 HFAStandard(4, &nValue);
2003 memcpy(pabyData + 18, &nValue, 4);
2004
2005 // Set each blockinfo.
2006 for( int iBlock = 0; iBlock < nBlocks; iBlock++ )
2007 {
2008 int nOffset = 22 + 14 * iBlock;
2009
2010 // fileCode.
2011 GInt16 nValue16 = 0;
2012 HFAStandard(2, &nValue16);
2013 memcpy(pabyData + nOffset, &nValue16, 2);
2014
2015 // Offset.
2016 if( bCreateCompressed )
2017 {
2018 // Flag it with zero offset. Allocate space when we compress it.
2019 nValue = 0;
2020 }
2021 else
2022 {
2023 nValue = HFAAllocateSpace(psInfo, nBytesPerBlock);
2024 }
2025 HFAStandard(4, &nValue);
2026 memcpy(pabyData + nOffset + 2, &nValue, 4);
2027
2028 // Size.
2029 if( bCreateCompressed )
2030 {
2031 // Flag with zero size. Don't know until we compress it.
2032 nValue = 0;
2033 }
2034 else
2035 {
2036 nValue = nBytesPerBlock;
2037 }
2038 HFAStandard(4, &nValue);
2039 memcpy(pabyData + nOffset + 6, &nValue, 4);
2040
2041 // logValid (false).
2042 nValue16 = 0;
2043 HFAStandard(2, &nValue16);
2044 memcpy(pabyData + nOffset + 10, &nValue16, 2);
2045
2046 // compressionType.
2047 if( bCreateCompressed )
2048 nValue16 = 1;
2049 else
2050 nValue16 = 0;
2051
2052 HFAStandard(2, &nValue16);
2053 memcpy(pabyData + nOffset + 12, &nValue16, 2);
2054 }
2055 }
2056
2057 // Create ExternalRasterDMS object.
2058 else if( bCreateLargeRaster )
2059 {
2060 HFAEntry *poEdms_State = HFAEntry::New(
2061 psInfo, "ExternalRasterDMS", "ImgExternalRaster", poEimg_Layer);
2062 poEdms_State->MakeData(
2063 static_cast<int>(8 + strlen(psInfo->pszIGEFilename) + 1 + 6 * 4));
2064
2065 poEdms_State->SetStringField("fileName.string", psInfo->pszIGEFilename);
2066
2067 poEdms_State->SetIntField(
2068 "layerStackValidFlagsOffset[0]",
2069 static_cast<int>(nStackValidFlagsOffset & 0xFFFFFFFF));
2070 poEdms_State->SetIntField(
2071 "layerStackValidFlagsOffset[1]",
2072 static_cast<int>(nStackValidFlagsOffset >> 32));
2073
2074 poEdms_State->SetIntField(
2075 "layerStackDataOffset[0]",
2076 static_cast<int>(nStackDataOffset & 0xFFFFFFFF));
2077 poEdms_State->SetIntField("layerStackDataOffset[1]",
2078 static_cast<int>(nStackDataOffset >> 32));
2079 poEdms_State->SetIntField("layerStackCount", nStackCount);
2080 poEdms_State->SetIntField("layerStackIndex", nStackIndex);
2081 }
2082 // Dependent...
2083 else if( bDependentLayer )
2084 {
2085 HFAEntry *poDepLayerName =
2086 HFAEntry::New(psInfo, "DependentLayerName",
2087 "Eimg_DependentLayerName", poEimg_Layer);
2088 poDepLayerName->MakeData(
2089 static_cast<int>(8 + strlen(pszLayerName) + 2));
2090
2091 poDepLayerName->SetStringField("ImageLayerName.string", pszLayerName);
2092 }
2093
2094 // Create the Ehfa_Layer.
2095 char chBandType = '\0';
2096
2097 if(eDataType == EPT_u1)
2098 chBandType = '1';
2099 else if(eDataType == EPT_u2)
2100 chBandType = '2';
2101 else if(eDataType == EPT_u4)
2102 chBandType = '4';
2103 else if(eDataType == EPT_u8)
2104 chBandType = 'c';
2105 else if(eDataType == EPT_s8)
2106 chBandType = 'C';
2107 else if(eDataType == EPT_u16)
2108 chBandType = 's';
2109 else if(eDataType == EPT_s16)
2110 chBandType = 'S';
2111 else if(eDataType == EPT_u32)
2112 // For some reason erdas imagine expects an L for unsigned 32 bit ints
2113 // otherwise it gives strange "out of memory errors".
2114 chBandType = 'L';
2115 else if(eDataType == EPT_s32)
2116 chBandType = 'L';
2117 else if(eDataType == EPT_f32)
2118 chBandType = 'f';
2119 else if(eDataType == EPT_f64)
2120 chBandType = 'd';
2121 else if(eDataType == EPT_c64)
2122 chBandType = 'm';
2123 else if(eDataType == EPT_c128)
2124 chBandType = 'M';
2125 else
2126 {
2127 CPLAssert(false);
2128 chBandType = 'c';
2129 }
2130
2131 // The first value in the entry below gives the number of pixels
2132 // within a block.
2133 char szLDict[128] = {};
2134 snprintf(szLDict, sizeof(szLDict), "{%d:%cdata,}RasterDMS,.",
2135 nBlockSize * nBlockSize, chBandType);
2136
2137 HFAEntry *poEhfa_Layer =
2138 HFAEntry::New(psInfo, "Ehfa_Layer", "Ehfa_Layer", poEimg_Layer);
2139 poEhfa_Layer->MakeData();
2140 poEhfa_Layer->SetPosition();
2141 const GUInt32 nLDict =
2142 HFAAllocateSpace(psInfo, static_cast<GUInt32>(strlen(szLDict) + 1));
2143
2144 poEhfa_Layer->SetStringField("type", "raster");
2145 poEhfa_Layer->SetIntField("dictionaryPtr", nLDict);
2146
2147 bool bRet = VSIFSeekL(psInfo->fp, nLDict, SEEK_SET) >= 0;
2148 bRet &= VSIFWriteL((void *)szLDict, strlen(szLDict) + 1, 1, psInfo->fp) > 0;
2149
2150 return bRet;
2151 }
2152
2153 /************************************************************************/
2154 /* HFACreate() */
2155 /************************************************************************/
2156
HFACreate(const char * pszFilename,int nXSize,int nYSize,int nBands,EPTType eDataType,char ** papszOptions)2157 HFAHandle HFACreate( const char *pszFilename,
2158 int nXSize, int nYSize, int nBands,
2159 EPTType eDataType, char **papszOptions )
2160
2161 {
2162 int nBlockSize = 64;
2163 const char *pszValue = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
2164
2165 if( pszValue != nullptr )
2166 {
2167 nBlockSize = atoi(pszValue);
2168 // Check for sane values.
2169 if( nBlockSize == 0 ||
2170 ((( nBlockSize < 32 ) || (nBlockSize > 2048)) &&
2171 !CPLTestBool(CPLGetConfigOption("FORCE_BLOCKSIZE", "NO"))) )
2172 {
2173 if( nBlockSize != 0 )
2174 CPLError(CE_Warning, CPLE_AppDefined,
2175 "Forcing BLOCKSIZE to %d", 64);
2176 nBlockSize = 64;
2177 }
2178 }
2179 bool bCreateLargeRaster = CPLFetchBool(papszOptions, "USE_SPILL", false);
2180 bool bCreateCompressed = CPLFetchBool(papszOptions, "COMPRESS", false) ||
2181 CPLFetchBool(papszOptions, "COMPRESSED", false);
2182 const bool bCreateAux = CPLFetchBool(papszOptions, "AUX", false);
2183
2184 char *pszFullFilename = nullptr;
2185 char *pszRawFilename = nullptr;
2186
2187 // Work out some details about the tiling scheme.
2188 const int nBlocksPerRow = DIV_ROUND_UP(nXSize, nBlockSize);
2189 const int nBlocksPerColumn = DIV_ROUND_UP(nYSize, nBlockSize);
2190 if( nBlocksPerRow > INT_MAX / nBlocksPerColumn )
2191 {
2192 CPLError(CE_Failure, CPLE_NotSupported, "Too many blocks");
2193 return nullptr;
2194 }
2195 const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
2196 const GInt64 nBytesPerBlock64 =
2197 (static_cast<GInt64>(nBlockSize) * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
2198 if( nBytesPerBlock64 > INT_MAX )
2199 {
2200 CPLError(CE_Failure, CPLE_NotSupported, "Too large block");
2201 return nullptr;
2202 }
2203 const int nBytesPerBlock = static_cast<int>(nBytesPerBlock64);
2204
2205 // Create the low level structure.
2206 HFAHandle psInfo = HFACreateLL(pszFilename);
2207 if( psInfo == nullptr )
2208 return nullptr;
2209
2210 // Create the DependentFile node if requested.
2211 const char *pszDependentFile =
2212 CSLFetchNameValue(papszOptions, "DEPENDENT_FILE");
2213
2214 if( pszDependentFile != nullptr )
2215 {
2216 HFAEntry *poDF = HFAEntry::New(psInfo, "DependentFile",
2217 "Eimg_DependentFile", psInfo->poRoot);
2218
2219 poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
2220 poDF->SetPosition();
2221 poDF->SetStringField("dependent.string", pszDependentFile);
2222 }
2223
2224 CPLDebug("HFACreate", "Blocks per row %d, blocks per column %d, "
2225 "total number of blocks %d, bytes per block %d.",
2226 nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock);
2227
2228 // Check whether we should create external large file with
2229 // image. We create a spill file if the amount of imagery is
2230 // close to 2GB. We don't check the amount of auxiliary
2231 // information, so in theory if there were an awful lot of
2232 // non-imagery data our approximate size could be smaller than
2233 // the file will actually we be. We leave room for 10MB of
2234 // auxiliary data.
2235 // We can also force spill file creation using option
2236 // SPILL_FILE=YES.
2237 const double dfApproxSize =
2238 static_cast<double>(nBytesPerBlock) * static_cast<double>(nBlocks) *
2239 static_cast<double>(nBands) + 10000000.0;
2240
2241 if( dfApproxSize > 2147483648.0 && !bCreateAux )
2242 bCreateLargeRaster = true;
2243
2244 // Erdas Imagine creates this entry even if an external spill file is used.
2245 if( !bCreateAux )
2246 {
2247 HFAEntry *poImgFormat = HFAEntry::New(
2248 psInfo, "IMGFormatInfo", "ImgFormatInfo831", psInfo->poRoot);
2249 poImgFormat->MakeData();
2250 if( bCreateLargeRaster )
2251 {
2252 poImgFormat->SetIntField("spaceUsedForRasterData", 0);
2253 // Can't be compressed if we are creating a spillfile.
2254 bCreateCompressed = false;
2255 }
2256 else
2257 {
2258 poImgFormat->SetIntField("spaceUsedForRasterData",
2259 nBytesPerBlock * nBlocks * nBands);
2260 }
2261 }
2262
2263 // Create external file and write its header.
2264 GIntBig nValidFlagsOffset = 0;
2265 GIntBig nDataOffset = 0;
2266
2267 if( bCreateLargeRaster )
2268 {
2269 if( !HFACreateSpillStack(psInfo, nXSize, nYSize, nBands,
2270 nBlockSize, eDataType,
2271 &nValidFlagsOffset, &nDataOffset) )
2272 {
2273 CPLFree(pszRawFilename);
2274 CPLFree(pszFullFilename);
2275 return nullptr;
2276 }
2277 }
2278
2279 // Create each band (layer).
2280 for( int iBand = 0; iBand < nBands; iBand++ )
2281 {
2282 char szName[128] = {};
2283
2284 snprintf(szName, sizeof(szName), "Layer_%d", iBand + 1);
2285
2286 if( !HFACreateLayer(psInfo, psInfo->poRoot, szName, FALSE, nBlockSize,
2287 bCreateCompressed, bCreateLargeRaster, bCreateAux,
2288 nXSize, nYSize, eDataType, papszOptions,
2289 nValidFlagsOffset, nDataOffset,
2290 nBands, iBand) )
2291 {
2292 CPL_IGNORE_RET_VAL(HFAClose(psInfo));
2293 return nullptr;
2294 }
2295 }
2296
2297 // Initialize the band information.
2298 HFAParseBandInfo(psInfo);
2299
2300 return psInfo;
2301 }
2302
2303 /************************************************************************/
2304 /* HFACreateOverview() */
2305 /* */
2306 /* Create an overview layer object for a band. */
2307 /************************************************************************/
2308
HFACreateOverview(HFAHandle hHFA,int nBand,int nOverviewLevel,const char * pszResampling)2309 int HFACreateOverview( HFAHandle hHFA, int nBand, int nOverviewLevel,
2310 const char *pszResampling )
2311
2312 {
2313 if( nBand < 1 || nBand > hHFA->nBands )
2314 return -1;
2315
2316 HFABand *poBand = hHFA->papoBand[nBand - 1];
2317 return poBand->CreateOverview(nOverviewLevel, pszResampling);
2318 }
2319
2320 /************************************************************************/
2321 /* HFAGetMetadata() */
2322 /* */
2323 /* Read metadata structured in a table called GDAL_MetaData. */
2324 /************************************************************************/
2325
HFAGetMetadata(HFAHandle hHFA,int nBand)2326 char **HFAGetMetadata( HFAHandle hHFA, int nBand )
2327
2328 {
2329 HFAEntry *poTable = nullptr;
2330
2331 if( nBand > 0 && nBand <= hHFA->nBands )
2332 poTable = hHFA->papoBand[nBand - 1]->poNode->GetChild();
2333 else if( nBand == 0 )
2334 poTable = hHFA->poRoot->GetChild();
2335 else
2336 return nullptr;
2337
2338 for( ;
2339 poTable != nullptr && !EQUAL(poTable->GetName(), "GDAL_MetaData");
2340 poTable = poTable->GetNext() ) {}
2341
2342 if( poTable == nullptr || !EQUAL(poTable->GetType(), "Edsc_Table") )
2343 return nullptr;
2344
2345 if( poTable->GetIntField("numRows") != 1 )
2346 {
2347 CPLDebug("HFADataset", "GDAL_MetaData.numRows = %d, expected 1!",
2348 poTable->GetIntField("numRows"));
2349 return nullptr;
2350 }
2351
2352 // Loop over each column. Each column will be one metadata
2353 // entry, with the title being the key, and the row value being
2354 // the value. There is only ever one row in GDAL_MetaData tables.
2355 char **papszMD = nullptr;
2356
2357 for( HFAEntry *poColumn = poTable->GetChild();
2358 poColumn != nullptr;
2359 poColumn = poColumn->GetNext() )
2360 {
2361 // Skip the #Bin_Function# entry.
2362 if( STARTS_WITH_CI(poColumn->GetName(), "#") )
2363 continue;
2364
2365 const char *pszValue = poColumn->GetStringField("dataType");
2366 if( pszValue == nullptr || !EQUAL(pszValue, "string") )
2367 continue;
2368
2369 const int columnDataPtr = poColumn->GetIntField("columnDataPtr");
2370 if( columnDataPtr <= 0 )
2371 continue;
2372
2373 // Read up to nMaxNumChars bytes from the indicated location.
2374 // allocate required space temporarily nMaxNumChars should have been
2375 // set by GDAL originally so we should trust it, but who knows.
2376 const int nMaxNumChars = poColumn->GetIntField("maxNumChars");
2377
2378 if( nMaxNumChars <= 0 )
2379 {
2380 papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), "");
2381 }
2382 else
2383 {
2384 char *pszMDValue =
2385 static_cast<char *>(VSI_MALLOC_VERBOSE(nMaxNumChars));
2386 if( pszMDValue == nullptr )
2387 {
2388 continue;
2389 }
2390
2391 if( VSIFSeekL(hHFA->fp, columnDataPtr, SEEK_SET) != 0 )
2392 {
2393 CPLFree(pszMDValue);
2394 continue;
2395 }
2396
2397 const int nMDBytes = static_cast<int>(
2398 VSIFReadL(pszMDValue, 1, nMaxNumChars, hHFA->fp));
2399 if( nMDBytes == 0 )
2400 {
2401 CPLFree(pszMDValue);
2402 continue;
2403 }
2404
2405 pszMDValue[nMaxNumChars - 1] = '\0';
2406
2407 papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), pszMDValue);
2408 CPLFree(pszMDValue);
2409 }
2410 }
2411
2412 return papszMD;
2413 }
2414
2415 /************************************************************************/
2416 /* HFASetGDALMetadata() */
2417 /* */
2418 /* This function is used to set metadata in a table called */
2419 /* GDAL_MetaData. It is called by HFASetMetadata() for all */
2420 /* metadata items that aren't some specific supported */
2421 /* information (like histogram or stats info). */
2422 /************************************************************************/
2423
2424 static CPLErr
HFASetGDALMetadata(HFAHandle hHFA,int nBand,char ** papszMD)2425 HFASetGDALMetadata( HFAHandle hHFA, int nBand, char **papszMD )
2426
2427 {
2428 if( papszMD == nullptr )
2429 return CE_None;
2430
2431 HFAEntry *poNode = nullptr;
2432
2433 if( nBand > 0 && nBand <= hHFA->nBands )
2434 poNode = hHFA->papoBand[nBand - 1]->poNode;
2435 else if( nBand == 0 )
2436 poNode = hHFA->poRoot;
2437 else
2438 return CE_Failure;
2439
2440 // Create the Descriptor table.
2441 // Check we have no table with this name already.
2442 HFAEntry *poEdsc_Table = poNode->GetNamedChild("GDAL_MetaData");
2443
2444 if( poEdsc_Table == nullptr || !EQUAL(poEdsc_Table->GetType(), "Edsc_Table") )
2445 poEdsc_Table =
2446 HFAEntry::New(hHFA, "GDAL_MetaData", "Edsc_Table", poNode);
2447
2448 poEdsc_Table->SetIntField("numrows", 1);
2449
2450 // Create the Binning function node. Do we really need this though?
2451 // Check it doesn't exist already.
2452 HFAEntry *poEdsc_BinFunction =
2453 poEdsc_Table->GetNamedChild("#Bin_Function#");
2454
2455 if( poEdsc_BinFunction == nullptr ||
2456 !EQUAL(poEdsc_BinFunction->GetType(), "Edsc_BinFunction") )
2457 poEdsc_BinFunction = HFAEntry::New(hHFA, "#Bin_Function#",
2458 "Edsc_BinFunction", poEdsc_Table);
2459
2460 // Because of the BaseData we have to hardcode the size.
2461 poEdsc_BinFunction->MakeData( 30 );
2462
2463 poEdsc_BinFunction->SetIntField("numBins", 1);
2464 poEdsc_BinFunction->SetStringField("binFunction", "direct");
2465 poEdsc_BinFunction->SetDoubleField("minLimit", 0.0);
2466 poEdsc_BinFunction->SetDoubleField("maxLimit", 0.0);
2467
2468 // Process each metadata item as a separate column.
2469 bool bRet = true;
2470 for( int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++ )
2471 {
2472 char *pszKey = nullptr;
2473 const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
2474 if( pszValue == nullptr )
2475 continue;
2476
2477 // Create the Edsc_Column.
2478 // Check it doesn't exist already.
2479 HFAEntry *poEdsc_Column = poEdsc_Table->GetNamedChild(pszKey);
2480
2481 if( poEdsc_Column == nullptr ||
2482 !EQUAL(poEdsc_Column->GetType(), "Edsc_Column") )
2483 poEdsc_Column =
2484 HFAEntry::New(hHFA, pszKey, "Edsc_Column", poEdsc_Table);
2485
2486 poEdsc_Column->SetIntField("numRows", 1);
2487 poEdsc_Column->SetStringField("dataType", "string");
2488 poEdsc_Column->SetIntField("maxNumChars",
2489 static_cast<GUInt32>(strlen(pszValue) + 1));
2490
2491 // Write the data out.
2492 const int nOffset =
2493 HFAAllocateSpace(hHFA, static_cast<GUInt32>(strlen(pszValue) + 1));
2494
2495 poEdsc_Column->SetIntField("columnDataPtr", nOffset);
2496
2497 bRet &= VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) >= 0;
2498 bRet &=
2499 VSIFWriteL((void *)pszValue, strlen(pszValue) + 1, 1, hHFA->fp) > 0;
2500
2501 CPLFree(pszKey);
2502 }
2503
2504 return bRet ? CE_None : CE_Failure;
2505 }
2506
2507 /************************************************************************/
2508 /* HFASetMetadata() */
2509 /************************************************************************/
2510
HFASetMetadata(HFAHandle hHFA,int nBand,char ** papszMD)2511 CPLErr HFASetMetadata( HFAHandle hHFA, int nBand, char **papszMD )
2512
2513 {
2514 char **papszGDALMD = nullptr;
2515
2516 if( CSLCount(papszMD) == 0 )
2517 return CE_None;
2518
2519 HFAEntry *poNode = nullptr;
2520
2521 if( nBand > 0 && nBand <= hHFA->nBands )
2522 poNode = hHFA->papoBand[nBand - 1]->poNode;
2523 else if( nBand == 0 )
2524 poNode = hHFA->poRoot;
2525 else
2526 return CE_Failure;
2527 #ifdef DEBUG
2528 // To please Clang Static Analyzer (CSA).
2529 if( poNode == nullptr )
2530 {
2531 CPLAssert(false);
2532 return CE_Failure;
2533 }
2534 #endif
2535
2536 // Check if the Metadata is an "known" entity which should be
2537 // stored in a better place.
2538 char *pszBinValues = nullptr;
2539 bool bCreatedHistogramParameters = false;
2540 bool bCreatedStatistics = false;
2541 const char *const *pszAuxMetaData = GetHFAAuxMetaDataList();
2542 // Check each metadata item.
2543 for( int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++ )
2544 {
2545 char *pszKey = nullptr;
2546 const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
2547 if( pszValue == nullptr )
2548 continue;
2549
2550 // Know look if its known.
2551 int i = 0; // Used after for.
2552 for( ; pszAuxMetaData[i] != nullptr; i += 4 )
2553 {
2554 if( EQUALN(pszAuxMetaData[i + 2], pszKey, strlen(pszKey)) )
2555 break;
2556 }
2557 if( pszAuxMetaData[i] != nullptr )
2558 {
2559 // Found one, get the right entry.
2560 HFAEntry *poEntry = nullptr;
2561
2562 if( strlen(pszAuxMetaData[i]) > 0 )
2563 poEntry = poNode->GetNamedChild(pszAuxMetaData[i]);
2564 else
2565 poEntry = poNode;
2566
2567 if( poEntry == nullptr && strlen(pszAuxMetaData[i + 3]) > 0 )
2568 {
2569 // Child does not yet exist --> create it,
2570 poEntry = HFAEntry::New(hHFA, pszAuxMetaData[i],
2571 pszAuxMetaData[i + 3], poNode);
2572
2573 if( STARTS_WITH_CI(pszAuxMetaData[i], "Statistics") )
2574 bCreatedStatistics = true;
2575
2576 if( STARTS_WITH_CI(pszAuxMetaData[i], "HistogramParameters") )
2577 {
2578 // A bit nasty. Need to set the string field for the object
2579 // first because the SetStringField sets the count for the
2580 // object BinFunction to the length of the string.
2581 poEntry->MakeData(70);
2582 poEntry->SetStringField("BinFunction.binFunctionType",
2583 "direct");
2584
2585 bCreatedHistogramParameters = true;
2586 }
2587 }
2588 if( poEntry == nullptr )
2589 {
2590 CPLFree(pszKey);
2591 continue;
2592 }
2593
2594 const char *pszFieldName = pszAuxMetaData[i + 1] + 1;
2595 switch( pszAuxMetaData[i + 1][0] )
2596 {
2597 case 'd':
2598 {
2599 double dfValue = CPLAtof(pszValue);
2600 poEntry->SetDoubleField(pszFieldName, dfValue);
2601 }
2602 break;
2603 case 'i':
2604 case 'l':
2605 {
2606 int nValue = atoi(pszValue);
2607 poEntry->SetIntField(pszFieldName, nValue);
2608 }
2609 break;
2610 case 's':
2611 case 'e':
2612 {
2613 poEntry->SetStringField(pszFieldName, pszValue);
2614 }
2615 break;
2616 default:
2617 CPLAssert(false);
2618 }
2619 }
2620 else if( STARTS_WITH_CI(pszKey, "STATISTICS_HISTOBINVALUES") )
2621 {
2622 CPLFree(pszBinValues);
2623 pszBinValues = CPLStrdup(pszValue);
2624 }
2625 else
2626 {
2627 papszGDALMD = CSLAddString(papszGDALMD, papszMD[iColumn]);
2628 }
2629
2630 CPLFree(pszKey);
2631 }
2632
2633 // Special case to write out the histogram.
2634 bool bRet = true;
2635 if( pszBinValues != nullptr )
2636 {
2637 HFAEntry *poEntry = poNode->GetNamedChild("HistogramParameters");
2638 if( poEntry != nullptr && bCreatedHistogramParameters )
2639 {
2640 // If this node exists we have added Histogram data -- complete with
2641 // some defaults.
2642 poEntry->SetIntField("SkipFactorX", 1);
2643 poEntry->SetIntField("SkipFactorY", 1);
2644
2645 const int nNumBins = poEntry->GetIntField("BinFunction.numBins");
2646 const double dMinLimit =
2647 poEntry->GetDoubleField("BinFunction.minLimit");
2648 const double dMaxLimit =
2649 poEntry->GetDoubleField("BinFunction.maxLimit");
2650
2651 // Fill the descriptor table - check it isn't there already.
2652 poEntry = poNode->GetNamedChild("Descriptor_Table");
2653 if( poEntry == nullptr || !EQUAL(poEntry->GetType(), "Edsc_Table") )
2654 poEntry =
2655 HFAEntry::New(hHFA, "Descriptor_Table",
2656 "Edsc_Table", poNode);
2657
2658 poEntry->SetIntField("numRows", nNumBins);
2659
2660 // Bin function.
2661 HFAEntry *poBinFunc = poEntry->GetNamedChild("#Bin_Function#");
2662 if( poBinFunc == nullptr ||
2663 !EQUAL(poBinFunc->GetType(), "Edsc_BinFunction") )
2664 poBinFunc =
2665 HFAEntry::New(hHFA, "#Bin_Function#",
2666 "Edsc_BinFunction", poEntry);
2667
2668 poBinFunc->MakeData(30);
2669 poBinFunc->SetIntField("numBins", nNumBins);
2670 poBinFunc->SetDoubleField("minLimit", dMinLimit);
2671 poBinFunc->SetDoubleField("maxLimit", dMaxLimit);
2672 // Direct for thematic layers, linear otherwise.
2673 if( STARTS_WITH_CI(poNode->GetStringField("layerType"),
2674 "thematic") )
2675 poBinFunc->SetStringField("binFunctionType", "direct");
2676 else
2677 poBinFunc->SetStringField("binFunctionType", "linear");
2678
2679 // We need a child named histogram.
2680 HFAEntry *poHisto = poEntry->GetNamedChild("Histogram");
2681 if( poHisto == nullptr || !EQUAL(poHisto->GetType(), "Edsc_Column") )
2682 poHisto =
2683 HFAEntry::New(hHFA, "Histogram", "Edsc_Column", poEntry);
2684
2685 poHisto->SetIntField("numRows", nNumBins);
2686 // Allocate space for the bin values.
2687 GUInt32 nOffset = HFAAllocateSpace(hHFA, nNumBins * 8);
2688 poHisto->SetIntField("columnDataPtr", nOffset);
2689 poHisto->SetStringField("dataType", "real");
2690 poHisto->SetIntField("maxNumChars", 0);
2691 // Write out histogram data.
2692 char *pszWork = pszBinValues;
2693 for( int nBin = 0; nBin < nNumBins; ++nBin )
2694 {
2695 char *pszEnd = strchr(pszWork, '|');
2696 if( pszEnd != nullptr )
2697 {
2698 *pszEnd = 0;
2699 bRet &=
2700 VSIFSeekL(hHFA->fp, nOffset + 8 * nBin, SEEK_SET) >= 0;
2701 double nValue = CPLAtof(pszWork);
2702 HFAStandard(8, &nValue);
2703
2704 bRet &= VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
2705 pszWork = pszEnd + 1;
2706 }
2707 }
2708 }
2709 else if( poEntry != nullptr )
2710 {
2711 // In this case, there are HistogramParameters present, but we did
2712 // not create them. However, we might be modifying them, in the case
2713 // where the data has changed and the histogram counts need to be
2714 // updated. It could be worse than that, but that is all we are
2715 // going to cope with for now. We are assuming that we did not
2716 // change any of the other stuff, like skip factors and so
2717 // forth. The main need for this case is for programs (such as
2718 // Imagine itself) which will happily modify the pixel values
2719 // without re-calculating the histogram counts.
2720 int nNumBins = poEntry->GetIntField("BinFunction.numBins");
2721 HFAEntry *poEntryDescrTbl =
2722 poNode->GetNamedChild("Descriptor_Table");
2723 HFAEntry *poHisto = nullptr;
2724 if( poEntryDescrTbl != nullptr)
2725 {
2726 poHisto = poEntryDescrTbl->GetNamedChild("Histogram");
2727 }
2728 if( poHisto != nullptr )
2729 {
2730 int nOffset = poHisto->GetIntField("columnDataPtr");
2731 // Write out histogram data.
2732 char *pszWork = pszBinValues;
2733
2734 // Check whether histogram counts were written as int or double
2735 bool bCountIsInt = true;
2736 const char *pszDataType = poHisto->GetStringField("dataType");
2737 if( STARTS_WITH_CI(pszDataType, "real") )
2738 {
2739 bCountIsInt = false;
2740 }
2741 for( int nBin = 0; nBin < nNumBins; ++nBin )
2742 {
2743 char *pszEnd = strchr(pszWork, '|');
2744 if( pszEnd != nullptr )
2745 {
2746 *pszEnd = 0;
2747 if( bCountIsInt )
2748 {
2749 // Histogram counts were written as ints, so
2750 // re-write them the same way.
2751 bRet &= VSIFSeekL(hHFA->fp, nOffset + 4 * nBin,
2752 SEEK_SET) >= 0;
2753 int nValue = atoi(pszWork);
2754 HFAStandard(4, &nValue);
2755 bRet &=
2756 VSIFWriteL((void *)&nValue, 4, 1, hHFA->fp) > 0;
2757 }
2758 else
2759 {
2760 // Histogram were written as doubles, as is now the
2761 // default behavior.
2762 bRet &= VSIFSeekL(hHFA->fp, nOffset + 8 * nBin,
2763 SEEK_SET) >= 0;
2764 double nValue = CPLAtof(pszWork);
2765 HFAStandard(8, &nValue);
2766 bRet &=
2767 VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
2768 }
2769 pszWork = pszEnd + 1;
2770 }
2771 }
2772 }
2773 }
2774 CPLFree(pszBinValues);
2775 }
2776
2777 // If we created a statistics node then try to create a
2778 // StatisticsParameters node too.
2779 if( bCreatedStatistics )
2780 {
2781 HFAEntry *poEntry =
2782 HFAEntry::New(hHFA, "StatisticsParameters",
2783 "Eimg_StatisticsParameters830", poNode);
2784
2785 poEntry->MakeData(70);
2786 // poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
2787
2788 poEntry->SetIntField("SkipFactorX", 1);
2789 poEntry->SetIntField("SkipFactorY", 1);
2790 }
2791
2792 // Write out metadata items without a special place.
2793 if( bRet && CSLCount(papszGDALMD) != 0 )
2794 {
2795 CPLErr eErr = HFASetGDALMetadata(hHFA, nBand, papszGDALMD);
2796
2797 CSLDestroy(papszGDALMD);
2798 return eErr;
2799 }
2800 else
2801 {
2802 CSLDestroy(papszGDALMD);
2803 return CE_Failure;
2804 }
2805 }
2806
2807 /************************************************************************/
2808 /* HFAGetIGEFilename() */
2809 /* */
2810 /* Returns the .ige filename if one is associated with this */
2811 /* object. For files not newly created we need to scan the */
2812 /* bands for spill files. Presumably there will only be one. */
2813 /* */
2814 /* NOTE: Returns full path, not just the filename portion. */
2815 /************************************************************************/
2816
HFAGetIGEFilename(HFAHandle hHFA)2817 const char *HFAGetIGEFilename( HFAHandle hHFA )
2818
2819 {
2820 if( hHFA->pszIGEFilename == nullptr )
2821 {
2822 std::vector<HFAEntry *> apoDMSList =
2823 hHFA->poRoot->FindChildren(nullptr, "ImgExternalRaster");
2824
2825 HFAEntry *poDMS = apoDMSList.empty() ? nullptr : apoDMSList[0];
2826
2827 // Get the IGE filename from if we have an ExternalRasterDMS.
2828 if( poDMS )
2829 {
2830 const char *pszRawFilename =
2831 poDMS->GetStringField("fileName.string");
2832
2833 if( pszRawFilename != nullptr )
2834 {
2835 VSIStatBufL sStatBuf;
2836 CPLString osFullFilename =
2837 CPLFormFilename(hHFA->pszPath, pszRawFilename, nullptr);
2838
2839 if( VSIStatL(osFullFilename, &sStatBuf) != 0 )
2840 {
2841 const CPLString osExtension =
2842 CPLGetExtension(pszRawFilename);
2843 const CPLString osBasename =
2844 CPLGetBasename(hHFA->pszFilename);
2845 osFullFilename =
2846 CPLFormFilename(hHFA->pszPath, osBasename, osExtension);
2847
2848 if( VSIStatL(osFullFilename, &sStatBuf) == 0 )
2849 hHFA->pszIGEFilename = CPLStrdup(
2850 CPLFormFilename(nullptr, osBasename, osExtension));
2851 else
2852 hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
2853 }
2854 else
2855 {
2856 hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
2857 }
2858 }
2859 }
2860 }
2861
2862 // Return the full filename.
2863 if( hHFA->pszIGEFilename )
2864 return CPLFormFilename(hHFA->pszPath, hHFA->pszIGEFilename, nullptr);
2865
2866 return nullptr;
2867 }
2868
2869 /************************************************************************/
2870 /* HFACreateSpillStack() */
2871 /* */
2872 /* Create a new stack of raster layers in the spill (.ige) */
2873 /* file. Create the spill file if it didn't exist before. */
2874 /************************************************************************/
2875
HFACreateSpillStack(HFAInfo_t * psInfo,int nXSize,int nYSize,int nLayers,int nBlockSize,EPTType eDataType,GIntBig * pnValidFlagsOffset,GIntBig * pnDataOffset)2876 bool HFACreateSpillStack( HFAInfo_t *psInfo, int nXSize, int nYSize,
2877 int nLayers, int nBlockSize, EPTType eDataType,
2878 GIntBig *pnValidFlagsOffset,
2879 GIntBig *pnDataOffset )
2880
2881 {
2882 // Form .ige filename.
2883 if( nBlockSize <= 0 )
2884 {
2885 CPLError(CE_Failure, CPLE_IllegalArg,
2886 "HFACreateSpillStack: nBlockXSize < 0");
2887 return false;
2888 }
2889
2890 if( psInfo->pszIGEFilename == nullptr )
2891 {
2892 if( EQUAL(CPLGetExtension(psInfo->pszFilename), "rrd") )
2893 psInfo->pszIGEFilename =
2894 CPLStrdup( CPLResetExtension(psInfo->pszFilename, "rde") );
2895 else if( EQUAL(CPLGetExtension(psInfo->pszFilename), "aux") )
2896 psInfo->pszIGEFilename =
2897 CPLStrdup(CPLResetExtension(psInfo->pszFilename, "axe"));
2898 else
2899 psInfo->pszIGEFilename =
2900 CPLStrdup(CPLResetExtension(psInfo->pszFilename, "ige"));
2901 }
2902
2903 char *pszFullFilename = CPLStrdup(
2904 CPLFormFilename(psInfo->pszPath, psInfo->pszIGEFilename, nullptr));
2905
2906 // Try and open it. If we fail, create it and write the magic header.
2907 static const char *const pszMagick = "ERDAS_IMG_EXTERNAL_RASTER";
2908
2909 bool bRet = true;
2910 VSILFILE *fpVSIL = VSIFOpenL(pszFullFilename, "r+b");
2911 if( fpVSIL == nullptr )
2912 {
2913 fpVSIL = VSIFOpenL(pszFullFilename, "w+");
2914 if( fpVSIL == nullptr )
2915 {
2916 CPLError(CE_Failure, CPLE_OpenFailed,
2917 "Failed to create spill file %s.\n%s",
2918 psInfo->pszIGEFilename, VSIStrerror(errno));
2919 return false;
2920 }
2921
2922 bRet &=
2923 VSIFWriteL((void *)pszMagick, strlen(pszMagick) + 1, 1, fpVSIL) > 0;
2924 }
2925
2926 CPLFree(pszFullFilename);
2927
2928 // Work out some details about the tiling scheme.
2929 const int nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2930 const int nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2931 // const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
2932 const int nBytesPerBlock =
2933 (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
2934
2935 const int nBytesPerRow = (nBlocksPerRow + 7) / 8;
2936 const int nBlockMapSize = nBytesPerRow * nBlocksPerColumn;
2937 // const int iFlagsSize = nBlockMapSize + 20;
2938
2939 // Write stack prefix information.
2940 bRet &= VSIFSeekL(fpVSIL, 0, SEEK_END) >= 0;
2941
2942 GByte bUnknown = 1;
2943 bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
2944
2945 GInt32 nValue32 = nLayers;
2946 HFAStandard(4, &nValue32);
2947 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2948 nValue32 = nXSize;
2949 HFAStandard(4, &nValue32);
2950 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2951 nValue32 = nYSize;
2952 HFAStandard(4, &nValue32);
2953 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2954 nValue32 = nBlockSize;
2955 HFAStandard(4, &nValue32);
2956 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2957 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2958 bUnknown = 3;
2959 bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
2960 bUnknown = 0;
2961 bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
2962
2963 // Write out ValidFlags section(s).
2964 *pnValidFlagsOffset = VSIFTellL(fpVSIL);
2965
2966 unsigned char *pabyBlockMap =
2967 static_cast<unsigned char *>(VSI_MALLOC_VERBOSE(nBlockMapSize));
2968 if( pabyBlockMap == nullptr )
2969 {
2970 CPL_IGNORE_RET_VAL(VSIFCloseL(fpVSIL));
2971 return false;
2972 }
2973
2974 memset(pabyBlockMap, 0xff, nBlockMapSize);
2975 for( int iBand = 0; iBand < nLayers; iBand++ )
2976 {
2977 nValue32 = 1; // Unknown
2978 HFAStandard(4, &nValue32);
2979 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2980 nValue32 = 0; // Unknown
2981 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2982 nValue32 = nBlocksPerColumn;
2983 HFAStandard(4, &nValue32);
2984 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2985 nValue32 = nBlocksPerRow;
2986 HFAStandard(4, &nValue32);
2987 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2988 nValue32 = 0x30000; // Unknown
2989 HFAStandard(4, &nValue32);
2990 bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
2991
2992 const int iRemainder = nBlocksPerRow % 8;
2993 CPLDebug("HFACreate",
2994 "Block map size %d, bytes per row %d, remainder %d.",
2995 nBlockMapSize, nBytesPerRow, iRemainder);
2996 if( iRemainder )
2997 {
2998 for( int i = nBytesPerRow - 1; i < nBlockMapSize; i += nBytesPerRow )
2999 pabyBlockMap[i] = static_cast<GByte>((1 << iRemainder) - 1);
3000 }
3001
3002 bRet &= VSIFWriteL(pabyBlockMap, nBlockMapSize, 1, fpVSIL) > 0;
3003 }
3004 CPLFree(pabyBlockMap);
3005 pabyBlockMap = nullptr;
3006
3007 // Extend the file to account for all the imagery space.
3008 const GIntBig nTileDataSize = static_cast<GIntBig>(nBytesPerBlock) *
3009 nBlocksPerRow * nBlocksPerColumn * nLayers;
3010
3011 *pnDataOffset = VSIFTellL(fpVSIL);
3012
3013 if( !bRet ||
3014 VSIFTruncateL(fpVSIL, nTileDataSize + *pnDataOffset) != 0 )
3015 {
3016 CPLError(CE_Failure, CPLE_FileIO,
3017 "Failed to extend %s to full size (" CPL_FRMT_GIB " bytes), "
3018 "likely out of disk space.\n%s",
3019 psInfo->pszIGEFilename,
3020 nTileDataSize + *pnDataOffset,
3021 VSIStrerror(errno));
3022
3023 CPL_IGNORE_RET_VAL(VSIFCloseL( fpVSIL ));
3024 return false;
3025 }
3026
3027 if( VSIFCloseL(fpVSIL) != 0 )
3028 return false;
3029
3030 return true;
3031 }
3032
3033 /************************************************************************/
3034 /* HFAReadAndValidatePoly() */
3035 /************************************************************************/
3036
HFAReadAndValidatePoly(HFAEntry * poTarget,const char * pszName,Efga_Polynomial * psRetPoly)3037 static bool HFAReadAndValidatePoly( HFAEntry *poTarget,
3038 const char *pszName,
3039 Efga_Polynomial *psRetPoly )
3040
3041 {
3042 memset(psRetPoly, 0, sizeof(Efga_Polynomial));
3043
3044 CPLString osFldName;
3045 osFldName.Printf("%sorder", pszName);
3046 psRetPoly->order = poTarget->GetIntField(osFldName);
3047
3048 if( psRetPoly->order < 1 || psRetPoly->order > 3 )
3049 return false;
3050
3051 // Validate that things are in a "well known" form.
3052 osFldName.Printf("%snumdimtransform", pszName);
3053 const int numdimtransform = poTarget->GetIntField(osFldName);
3054
3055 osFldName.Printf("%snumdimpolynomial", pszName);
3056 const int numdimpolynomial = poTarget->GetIntField(osFldName);
3057
3058 osFldName.Printf("%stermcount", pszName);
3059 const int termcount = poTarget->GetIntField(osFldName);
3060
3061 if( numdimtransform != 2 || numdimpolynomial != 2 )
3062 return false;
3063
3064 if( (psRetPoly->order == 1 && termcount != 3) ||
3065 (psRetPoly->order == 2 && termcount != 6) ||
3066 (psRetPoly->order == 3 && termcount != 10) )
3067 return false;
3068
3069 // We don't check the exponent organization for now. Hopefully
3070 // it is always standard.
3071
3072 // Get coefficients.
3073 for( int i = 0; i < termcount*2 - 2; i++ )
3074 {
3075 osFldName.Printf("%spolycoefmtx[%d]", pszName, i);
3076 psRetPoly->polycoefmtx[i] = poTarget->GetDoubleField(osFldName);
3077 }
3078
3079 for( int i = 0; i < 2; i++ )
3080 {
3081 osFldName.Printf("%spolycoefvector[%d]", pszName, i);
3082 psRetPoly->polycoefvector[i] = poTarget->GetDoubleField(osFldName);
3083 }
3084
3085 return true;
3086 }
3087
3088 /************************************************************************/
3089 /* HFAReadXFormStack() */
3090 /************************************************************************/
3091
HFAReadXFormStack(HFAHandle hHFA,Efga_Polynomial ** ppasPolyListForward,Efga_Polynomial ** ppasPolyListReverse)3092 int HFAReadXFormStack( HFAHandle hHFA,
3093 Efga_Polynomial **ppasPolyListForward,
3094 Efga_Polynomial **ppasPolyListReverse )
3095
3096 {
3097 if( hHFA->nBands == 0 )
3098 return 0;
3099
3100 // Get the HFA node.
3101 HFAEntry *poXFormHeader =
3102 hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm");
3103 if( poXFormHeader == nullptr )
3104 return 0;
3105
3106 // Loop over children, collecting XForms.
3107 int nStepCount = 0;
3108 *ppasPolyListForward = nullptr;
3109 *ppasPolyListReverse = nullptr;
3110
3111 for( HFAEntry *poXForm = poXFormHeader->GetChild();
3112 poXForm != nullptr;
3113 poXForm = poXForm->GetNext() )
3114 {
3115 bool bSuccess = false;
3116 Efga_Polynomial sForward;
3117 Efga_Polynomial sReverse;
3118 memset(&sForward, 0, sizeof(sForward));
3119 memset(&sReverse, 0, sizeof(sReverse));
3120
3121 if( EQUAL(poXForm->GetType(), "Efga_Polynomial") )
3122 {
3123 bSuccess = HFAReadAndValidatePoly(poXForm, "", &sForward);
3124
3125 if( bSuccess )
3126 {
3127 double adfGT[6] = {
3128 sForward.polycoefvector[0],
3129 sForward.polycoefmtx[0],
3130 sForward.polycoefmtx[2],
3131 sForward.polycoefvector[1],
3132 sForward.polycoefmtx[1],
3133 sForward.polycoefmtx[3]
3134 };
3135
3136 double adfInvGT[6] = {};
3137 bSuccess = HFAInvGeoTransform(adfGT, adfInvGT);
3138 if( !bSuccess )
3139 memset(adfInvGT, 0, sizeof(adfInvGT));
3140
3141 sReverse.order = sForward.order;
3142 sReverse.polycoefvector[0] = adfInvGT[0];
3143 sReverse.polycoefmtx[0] = adfInvGT[1];
3144 sReverse.polycoefmtx[2] = adfInvGT[2];
3145 sReverse.polycoefvector[1] = adfInvGT[3];
3146 sReverse.polycoefmtx[1] = adfInvGT[4];
3147 sReverse.polycoefmtx[3] = adfInvGT[5];
3148 }
3149 }
3150 else if( EQUAL(poXForm->GetType(), "GM_PolyPair") )
3151 {
3152 bSuccess =
3153 HFAReadAndValidatePoly(poXForm, "forward.", &sForward) &&
3154 HFAReadAndValidatePoly(poXForm, "reverse.", &sReverse);
3155 }
3156
3157 if( bSuccess )
3158 {
3159 nStepCount++;
3160 *ppasPolyListForward = static_cast<Efga_Polynomial *>(CPLRealloc(
3161 *ppasPolyListForward, sizeof(Efga_Polynomial) * nStepCount));
3162 memcpy(*ppasPolyListForward + nStepCount - 1, &sForward,
3163 sizeof(sForward));
3164
3165 *ppasPolyListReverse = static_cast<Efga_Polynomial *>(CPLRealloc(
3166 *ppasPolyListReverse, sizeof(Efga_Polynomial) * nStepCount));
3167 memcpy(*ppasPolyListReverse + nStepCount - 1, &sReverse,
3168 sizeof(sReverse));
3169 }
3170 }
3171
3172 return nStepCount;
3173 }
3174
3175 /************************************************************************/
3176 /* HFAEvaluateXFormStack() */
3177 /************************************************************************/
3178
HFAEvaluateXFormStack(int nStepCount,int bForward,Efga_Polynomial * pasPolyList,double * pdfX,double * pdfY)3179 int HFAEvaluateXFormStack( int nStepCount, int bForward,
3180 Efga_Polynomial *pasPolyList,
3181 double *pdfX, double *pdfY )
3182
3183 {
3184 for( int iStep = 0; iStep < nStepCount; iStep++ )
3185 {
3186 const Efga_Polynomial *psStep =
3187 bForward ? pasPolyList + iStep
3188 : pasPolyList + nStepCount - iStep - 1;
3189
3190 double dfXOut = 0.0;
3191 double dfYOut = 0.0;
3192 if( psStep->order == 1 )
3193 {
3194 dfXOut = psStep->polycoefvector[0] +
3195 psStep->polycoefmtx[0] * *pdfX +
3196 psStep->polycoefmtx[2] * *pdfY;
3197
3198 dfYOut = psStep->polycoefvector[1] +
3199 psStep->polycoefmtx[1] * *pdfX +
3200 psStep->polycoefmtx[3] * *pdfY;
3201
3202 *pdfX = dfXOut;
3203 *pdfY = dfYOut;
3204 }
3205 else if( psStep->order == 2 )
3206 {
3207 dfXOut = psStep->polycoefvector[0] +
3208 psStep->polycoefmtx[0] * *pdfX +
3209 psStep->polycoefmtx[2] * *pdfY +
3210 psStep->polycoefmtx[4] * *pdfX * *pdfX +
3211 psStep->polycoefmtx[6] * *pdfX * *pdfY +
3212 psStep->polycoefmtx[8] * *pdfY * *pdfY;
3213 dfYOut = psStep->polycoefvector[1] +
3214 psStep->polycoefmtx[1] * *pdfX +
3215 psStep->polycoefmtx[3] * *pdfY +
3216 psStep->polycoefmtx[5] * *pdfX * *pdfX +
3217 psStep->polycoefmtx[7] * *pdfX * *pdfY +
3218 psStep->polycoefmtx[9] * *pdfY * *pdfY;
3219
3220 *pdfX = dfXOut;
3221 *pdfY = dfYOut;
3222 }
3223 else if( psStep->order == 3 )
3224 {
3225 dfXOut = psStep->polycoefvector[0] +
3226 psStep->polycoefmtx[0] * *pdfX +
3227 psStep->polycoefmtx[2] * *pdfY +
3228 psStep->polycoefmtx[4] * *pdfX * *pdfX +
3229 psStep->polycoefmtx[6] * *pdfX * *pdfY +
3230 psStep->polycoefmtx[8] * *pdfY * *pdfY +
3231 psStep->polycoefmtx[10] * *pdfX * *pdfX * *pdfX +
3232 psStep->polycoefmtx[12] * *pdfX * *pdfX * *pdfY +
3233 psStep->polycoefmtx[14] * *pdfX * *pdfY * *pdfY +
3234 psStep->polycoefmtx[16] * *pdfY * *pdfY * *pdfY;
3235 dfYOut = psStep->polycoefvector[1] +
3236 psStep->polycoefmtx[1] * *pdfX +
3237 psStep->polycoefmtx[3] * *pdfY +
3238 psStep->polycoefmtx[5] * *pdfX * *pdfX +
3239 psStep->polycoefmtx[7] * *pdfX * *pdfY +
3240 psStep->polycoefmtx[9] * *pdfY * *pdfY +
3241 psStep->polycoefmtx[11] * *pdfX * *pdfX * *pdfX +
3242 psStep->polycoefmtx[13] * *pdfX * *pdfX * *pdfY +
3243 psStep->polycoefmtx[15] * *pdfX * *pdfY * *pdfY +
3244 psStep->polycoefmtx[17] * *pdfY * *pdfY * *pdfY;
3245
3246 *pdfX = dfXOut;
3247 *pdfY = dfYOut;
3248 }
3249 else
3250 return FALSE;
3251 }
3252
3253 return TRUE;
3254 }
3255
3256 /************************************************************************/
3257 /* HFAWriteXFormStack() */
3258 /************************************************************************/
3259
HFAWriteXFormStack(HFAHandle hHFA,int nBand,int nXFormCount,Efga_Polynomial ** ppasPolyListForward,Efga_Polynomial ** ppasPolyListReverse)3260 CPLErr HFAWriteXFormStack( HFAHandle hHFA, int nBand, int nXFormCount,
3261 Efga_Polynomial **ppasPolyListForward,
3262 Efga_Polynomial **ppasPolyListReverse )
3263
3264 {
3265 if( nXFormCount == 0 )
3266 return CE_None;
3267
3268 if( ppasPolyListForward[0]->order != 1 )
3269 {
3270 CPLError(
3271 CE_Failure, CPLE_AppDefined,
3272 "For now HFAWriteXFormStack() only supports order 1 polynomials");
3273 return CE_Failure;
3274 }
3275
3276 if( nBand < 0 || nBand > hHFA->nBands )
3277 return CE_Failure;
3278
3279 // If no band number is provided, operate on all bands.
3280 if( nBand == 0 )
3281 {
3282 for( nBand = 1; nBand <= hHFA->nBands; nBand++ )
3283 {
3284 CPLErr eErr =
3285 HFAWriteXFormStack(hHFA, nBand, nXFormCount,
3286 ppasPolyListForward,
3287 ppasPolyListReverse);
3288 if( eErr != CE_None )
3289 return eErr;
3290 }
3291
3292 return CE_None;
3293 }
3294
3295 // Fetch our band node.
3296 HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
3297 HFAEntry *poXFormHeader = poBandNode->GetNamedChild("MapToPixelXForm");
3298 if( poXFormHeader == nullptr )
3299 {
3300 poXFormHeader = HFAEntry::New(hHFA, "MapToPixelXForm",
3301 "Exfr_GenericXFormHeader", poBandNode);
3302 poXFormHeader->MakeData(23);
3303 poXFormHeader->SetPosition();
3304 poXFormHeader->SetStringField("titleList.string", "Affine");
3305 }
3306
3307 // Loop over XForms.
3308 for( int iXForm = 0; iXForm < nXFormCount; iXForm++ )
3309 {
3310 Efga_Polynomial *psForward = *ppasPolyListForward + iXForm;
3311 CPLString osXFormName;
3312 osXFormName.Printf("XForm%d", iXForm);
3313
3314 HFAEntry *poXForm = poXFormHeader->GetNamedChild(osXFormName);
3315
3316 if( poXForm == nullptr )
3317 {
3318 poXForm = HFAEntry::New(hHFA, osXFormName, "Efga_Polynomial",
3319 poXFormHeader);
3320 poXForm->MakeData(136);
3321 poXForm->SetPosition();
3322 }
3323
3324 poXForm->SetIntField("order", 1);
3325 poXForm->SetIntField("numdimtransform", 2);
3326 poXForm->SetIntField("numdimpolynomial", 2);
3327 poXForm->SetIntField("termcount", 3);
3328 poXForm->SetIntField("exponentlist[0]", 0);
3329 poXForm->SetIntField("exponentlist[1]", 0);
3330 poXForm->SetIntField("exponentlist[2]", 1);
3331 poXForm->SetIntField("exponentlist[3]", 0);
3332 poXForm->SetIntField("exponentlist[4]", 0);
3333 poXForm->SetIntField("exponentlist[5]", 1);
3334
3335 poXForm->SetIntField("polycoefmtx[-3]", EPT_f64);
3336 poXForm->SetIntField("polycoefmtx[-2]", 2);
3337 poXForm->SetIntField("polycoefmtx[-1]", 2);
3338 poXForm->SetDoubleField("polycoefmtx[0]", psForward->polycoefmtx[0]);
3339 poXForm->SetDoubleField("polycoefmtx[1]", psForward->polycoefmtx[1]);
3340 poXForm->SetDoubleField("polycoefmtx[2]", psForward->polycoefmtx[2]);
3341 poXForm->SetDoubleField("polycoefmtx[3]", psForward->polycoefmtx[3]);
3342
3343 poXForm->SetIntField("polycoefvector[-3]", EPT_f64);
3344 poXForm->SetIntField("polycoefvector[-2]", 1);
3345 poXForm->SetIntField("polycoefvector[-1]", 2);
3346 poXForm->SetDoubleField("polycoefvector[0]",
3347 psForward->polycoefvector[0]);
3348 poXForm->SetDoubleField("polycoefvector[1]",
3349 psForward->polycoefvector[1]);
3350 }
3351
3352 return CE_None;
3353 }
3354
3355 /************************************************************************/
3356 /* HFAReadCameraModel() */
3357 /************************************************************************/
3358
HFAReadCameraModel(HFAHandle hHFA)3359 char **HFAReadCameraModel( HFAHandle hHFA )
3360
3361 {
3362 if( hHFA->nBands == 0 )
3363 return nullptr;
3364
3365 // Get the camera model node, and confirm its type.
3366 HFAEntry *poXForm =
3367 hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
3368 if( poXForm == nullptr )
3369 return nullptr;
3370
3371 if( !EQUAL(poXForm->GetType(), "Camera_ModelX") )
3372 return nullptr;
3373
3374 // Convert the values to metadata.
3375 char **papszMD = nullptr;
3376 static const char * const apszFields[] = {
3377 "direction", "refType", "demsource", "PhotoDirection", "RotationSystem",
3378 "demfilename", "demzunits",
3379 "forSrcAffine[0]", "forSrcAffine[1]", "forSrcAffine[2]",
3380 "forSrcAffine[3]", "forSrcAffine[4]", "forSrcAffine[5]",
3381 "forDstAffine[0]", "forDstAffine[1]", "forDstAffine[2]",
3382 "forDstAffine[3]", "forDstAffine[4]", "forDstAffine[5]",
3383 "invSrcAffine[0]", "invSrcAffine[1]", "invSrcAffine[2]",
3384 "invSrcAffine[3]", "invSrcAffine[4]", "invSrcAffine[5]",
3385 "invDstAffine[0]", "invDstAffine[1]", "invDstAffine[2]",
3386 "invDstAffine[3]", "invDstAffine[4]", "invDstAffine[5]",
3387 "z_mean", "lat0", "lon0",
3388 "coeffs[0]", "coeffs[1]", "coeffs[2]",
3389 "coeffs[3]", "coeffs[4]", "coeffs[5]",
3390 "coeffs[6]", "coeffs[7]", "coeffs[8]",
3391 "LensDistortion[0]", "LensDistortion[1]", "LensDistortion[2]",
3392 nullptr };
3393
3394 const char *pszValue = nullptr;
3395 for( int i = 0; apszFields[i] != nullptr; i++ )
3396 {
3397 pszValue = poXForm->GetStringField(apszFields[i]);
3398 if( pszValue == nullptr )
3399 pszValue = "";
3400
3401 papszMD = CSLSetNameValue(papszMD, apszFields[i], pszValue);
3402 }
3403
3404 // Create a pseudo-entry for the MIFObject with the outputProjection.
3405 HFAEntry *poProjInfo =
3406 HFAEntry::BuildEntryFromMIFObject(poXForm, "outputProjection");
3407 if( poProjInfo )
3408 {
3409 // Fetch the datum.
3410 Eprj_Datum sDatum;
3411
3412 memset(&sDatum, 0, sizeof(sDatum));
3413
3414 sDatum.datumname =
3415 (char *)poProjInfo->GetStringField("earthModel.datum.datumname");
3416
3417 const int nDatumType = poProjInfo->GetIntField("earthModel.datum.type");
3418 if( nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE )
3419 {
3420 CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
3421 sDatum.type = EPRJ_DATUM_NONE;
3422 }
3423 else
3424 {
3425 sDatum.type = static_cast<Eprj_DatumType>(nDatumType);
3426 }
3427
3428 for( int i = 0; i < 7; i++ )
3429 {
3430 char szFieldName[60] = {};
3431
3432 snprintf(szFieldName, sizeof(szFieldName),
3433 "earthModel.datum.params[%d]", i);
3434 sDatum.params[i] = poProjInfo->GetDoubleField(szFieldName);
3435 }
3436
3437 sDatum.gridname =
3438 (char *)poProjInfo->GetStringField("earthModel.datum.gridname");
3439
3440 // Fetch the projection parameters.
3441 Eprj_ProParameters sPro;
3442
3443 memset(&sPro, 0, sizeof(sPro));
3444
3445 sPro.proType =
3446 (Eprj_ProType)poProjInfo->GetIntField("projectionObject.proType");
3447 sPro.proNumber = poProjInfo->GetIntField("projectionObject.proNumber");
3448 sPro.proExeName =
3449 (char *)poProjInfo->GetStringField("projectionObject.proExeName");
3450 sPro.proName =
3451 (char *)poProjInfo->GetStringField("projectionObject.proName");
3452 sPro.proZone = poProjInfo->GetIntField("projectionObject.proZone");
3453
3454 for( int i = 0; i < 15; i++ )
3455 {
3456 char szFieldName[40] = {};
3457
3458 snprintf(szFieldName, sizeof(szFieldName),
3459 "projectionObject.proParams[%d]", i);
3460 sPro.proParams[i] = poProjInfo->GetDoubleField(szFieldName);
3461 }
3462
3463 // Fetch the spheroid.
3464 sPro.proSpheroid.sphereName = (char *)poProjInfo->GetStringField(
3465 "earthModel.proSpheroid.sphereName");
3466 sPro.proSpheroid.a =
3467 poProjInfo->GetDoubleField("earthModel.proSpheroid.a");
3468 sPro.proSpheroid.b =
3469 poProjInfo->GetDoubleField("earthModel.proSpheroid.b");
3470 sPro.proSpheroid.eSquared =
3471 poProjInfo->GetDoubleField("earthModel.proSpheroid.eSquared");
3472 sPro.proSpheroid.radius =
3473 poProjInfo->GetDoubleField("earthModel.proSpheroid.radius");
3474
3475 // Fetch the projection info.
3476 // poProjInfo->DumpFieldValues( stdout, "" );
3477
3478 char *pszProjection = HFAPCSStructToWKT(&sDatum, &sPro, nullptr, nullptr);
3479
3480 if( pszProjection )
3481 {
3482 papszMD =
3483 CSLSetNameValue(papszMD, "outputProjection", pszProjection);
3484 CPLFree(pszProjection);
3485 }
3486
3487 delete poProjInfo;
3488 }
3489
3490 // Fetch the horizontal units.
3491 pszValue = poXForm->GetStringField("outputHorizontalUnits.string");
3492 if( pszValue == nullptr )
3493 pszValue = "";
3494
3495 papszMD = CSLSetNameValue(papszMD, "outputHorizontalUnits", pszValue);
3496
3497 // Fetch the elevationinfo.
3498 HFAEntry *poElevInfo =
3499 HFAEntry::BuildEntryFromMIFObject(poXForm, "outputElevationInfo");
3500 if( poElevInfo )
3501 {
3502 // poElevInfo->DumpFieldValues( stdout, "" );
3503
3504 if( poElevInfo->GetDataSize() != 0 )
3505 {
3506 static const char *const apszEFields[] = {
3507 "verticalDatum.datumname",
3508 "verticalDatum.type",
3509 "elevationUnit",
3510 "elevationType",
3511 nullptr };
3512
3513 for( int i = 0; apszEFields[i] != nullptr; i++ )
3514 {
3515 pszValue = poElevInfo->GetStringField(apszEFields[i]);
3516 if( pszValue == nullptr )
3517 pszValue = "";
3518
3519 papszMD = CSLSetNameValue(papszMD, apszEFields[i], pszValue);
3520 }
3521 }
3522
3523 delete poElevInfo;
3524 }
3525
3526 return papszMD;
3527 }
3528
3529 /************************************************************************/
3530 /* HFAReadElevationUnit() */
3531 /************************************************************************/
3532
HFAReadElevationUnit(HFAHandle hHFA,int iBand)3533 const char *HFAReadElevationUnit( HFAHandle hHFA, int iBand )
3534 {
3535 if( hHFA->nBands <= iBand )
3536 return nullptr;
3537
3538 HFABand *poBand(hHFA->papoBand[iBand]);
3539 if( poBand == nullptr || poBand->poNode == nullptr )
3540 {
3541 return nullptr;
3542 }
3543 HFAEntry *poElevInfo = poBand->poNode->GetNamedChild("Elevation_Info");
3544 if( poElevInfo == nullptr )
3545 {
3546 return nullptr;
3547 }
3548 return poElevInfo->GetStringField("elevationUnit");
3549 }
3550
3551 /************************************************************************/
3552 /* HFASetGeoTransform() */
3553 /* */
3554 /* Set a MapInformation and XForm block. Allows for rotated */
3555 /* and shared geotransforms. */
3556 /************************************************************************/
3557
HFASetGeoTransform(HFAHandle hHFA,const char * pszProName,const char * pszUnits,double * padfGeoTransform)3558 CPLErr HFASetGeoTransform( HFAHandle hHFA,
3559 const char *pszProName,
3560 const char *pszUnits,
3561 double *padfGeoTransform )
3562
3563 {
3564 // Write MapInformation.
3565 for( int nBand = 1; nBand <= hHFA->nBands; nBand++ )
3566 {
3567 HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
3568
3569 HFAEntry *poMI = poBandNode->GetNamedChild("MapInformation");
3570 if( poMI == nullptr )
3571 {
3572 poMI = HFAEntry::New(hHFA, "MapInformation", "Eimg_MapInformation",
3573 poBandNode);
3574 poMI->MakeData(
3575 static_cast<int>(18 + strlen(pszProName) + strlen(pszUnits)));
3576 poMI->SetPosition();
3577 }
3578
3579 poMI->SetStringField("projection.string", pszProName);
3580 poMI->SetStringField("units.string", pszUnits);
3581 }
3582
3583 // Write XForm.
3584 double adfAdjTransform[6] = {};
3585
3586 // Offset by half pixel.
3587
3588 memcpy(adfAdjTransform, padfGeoTransform, sizeof(double) * 6);
3589 adfAdjTransform[0] += adfAdjTransform[1] * 0.5;
3590 adfAdjTransform[0] += adfAdjTransform[2] * 0.5;
3591 adfAdjTransform[3] += adfAdjTransform[4] * 0.5;
3592 adfAdjTransform[3] += adfAdjTransform[5] * 0.5;
3593
3594 // Invert.
3595 double adfRevTransform[6] = {};
3596 if( !HFAInvGeoTransform( adfAdjTransform, adfRevTransform ) )
3597 memset(adfRevTransform, 0, sizeof(adfRevTransform));
3598
3599 // Assign to polynomial object.
3600
3601 Efga_Polynomial sForward;
3602 memset(&sForward, 0, sizeof(sForward));
3603 Efga_Polynomial *psForward = &sForward;
3604 sForward.order = 1;
3605 sForward.polycoefvector[0] = adfRevTransform[0];
3606 sForward.polycoefmtx[0] = adfRevTransform[1];
3607 sForward.polycoefmtx[1] = adfRevTransform[4];
3608 sForward.polycoefvector[1] = adfRevTransform[3];
3609 sForward.polycoefmtx[2] = adfRevTransform[2];
3610 sForward.polycoefmtx[3] = adfRevTransform[5];
3611
3612 Efga_Polynomial sReverse = sForward;
3613 Efga_Polynomial *psReverse = &sReverse;
3614
3615 return HFAWriteXFormStack(hHFA, 0, 1, &psForward, &psReverse);
3616 }
3617
3618 /************************************************************************/
3619 /* HFARenameReferences() */
3620 /* */
3621 /* Rename references in this .img file from the old basename to */
3622 /* a new basename. This should be passed on to .aux and .rrd */
3623 /* files and should include references to .aux, .rrd and .ige. */
3624 /************************************************************************/
3625
HFARenameReferences(HFAHandle hHFA,const char * pszNewBase,const char * pszOldBase)3626 CPLErr HFARenameReferences( HFAHandle hHFA,
3627 const char *pszNewBase,
3628 const char *pszOldBase )
3629
3630 {
3631 // Handle RRDNamesList updates.
3632 std::vector<HFAEntry *> apoNodeList =
3633 hHFA->poRoot->FindChildren("RRDNamesList", nullptr);
3634
3635 for( size_t iNode = 0; iNode < apoNodeList.size(); iNode++ )
3636 {
3637 HFAEntry *poRRDNL = apoNodeList[iNode];
3638 std::vector<CPLString> aosNL;
3639
3640 // Collect all the existing names.
3641 const int nNameCount = poRRDNL->GetFieldCount("nameList");
3642
3643 CPLString osAlgorithm = poRRDNL->GetStringField("algorithm.string");
3644 for( int i = 0; i < nNameCount; i++ )
3645 {
3646 CPLString osFN;
3647 osFN.Printf("nameList[%d].string", i);
3648 aosNL.push_back(poRRDNL->GetStringField(osFN));
3649 }
3650
3651 // Adjust the names to the new form.
3652 for( int i = 0; i < nNameCount; i++ )
3653 {
3654 if( strncmp(aosNL[i], pszOldBase, strlen(pszOldBase)) == 0 )
3655 {
3656 CPLString osNew = pszNewBase;
3657 osNew += aosNL[i].c_str() + strlen(pszOldBase);
3658 aosNL[i] = osNew;
3659 }
3660 }
3661
3662 // Try to make sure the RRDNamesList is big enough to hold the
3663 // adjusted name list.
3664 if( strlen(pszNewBase) > strlen(pszOldBase) )
3665 {
3666 CPLDebug("HFA", "Growing RRDNamesList to hold new names");
3667 poRRDNL->MakeData(static_cast<int>(
3668 poRRDNL->GetDataSize() +
3669 nNameCount * (strlen(pszNewBase) - strlen(pszOldBase))));
3670 }
3671
3672 // Initialize the whole thing to zeros for a clean start.
3673 memset(poRRDNL->GetData(), 0, poRRDNL->GetDataSize());
3674
3675 // Write the updates back to the file.
3676 poRRDNL->SetStringField("algorithm.string", osAlgorithm);
3677 for( int i = 0; i < nNameCount; i++ )
3678 {
3679 CPLString osFN;
3680 osFN.Printf("nameList[%d].string", i);
3681 poRRDNL->SetStringField(osFN, aosNL[i]);
3682 }
3683 }
3684
3685 // Spill file references.
3686 apoNodeList =
3687 hHFA->poRoot->FindChildren("ExternalRasterDMS", "ImgExternalRaster");
3688
3689 for( size_t iNode = 0; iNode < apoNodeList.size(); iNode++ )
3690 {
3691 HFAEntry *poERDMS = apoNodeList[iNode];
3692
3693 if( poERDMS == nullptr )
3694 continue;
3695
3696 // Fetch all existing values.
3697 CPLString osFileName = poERDMS->GetStringField("fileName.string");
3698
3699 GInt32 anValidFlagsOffset[2] = {
3700 poERDMS->GetIntField("layerStackValidFlagsOffset[0]"),
3701 poERDMS->GetIntField("layerStackValidFlagsOffset[1]")
3702 };
3703
3704 GInt32 anStackDataOffset[2] = {
3705 poERDMS->GetIntField("layerStackDataOffset[0]"),
3706 poERDMS->GetIntField("layerStackDataOffset[1]")
3707 };
3708
3709 const GInt32 nStackCount = poERDMS->GetIntField("layerStackCount");
3710 const GInt32 nStackIndex = poERDMS->GetIntField("layerStackIndex");
3711
3712 // Update the filename.
3713 if( strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0 )
3714 {
3715 CPLString osNew = pszNewBase;
3716 osNew += osFileName.c_str() + strlen(pszOldBase);
3717 osFileName = osNew;
3718 }
3719
3720 // Grow the node if needed.
3721 if( strlen(pszNewBase) > strlen(pszOldBase) )
3722 {
3723 CPLDebug("HFA", "Growing ExternalRasterDMS to hold new names");
3724 poERDMS->MakeData(
3725 static_cast<int>(poERDMS->GetDataSize() +
3726 (strlen(pszNewBase) - strlen(pszOldBase))));
3727 }
3728
3729 // Initialize the whole thing to zeros for a clean start.
3730 memset(poERDMS->GetData(), 0, poERDMS->GetDataSize());
3731
3732 // Write it all out again, this may change the size of the node.
3733 poERDMS->SetStringField("fileName.string", osFileName);
3734 poERDMS->SetIntField("layerStackValidFlagsOffset[0]",
3735 anValidFlagsOffset[0]);
3736 poERDMS->SetIntField("layerStackValidFlagsOffset[1]",
3737 anValidFlagsOffset[1]);
3738
3739 poERDMS->SetIntField("layerStackDataOffset[0]", anStackDataOffset[0]);
3740 poERDMS->SetIntField("layerStackDataOffset[1]", anStackDataOffset[1]);
3741
3742 poERDMS->SetIntField("layerStackCount", nStackCount);
3743 poERDMS->SetIntField("layerStackIndex", nStackIndex);
3744 }
3745
3746 // DependentFile.
3747 apoNodeList =
3748 hHFA->poRoot->FindChildren("DependentFile", "Eimg_DependentFile");
3749
3750 for( size_t iNode = 0; iNode < apoNodeList.size(); iNode++ )
3751 {
3752 CPLString osFileName =
3753 apoNodeList[iNode]->GetStringField("dependent.string");
3754
3755 // Grow the node if needed.
3756 if( strlen(pszNewBase) > strlen(pszOldBase) )
3757 {
3758 CPLDebug("HFA", "Growing DependentFile to hold new names");
3759 apoNodeList[iNode]->MakeData(
3760 static_cast<int>(apoNodeList[iNode]->GetDataSize() +
3761 (strlen(pszNewBase) - strlen(pszOldBase))));
3762 }
3763
3764 // Update the filename.
3765 if( strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0 )
3766 {
3767 CPLString osNew = pszNewBase;
3768 osNew += osFileName.c_str() + strlen(pszOldBase);
3769 osFileName = osNew;
3770 }
3771
3772 apoNodeList[iNode]->SetStringField("dependent.string", osFileName);
3773 }
3774
3775 return CE_None;
3776 }
3777