Line data Source code
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 : * SPDX-License-Identifier: MIT
13 : ******************************************************************************
14 : *
15 : * hfaopen.cpp
16 : *
17 : * Supporting routines for reading Erdas Imagine (.imf) Hierarchical
18 : * File Architecture files. This is intended to be a library independent
19 : * of the GDAL core, but dependent on the Common Portability Library.
20 : *
21 : */
22 :
23 : #include "cpl_port.h"
24 : #include "hfa_p.h"
25 :
26 : #include <cerrno>
27 : #include <climits>
28 : #include <cmath>
29 : #include <cstddef>
30 : #include <cstdio>
31 : #include <cstdlib>
32 : #include <cstring>
33 : #include <algorithm>
34 : #include <memory>
35 : #include <string>
36 : #include <vector>
37 :
38 : #include "cpl_conv.h"
39 : #include "cpl_error.h"
40 : #include "cpl_string.h"
41 : #include "cpl_vsi.h"
42 : #include "gdal_priv.h"
43 : #include "hfa.h"
44 : #include "ogr_proj_p.h"
45 : #include "proj.h"
46 :
47 : constexpr double R2D = 180.0 / M_PI;
48 :
49 : constexpr double RAD2ARCSEC = 648000.0 / M_PI;
50 :
51 : static const char *const apszAuxMetadataItems[] = {
52 : // node/entry field_name metadata_key type
53 : "Statistics",
54 : "dminimum",
55 : "STATISTICS_MINIMUM",
56 : "Esta_Statistics",
57 : "Statistics",
58 : "dmaximum",
59 : "STATISTICS_MAXIMUM",
60 : "Esta_Statistics",
61 : "Statistics",
62 : "dmean",
63 : "STATISTICS_MEAN",
64 : "Esta_Statistics",
65 : "Statistics",
66 : "dmedian",
67 : "STATISTICS_MEDIAN",
68 : "Esta_Statistics",
69 : "Statistics",
70 : "dmode",
71 : "STATISTICS_MODE",
72 : "Esta_Statistics",
73 : "Statistics",
74 : "dstddev",
75 : "STATISTICS_STDDEV",
76 : "Esta_Statistics",
77 : "HistogramParameters",
78 : "lBinFunction.numBins",
79 : "STATISTICS_HISTONUMBINS",
80 : "Eimg_StatisticsParameters830",
81 : "HistogramParameters",
82 : "dBinFunction.minLimit",
83 : "STATISTICS_HISTOMIN",
84 : "Eimg_StatisticsParameters830",
85 : "HistogramParameters",
86 : "dBinFunction.maxLimit",
87 : "STATISTICS_HISTOMAX",
88 : "Eimg_StatisticsParameters830",
89 : "StatisticsParameters",
90 : "lSkipFactorX",
91 : "STATISTICS_SKIPFACTORX",
92 : "",
93 : "StatisticsParameters",
94 : "lSkipFactorY",
95 : "STATISTICS_SKIPFACTORY",
96 : "",
97 : "StatisticsParameters",
98 : "dExcludedValues",
99 : "STATISTICS_EXCLUDEDVALUES",
100 : "",
101 : "",
102 : "elayerType",
103 : "LAYER_TYPE",
104 : "",
105 : "RRDInfoList",
106 : "salgorithm.string",
107 : "OVERVIEWS_ALGORITHM",
108 : "Emif_String",
109 : nullptr};
110 :
111 688 : const char *const *GetHFAAuxMetaDataList()
112 : {
113 688 : return apszAuxMetadataItems;
114 : }
115 :
116 : /************************************************************************/
117 : /* HFAGetDictionary() */
118 : /************************************************************************/
119 :
120 562 : static char *HFAGetDictionary(HFAHandle hHFA)
121 :
122 : {
123 562 : int nDictMax = 100;
124 562 : char *pszDictionary = static_cast<char *>(CPLMalloc(nDictMax));
125 562 : int nDictSize = 0;
126 :
127 562 : if (VSIFSeekL(hHFA->fp, static_cast<vsi_l_offset>(hHFA->nDictionaryPos),
128 562 : SEEK_SET) < 0)
129 : {
130 0 : pszDictionary[nDictSize] = '\0';
131 0 : return pszDictionary;
132 : }
133 :
134 : while (true)
135 : {
136 2075220 : if (nDictSize >= nDictMax - 1)
137 : {
138 2791 : nDictMax = nDictSize * 2 + 100;
139 : pszDictionary =
140 2791 : static_cast<char *>(CPLRealloc(pszDictionary, nDictMax));
141 : }
142 :
143 2075220 : if (VSIFReadL(pszDictionary + nDictSize, 1, 1, hHFA->fp) < 1 ||
144 4148190 : pszDictionary[nDictSize] == '\0' ||
145 2072970 : (nDictSize > 2 && pszDictionary[nDictSize - 2] == ',' &&
146 161989 : pszDictionary[nDictSize - 1] == '.'))
147 562 : break;
148 :
149 2074650 : nDictSize++;
150 : }
151 :
152 562 : pszDictionary[nDictSize] = '\0';
153 :
154 562 : return pszDictionary;
155 : }
156 :
157 : /************************************************************************/
158 : /* HFAOpen() */
159 : /************************************************************************/
160 :
161 562 : HFAHandle HFAOpen(const char *pszFilename, const char *pszAccess)
162 :
163 : {
164 562 : VSILFILE *fp = VSIFOpenL(
165 : pszFilename,
166 562 : (EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb")) ? "rb" : "r+b");
167 :
168 : // Should this be changed to use some sort of CPLFOpen() which will
169 : // set the error?
170 562 : if (fp == nullptr)
171 : {
172 0 : CPLError(CE_Failure, CPLE_OpenFailed, "File open of %s failed.",
173 : pszFilename);
174 :
175 0 : return nullptr;
176 : }
177 :
178 : // Read and verify the header.
179 562 : char szHeader[16] = {};
180 562 : if (VSIFReadL(szHeader, 16, 1, fp) < 1)
181 : {
182 0 : CPLError(CE_Failure, CPLE_AppDefined,
183 : "Attempt to read 16 byte header failed for\n%s.", pszFilename);
184 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
185 0 : return nullptr;
186 : }
187 :
188 562 : if (!STARTS_WITH_CI(szHeader, "EHFA_HEADER_TAG"))
189 : {
190 0 : CPLError(CE_Failure, CPLE_AppDefined,
191 : "File %s is not an Imagine HFA file ... header wrong.",
192 : pszFilename);
193 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
194 0 : return nullptr;
195 : }
196 :
197 : // Create the HFAInfo_t.
198 : HFAInfo_t *psInfo =
199 562 : static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
200 :
201 562 : psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
202 562 : psInfo->pszPath = CPLStrdup(CPLGetPathSafe(pszFilename).c_str());
203 562 : psInfo->fp = fp;
204 562 : if (EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb"))
205 343 : psInfo->eAccess = HFA_ReadOnly;
206 : else
207 219 : psInfo->eAccess = HFA_Update;
208 562 : psInfo->bTreeDirty = false;
209 :
210 : // Where is the header?
211 562 : GUInt32 nHeaderPos = 0;
212 562 : bool bRet = VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, fp) > 0;
213 : HFAStandard(4, &nHeaderPos);
214 :
215 : // Read the header.
216 562 : bRet &= VSIFSeekL(fp, static_cast<vsi_l_offset>(nHeaderPos), SEEK_SET) >= 0;
217 :
218 562 : bRet &= VSIFReadL(&(psInfo->nVersion), sizeof(GInt32), 1, fp) > 0;
219 : HFAStandard(4, &(psInfo->nVersion));
220 :
221 562 : bRet &= VSIFReadL(szHeader, 4, 1, fp) > 0; // Skip freeList.
222 :
223 562 : bRet &= VSIFReadL(&(psInfo->nRootPos), sizeof(GInt32), 1, fp) > 0;
224 : HFAStandard(4, &(psInfo->nRootPos));
225 :
226 562 : bRet &= VSIFReadL(&(psInfo->nEntryHeaderLength), sizeof(GInt16), 1, fp) > 0;
227 : HFAStandard(2, &(psInfo->nEntryHeaderLength));
228 :
229 562 : bRet &= VSIFReadL(&(psInfo->nDictionaryPos), sizeof(GInt32), 1, fp) > 0;
230 : HFAStandard(4, &(psInfo->nDictionaryPos));
231 :
232 : // Collect file size.
233 562 : bRet &= VSIFSeekL(fp, 0, SEEK_END) >= 0;
234 562 : if (!bRet)
235 : {
236 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
237 0 : CPLFree(psInfo->pszFilename);
238 0 : CPLFree(psInfo->pszPath);
239 0 : CPLFree(psInfo);
240 0 : return nullptr;
241 : }
242 562 : psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
243 :
244 : // Instantiate the root entry.
245 562 : psInfo->poRoot = HFAEntry::New(psInfo, psInfo->nRootPos, nullptr, nullptr);
246 562 : if (psInfo->poRoot == nullptr)
247 : {
248 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
249 0 : CPLFree(psInfo->pszFilename);
250 0 : CPLFree(psInfo->pszPath);
251 0 : CPLFree(psInfo);
252 0 : return nullptr;
253 : }
254 :
255 : // Read the dictionary.
256 562 : psInfo->pszDictionary = HFAGetDictionary(psInfo);
257 562 : psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
258 :
259 : // Collect band definitions.
260 562 : HFAParseBandInfo(psInfo);
261 :
262 562 : return psInfo;
263 : }
264 :
265 : /************************************************************************/
266 : /* HFACreateDependent() */
267 : /* */
268 : /* Create a .rrd file for the named file if it does not exist, */
269 : /* or return the existing dependent if it already exists. */
270 : /************************************************************************/
271 :
272 2 : HFAInfo_t *HFACreateDependent(HFAInfo_t *psBase)
273 :
274 : {
275 2 : if (psBase->psDependent != nullptr)
276 1 : return psBase->psDependent;
277 :
278 : // Create desired RRD filename.
279 2 : const CPLString oBasename = CPLGetBasenameSafe(psBase->pszFilename);
280 : const CPLString oRRDFilename =
281 2 : CPLFormFilenameSafe(psBase->pszPath, oBasename, "rrd");
282 :
283 : // Does this file already exist? If so, re-use it.
284 1 : VSILFILE *fp = VSIFOpenL(oRRDFilename, "rb");
285 1 : if (fp != nullptr)
286 : {
287 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
288 0 : psBase->psDependent = HFAOpen(oRRDFilename, "rb");
289 : // FIXME? this is not going to be reused but recreated...
290 : }
291 :
292 : // Otherwise create it now.
293 1 : HFAInfo_t *psDep = HFACreateLL(oRRDFilename);
294 1 : psBase->psDependent = psDep;
295 1 : if (psDep == nullptr)
296 0 : return nullptr;
297 :
298 : /* -------------------------------------------------------------------- */
299 : /* Add the DependentFile node with the pointer back to the */
300 : /* parent. When working from an .aux file we really want the */
301 : /* .rrd to point back to the original file, not the .aux file. */
302 : /* -------------------------------------------------------------------- */
303 1 : HFAEntry *poEntry = psBase->poRoot->GetNamedChild("DependentFile");
304 1 : const char *pszDependentFile = nullptr;
305 1 : if (poEntry != nullptr)
306 0 : pszDependentFile = poEntry->GetStringField("dependent.string");
307 1 : if (pszDependentFile == nullptr)
308 1 : pszDependentFile = psBase->pszFilename;
309 :
310 1 : HFAEntry *poDF = HFAEntry::New(psDep, "DependentFile", "Eimg_DependentFile",
311 : psDep->poRoot);
312 :
313 1 : poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
314 1 : poDF->SetPosition();
315 1 : poDF->SetStringField("dependent.string", pszDependentFile);
316 :
317 1 : return psDep;
318 : }
319 :
320 : /************************************************************************/
321 : /* HFAGetDependent() */
322 : /************************************************************************/
323 :
324 50 : HFAInfo_t *HFAGetDependent(HFAInfo_t *psBase, const char *pszFilename)
325 :
326 : {
327 50 : if (EQUAL(pszFilename, psBase->pszFilename))
328 31 : return psBase;
329 :
330 19 : if (psBase->psDependent != nullptr)
331 : {
332 2 : if (EQUAL(pszFilename, psBase->psDependent->pszFilename))
333 2 : return psBase->psDependent;
334 : else
335 0 : return nullptr;
336 : }
337 :
338 : // Try to open the dependent file.
339 17 : const char *pszMode = psBase->eAccess == HFA_Update ? "r+b" : "rb";
340 :
341 17 : char *pszDependent = CPLStrdup(
342 34 : CPLFormFilenameSafe(psBase->pszPath, pszFilename, nullptr).c_str());
343 :
344 17 : VSILFILE *fp = VSIFOpenL(pszDependent, pszMode);
345 17 : if (fp != nullptr)
346 : {
347 13 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
348 13 : psBase->psDependent = HFAOpen(pszDependent, pszMode);
349 : }
350 :
351 17 : CPLFree(pszDependent);
352 :
353 17 : return psBase->psDependent;
354 : }
355 :
356 : /************************************************************************/
357 : /* HFAParseBandInfo() */
358 : /* */
359 : /* This is used by HFAOpen() and HFACreate() to initialize the */
360 : /* band structures. */
361 : /************************************************************************/
362 :
363 754 : CPLErr HFAParseBandInfo(HFAInfo_t *psInfo)
364 :
365 : {
366 : // Find the first band node.
367 754 : psInfo->nBands = 0;
368 754 : HFAEntry *poNode = psInfo->poRoot->GetChild();
369 2463 : while (poNode != nullptr)
370 : {
371 1709 : if (EQUAL(poNode->GetType(), "Eimg_Layer") &&
372 2563 : poNode->GetIntField("width") > 0 &&
373 854 : poNode->GetIntField("height") > 0)
374 : {
375 854 : if (psInfo->nBands == 0)
376 : {
377 736 : psInfo->nXSize = poNode->GetIntField("width");
378 736 : psInfo->nYSize = poNode->GetIntField("height");
379 : }
380 236 : else if (poNode->GetIntField("width") != psInfo->nXSize ||
381 118 : poNode->GetIntField("height") != psInfo->nYSize)
382 : {
383 0 : return CE_Failure;
384 : }
385 :
386 1708 : psInfo->papoBand = static_cast<HFABand **>(CPLRealloc(
387 854 : psInfo->papoBand, sizeof(HFABand *) * (psInfo->nBands + 1)));
388 854 : psInfo->papoBand[psInfo->nBands] = new HFABand(psInfo, poNode);
389 854 : if (psInfo->papoBand[psInfo->nBands]->nWidth == 0)
390 : {
391 0 : delete psInfo->papoBand[psInfo->nBands];
392 0 : return CE_Failure;
393 : }
394 854 : psInfo->nBands++;
395 : }
396 :
397 1709 : poNode = poNode->GetNext();
398 : }
399 :
400 754 : return CE_None;
401 : }
402 :
403 : /************************************************************************/
404 : /* HFAClose() */
405 : /************************************************************************/
406 :
407 763 : int HFAClose(HFAHandle hHFA)
408 :
409 : {
410 763 : if (hHFA->eAccess == HFA_Update &&
411 420 : (hHFA->bTreeDirty || (hHFA->poDictionary != nullptr &&
412 16 : hHFA->poDictionary->bDictionaryTextDirty)))
413 400 : HFAFlush(hHFA);
414 :
415 763 : int nRet = 0;
416 763 : if (hHFA->psDependent != nullptr)
417 : {
418 13 : if (HFAClose(hHFA->psDependent) != 0)
419 0 : nRet = -1;
420 : }
421 :
422 763 : delete hHFA->poRoot;
423 :
424 763 : if (VSIFCloseL(hHFA->fp) != 0)
425 0 : nRet = -1;
426 :
427 763 : if (hHFA->poDictionary != nullptr)
428 759 : delete hHFA->poDictionary;
429 :
430 763 : CPLFree(hHFA->pszDictionary);
431 763 : CPLFree(hHFA->pszFilename);
432 763 : CPLFree(hHFA->pszIGEFilename);
433 763 : CPLFree(hHFA->pszPath);
434 :
435 1617 : for (int i = 0; i < hHFA->nBands; i++)
436 : {
437 854 : delete hHFA->papoBand[i];
438 : }
439 :
440 763 : CPLFree(hHFA->papoBand);
441 :
442 763 : if (hHFA->pProParameters != nullptr)
443 : {
444 203 : Eprj_ProParameters *psProParams =
445 : (Eprj_ProParameters *)hHFA->pProParameters;
446 :
447 203 : CPLFree(psProParams->proExeName);
448 203 : CPLFree(psProParams->proName);
449 203 : CPLFree(psProParams->proSpheroid.sphereName);
450 :
451 203 : CPLFree(psProParams);
452 : }
453 :
454 763 : if (hHFA->pDatum != nullptr)
455 : {
456 203 : CPLFree(((Eprj_Datum *)hHFA->pDatum)->datumname);
457 203 : CPLFree(((Eprj_Datum *)hHFA->pDatum)->gridname);
458 203 : CPLFree(hHFA->pDatum);
459 : }
460 :
461 763 : if (hHFA->pMapInfo != nullptr)
462 : {
463 246 : CPLFree(((Eprj_MapInfo *)hHFA->pMapInfo)->proName);
464 246 : CPLFree(((Eprj_MapInfo *)hHFA->pMapInfo)->units);
465 246 : CPLFree(hHFA->pMapInfo);
466 : }
467 :
468 763 : CPLFree(hHFA);
469 763 : return nRet;
470 : }
471 :
472 : /************************************************************************/
473 : /* HFARemove() */
474 : /* Used from HFADelete() function. */
475 : /************************************************************************/
476 :
477 0 : static CPLErr HFARemove(const char *pszFilename)
478 :
479 : {
480 : VSIStatBufL sStat;
481 :
482 0 : if (VSIStatL(pszFilename, &sStat) == 0 && VSI_ISREG(sStat.st_mode))
483 : {
484 0 : if (VSIUnlink(pszFilename) == 0)
485 0 : return CE_None;
486 : else
487 : {
488 0 : CPLError(CE_Failure, CPLE_AppDefined,
489 : "Attempt to unlink %s failed.", pszFilename);
490 0 : return CE_Failure;
491 : }
492 : }
493 :
494 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unable to delete %s, not a file.",
495 : pszFilename);
496 0 : return CE_Failure;
497 : }
498 :
499 : /************************************************************************/
500 : /* HFADelete() */
501 : /************************************************************************/
502 :
503 0 : CPLErr HFADelete(const char *pszFilename)
504 :
505 : {
506 0 : HFAInfo_t *psInfo = HFAOpen(pszFilename, "rb");
507 0 : HFAEntry *poDMS = nullptr;
508 0 : HFAEntry *poLayer = nullptr;
509 0 : HFAEntry *poNode = nullptr;
510 :
511 0 : if (psInfo != nullptr)
512 : {
513 0 : poNode = psInfo->poRoot->GetChild();
514 0 : while ((poNode != nullptr) && (poLayer == nullptr))
515 : {
516 0 : if (EQUAL(poNode->GetType(), "Eimg_Layer"))
517 : {
518 0 : poLayer = poNode;
519 : }
520 0 : poNode = poNode->GetNext();
521 : }
522 :
523 0 : if (poLayer != nullptr)
524 0 : poDMS = poLayer->GetNamedChild("ExternalRasterDMS");
525 :
526 0 : if (poDMS)
527 : {
528 : const char *pszRawFilename =
529 0 : poDMS->GetStringField("fileName.string");
530 :
531 0 : if (pszRawFilename != nullptr)
532 : {
533 0 : if (CPLHasPathTraversal(pszRawFilename))
534 : {
535 0 : CPLError(CE_Warning, CPLE_AppDefined,
536 : "Path traversal detected in %s", pszRawFilename);
537 : }
538 : else
539 : {
540 0 : HFARemove(CPLFormFilenameSafe(psInfo->pszPath,
541 : pszRawFilename, nullptr)
542 : .c_str());
543 : }
544 : }
545 : }
546 :
547 0 : CPL_IGNORE_RET_VAL(HFAClose(psInfo));
548 : }
549 0 : return HFARemove(pszFilename);
550 : }
551 :
552 : /************************************************************************/
553 : /* HFAGetRasterInfo() */
554 : /************************************************************************/
555 :
556 547 : CPLErr HFAGetRasterInfo(HFAHandle hHFA, int *pnXSize, int *pnYSize,
557 : int *pnBands)
558 :
559 : {
560 547 : if (pnXSize != nullptr)
561 547 : *pnXSize = hHFA->nXSize;
562 547 : if (pnYSize != nullptr)
563 547 : *pnYSize = hHFA->nYSize;
564 547 : if (pnBands != nullptr)
565 547 : *pnBands = hHFA->nBands;
566 547 : return CE_None;
567 : }
568 :
569 : /************************************************************************/
570 : /* HFAGetBandInfo() */
571 : /************************************************************************/
572 :
573 679 : CPLErr HFAGetBandInfo(HFAHandle hHFA, int nBand, EPTType *peDataType,
574 : int *pnBlockXSize, int *pnBlockYSize,
575 : int *pnCompressionType)
576 :
577 : {
578 679 : if (nBand < 0 || nBand > hHFA->nBands)
579 : {
580 0 : CPLAssert(false);
581 : return CE_Failure;
582 : }
583 :
584 679 : HFABand *poBand = hHFA->papoBand[nBand - 1];
585 :
586 679 : if (peDataType != nullptr)
587 679 : *peDataType = poBand->eDataType;
588 :
589 679 : if (pnBlockXSize != nullptr)
590 679 : *pnBlockXSize = poBand->nBlockXSize;
591 :
592 679 : if (pnBlockYSize != nullptr)
593 679 : *pnBlockYSize = poBand->nBlockYSize;
594 :
595 : // Get compression code from RasterDMS.
596 679 : if (pnCompressionType != nullptr)
597 : {
598 679 : *pnCompressionType = 0;
599 :
600 679 : HFAEntry *poDMS = poBand->poNode->GetNamedChild("RasterDMS");
601 :
602 679 : if (poDMS != nullptr)
603 585 : *pnCompressionType = poDMS->GetIntField("compressionType");
604 : }
605 :
606 679 : return CE_None;
607 : }
608 :
609 : /************************************************************************/
610 : /* HFAGetBandNoData() */
611 : /* */
612 : /* returns TRUE if value is set, otherwise FALSE. */
613 : /************************************************************************/
614 :
615 223 : int HFAGetBandNoData(HFAHandle hHFA, int nBand, double *pdfNoData)
616 :
617 : {
618 223 : if (nBand < 0 || nBand > hHFA->nBands)
619 : {
620 0 : CPLAssert(false);
621 : return CE_Failure;
622 : }
623 :
624 223 : HFABand *poBand = hHFA->papoBand[nBand - 1];
625 :
626 223 : if (!poBand->bNoDataSet && poBand->nOverviews > 0)
627 : {
628 13 : poBand = poBand->papoOverviews[0];
629 13 : if (poBand == nullptr)
630 0 : return FALSE;
631 : }
632 :
633 223 : *pdfNoData = poBand->dfNoData;
634 223 : return poBand->bNoDataSet;
635 : }
636 :
637 : /************************************************************************/
638 : /* HFASetBandNoData() */
639 : /* */
640 : /* attempts to set a no-data value on the given band */
641 : /************************************************************************/
642 :
643 9 : CPLErr HFASetBandNoData(HFAHandle hHFA, int nBand, double dfValue)
644 :
645 : {
646 9 : if (nBand < 0 || nBand > hHFA->nBands)
647 : {
648 0 : CPLAssert(false);
649 : return CE_Failure;
650 : }
651 :
652 9 : HFABand *poBand = hHFA->papoBand[nBand - 1];
653 :
654 9 : return poBand->SetNoDataValue(dfValue);
655 : }
656 :
657 : /************************************************************************/
658 : /* HFAGetOverviewCount() */
659 : /************************************************************************/
660 :
661 127 : int HFAGetOverviewCount(HFAHandle hHFA, int nBand)
662 :
663 : {
664 127 : if (nBand < 0 || nBand > hHFA->nBands)
665 : {
666 0 : CPLAssert(false);
667 : return CE_Failure;
668 : }
669 :
670 127 : HFABand *poBand = hHFA->papoBand[nBand - 1];
671 127 : poBand->LoadOverviews();
672 :
673 127 : return poBand->nOverviews;
674 : }
675 :
676 : /************************************************************************/
677 : /* HFAGetOverviewInfo() */
678 : /************************************************************************/
679 :
680 56 : CPLErr HFAGetOverviewInfo(HFAHandle hHFA, int nBand, int iOverview,
681 : int *pnXSize, int *pnYSize, int *pnBlockXSize,
682 : int *pnBlockYSize, EPTType *peHFADataType)
683 :
684 : {
685 56 : if (nBand < 0 || nBand > hHFA->nBands)
686 : {
687 0 : CPLAssert(false);
688 : return CE_Failure;
689 : }
690 :
691 56 : HFABand *poBand = hHFA->papoBand[nBand - 1];
692 56 : poBand->LoadOverviews();
693 :
694 56 : if (iOverview < 0 || iOverview >= poBand->nOverviews)
695 : {
696 0 : CPLAssert(false);
697 : return CE_Failure;
698 : }
699 56 : poBand = poBand->papoOverviews[iOverview];
700 56 : if (poBand == nullptr)
701 : {
702 0 : return CE_Failure;
703 : }
704 :
705 56 : if (pnXSize != nullptr)
706 56 : *pnXSize = poBand->nWidth;
707 :
708 56 : if (pnYSize != nullptr)
709 56 : *pnYSize = poBand->nHeight;
710 :
711 56 : if (pnBlockXSize != nullptr)
712 56 : *pnBlockXSize = poBand->nBlockXSize;
713 :
714 56 : if (pnBlockYSize != nullptr)
715 56 : *pnBlockYSize = poBand->nBlockYSize;
716 :
717 56 : if (peHFADataType != nullptr)
718 56 : *peHFADataType = poBand->eDataType;
719 :
720 56 : return CE_None;
721 : }
722 :
723 : /************************************************************************/
724 : /* HFAGetRasterBlock() */
725 : /************************************************************************/
726 :
727 0 : CPLErr HFAGetRasterBlock(HFAHandle hHFA, int nBand, int nXBlock, int nYBlock,
728 : void *pData)
729 :
730 : {
731 0 : return HFAGetRasterBlockEx(hHFA, nBand, nXBlock, nYBlock, pData, -1);
732 : }
733 :
734 : /************************************************************************/
735 : /* HFAGetRasterBlockEx() */
736 : /************************************************************************/
737 :
738 1330 : CPLErr HFAGetRasterBlockEx(HFAHandle hHFA, int nBand, int nXBlock, int nYBlock,
739 : void *pData, int nDataSize)
740 :
741 : {
742 1330 : if (nBand < 1 || nBand > hHFA->nBands)
743 0 : return CE_Failure;
744 :
745 1330 : return hHFA->papoBand[nBand - 1]->GetRasterBlock(nXBlock, nYBlock, pData,
746 1330 : nDataSize);
747 : }
748 :
749 : /************************************************************************/
750 : /* HFAGetOverviewRasterBlock() */
751 : /************************************************************************/
752 :
753 0 : CPLErr HFAGetOverviewRasterBlock(HFAHandle hHFA, int nBand, int iOverview,
754 : int nXBlock, int nYBlock, void *pData)
755 :
756 : {
757 0 : return HFAGetOverviewRasterBlockEx(hHFA, nBand, iOverview, nXBlock, nYBlock,
758 0 : pData, -1);
759 : }
760 :
761 : /************************************************************************/
762 : /* HFAGetOverviewRasterBlockEx() */
763 : /************************************************************************/
764 :
765 24 : CPLErr HFAGetOverviewRasterBlockEx(HFAHandle hHFA, int nBand, int iOverview,
766 : int nXBlock, int nYBlock, void *pData,
767 : int nDataSize)
768 :
769 : {
770 24 : if (nBand < 1 || nBand > hHFA->nBands)
771 0 : return CE_Failure;
772 :
773 24 : if (iOverview < 0 || iOverview >= hHFA->papoBand[nBand - 1]->nOverviews)
774 0 : return CE_Failure;
775 :
776 24 : return hHFA->papoBand[nBand - 1]->papoOverviews[iOverview]->GetRasterBlock(
777 24 : nXBlock, nYBlock, pData, nDataSize);
778 : }
779 :
780 : /************************************************************************/
781 : /* HFASetRasterBlock() */
782 : /************************************************************************/
783 :
784 105 : CPLErr HFASetRasterBlock(HFAHandle hHFA, int nBand, int nXBlock, int nYBlock,
785 : void *pData)
786 :
787 : {
788 105 : if (nBand < 1 || nBand > hHFA->nBands)
789 0 : return CE_Failure;
790 :
791 105 : return hHFA->papoBand[nBand - 1]->SetRasterBlock(nXBlock, nYBlock, pData);
792 : }
793 :
794 : /************************************************************************/
795 : /* HFASetRasterBlock() */
796 : /************************************************************************/
797 :
798 29 : CPLErr HFASetOverviewRasterBlock(HFAHandle hHFA, int nBand, int iOverview,
799 : int nXBlock, int nYBlock, void *pData)
800 :
801 : {
802 29 : if (nBand < 1 || nBand > hHFA->nBands)
803 0 : return CE_Failure;
804 :
805 29 : if (iOverview < 0 || iOverview >= hHFA->papoBand[nBand - 1]->nOverviews)
806 0 : return CE_Failure;
807 :
808 29 : return hHFA->papoBand[nBand - 1]->papoOverviews[iOverview]->SetRasterBlock(
809 29 : nXBlock, nYBlock, pData);
810 : }
811 :
812 : /************************************************************************/
813 : /* HFAGetBandName() */
814 : /************************************************************************/
815 :
816 214 : const char *HFAGetBandName(HFAHandle hHFA, int nBand)
817 : {
818 214 : if (nBand < 1 || nBand > hHFA->nBands)
819 0 : return "";
820 :
821 214 : return hHFA->papoBand[nBand - 1]->GetBandName();
822 : }
823 :
824 : /************************************************************************/
825 : /* HFASetBandName() */
826 : /************************************************************************/
827 :
828 7 : void HFASetBandName(HFAHandle hHFA, int nBand, const char *pszName)
829 : {
830 7 : if (nBand < 1 || nBand > hHFA->nBands)
831 0 : return;
832 :
833 7 : hHFA->papoBand[nBand - 1]->SetBandName(pszName);
834 : }
835 :
836 : /************************************************************************/
837 : /* HFAGetDataTypeBits() */
838 : /************************************************************************/
839 :
840 4730 : int HFAGetDataTypeBits(EPTType eDataType)
841 :
842 : {
843 4730 : switch (eDataType)
844 : {
845 76 : case EPT_u1:
846 76 : return 1;
847 :
848 42 : case EPT_u2:
849 42 : return 2;
850 :
851 18 : case EPT_u4:
852 18 : return 4;
853 :
854 1345 : case EPT_u8:
855 : case EPT_s8:
856 1345 : return 8;
857 :
858 1149 : case EPT_u16:
859 : case EPT_s16:
860 1149 : return 16;
861 :
862 268 : case EPT_u32:
863 : case EPT_s32:
864 : case EPT_f32:
865 268 : return 32;
866 :
867 1793 : case EPT_f64:
868 : case EPT_c64:
869 1793 : return 64;
870 :
871 39 : case EPT_c128:
872 39 : return 128;
873 : }
874 :
875 0 : CPLAssert(false);
876 : return 1;
877 : }
878 :
879 : /************************************************************************/
880 : /* HFAGetDataTypeName() */
881 : /************************************************************************/
882 :
883 0 : const char *HFAGetDataTypeName(EPTType eDataType)
884 :
885 : {
886 0 : switch (eDataType)
887 : {
888 0 : case EPT_u1:
889 0 : return "u1";
890 :
891 0 : case EPT_u2:
892 0 : return "u2";
893 :
894 0 : case EPT_u4:
895 0 : return "u4";
896 :
897 0 : case EPT_u8:
898 0 : return "u8";
899 :
900 0 : case EPT_s8:
901 0 : return "s8";
902 :
903 0 : case EPT_u16:
904 0 : return "u16";
905 :
906 0 : case EPT_s16:
907 0 : return "s16";
908 :
909 0 : case EPT_u32:
910 0 : return "u32";
911 :
912 0 : case EPT_s32:
913 0 : return "s32";
914 :
915 0 : case EPT_f32:
916 0 : return "f32";
917 :
918 0 : case EPT_f64:
919 0 : return "f64";
920 :
921 0 : case EPT_c64:
922 0 : return "c64";
923 :
924 0 : case EPT_c128:
925 0 : return "c128";
926 :
927 0 : default:
928 0 : CPLAssert(false);
929 : return "unknown";
930 : }
931 : }
932 :
933 : /************************************************************************/
934 : /* HFAGetMapInfo() */
935 : /************************************************************************/
936 :
937 1086 : const Eprj_MapInfo *HFAGetMapInfo(HFAHandle hHFA)
938 :
939 : {
940 1086 : if (hHFA->nBands < 1)
941 0 : return nullptr;
942 :
943 : // Do we already have it?
944 1086 : if (hHFA->pMapInfo != nullptr)
945 246 : return (Eprj_MapInfo *)hHFA->pMapInfo;
946 :
947 : // Get the HFA node. If we don't find it under the usual name
948 : // we search for any node of the right type (#3338).
949 840 : HFAEntry *poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild("Map_Info");
950 840 : if (poMIEntry == nullptr)
951 : {
952 594 : for (HFAEntry *poChild = hHFA->papoBand[0]->poNode->GetChild();
953 2060 : poChild != nullptr && poMIEntry == nullptr;
954 1466 : poChild = poChild->GetNext())
955 : {
956 1466 : if (EQUAL(poChild->GetType(), "Eprj_MapInfo"))
957 0 : poMIEntry = poChild;
958 : }
959 : }
960 :
961 840 : if (poMIEntry == nullptr)
962 : {
963 594 : return nullptr;
964 : }
965 :
966 : // Allocate the structure.
967 : Eprj_MapInfo *psMapInfo =
968 246 : static_cast<Eprj_MapInfo *>(CPLCalloc(sizeof(Eprj_MapInfo), 1));
969 :
970 : // Fetch the fields.
971 246 : psMapInfo->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
972 :
973 246 : psMapInfo->upperLeftCenter.x =
974 246 : poMIEntry->GetDoubleField("upperLeftCenter.x");
975 246 : psMapInfo->upperLeftCenter.y =
976 246 : poMIEntry->GetDoubleField("upperLeftCenter.y");
977 :
978 246 : psMapInfo->lowerRightCenter.x =
979 246 : poMIEntry->GetDoubleField("lowerRightCenter.x");
980 246 : psMapInfo->lowerRightCenter.y =
981 246 : poMIEntry->GetDoubleField("lowerRightCenter.y");
982 :
983 246 : CPLErr eErr = CE_None;
984 246 : psMapInfo->pixelSize.width =
985 246 : poMIEntry->GetDoubleField("pixelSize.width", &eErr);
986 246 : psMapInfo->pixelSize.height =
987 246 : poMIEntry->GetDoubleField("pixelSize.height", &eErr);
988 :
989 : // The following is basically a hack to get files with
990 : // non-standard MapInfo's that misname the pixelSize fields. (#3338)
991 246 : if (eErr != CE_None)
992 : {
993 0 : psMapInfo->pixelSize.width = poMIEntry->GetDoubleField("pixelSize.x");
994 0 : psMapInfo->pixelSize.height = poMIEntry->GetDoubleField("pixelSize.y");
995 : }
996 :
997 246 : psMapInfo->units = CPLStrdup(poMIEntry->GetStringField("units"));
998 :
999 246 : hHFA->pMapInfo = (void *)psMapInfo;
1000 :
1001 246 : return psMapInfo;
1002 : }
1003 :
1004 : /************************************************************************/
1005 : /* HFAInvGeoTransform() */
1006 : /************************************************************************/
1007 :
1008 6 : static bool HFAInvGeoTransform(const double *gt_in, double *gt_out)
1009 :
1010 : {
1011 : // Assume a 3rd row that is [1 0 0].
1012 : // Compute determinate.
1013 6 : const double det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
1014 :
1015 6 : if (fabs(det) < 1.0e-15)
1016 0 : return false;
1017 :
1018 6 : const double inv_det = 1.0 / det;
1019 :
1020 : // Compute adjoint, and divide by determinate.
1021 6 : gt_out[1] = gt_in[5] * inv_det;
1022 6 : gt_out[4] = -gt_in[4] * inv_det;
1023 :
1024 6 : gt_out[2] = -gt_in[2] * inv_det;
1025 6 : gt_out[5] = gt_in[1] * inv_det;
1026 :
1027 6 : gt_out[0] = (gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
1028 6 : gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
1029 :
1030 6 : return true;
1031 : }
1032 :
1033 : /************************************************************************/
1034 : /* HFAGetGeoTransform() */
1035 : /************************************************************************/
1036 :
1037 543 : int HFAGetGeoTransform(HFAHandle hHFA, double *padfGeoTransform)
1038 :
1039 : {
1040 543 : const Eprj_MapInfo *psMapInfo = HFAGetMapInfo(hHFA);
1041 :
1042 543 : padfGeoTransform[0] = 0.0;
1043 543 : padfGeoTransform[1] = 1.0;
1044 543 : padfGeoTransform[2] = 0.0;
1045 543 : padfGeoTransform[3] = 0.0;
1046 543 : padfGeoTransform[4] = 0.0;
1047 543 : padfGeoTransform[5] = 1.0;
1048 :
1049 : // Simple (north up) MapInfo approach.
1050 543 : if (psMapInfo != nullptr)
1051 : {
1052 246 : padfGeoTransform[0] =
1053 246 : psMapInfo->upperLeftCenter.x - psMapInfo->pixelSize.width * 0.5;
1054 246 : padfGeoTransform[1] = psMapInfo->pixelSize.width;
1055 246 : if (padfGeoTransform[1] == 0.0)
1056 0 : padfGeoTransform[1] = 1.0;
1057 246 : padfGeoTransform[2] = 0.0;
1058 246 : if (psMapInfo->upperLeftCenter.y >= psMapInfo->lowerRightCenter.y)
1059 246 : padfGeoTransform[5] = -psMapInfo->pixelSize.height;
1060 : else
1061 0 : padfGeoTransform[5] = psMapInfo->pixelSize.height;
1062 246 : if (padfGeoTransform[5] == 0.0)
1063 0 : padfGeoTransform[5] = 1.0;
1064 :
1065 246 : padfGeoTransform[3] =
1066 246 : psMapInfo->upperLeftCenter.y - padfGeoTransform[5] * 0.5;
1067 246 : padfGeoTransform[4] = 0.0;
1068 :
1069 : // Special logic to fixup odd angular units.
1070 246 : if (EQUAL(psMapInfo->units, "ds"))
1071 : {
1072 0 : padfGeoTransform[0] /= 3600.0;
1073 0 : padfGeoTransform[1] /= 3600.0;
1074 0 : padfGeoTransform[2] /= 3600.0;
1075 0 : padfGeoTransform[3] /= 3600.0;
1076 0 : padfGeoTransform[4] /= 3600.0;
1077 0 : padfGeoTransform[5] /= 3600.0;
1078 : }
1079 :
1080 246 : return TRUE;
1081 : }
1082 :
1083 : // Try for a MapToPixelXForm affine polynomial supporting
1084 : // rotated and sheared affine transformations.
1085 297 : if (hHFA->nBands == 0)
1086 0 : return FALSE;
1087 :
1088 : HFAEntry *poXForm0 =
1089 297 : hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
1090 :
1091 297 : if (poXForm0 == nullptr)
1092 291 : return FALSE;
1093 :
1094 6 : if (poXForm0->GetIntField("order") != 1 ||
1095 5 : poXForm0->GetIntField("numdimtransform") != 2 ||
1096 16 : poXForm0->GetIntField("numdimpolynomial") != 2 ||
1097 5 : poXForm0->GetIntField("termcount") != 3)
1098 1 : return FALSE;
1099 :
1100 : // Verify that there aren't any further xform steps.
1101 5 : if (hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm1") !=
1102 : nullptr)
1103 1 : return FALSE;
1104 :
1105 : // We should check that the exponent list is 0 0 1 0 0 1, but
1106 : // we don't because we are lazy.
1107 :
1108 : // Fetch geotransform values.
1109 4 : double adfXForm[6] = {poXForm0->GetDoubleField("polycoefvector[0]"),
1110 8 : poXForm0->GetDoubleField("polycoefmtx[0]"),
1111 8 : poXForm0->GetDoubleField("polycoefmtx[2]"),
1112 8 : poXForm0->GetDoubleField("polycoefvector[1]"),
1113 8 : poXForm0->GetDoubleField("polycoefmtx[1]"),
1114 4 : poXForm0->GetDoubleField("polycoefmtx[3]")};
1115 :
1116 : // Invert.
1117 :
1118 4 : if (!HFAInvGeoTransform(adfXForm, padfGeoTransform))
1119 0 : memset(padfGeoTransform, 0, 6 * sizeof(double));
1120 :
1121 : // Adjust origin from center of top left pixel to top left corner
1122 : // of top left pixel.
1123 4 : padfGeoTransform[0] -= padfGeoTransform[1] * 0.5;
1124 4 : padfGeoTransform[0] -= padfGeoTransform[2] * 0.5;
1125 4 : padfGeoTransform[3] -= padfGeoTransform[4] * 0.5;
1126 4 : padfGeoTransform[3] -= padfGeoTransform[5] * 0.5;
1127 :
1128 4 : return TRUE;
1129 : }
1130 :
1131 : /************************************************************************/
1132 : /* HFASetMapInfo() */
1133 : /************************************************************************/
1134 :
1135 145 : CPLErr HFASetMapInfo(HFAHandle hHFA, const Eprj_MapInfo *poMapInfo)
1136 :
1137 : {
1138 : // Loop over bands, setting information on each one.
1139 328 : for (int iBand = 0; iBand < hHFA->nBands; iBand++)
1140 : {
1141 : // Create a new Map_Info if there isn't one present already.
1142 : HFAEntry *poMIEntry =
1143 183 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Map_Info");
1144 183 : if (poMIEntry == nullptr)
1145 : {
1146 182 : poMIEntry = HFAEntry::New(hHFA, "Map_Info", "Eprj_MapInfo",
1147 182 : hHFA->papoBand[iBand]->poNode);
1148 : }
1149 :
1150 183 : poMIEntry->MarkDirty();
1151 :
1152 : // Ensure we have enough space for all the data.
1153 : // TODO(schwehr): Explain 48 and 40 constants.
1154 183 : const int nSize =
1155 183 : static_cast<int>(48 + 40 + strlen(poMapInfo->proName) + 1 +
1156 183 : strlen(poMapInfo->units) + 1);
1157 :
1158 183 : GByte *pabyData = poMIEntry->MakeData(nSize);
1159 183 : memset(pabyData, 0, nSize);
1160 :
1161 183 : poMIEntry->SetPosition();
1162 :
1163 : // Write the various fields.
1164 183 : poMIEntry->SetStringField("proName", poMapInfo->proName);
1165 :
1166 183 : poMIEntry->SetDoubleField("upperLeftCenter.x",
1167 183 : poMapInfo->upperLeftCenter.x);
1168 183 : poMIEntry->SetDoubleField("upperLeftCenter.y",
1169 183 : poMapInfo->upperLeftCenter.y);
1170 :
1171 183 : poMIEntry->SetDoubleField("lowerRightCenter.x",
1172 183 : poMapInfo->lowerRightCenter.x);
1173 183 : poMIEntry->SetDoubleField("lowerRightCenter.y",
1174 183 : poMapInfo->lowerRightCenter.y);
1175 :
1176 183 : poMIEntry->SetDoubleField("pixelSize.width",
1177 183 : poMapInfo->pixelSize.width);
1178 183 : poMIEntry->SetDoubleField("pixelSize.height",
1179 183 : poMapInfo->pixelSize.height);
1180 :
1181 183 : poMIEntry->SetStringField("units", poMapInfo->units);
1182 : }
1183 :
1184 145 : return CE_None;
1185 : }
1186 :
1187 : /************************************************************************/
1188 : /* HFAGetPEString() */
1189 : /* */
1190 : /* Some files have a ProjectionX node containing the ESRI style */
1191 : /* PE_STRING. This function allows fetching from it. */
1192 : /************************************************************************/
1193 :
1194 81 : char *HFAGetPEString(HFAHandle hHFA)
1195 :
1196 : {
1197 81 : if (hHFA->nBands == 0)
1198 0 : return nullptr;
1199 :
1200 : // Get the HFA node.
1201 81 : HFAEntry *poProX = hHFA->papoBand[0]->poNode->GetNamedChild("ProjectionX");
1202 81 : if (poProX == nullptr)
1203 62 : return nullptr;
1204 :
1205 19 : const char *pszType = poProX->GetStringField("projection.type.string");
1206 19 : if (pszType == nullptr || !EQUAL(pszType, "PE_COORDSYS"))
1207 0 : return nullptr;
1208 :
1209 : // Use a gross hack to scan ahead to the actual projection
1210 : // string. We do it this way because we don't have general
1211 : // handling for MIFObjects.
1212 19 : GByte *pabyData = poProX->GetData();
1213 19 : int nDataSize = poProX->GetDataSize();
1214 :
1215 1767 : while (nDataSize > 10 &&
1216 1767 : !STARTS_WITH_CI((const char *)pabyData, "PE_COORDSYS,."))
1217 : {
1218 1748 : pabyData++;
1219 1748 : nDataSize--;
1220 : }
1221 :
1222 19 : if (nDataSize < 31)
1223 0 : return nullptr;
1224 :
1225 : // Skip ahead to the actual string.
1226 19 : pabyData += 30;
1227 : // nDataSize -= 30;
1228 :
1229 19 : return CPLStrdup((const char *)pabyData);
1230 : }
1231 :
1232 : /************************************************************************/
1233 : /* HFASetPEString() */
1234 : /************************************************************************/
1235 :
1236 133 : CPLErr HFASetPEString(HFAHandle hHFA, const char *pszPEString)
1237 :
1238 : {
1239 133 : if (!CPLTestBool(CPLGetConfigOption("HFA_WRITE_PE_STRING", "YES")))
1240 0 : return CE_None;
1241 :
1242 : // Loop over bands, setting information on each one.
1243 304 : for (int iBand = 0; iBand < hHFA->nBands; iBand++)
1244 : {
1245 : // Verify we don't already have the node, since update-in-place
1246 : // is likely to be more complicated.
1247 : HFAEntry *poProX =
1248 171 : hHFA->papoBand[iBand]->poNode->GetNamedChild("ProjectionX");
1249 :
1250 : // If we are setting an empty string then a missing entry is equivalent.
1251 171 : if (strlen(pszPEString) == 0 && poProX == nullptr)
1252 1 : continue;
1253 :
1254 : // Create the node.
1255 170 : if (poProX == nullptr)
1256 : {
1257 338 : poProX = HFAEntry::New(hHFA, "ProjectionX", "Eprj_MapProjection842",
1258 169 : hHFA->papoBand[iBand]->poNode);
1259 169 : if (poProX->GetTypeObject() == nullptr)
1260 0 : return CE_Failure;
1261 : }
1262 :
1263 : // Prepare the data area with some extra space just in case.
1264 : GByte *pabyData =
1265 170 : poProX->MakeData(static_cast<int>(700 + strlen(pszPEString)));
1266 170 : if (!pabyData)
1267 0 : return CE_Failure;
1268 :
1269 170 : memset(pabyData, 0, 250 + strlen(pszPEString));
1270 :
1271 170 : poProX->SetPosition();
1272 :
1273 170 : poProX->SetStringField("projection.type.string", "PE_COORDSYS");
1274 170 : poProX->SetStringField("projection.MIFDictionary.string",
1275 : "{0:pcstring,}Emif_String,{1:x{0:pcstring,}"
1276 : "Emif_String,coordSys,}PE_COORDSYS,.");
1277 :
1278 : // Use a gross hack to scan ahead to the actual projection
1279 : // string. We do it this way because we don't have general
1280 : // handling for MIFObjects
1281 170 : pabyData = poProX->GetData();
1282 170 : int nDataSize = poProX->GetDataSize();
1283 170 : GUInt32 iOffset = poProX->GetDataPos();
1284 :
1285 15810 : while (nDataSize > 10 &&
1286 15810 : !STARTS_WITH_CI((const char *)pabyData, "PE_COORDSYS,."))
1287 : {
1288 15640 : pabyData++;
1289 15640 : nDataSize--;
1290 15640 : iOffset++;
1291 : }
1292 :
1293 170 : CPLAssert(nDataSize > static_cast<int>(strlen(pszPEString)) + 10);
1294 :
1295 170 : pabyData += 14;
1296 170 : iOffset += 14;
1297 :
1298 : // Set the size and offset of the mifobject.
1299 170 : iOffset += 8;
1300 :
1301 170 : GUInt32 nSize = static_cast<GUInt32>(strlen(pszPEString) + 9);
1302 :
1303 : HFAStandard(4, &nSize);
1304 170 : memcpy(pabyData, &nSize, 4);
1305 170 : pabyData += 4;
1306 :
1307 : HFAStandard(4, &iOffset);
1308 170 : memcpy(pabyData, &iOffset, 4);
1309 170 : pabyData += 4;
1310 :
1311 : // Set the size and offset of the string value.
1312 170 : nSize = static_cast<GUInt32>(strlen(pszPEString) + 1);
1313 :
1314 : HFAStandard(4, &nSize);
1315 170 : memcpy(pabyData, &nSize, 4);
1316 170 : pabyData += 4;
1317 :
1318 170 : iOffset = 8;
1319 : HFAStandard(4, &iOffset);
1320 170 : memcpy(pabyData, &iOffset, 4);
1321 170 : pabyData += 4;
1322 :
1323 : // Place the string itself.
1324 170 : memcpy(pabyData, pszPEString, strlen(pszPEString) + 1);
1325 :
1326 170 : poProX->SetStringField("title.string", "PE");
1327 : }
1328 :
1329 133 : return CE_None;
1330 : }
1331 :
1332 : /************************************************************************/
1333 : /* HFAGetProParameters() */
1334 : /************************************************************************/
1335 :
1336 543 : const Eprj_ProParameters *HFAGetProParameters(HFAHandle hHFA)
1337 :
1338 : {
1339 543 : if (hHFA->nBands < 1)
1340 0 : return nullptr;
1341 :
1342 : // Do we already have it?
1343 543 : if (hHFA->pProParameters != nullptr)
1344 0 : return (Eprj_ProParameters *)hHFA->pProParameters;
1345 :
1346 : // Get the HFA node.
1347 : HFAEntry *poMIEntry =
1348 543 : hHFA->papoBand[0]->poNode->GetNamedChild("Projection");
1349 543 : if (poMIEntry == nullptr)
1350 340 : return nullptr;
1351 :
1352 : // Allocate the structure.
1353 : Eprj_ProParameters *psProParams = static_cast<Eprj_ProParameters *>(
1354 203 : CPLCalloc(sizeof(Eprj_ProParameters), 1));
1355 :
1356 : // Fetch the fields.
1357 203 : const int proType = poMIEntry->GetIntField("proType");
1358 203 : if (proType != EPRJ_INTERNAL && proType != EPRJ_EXTERNAL)
1359 : {
1360 0 : CPLError(CE_Failure, CPLE_AppDefined, "Wrong value for proType");
1361 0 : CPLFree(psProParams);
1362 0 : return nullptr;
1363 : }
1364 203 : psProParams->proType = static_cast<Eprj_ProType>(proType);
1365 203 : psProParams->proNumber = poMIEntry->GetIntField("proNumber");
1366 203 : psProParams->proExeName =
1367 203 : CPLStrdup(poMIEntry->GetStringField("proExeName"));
1368 203 : psProParams->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
1369 203 : psProParams->proZone = poMIEntry->GetIntField("proZone");
1370 :
1371 3248 : for (int i = 0; i < 15; i++)
1372 : {
1373 3045 : char szFieldName[40] = {};
1374 :
1375 3045 : snprintf(szFieldName, sizeof(szFieldName), "proParams[%d]", i);
1376 3045 : psProParams->proParams[i] = poMIEntry->GetDoubleField(szFieldName);
1377 : }
1378 :
1379 203 : psProParams->proSpheroid.sphereName =
1380 203 : CPLStrdup(poMIEntry->GetStringField("proSpheroid.sphereName"));
1381 203 : psProParams->proSpheroid.a = poMIEntry->GetDoubleField("proSpheroid.a");
1382 203 : psProParams->proSpheroid.b = poMIEntry->GetDoubleField("proSpheroid.b");
1383 203 : psProParams->proSpheroid.eSquared =
1384 203 : poMIEntry->GetDoubleField("proSpheroid.eSquared");
1385 203 : psProParams->proSpheroid.radius =
1386 203 : poMIEntry->GetDoubleField("proSpheroid.radius");
1387 :
1388 203 : hHFA->pProParameters = (void *)psProParams;
1389 :
1390 203 : return psProParams;
1391 : }
1392 :
1393 : /************************************************************************/
1394 : /* HFASetProParameters() */
1395 : /************************************************************************/
1396 :
1397 131 : CPLErr HFASetProParameters(HFAHandle hHFA, const Eprj_ProParameters *poPro)
1398 :
1399 : {
1400 : // Loop over bands, setting information on each one.
1401 300 : for (int iBand = 0; iBand < hHFA->nBands; iBand++)
1402 : {
1403 : // Create a new Projection if there isn't one present already.
1404 : HFAEntry *poMIEntry =
1405 169 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1406 169 : if (poMIEntry == nullptr)
1407 : {
1408 168 : poMIEntry = HFAEntry::New(hHFA, "Projection", "Eprj_ProParameters",
1409 168 : hHFA->papoBand[iBand]->poNode);
1410 : }
1411 :
1412 169 : poMIEntry->MarkDirty();
1413 :
1414 : // Ensure we have enough space for all the data.
1415 : // TODO(schwehr): Explain all these constants.
1416 169 : int nSize =
1417 169 : static_cast<int>(34 + 15 * 8 + 8 + strlen(poPro->proName) + 1 + 32 +
1418 169 : 8 + strlen(poPro->proSpheroid.sphereName) + 1);
1419 :
1420 169 : if (poPro->proExeName != nullptr)
1421 1 : nSize += static_cast<int>(strlen(poPro->proExeName) + 1);
1422 :
1423 169 : GByte *pabyData = poMIEntry->MakeData(nSize);
1424 169 : if (!pabyData)
1425 0 : return CE_Failure;
1426 :
1427 169 : poMIEntry->SetPosition();
1428 :
1429 : // Initialize the whole thing to zeros for a clean start.
1430 169 : memset(poMIEntry->GetData(), 0, poMIEntry->GetDataSize());
1431 :
1432 : // Write the various fields.
1433 169 : poMIEntry->SetIntField("proType", poPro->proType);
1434 :
1435 169 : poMIEntry->SetIntField("proNumber", poPro->proNumber);
1436 :
1437 169 : poMIEntry->SetStringField("proExeName", poPro->proExeName);
1438 169 : poMIEntry->SetStringField("proName", poPro->proName);
1439 169 : poMIEntry->SetIntField("proZone", poPro->proZone);
1440 169 : poMIEntry->SetDoubleField("proParams[0]", poPro->proParams[0]);
1441 169 : poMIEntry->SetDoubleField("proParams[1]", poPro->proParams[1]);
1442 169 : poMIEntry->SetDoubleField("proParams[2]", poPro->proParams[2]);
1443 169 : poMIEntry->SetDoubleField("proParams[3]", poPro->proParams[3]);
1444 169 : poMIEntry->SetDoubleField("proParams[4]", poPro->proParams[4]);
1445 169 : poMIEntry->SetDoubleField("proParams[5]", poPro->proParams[5]);
1446 169 : poMIEntry->SetDoubleField("proParams[6]", poPro->proParams[6]);
1447 169 : poMIEntry->SetDoubleField("proParams[7]", poPro->proParams[7]);
1448 169 : poMIEntry->SetDoubleField("proParams[8]", poPro->proParams[8]);
1449 169 : poMIEntry->SetDoubleField("proParams[9]", poPro->proParams[9]);
1450 169 : poMIEntry->SetDoubleField("proParams[10]", poPro->proParams[10]);
1451 169 : poMIEntry->SetDoubleField("proParams[11]", poPro->proParams[11]);
1452 169 : poMIEntry->SetDoubleField("proParams[12]", poPro->proParams[12]);
1453 169 : poMIEntry->SetDoubleField("proParams[13]", poPro->proParams[13]);
1454 169 : poMIEntry->SetDoubleField("proParams[14]", poPro->proParams[14]);
1455 169 : poMIEntry->SetStringField("proSpheroid.sphereName",
1456 169 : poPro->proSpheroid.sphereName);
1457 169 : poMIEntry->SetDoubleField("proSpheroid.a", poPro->proSpheroid.a);
1458 169 : poMIEntry->SetDoubleField("proSpheroid.b", poPro->proSpheroid.b);
1459 169 : poMIEntry->SetDoubleField("proSpheroid.eSquared",
1460 169 : poPro->proSpheroid.eSquared);
1461 169 : poMIEntry->SetDoubleField("proSpheroid.radius",
1462 169 : poPro->proSpheroid.radius);
1463 : }
1464 :
1465 131 : return CE_None;
1466 : }
1467 :
1468 : /************************************************************************/
1469 : /* HFAGetDatum() */
1470 : /************************************************************************/
1471 :
1472 543 : const Eprj_Datum *HFAGetDatum(HFAHandle hHFA)
1473 :
1474 : {
1475 543 : if (hHFA->nBands < 1)
1476 0 : return nullptr;
1477 :
1478 : // Do we already have it?
1479 543 : if (hHFA->pDatum != nullptr)
1480 0 : return (Eprj_Datum *)hHFA->pDatum;
1481 :
1482 : // Get the HFA node.
1483 : HFAEntry *poMIEntry =
1484 543 : hHFA->papoBand[0]->poNode->GetNamedChild("Projection.Datum");
1485 543 : if (poMIEntry == nullptr)
1486 340 : return nullptr;
1487 :
1488 : // Allocate the structure.
1489 : Eprj_Datum *psDatum =
1490 203 : static_cast<Eprj_Datum *>(CPLCalloc(sizeof(Eprj_Datum), 1));
1491 :
1492 : // Fetch the fields.
1493 203 : psDatum->datumname = CPLStrdup(poMIEntry->GetStringField("datumname"));
1494 203 : const int nDatumType = poMIEntry->GetIntField("type");
1495 203 : if (nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE)
1496 : {
1497 1 : CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
1498 1 : psDatum->type = EPRJ_DATUM_NONE;
1499 : }
1500 : else
1501 202 : psDatum->type = static_cast<Eprj_DatumType>(nDatumType);
1502 :
1503 1624 : for (int i = 0; i < 7; i++)
1504 : {
1505 1421 : char szFieldName[30] = {};
1506 1421 : snprintf(szFieldName, sizeof(szFieldName), "params[%d]", i);
1507 1421 : psDatum->params[i] = poMIEntry->GetDoubleField(szFieldName);
1508 : }
1509 :
1510 203 : psDatum->gridname = CPLStrdup(poMIEntry->GetStringField("gridname"));
1511 :
1512 203 : hHFA->pDatum = (void *)psDatum;
1513 :
1514 203 : return psDatum;
1515 : }
1516 :
1517 : /************************************************************************/
1518 : /* HFASetDatum() */
1519 : /************************************************************************/
1520 :
1521 131 : CPLErr HFASetDatum(HFAHandle hHFA, const Eprj_Datum *poDatum)
1522 :
1523 : {
1524 : // Loop over bands, setting information on each one.
1525 300 : for (int iBand = 0; iBand < hHFA->nBands; iBand++)
1526 : {
1527 : // Create a new Projection if there isn't one present already.
1528 : HFAEntry *poProParams =
1529 169 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1530 169 : if (poProParams == nullptr)
1531 : {
1532 0 : CPLError(CE_Failure, CPLE_AppDefined,
1533 : "Can't add Eprj_Datum with no Eprj_ProjParameters.");
1534 0 : return CE_Failure;
1535 : }
1536 :
1537 169 : HFAEntry *poDatumEntry = poProParams->GetNamedChild("Datum");
1538 169 : if (poDatumEntry == nullptr)
1539 : {
1540 : poDatumEntry =
1541 168 : HFAEntry::New(hHFA, "Datum", "Eprj_Datum", poProParams);
1542 : }
1543 :
1544 169 : poDatumEntry->MarkDirty();
1545 :
1546 : // Ensure we have enough space for all the data.
1547 : // TODO(schwehr): Explain constants.
1548 169 : int nSize =
1549 169 : static_cast<int>(26 + strlen(poDatum->datumname) + 1 + 7 * 8);
1550 :
1551 169 : if (poDatum->gridname != nullptr)
1552 30 : nSize += static_cast<int>(strlen(poDatum->gridname) + 1);
1553 :
1554 169 : GByte *pabyData = poDatumEntry->MakeData(nSize);
1555 169 : if (!pabyData)
1556 0 : return CE_Failure;
1557 :
1558 169 : poDatumEntry->SetPosition();
1559 :
1560 : // Initialize the whole thing to zeros for a clean start.
1561 169 : memset(poDatumEntry->GetData(), 0, poDatumEntry->GetDataSize());
1562 :
1563 : // Write the various fields.
1564 169 : poDatumEntry->SetStringField("datumname", poDatum->datumname);
1565 169 : poDatumEntry->SetIntField("type", poDatum->type);
1566 :
1567 169 : poDatumEntry->SetDoubleField("params[0]", poDatum->params[0]);
1568 169 : poDatumEntry->SetDoubleField("params[1]", poDatum->params[1]);
1569 169 : poDatumEntry->SetDoubleField("params[2]", poDatum->params[2]);
1570 169 : poDatumEntry->SetDoubleField("params[3]", poDatum->params[3]);
1571 169 : poDatumEntry->SetDoubleField("params[4]", poDatum->params[4]);
1572 169 : poDatumEntry->SetDoubleField("params[5]", poDatum->params[5]);
1573 169 : poDatumEntry->SetDoubleField("params[6]", poDatum->params[6]);
1574 :
1575 169 : poDatumEntry->SetStringField("gridname", poDatum->gridname);
1576 : }
1577 :
1578 131 : return CE_None;
1579 : }
1580 :
1581 : /************************************************************************/
1582 : /* HFAGetPCT() */
1583 : /* */
1584 : /* Read the PCT from a band, if it has one. */
1585 : /************************************************************************/
1586 :
1587 623 : CPLErr HFAGetPCT(HFAHandle hHFA, int nBand, int *pnColors, double **ppadfRed,
1588 : double **ppadfGreen, double **ppadfBlue, double **ppadfAlpha,
1589 : double **ppadfBins)
1590 :
1591 : {
1592 623 : if (nBand < 1 || nBand > hHFA->nBands)
1593 0 : return CE_Failure;
1594 :
1595 623 : return hHFA->papoBand[nBand - 1]->GetPCT(pnColors, ppadfRed, ppadfGreen,
1596 623 : ppadfBlue, ppadfAlpha, ppadfBins);
1597 : }
1598 :
1599 : /************************************************************************/
1600 : /* HFASetPCT() */
1601 : /* */
1602 : /* Set the PCT on a band. */
1603 : /************************************************************************/
1604 :
1605 3 : CPLErr HFASetPCT(HFAHandle hHFA, int nBand, int nColors, double *padfRed,
1606 : double *padfGreen, double *padfBlue, double *padfAlpha)
1607 :
1608 : {
1609 3 : if (nBand < 1 || nBand > hHFA->nBands)
1610 0 : return CE_Failure;
1611 :
1612 3 : return hHFA->papoBand[nBand - 1]->SetPCT(nColors, padfRed, padfGreen,
1613 3 : padfBlue, padfAlpha);
1614 : }
1615 :
1616 : /************************************************************************/
1617 : /* HFAGetDataRange() */
1618 : /************************************************************************/
1619 :
1620 0 : CPLErr HFAGetDataRange(HFAHandle hHFA, int nBand, double *pdfMin,
1621 : double *pdfMax)
1622 :
1623 : {
1624 0 : if (nBand < 1 || nBand > hHFA->nBands)
1625 0 : return CE_Failure;
1626 :
1627 : HFAEntry *poBinInfo =
1628 0 : hHFA->papoBand[nBand - 1]->poNode->GetNamedChild("Statistics");
1629 :
1630 0 : if (poBinInfo == nullptr)
1631 0 : return CE_Failure;
1632 :
1633 0 : *pdfMin = poBinInfo->GetDoubleField("minimum");
1634 0 : *pdfMax = poBinInfo->GetDoubleField("maximum");
1635 :
1636 0 : if (*pdfMax > *pdfMin)
1637 0 : return CE_None;
1638 : else
1639 0 : return CE_Failure;
1640 : }
1641 :
1642 : /************************************************************************/
1643 : /* HFADumpNode() */
1644 : /************************************************************************/
1645 :
1646 0 : static void HFADumpNode(HFAEntry *poEntry, int nIndent, bool bVerbose, FILE *fp)
1647 :
1648 : {
1649 0 : std::string osSpaces(nIndent * 2, ' ');
1650 :
1651 0 : fprintf(fp, "%s%s(%s) @ %u + %u @ %u\n", osSpaces.c_str(),
1652 : poEntry->GetName(), poEntry->GetType(), poEntry->GetFilePos(),
1653 : poEntry->GetDataSize(), poEntry->GetDataPos());
1654 :
1655 0 : if (bVerbose)
1656 : {
1657 0 : osSpaces += "+ ";
1658 0 : poEntry->DumpFieldValues(fp, osSpaces.c_str());
1659 0 : fprintf(fp, "\n");
1660 : }
1661 :
1662 0 : if (poEntry->GetChild() != nullptr)
1663 0 : HFADumpNode(poEntry->GetChild(), nIndent + 1, bVerbose, fp);
1664 :
1665 0 : if (poEntry->GetNext() != nullptr)
1666 0 : HFADumpNode(poEntry->GetNext(), nIndent, bVerbose, fp);
1667 0 : }
1668 :
1669 : /************************************************************************/
1670 : /* HFADumpTree() */
1671 : /* */
1672 : /* Dump the tree of information in a HFA file. */
1673 : /************************************************************************/
1674 :
1675 0 : void HFADumpTree(HFAHandle hHFA, FILE *fpOut)
1676 :
1677 : {
1678 0 : HFADumpNode(hHFA->poRoot, 0, true, fpOut);
1679 0 : }
1680 :
1681 : /************************************************************************/
1682 : /* HFADumpDictionary() */
1683 : /* */
1684 : /* Dump the dictionary (in raw, and parsed form) to the named */
1685 : /* device. */
1686 : /************************************************************************/
1687 :
1688 0 : void HFADumpDictionary(HFAHandle hHFA, FILE *fpOut)
1689 :
1690 : {
1691 0 : fprintf(fpOut, "%s\n", hHFA->pszDictionary);
1692 :
1693 0 : hHFA->poDictionary->Dump(fpOut);
1694 0 : }
1695 :
1696 : /************************************************************************/
1697 : /* HFAStandard() */
1698 : /* */
1699 : /* Swap byte order on MSB systems. */
1700 : /************************************************************************/
1701 :
1702 : #ifdef CPL_MSB
1703 : void HFAStandard(int nBytes, void *pData)
1704 :
1705 : {
1706 : GByte *pabyData = static_cast<GByte *>(pData);
1707 :
1708 : for (int i = nBytes / 2 - 1; i >= 0; i--)
1709 : {
1710 : GByte byTemp = pabyData[i];
1711 : pabyData[i] = pabyData[nBytes - i - 1];
1712 : pabyData[nBytes - i - 1] = byTemp;
1713 : }
1714 : }
1715 : #endif
1716 :
1717 : /* ==================================================================== */
1718 : /* Default data dictionary. Emitted verbatim into the imagine */
1719 : /* file. */
1720 : /* ==================================================================== */
1721 :
1722 : static const char *const aszDefaultDD[] = {
1723 : "{1:lversion,1:LfreeList,1:LrootEntryPtr,1:sentryHeaderLength,1:"
1724 : "LdictionaryPtr,}Ehfa_File,{1:Lnext,1:Lprev,1:Lparent,1:Lchild,1:Ldata,1:"
1725 : "ldataSize,64:cname,32:ctype,1:tmodTime,}Ehfa_Entry,{16:clabel,1:"
1726 : "LheaderPtr,}Ehfa_HeaderTag,{1:LfreeList,1:lfreeSize,}Ehfa_FreeListNode,{1:"
1727 : "lsize,1:Lptr,}Ehfa_Data,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft "
1728 : "of real-valued data,layerType,",
1729 : "1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:"
1730 : "lblockWidth,1:lblockHeight,}Eimg_Layer,{1:lwidth,1:lheight,1:e3:thematic,"
1731 : "athematic,fft of real-valued "
1732 : "data,layerType,1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,"
1733 : "pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer_SubSample,{1:e2:raster,"
1734 : "vector,type,1:LdictionaryPtr,}Ehfa_Layer,{1:LspaceUsedForRasterData,}"
1735 : "ImgFormatInfo831,{1:sfileCode,1:Loffset,1:lsize,1:e2:false,true,logvalid,",
1736 : "1:e2:no compression,ESRI GRID "
1737 : "compression,compressionType,}Edms_VirtualBlockInfo,{1:lmin,1:lmax,}Edms_"
1738 : "FreeIDList,{1:lnumvirtualblocks,1:lnumobjectsperblock,1:lnextobjectnum,1:"
1739 : "e2:no compression,RLC "
1740 : "compression,compressionType,0:poEdms_VirtualBlockInfo,blockinfo,0:poEdms_"
1741 : "FreeIDList,freelist,1:tmodTime,}Edms_State,{0:pcstring,}Emif_String,{1:"
1742 : "oEmif_String,fileName,2:LlayerStackValidFlagsOffset,2:"
1743 : "LlayerStackDataOffset,1:LlayerStackCount,1:LlayerStackIndex,}"
1744 : "ImgExternalRaster,{1:oEmif_String,algorithm,0:poEmif_String,nameList,}"
1745 : "Eimg_RRDNamesList,{1:oEmif_String,projection,1:oEmif_String,units,}Eimg_"
1746 : "MapInformation,",
1747 : "{1:oEmif_String,dependent,}Eimg_DependentFile,{1:oEmif_String,"
1748 : "ImageLayerName,}Eimg_DependentLayerName,{1:lnumrows,1:lnumcolumns,1:e13:"
1749 : "EGDA_TYPE_U1,EGDA_TYPE_U2,EGDA_TYPE_U4,EGDA_TYPE_U8,EGDA_TYPE_S8,EGDA_"
1750 : "TYPE_U16,EGDA_TYPE_S16,EGDA_TYPE_U32,EGDA_TYPE_S32,EGDA_TYPE_F32,EGDA_"
1751 : "TYPE_F64,EGDA_TYPE_C64,EGDA_TYPE_C128,datatype,1:e4:EGDA_SCALAR_OBJECT,"
1752 : "EGDA_TABLE_OBJECT,EGDA_MATRIX_OBJECT,EGDA_RASTER_OBJECT,objecttype,}Egda_"
1753 : "BaseData,{1:*bvalueBD,}Eimg_NonInitializedValue,{1:dx,1:dy,}Eprj_"
1754 : "Coordinate,{1:dwidth,1:dheight,}Eprj_Size,{0:pcproName,1:*oEprj_"
1755 : "Coordinate,upperLeftCenter,",
1756 : "1:*oEprj_Coordinate,lowerRightCenter,1:*oEprj_Size,pixelSize,0:pcunits,}"
1757 : "Eprj_MapInfo,{0:pcdatumname,1:e3:EPRJ_DATUM_PARAMETRIC,EPRJ_DATUM_GRID,"
1758 : "EPRJ_DATUM_REGRESSION,type,0:pdparams,0:pcgridname,}Eprj_Datum,{0:"
1759 : "pcsphereName,1:da,1:db,1:deSquared,1:dradius,}Eprj_Spheroid,{1:e2:EPRJ_"
1760 : "INTERNAL,EPRJ_EXTERNAL,proType,1:lproNumber,0:pcproExeName,0:pcproName,1:"
1761 : "lproZone,0:pdproParams,1:*oEprj_Spheroid,proSpheroid,}Eprj_ProParameters,{"
1762 : "1:dminimum,1:dmaximum,1:dmean,1:dmedian,1:dmode,1:dstddev,}Esta_"
1763 : "Statistics,{1:lnumBins,1:e4:direct,linear,logarithmic,explicit,"
1764 : "binFunctionType,1:dminLimit,1:dmaxLimit,1:*bbinLimits,}Edsc_BinFunction,{"
1765 : "0:poEmif_String,LayerNames,1:*bExcludedValues,1:oEmif_String,AOIname,",
1766 : "1:lSkipFactorX,1:lSkipFactorY,1:*oEdsc_BinFunction,BinFunction,}Eimg_"
1767 : "StatisticsParameters830,{1:lnumrows,}Edsc_Table,{1:lnumRows,1:"
1768 : "LcolumnDataPtr,1:e4:integer,real,complex,string,dataType,1:lmaxNumChars,}"
1769 : "Edsc_Column,{1:lposition,0:pcname,1:e2:EMSC_FALSE,EMSC_TRUE,editable,1:e3:"
1770 : "LEFT,CENTER,RIGHT,alignment,0:pcformat,1:e3:DEFAULT,APPLY,AUTO-APPLY,"
1771 : "formulamode,0:pcformula,1:dcolumnwidth,0:pcunits,1:e5:NO_COLOR,RED,GREEN,"
1772 : "BLUE,COLOR,colorflag,0:pcgreenname,0:pcbluename,}Eded_ColumnAttributes_1,{"
1773 : "1:lversion,1:lnumobjects,1:e2:EAOI_UNION,EAOI_INTERSECTION,operation,}"
1774 : "Eaoi_AreaOfInterest,",
1775 : "{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,"
1776 : "MIFDictionary,0:pCMIFObject,}Emif_MIFObject,",
1777 : "{1:x{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,"
1778 : "MIFDictionary,0:pCMIFObject,}Emif_MIFObject,projection,1:x{0:pcstring,}"
1779 : "Emif_String,title,}Eprj_MapProjection842,",
1780 : "{0:poEmif_String,titleList,}Exfr_GenericXFormHeader,{1:lorder,1:"
1781 : "lnumdimtransform,1:lnumdimpolynomial,1:ltermcount,0:plexponentlist,1:*"
1782 : "bpolycoefmtx,1:*bpolycoefvector,}Efga_Polynomial,",
1783 : ".",
1784 : nullptr};
1785 :
1786 : /************************************************************************/
1787 : /* HFACreateLL() */
1788 : /* */
1789 : /* Low level creation of an Imagine file. Writes out the */
1790 : /* Ehfa_HeaderTag, dictionary and Ehfa_File. */
1791 : /************************************************************************/
1792 :
1793 204 : HFAHandle HFACreateLL(const char *pszFilename)
1794 :
1795 : {
1796 : // Create the file in the file system.
1797 204 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
1798 204 : if (fp == nullptr)
1799 : {
1800 3 : CPLError(CE_Failure, CPLE_OpenFailed, "Creation of file %s failed.",
1801 : pszFilename);
1802 3 : return nullptr;
1803 : }
1804 :
1805 : // Create the HFAInfo_t.
1806 : HFAInfo_t *psInfo =
1807 201 : static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
1808 :
1809 201 : psInfo->fp = fp;
1810 201 : psInfo->eAccess = HFA_Update;
1811 201 : psInfo->nXSize = 0;
1812 201 : psInfo->nYSize = 0;
1813 201 : psInfo->nBands = 0;
1814 201 : psInfo->papoBand = nullptr;
1815 201 : psInfo->pMapInfo = nullptr;
1816 201 : psInfo->pDatum = nullptr;
1817 201 : psInfo->pProParameters = nullptr;
1818 201 : psInfo->bTreeDirty = false;
1819 201 : psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
1820 201 : psInfo->pszPath = CPLStrdup(CPLGetPathSafe(pszFilename).c_str());
1821 :
1822 : // Write out the Ehfa_HeaderTag.
1823 201 : bool bRet = VSIFWriteL((void *)"EHFA_HEADER_TAG", 1, 16, fp) > 0;
1824 :
1825 201 : GInt32 nHeaderPos = 20;
1826 : HFAStandard(4, &nHeaderPos);
1827 201 : bRet &= VSIFWriteL(&nHeaderPos, 4, 1, fp) > 0;
1828 :
1829 : // Write the Ehfa_File node, locked in at offset 20.
1830 201 : GInt32 nVersion = 1;
1831 201 : GInt32 nFreeList = 0;
1832 201 : GInt32 nRootEntry = 0;
1833 201 : GInt16 nEntryHeaderLength = 128;
1834 201 : GInt32 nDictionaryPtr = 38;
1835 :
1836 201 : psInfo->nEntryHeaderLength = nEntryHeaderLength;
1837 201 : psInfo->nRootPos = 0;
1838 201 : psInfo->nDictionaryPos = nDictionaryPtr;
1839 201 : psInfo->nVersion = nVersion;
1840 :
1841 : HFAStandard(4, &nVersion);
1842 : HFAStandard(4, &nFreeList);
1843 : HFAStandard(4, &nRootEntry);
1844 : HFAStandard(2, &nEntryHeaderLength);
1845 : HFAStandard(4, &nDictionaryPtr);
1846 :
1847 201 : bRet &= VSIFWriteL(&nVersion, 4, 1, fp) > 0;
1848 201 : bRet &= VSIFWriteL(&nFreeList, 4, 1, fp) > 0;
1849 201 : bRet &= VSIFWriteL(&nRootEntry, 4, 1, fp) > 0;
1850 201 : bRet &= VSIFWriteL(&nEntryHeaderLength, 2, 1, fp) > 0;
1851 201 : bRet &= VSIFWriteL(&nDictionaryPtr, 4, 1, fp) > 0;
1852 :
1853 : // Write the dictionary, locked in at location 38. Note that
1854 : // we jump through a bunch of hoops to operate on the
1855 : // dictionary in chunks because some compiles (such as VC++)
1856 : // don't allow particularly large static strings.
1857 201 : int nDictLen = 0;
1858 :
1859 2211 : for (int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++)
1860 2010 : nDictLen += static_cast<int>(strlen(aszDefaultDD[iChunk]));
1861 :
1862 201 : psInfo->pszDictionary = static_cast<char *>(CPLMalloc(nDictLen + 1));
1863 201 : psInfo->pszDictionary[0] = '\0';
1864 :
1865 2211 : for (int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++)
1866 2010 : strcat(psInfo->pszDictionary, aszDefaultDD[iChunk]);
1867 :
1868 402 : bRet &= VSIFWriteL((void *)psInfo->pszDictionary,
1869 201 : strlen(psInfo->pszDictionary) + 1, 1, fp) > 0;
1870 201 : if (!bRet)
1871 : {
1872 4 : CPL_IGNORE_RET_VAL(HFAClose(psInfo));
1873 4 : return nullptr;
1874 : }
1875 :
1876 197 : psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
1877 :
1878 197 : psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
1879 :
1880 : // Create a root entry.
1881 197 : psInfo->poRoot = new HFAEntry(psInfo, "root", "root", nullptr);
1882 :
1883 : // If an .ige or .rrd file exists with the same base name,
1884 : // delete them. (#1784)
1885 197 : CPLString osExtension = CPLGetExtensionSafe(pszFilename);
1886 197 : if (!EQUAL(osExtension, "rrd") && !EQUAL(osExtension, "aux"))
1887 : {
1888 382 : CPLString osPath = CPLGetPathSafe(pszFilename);
1889 382 : CPLString osBasename = CPLGetBasenameSafe(pszFilename);
1890 : VSIStatBufL sStatBuf;
1891 382 : CPLString osSupFile = CPLFormCIFilenameSafe(osPath, osBasename, "rrd");
1892 :
1893 191 : if (VSIStatL(osSupFile, &sStatBuf) == 0)
1894 0 : VSIUnlink(osSupFile);
1895 :
1896 191 : osSupFile = CPLFormCIFilenameSafe(osPath, osBasename, "ige");
1897 :
1898 191 : if (VSIStatL(osSupFile, &sStatBuf) == 0)
1899 0 : VSIUnlink(osSupFile);
1900 : }
1901 :
1902 197 : return psInfo;
1903 : }
1904 :
1905 : /************************************************************************/
1906 : /* HFAAllocateSpace() */
1907 : /* */
1908 : /* Return an area in the file to the caller to write the */
1909 : /* requested number of bytes. Currently this is always at the */
1910 : /* end of the file, but eventually we might actually keep track */
1911 : /* of free space. The HFAInfo_t's concept of file size is */
1912 : /* updated, even if nothing ever gets written to this region. */
1913 : /* */
1914 : /* Returns the offset to the requested space, or zero one */
1915 : /* failure. */
1916 : /************************************************************************/
1917 :
1918 2880 : vsi_l_offset HFAAllocateSpace(HFAInfo_t *psInfo, GUInt32 nBytes)
1919 :
1920 : {
1921 2880 : psInfo->nEndOfFile += nBytes;
1922 2880 : return psInfo->nEndOfFile - nBytes;
1923 : }
1924 :
1925 : /************************************************************************/
1926 : /* HFAFlush() */
1927 : /* */
1928 : /* Write out any dirty tree information to disk, putting the */
1929 : /* disk file in a consistent state. */
1930 : /************************************************************************/
1931 :
1932 400 : CPLErr HFAFlush(HFAHandle hHFA)
1933 :
1934 : {
1935 400 : if (!hHFA->bTreeDirty && !hHFA->poDictionary->bDictionaryTextDirty)
1936 0 : return CE_None;
1937 :
1938 400 : CPLAssert(hHFA->poRoot != nullptr);
1939 :
1940 : // Flush HFAEntry tree to disk.
1941 400 : if (hHFA->bTreeDirty)
1942 : {
1943 400 : const CPLErr eErr = hHFA->poRoot->FlushToDisk();
1944 400 : if (eErr != CE_None)
1945 6 : return eErr;
1946 :
1947 394 : hHFA->bTreeDirty = false;
1948 : }
1949 :
1950 : // Flush Dictionary to disk.
1951 394 : GUInt32 nNewDictionaryPos = hHFA->nDictionaryPos;
1952 394 : bool bRet = true;
1953 394 : if (hHFA->poDictionary->bDictionaryTextDirty)
1954 : {
1955 0 : bRet &= VSIFSeekL(hHFA->fp, 0, SEEK_END) >= 0;
1956 0 : nNewDictionaryPos = static_cast<GUInt32>(VSIFTellL(hHFA->fp));
1957 0 : bRet &=
1958 0 : VSIFWriteL(hHFA->poDictionary->osDictionaryText.c_str(),
1959 0 : strlen(hHFA->poDictionary->osDictionaryText.c_str()) + 1,
1960 0 : 1, hHFA->fp) > 0;
1961 0 : hHFA->poDictionary->bDictionaryTextDirty = false;
1962 : }
1963 :
1964 : // Do we need to update the Ehfa_File pointer to the root node?
1965 596 : if (hHFA->nRootPos != hHFA->poRoot->GetFilePos() ||
1966 202 : nNewDictionaryPos != hHFA->nDictionaryPos)
1967 : {
1968 192 : GUInt32 nHeaderPos = 0;
1969 :
1970 192 : bRet &= VSIFSeekL(hHFA->fp, 16, SEEK_SET) >= 0;
1971 192 : bRet &= VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, hHFA->fp) > 0;
1972 : HFAStandard(4, &nHeaderPos);
1973 :
1974 192 : GUInt32 nOffset = hHFA->poRoot->GetFilePos();
1975 192 : hHFA->nRootPos = nOffset;
1976 : HFAStandard(4, &nOffset);
1977 192 : bRet &= VSIFSeekL(hHFA->fp, static_cast<vsi_l_offset>(nHeaderPos) + 8,
1978 192 : SEEK_SET) >= 0;
1979 192 : bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
1980 :
1981 192 : nOffset = nNewDictionaryPos;
1982 192 : hHFA->nDictionaryPos = nNewDictionaryPos;
1983 : HFAStandard(4, &nOffset);
1984 192 : bRet &= VSIFSeekL(hHFA->fp, static_cast<vsi_l_offset>(nHeaderPos) + 14,
1985 192 : SEEK_SET) >= 0;
1986 192 : bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
1987 : }
1988 :
1989 394 : return bRet ? CE_None : CE_Failure;
1990 : }
1991 :
1992 : /************************************************************************/
1993 : /* HFACreateLayer() */
1994 : /* */
1995 : /* Create a layer object, and corresponding RasterDMS. */
1996 : /* Suitable for use with primary layers, and overviews. */
1997 : /************************************************************************/
1998 :
1999 245 : int HFACreateLayer(HFAHandle psInfo, HFAEntry *poParent,
2000 : const char *pszLayerName, int bOverview, int nBlockSize,
2001 : int bCreateCompressed, int bCreateLargeRaster,
2002 : int bDependentLayer, int nXSize, int nYSize,
2003 : EPTType eDataType, char ** /* papszOptions */,
2004 : // These are only related to external (large) files.
2005 : GIntBig nStackValidFlagsOffset, GIntBig nStackDataOffset,
2006 : int nStackCount, int nStackIndex)
2007 :
2008 : {
2009 245 : const char *pszLayerType =
2010 245 : bOverview ? "Eimg_Layer_SubSample" : "Eimg_Layer";
2011 :
2012 245 : if (nBlockSize <= 0)
2013 : {
2014 0 : CPLError(CE_Failure, CPLE_IllegalArg,
2015 : "HFACreateLayer: nBlockXSize < 0");
2016 0 : return FALSE;
2017 : }
2018 :
2019 : // Work out some details about the tiling scheme.
2020 245 : const int nBlocksPerRow = DIV_ROUND_UP(nXSize, nBlockSize);
2021 245 : const int nBlocksPerColumn = DIV_ROUND_UP(nYSize, nBlockSize);
2022 245 : const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
2023 : const int nBytesPerBlock =
2024 245 : (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
2025 :
2026 : // Create the Eimg_Layer for the band.
2027 : HFAEntry *poEimg_Layer =
2028 245 : HFAEntry::New(psInfo, pszLayerName, pszLayerType, poParent);
2029 :
2030 245 : poEimg_Layer->SetIntField("width", nXSize);
2031 245 : poEimg_Layer->SetIntField("height", nYSize);
2032 245 : poEimg_Layer->SetStringField("layerType", "athematic");
2033 245 : poEimg_Layer->SetIntField("pixelType", eDataType);
2034 245 : poEimg_Layer->SetIntField("blockWidth", nBlockSize);
2035 245 : poEimg_Layer->SetIntField("blockHeight", nBlockSize);
2036 :
2037 : // Create the RasterDMS (block list). This is a complex type
2038 : // with pointers, and variable size. We set the superstructure
2039 : // ourselves rather than trying to have the HFA type management
2040 : // system do it for us (since this would be hard to implement).
2041 245 : if (!bCreateLargeRaster && !bDependentLayer)
2042 : {
2043 : HFAEntry *poEdms_State =
2044 233 : HFAEntry::New(psInfo, "RasterDMS", "Edms_State", poEimg_Layer);
2045 :
2046 : // TODO(schwehr): Explain constants.
2047 233 : const int nDmsSize = 14 * nBlocks + 38;
2048 233 : GByte *pabyData = poEdms_State->MakeData(nDmsSize);
2049 :
2050 : // Set some simple values.
2051 233 : poEdms_State->SetIntField("numvirtualblocks", nBlocks);
2052 233 : poEdms_State->SetIntField("numobjectsperblock",
2053 : nBlockSize * nBlockSize);
2054 233 : poEdms_State->SetIntField("nextobjectnum",
2055 233 : nBlockSize * nBlockSize * nBlocks);
2056 :
2057 : // Is file compressed or not?
2058 233 : if (bCreateCompressed)
2059 : {
2060 10 : poEdms_State->SetStringField("compressionType", "RLC compression");
2061 : }
2062 : else
2063 : {
2064 223 : poEdms_State->SetStringField("compressionType", "no compression");
2065 : }
2066 :
2067 : // We need to hardcode file offset into the data, so locate it now.
2068 233 : poEdms_State->SetPosition();
2069 :
2070 : // Set block info headers.
2071 :
2072 : // Blockinfo count.
2073 233 : GUInt32 nValue = nBlocks;
2074 : HFAStandard(4, &nValue);
2075 233 : memcpy(pabyData + 14, &nValue, 4);
2076 :
2077 : // Blockinfo position.
2078 233 : nValue = poEdms_State->GetDataPos() + 22;
2079 : HFAStandard(4, &nValue);
2080 233 : memcpy(pabyData + 18, &nValue, 4);
2081 :
2082 : // Set each blockinfo.
2083 734 : for (int iBlock = 0; iBlock < nBlocks; iBlock++)
2084 : {
2085 501 : int nOffset = 22 + 14 * iBlock;
2086 :
2087 : // fileCode.
2088 501 : GInt16 nValue16 = 0;
2089 : HFAStandard(2, &nValue16);
2090 501 : memcpy(pabyData + nOffset, &nValue16, 2);
2091 :
2092 : // Offset.
2093 501 : if (bCreateCompressed)
2094 : {
2095 : // Flag it with zero offset. Allocate space when we compress it.
2096 18 : nValue = 0;
2097 : }
2098 : else
2099 : {
2100 483 : const auto nValue64 = HFAAllocateSpace(psInfo, nBytesPerBlock);
2101 483 : if (nValue64 > UINT32_MAX)
2102 0 : return FALSE;
2103 483 : nValue = static_cast<GUInt32>(nValue64);
2104 : }
2105 : HFAStandard(4, &nValue);
2106 501 : memcpy(pabyData + nOffset + 2, &nValue, 4);
2107 :
2108 : // Size.
2109 501 : if (bCreateCompressed)
2110 : {
2111 : // Flag with zero size. Don't know until we compress it.
2112 18 : nValue = 0;
2113 : }
2114 : else
2115 : {
2116 483 : nValue = nBytesPerBlock;
2117 : }
2118 : HFAStandard(4, &nValue);
2119 501 : memcpy(pabyData + nOffset + 6, &nValue, 4);
2120 :
2121 : // logValid (false).
2122 501 : nValue16 = 0;
2123 : HFAStandard(2, &nValue16);
2124 501 : memcpy(pabyData + nOffset + 10, &nValue16, 2);
2125 :
2126 : // compressionType.
2127 501 : if (bCreateCompressed)
2128 18 : nValue16 = 1;
2129 : else
2130 483 : nValue16 = 0;
2131 :
2132 : HFAStandard(2, &nValue16);
2133 501 : memcpy(pabyData + nOffset + 12, &nValue16, 2);
2134 233 : }
2135 : }
2136 :
2137 : // Create ExternalRasterDMS object.
2138 12 : else if (bCreateLargeRaster)
2139 : {
2140 7 : HFAEntry *poEdms_State = HFAEntry::New(
2141 : psInfo, "ExternalRasterDMS", "ImgExternalRaster", poEimg_Layer);
2142 7 : poEdms_State->MakeData(
2143 7 : static_cast<int>(8 + strlen(psInfo->pszIGEFilename) + 1 + 6 * 4));
2144 :
2145 7 : poEdms_State->SetStringField("fileName.string", psInfo->pszIGEFilename);
2146 :
2147 7 : poEdms_State->SetIntField(
2148 : "layerStackValidFlagsOffset[0]",
2149 : static_cast<int>(nStackValidFlagsOffset & 0xFFFFFFFF));
2150 7 : poEdms_State->SetIntField(
2151 : "layerStackValidFlagsOffset[1]",
2152 7 : static_cast<int>(nStackValidFlagsOffset >> 32));
2153 :
2154 7 : poEdms_State->SetIntField(
2155 : "layerStackDataOffset[0]",
2156 : static_cast<int>(nStackDataOffset & 0xFFFFFFFF));
2157 7 : poEdms_State->SetIntField("layerStackDataOffset[1]",
2158 7 : static_cast<int>(nStackDataOffset >> 32));
2159 7 : poEdms_State->SetIntField("layerStackCount", nStackCount);
2160 7 : poEdms_State->SetIntField("layerStackIndex", nStackIndex);
2161 : }
2162 : // Dependent...
2163 5 : else if (bDependentLayer)
2164 : {
2165 : HFAEntry *poDepLayerName =
2166 5 : HFAEntry::New(psInfo, "DependentLayerName",
2167 : "Eimg_DependentLayerName", poEimg_Layer);
2168 5 : poDepLayerName->MakeData(
2169 5 : static_cast<int>(8 + strlen(pszLayerName) + 2));
2170 :
2171 5 : poDepLayerName->SetStringField("ImageLayerName.string", pszLayerName);
2172 : }
2173 :
2174 : // Create the Ehfa_Layer.
2175 245 : char chBandType = '\0';
2176 :
2177 245 : if (eDataType == EPT_u1)
2178 2 : chBandType = '1';
2179 243 : else if (eDataType == EPT_u2)
2180 2 : chBandType = '2';
2181 241 : else if (eDataType == EPT_u4)
2182 2 : chBandType = '4';
2183 239 : else if (eDataType == EPT_u8)
2184 151 : chBandType = 'c';
2185 88 : else if (eDataType == EPT_s8)
2186 4 : chBandType = 'C';
2187 84 : else if (eDataType == EPT_u16)
2188 13 : chBandType = 's';
2189 71 : else if (eDataType == EPT_s16)
2190 9 : chBandType = 'S';
2191 62 : else if (eDataType == EPT_u32)
2192 : // For some reason erdas imagine expects an L for unsigned 32 bit ints
2193 : // otherwise it gives strange "out of memory errors".
2194 10 : chBandType = 'L';
2195 52 : else if (eDataType == EPT_s32)
2196 11 : chBandType = 'L';
2197 41 : else if (eDataType == EPT_f32)
2198 11 : chBandType = 'f';
2199 30 : else if (eDataType == EPT_f64)
2200 12 : chBandType = 'd';
2201 18 : else if (eDataType == EPT_c64)
2202 9 : chBandType = 'm';
2203 9 : else if (eDataType == EPT_c128)
2204 9 : chBandType = 'M';
2205 : else
2206 : {
2207 0 : CPLAssert(false);
2208 : chBandType = 'c';
2209 : }
2210 :
2211 : // The first value in the entry below gives the number of pixels
2212 : // within a block.
2213 245 : char szLDict[128] = {};
2214 245 : snprintf(szLDict, sizeof(szLDict), "{%d:%cdata,}RasterDMS,.",
2215 : nBlockSize * nBlockSize, chBandType);
2216 :
2217 : HFAEntry *poEhfa_Layer =
2218 245 : HFAEntry::New(psInfo, "Ehfa_Layer", "Ehfa_Layer", poEimg_Layer);
2219 245 : poEhfa_Layer->MakeData();
2220 245 : poEhfa_Layer->SetPosition();
2221 : const auto nLDict =
2222 245 : HFAAllocateSpace(psInfo, static_cast<GUInt32>(strlen(szLDict) + 1));
2223 245 : if (nLDict > static_cast<unsigned>(INT_MAX))
2224 0 : return FALSE;
2225 :
2226 245 : poEhfa_Layer->SetStringField("type", "raster");
2227 245 : poEhfa_Layer->SetIntField("dictionaryPtr", static_cast<int>(nLDict));
2228 :
2229 245 : bool bRet = VSIFSeekL(psInfo->fp, nLDict, SEEK_SET) >= 0;
2230 245 : bRet &= VSIFWriteL((void *)szLDict, strlen(szLDict) + 1, 1, psInfo->fp) > 0;
2231 :
2232 245 : return bRet;
2233 : }
2234 :
2235 : /************************************************************************/
2236 : /* HFACreate() */
2237 : /************************************************************************/
2238 :
2239 204 : HFAHandle HFACreate(const char *pszFilename, int nXSize, int nYSize, int nBands,
2240 : EPTType eDataType, char **papszOptions)
2241 :
2242 : {
2243 204 : if (nXSize == 0 || nYSize == 0)
2244 : {
2245 1 : CPLError(CE_Failure, CPLE_NotSupported,
2246 : "nXSize == 0 || nYSize == 0 not supported");
2247 1 : return nullptr;
2248 : }
2249 203 : int nBlockSize = 64;
2250 203 : const char *pszValue = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
2251 :
2252 203 : if (pszValue != nullptr)
2253 : {
2254 1 : nBlockSize = atoi(pszValue);
2255 : // Check for sane values.
2256 2 : if (nBlockSize == 0 ||
2257 1 : (((nBlockSize < 32) || (nBlockSize > 2048)) &&
2258 0 : !CPLTestBool(CPLGetConfigOption("FORCE_BLOCKSIZE", "NO"))))
2259 : {
2260 0 : if (nBlockSize != 0)
2261 0 : CPLError(CE_Warning, CPLE_AppDefined, "Forcing BLOCKSIZE to %d",
2262 : 64);
2263 0 : nBlockSize = 64;
2264 : }
2265 : }
2266 203 : bool bCreateLargeRaster = CPLFetchBool(papszOptions, "USE_SPILL", false);
2267 402 : bool bCreateCompressed = CPLFetchBool(papszOptions, "COMPRESS", false) ||
2268 199 : CPLFetchBool(papszOptions, "COMPRESSED", false);
2269 203 : const bool bCreateAux = CPLFetchBool(papszOptions, "AUX", false);
2270 :
2271 203 : char *pszFullFilename = nullptr;
2272 203 : char *pszRawFilename = nullptr;
2273 :
2274 : // Work out some details about the tiling scheme.
2275 203 : const int nBlocksPerRow = DIV_ROUND_UP(nXSize, nBlockSize);
2276 203 : const int nBlocksPerColumn = DIV_ROUND_UP(nYSize, nBlockSize);
2277 203 : if (nBlocksPerRow > INT_MAX / nBlocksPerColumn)
2278 : {
2279 0 : CPLError(CE_Failure, CPLE_NotSupported, "Too many blocks");
2280 0 : return nullptr;
2281 : }
2282 203 : const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
2283 : const GInt64 nBytesPerBlock64 =
2284 406 : (static_cast<GInt64>(nBlockSize) * nBlockSize *
2285 203 : HFAGetDataTypeBits(eDataType) +
2286 : 7) /
2287 203 : 8;
2288 203 : if (nBytesPerBlock64 > INT_MAX)
2289 : {
2290 0 : CPLError(CE_Failure, CPLE_NotSupported, "Too large block");
2291 0 : return nullptr;
2292 : }
2293 203 : const int nBytesPerBlock = static_cast<int>(nBytesPerBlock64);
2294 :
2295 : // Create the low level structure.
2296 203 : HFAHandle psInfo = HFACreateLL(pszFilename);
2297 203 : if (psInfo == nullptr)
2298 7 : return nullptr;
2299 :
2300 : // Create the DependentFile node if requested.
2301 : const char *pszDependentFile =
2302 196 : CSLFetchNameValue(papszOptions, "DEPENDENT_FILE");
2303 :
2304 196 : if (pszDependentFile != nullptr)
2305 : {
2306 5 : HFAEntry *poDF = HFAEntry::New(psInfo, "DependentFile",
2307 : "Eimg_DependentFile", psInfo->poRoot);
2308 :
2309 5 : poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
2310 5 : poDF->SetPosition();
2311 5 : poDF->SetStringField("dependent.string", pszDependentFile);
2312 : }
2313 :
2314 196 : CPLDebug("HFACreate",
2315 : "Blocks per row %d, blocks per column %d, "
2316 : "total number of blocks %d, bytes per block %d.",
2317 : nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock);
2318 :
2319 : // Check whether we should create external large file with
2320 : // image. We create a spill file if the amount of imagery is
2321 : // close to 2GB. We don't check the amount of auxiliary
2322 : // information, so in theory if there were an awful lot of
2323 : // non-imagery data our approximate size could be smaller than
2324 : // the file will actually we be. We leave room for 10MB of
2325 : // auxiliary data.
2326 : // We can also force spill file creation using option
2327 : // SPILL_FILE=YES.
2328 196 : const double dfApproxSize = static_cast<double>(nBytesPerBlock) *
2329 196 : static_cast<double>(nBlocks) *
2330 196 : static_cast<double>(nBands) +
2331 : 10000000.0;
2332 :
2333 196 : if (dfApproxSize > 2147483648.0 && !bCreateAux)
2334 0 : bCreateLargeRaster = true;
2335 :
2336 : // Erdas Imagine creates this entry even if an external spill file is used.
2337 196 : if (!bCreateAux)
2338 : {
2339 191 : HFAEntry *poImgFormat = HFAEntry::New(
2340 : psInfo, "IMGFormatInfo", "ImgFormatInfo831", psInfo->poRoot);
2341 191 : poImgFormat->MakeData();
2342 191 : if (bCreateLargeRaster)
2343 : {
2344 7 : poImgFormat->SetIntField("spaceUsedForRasterData", 0);
2345 : // Can't be compressed if we are creating a spillfile.
2346 7 : bCreateCompressed = false;
2347 : }
2348 : else
2349 : {
2350 184 : poImgFormat->SetIntField("spaceUsedForRasterData",
2351 184 : nBytesPerBlock * nBlocks * nBands);
2352 : }
2353 : }
2354 :
2355 : // Create external file and write its header.
2356 196 : GIntBig nValidFlagsOffset = 0;
2357 196 : GIntBig nDataOffset = 0;
2358 :
2359 196 : if (bCreateLargeRaster)
2360 : {
2361 7 : if (!HFACreateSpillStack(psInfo, nXSize, nYSize, nBands, nBlockSize,
2362 : eDataType, &nValidFlagsOffset, &nDataOffset))
2363 : {
2364 0 : CPLFree(pszRawFilename);
2365 0 : CPLFree(pszFullFilename);
2366 0 : return nullptr;
2367 : }
2368 : }
2369 :
2370 : // Create each band (layer).
2371 425 : for (int iBand = 0; iBand < nBands; iBand++)
2372 : {
2373 233 : char szName[128] = {};
2374 :
2375 233 : snprintf(szName, sizeof(szName), "Layer_%d", iBand + 1);
2376 :
2377 233 : if (!HFACreateLayer(psInfo, psInfo->poRoot, szName, FALSE, nBlockSize,
2378 : bCreateCompressed, bCreateLargeRaster, bCreateAux,
2379 : nXSize, nYSize, eDataType, papszOptions,
2380 : nValidFlagsOffset, nDataOffset, nBands, iBand))
2381 : {
2382 4 : CPL_IGNORE_RET_VAL(HFAClose(psInfo));
2383 4 : return nullptr;
2384 : }
2385 : }
2386 :
2387 : // Initialize the band information.
2388 192 : HFAParseBandInfo(psInfo);
2389 :
2390 192 : return psInfo;
2391 : }
2392 :
2393 : /************************************************************************/
2394 : /* HFACreateOverview() */
2395 : /* */
2396 : /* Create an overview layer object for a band. */
2397 : /************************************************************************/
2398 :
2399 12 : int HFACreateOverview(HFAHandle hHFA, int nBand, int nOverviewLevel,
2400 : const char *pszResampling)
2401 :
2402 : {
2403 12 : if (nBand < 1 || nBand > hHFA->nBands)
2404 0 : return -1;
2405 :
2406 12 : HFABand *poBand = hHFA->papoBand[nBand - 1];
2407 12 : return poBand->CreateOverview(nOverviewLevel, pszResampling);
2408 : }
2409 :
2410 : /************************************************************************/
2411 : /* HFAGetMetadata() */
2412 : /* */
2413 : /* Read metadata structured in a table called GDAL_MetaData. */
2414 : /************************************************************************/
2415 :
2416 1166 : char **HFAGetMetadata(HFAHandle hHFA, int nBand)
2417 :
2418 : {
2419 1166 : HFAEntry *poTable = nullptr;
2420 :
2421 1166 : if (nBand > 0 && nBand <= hHFA->nBands)
2422 623 : poTable = hHFA->papoBand[nBand - 1]->poNode->GetChild();
2423 543 : else if (nBand == 0)
2424 543 : poTable = hHFA->poRoot->GetChild();
2425 : else
2426 0 : return nullptr;
2427 :
2428 4777 : for (; poTable != nullptr && !EQUAL(poTable->GetName(), "GDAL_MetaData");
2429 3611 : poTable = poTable->GetNext())
2430 : {
2431 : }
2432 :
2433 1166 : if (poTable == nullptr || !EQUAL(poTable->GetType(), "Edsc_Table"))
2434 1071 : return nullptr;
2435 :
2436 95 : if (poTable->GetIntField("numRows") != 1)
2437 : {
2438 0 : CPLDebug("HFADataset", "GDAL_MetaData.numRows = %d, expected 1!",
2439 : poTable->GetIntField("numRows"));
2440 0 : return nullptr;
2441 : }
2442 :
2443 : // Loop over each column. Each column will be one metadata
2444 : // entry, with the title being the key, and the row value being
2445 : // the value. There is only ever one row in GDAL_MetaData tables.
2446 95 : char **papszMD = nullptr;
2447 :
2448 286 : for (HFAEntry *poColumn = poTable->GetChild(); poColumn != nullptr;
2449 191 : poColumn = poColumn->GetNext())
2450 : {
2451 : // Skip the #Bin_Function# entry.
2452 191 : if (STARTS_WITH_CI(poColumn->GetName(), "#"))
2453 95 : continue;
2454 :
2455 96 : const char *pszValue = poColumn->GetStringField("dataType");
2456 96 : if (pszValue == nullptr || !EQUAL(pszValue, "string"))
2457 0 : continue;
2458 :
2459 96 : const int columnDataPtr = poColumn->GetIntField("columnDataPtr");
2460 96 : if (columnDataPtr <= 0)
2461 0 : continue;
2462 :
2463 : // Read up to nMaxNumChars bytes from the indicated location.
2464 : // allocate required space temporarily nMaxNumChars should have been
2465 : // set by GDAL originally so we should trust it, but who knows.
2466 96 : const int nMaxNumChars = poColumn->GetIntField("maxNumChars");
2467 :
2468 96 : if (nMaxNumChars <= 0)
2469 : {
2470 0 : papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), "");
2471 : }
2472 : else
2473 : {
2474 : char *pszMDValue =
2475 96 : static_cast<char *>(VSI_MALLOC_VERBOSE(nMaxNumChars));
2476 96 : if (pszMDValue == nullptr)
2477 : {
2478 0 : continue;
2479 : }
2480 :
2481 96 : if (VSIFSeekL(hHFA->fp, static_cast<vsi_l_offset>(columnDataPtr),
2482 96 : SEEK_SET) != 0)
2483 : {
2484 0 : CPLFree(pszMDValue);
2485 0 : continue;
2486 : }
2487 :
2488 : const int nMDBytes = static_cast<int>(
2489 96 : VSIFReadL(pszMDValue, 1, nMaxNumChars, hHFA->fp));
2490 96 : if (nMDBytes == 0)
2491 : {
2492 0 : CPLFree(pszMDValue);
2493 0 : continue;
2494 : }
2495 :
2496 96 : pszMDValue[nMaxNumChars - 1] = '\0';
2497 :
2498 96 : papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), pszMDValue);
2499 96 : CPLFree(pszMDValue);
2500 : }
2501 : }
2502 :
2503 95 : return papszMD;
2504 : }
2505 :
2506 : /************************************************************************/
2507 : /* HFASetGDALMetadata() */
2508 : /* */
2509 : /* This function is used to set metadata in a table called */
2510 : /* GDAL_MetaData. It is called by HFASetMetadata() for all */
2511 : /* metadata items that aren't some specific supported */
2512 : /* information (like histogram or stats info). */
2513 : /************************************************************************/
2514 :
2515 42 : static CPLErr HFASetGDALMetadata(HFAHandle hHFA, int nBand, char **papszMD)
2516 :
2517 : {
2518 42 : if (papszMD == nullptr)
2519 0 : return CE_None;
2520 :
2521 42 : HFAEntry *poNode = nullptr;
2522 :
2523 42 : if (nBand > 0 && nBand <= hHFA->nBands)
2524 6 : poNode = hHFA->papoBand[nBand - 1]->poNode;
2525 36 : else if (nBand == 0)
2526 36 : poNode = hHFA->poRoot;
2527 : else
2528 0 : return CE_Failure;
2529 :
2530 : // Create the Descriptor table.
2531 : // Check we have no table with this name already.
2532 42 : HFAEntry *poEdsc_Table = poNode->GetNamedChild("GDAL_MetaData");
2533 :
2534 43 : if (poEdsc_Table == nullptr ||
2535 1 : !EQUAL(poEdsc_Table->GetType(), "Edsc_Table"))
2536 : poEdsc_Table =
2537 41 : HFAEntry::New(hHFA, "GDAL_MetaData", "Edsc_Table", poNode);
2538 :
2539 42 : poEdsc_Table->SetIntField("numrows", 1);
2540 :
2541 : // Create the Binning function node. Do we really need this though?
2542 : // Check it doesn't exist already.
2543 : HFAEntry *poEdsc_BinFunction =
2544 42 : poEdsc_Table->GetNamedChild("#Bin_Function#");
2545 :
2546 43 : if (poEdsc_BinFunction == nullptr ||
2547 1 : !EQUAL(poEdsc_BinFunction->GetType(), "Edsc_BinFunction"))
2548 41 : poEdsc_BinFunction = HFAEntry::New(hHFA, "#Bin_Function#",
2549 : "Edsc_BinFunction", poEdsc_Table);
2550 :
2551 : // Because of the BaseData we have to hardcode the size.
2552 42 : poEdsc_BinFunction->MakeData(30);
2553 :
2554 42 : poEdsc_BinFunction->SetIntField("numBins", 1);
2555 42 : poEdsc_BinFunction->SetStringField("binFunction", "direct");
2556 42 : poEdsc_BinFunction->SetDoubleField("minLimit", 0.0);
2557 42 : poEdsc_BinFunction->SetDoubleField("maxLimit", 0.0);
2558 :
2559 : // Process each metadata item as a separate column.
2560 42 : bool bRet = true;
2561 85 : for (int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++)
2562 : {
2563 43 : char *pszKey = nullptr;
2564 43 : const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
2565 43 : if (pszValue == nullptr)
2566 0 : continue;
2567 :
2568 : // Create the Edsc_Column.
2569 : // Check it doesn't exist already.
2570 43 : HFAEntry *poEdsc_Column = poEdsc_Table->GetNamedChild(pszKey);
2571 :
2572 44 : if (poEdsc_Column == nullptr ||
2573 1 : !EQUAL(poEdsc_Column->GetType(), "Edsc_Column"))
2574 : poEdsc_Column =
2575 42 : HFAEntry::New(hHFA, pszKey, "Edsc_Column", poEdsc_Table);
2576 :
2577 43 : CPLFree(pszKey);
2578 :
2579 43 : poEdsc_Column->SetIntField("numRows", 1);
2580 43 : poEdsc_Column->SetStringField("dataType", "string");
2581 43 : poEdsc_Column->SetIntField("maxNumChars",
2582 43 : static_cast<GUInt32>(strlen(pszValue) + 1));
2583 :
2584 : // Write the data out.
2585 : const auto nOffset =
2586 43 : HFAAllocateSpace(hHFA, static_cast<GUInt32>(strlen(pszValue) + 1));
2587 43 : if (nOffset > static_cast<unsigned>(INT_MAX))
2588 : {
2589 0 : return CE_Failure;
2590 : }
2591 :
2592 43 : poEdsc_Column->SetIntField("columnDataPtr", static_cast<int>(nOffset));
2593 :
2594 43 : bRet &= VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) >= 0;
2595 43 : bRet &=
2596 43 : VSIFWriteL((void *)pszValue, strlen(pszValue) + 1, 1, hHFA->fp) > 0;
2597 : }
2598 :
2599 42 : return bRet ? CE_None : CE_Failure;
2600 : }
2601 :
2602 : /************************************************************************/
2603 : /* HFASetMetadata() */
2604 : /************************************************************************/
2605 :
2606 65 : CPLErr HFASetMetadata(HFAHandle hHFA, int nBand, CSLConstList papszMD)
2607 :
2608 : {
2609 65 : char **papszGDALMD = nullptr;
2610 :
2611 65 : if (CSLCount(papszMD) == 0)
2612 0 : return CE_None;
2613 :
2614 65 : HFAEntry *poNode = nullptr;
2615 :
2616 65 : if (nBand > 0 && nBand <= hHFA->nBands)
2617 15 : poNode = hHFA->papoBand[nBand - 1]->poNode;
2618 50 : else if (nBand == 0)
2619 50 : poNode = hHFA->poRoot;
2620 : else
2621 0 : return CE_Failure;
2622 : #ifdef DEBUG
2623 : // To please Clang Static Analyzer (CSA).
2624 65 : if (poNode == nullptr)
2625 : {
2626 0 : CPLAssert(false);
2627 : return CE_Failure;
2628 : }
2629 : #endif
2630 :
2631 : // Check if the Metadata is an "known" entity which should be
2632 : // stored in a better place.
2633 65 : char *pszBinValues = nullptr;
2634 65 : bool bCreatedHistogramParameters = false;
2635 65 : bool bCreatedStatistics = false;
2636 65 : const char *const *pszAuxMetaData = GetHFAAuxMetaDataList();
2637 : // Check each metadata item.
2638 229 : for (int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++)
2639 : {
2640 164 : char *pszKey = nullptr;
2641 164 : const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
2642 164 : if (pszValue == nullptr)
2643 28 : continue;
2644 :
2645 : // Know look if its known.
2646 150 : int i = 0; // Used after for.
2647 1408 : for (; pszAuxMetaData[i] != nullptr; i += 4)
2648 : {
2649 1358 : if (EQUALN(pszAuxMetaData[i + 2], pszKey, strlen(pszKey)))
2650 100 : break;
2651 : }
2652 150 : if (pszAuxMetaData[i] != nullptr)
2653 : {
2654 : // Found one, get the right entry.
2655 100 : HFAEntry *poEntry = nullptr;
2656 :
2657 100 : if (strlen(pszAuxMetaData[i]) > 0)
2658 90 : poEntry = poNode->GetNamedChild(pszAuxMetaData[i]);
2659 : else
2660 10 : poEntry = poNode;
2661 :
2662 100 : if (poEntry == nullptr && strlen(pszAuxMetaData[i + 3]) > 0)
2663 : {
2664 : // Child does not yet exist --> create it,
2665 28 : poEntry = HFAEntry::New(hHFA, pszAuxMetaData[i],
2666 14 : pszAuxMetaData[i + 3], poNode);
2667 :
2668 14 : if (STARTS_WITH_CI(pszAuxMetaData[i], "Statistics"))
2669 8 : bCreatedStatistics = true;
2670 :
2671 14 : if (STARTS_WITH_CI(pszAuxMetaData[i], "HistogramParameters"))
2672 : {
2673 : // A bit nasty. Need to set the string field for the object
2674 : // first because the SetStringField sets the count for the
2675 : // object BinFunction to the length of the string.
2676 6 : poEntry->MakeData(70);
2677 6 : poEntry->SetStringField("BinFunction.binFunctionType",
2678 : "direct");
2679 :
2680 6 : bCreatedHistogramParameters = true;
2681 : }
2682 : }
2683 100 : if (poEntry == nullptr)
2684 : {
2685 14 : CPLFree(pszKey);
2686 14 : continue;
2687 : }
2688 :
2689 86 : const char *pszFieldName = pszAuxMetaData[i + 1] + 1;
2690 86 : switch (pszAuxMetaData[i + 1][0])
2691 : {
2692 67 : case 'd':
2693 : {
2694 67 : double dfValue = CPLAtof(pszValue);
2695 67 : poEntry->SetDoubleField(pszFieldName, dfValue);
2696 : }
2697 67 : break;
2698 9 : case 'i':
2699 : case 'l':
2700 : {
2701 9 : int nValue = atoi(pszValue);
2702 9 : poEntry->SetIntField(pszFieldName, nValue);
2703 : }
2704 9 : break;
2705 10 : case 's':
2706 : case 'e':
2707 : {
2708 10 : poEntry->SetStringField(pszFieldName, pszValue);
2709 : }
2710 10 : break;
2711 0 : default:
2712 0 : CPLAssert(false);
2713 : }
2714 : }
2715 50 : else if (STARTS_WITH_CI(pszKey, "STATISTICS_HISTOBINVALUES"))
2716 : {
2717 7 : CPLFree(pszBinValues);
2718 7 : pszBinValues = CPLStrdup(pszValue);
2719 : }
2720 : else
2721 : {
2722 43 : papszGDALMD = CSLAddString(papszGDALMD, papszMD[iColumn]);
2723 : }
2724 :
2725 136 : CPLFree(pszKey);
2726 : }
2727 :
2728 : // Special case to write out the histogram.
2729 65 : bool bRet = true;
2730 65 : if (pszBinValues != nullptr)
2731 : {
2732 7 : HFAEntry *poEntry = poNode->GetNamedChild("HistogramParameters");
2733 7 : if (poEntry != nullptr && bCreatedHistogramParameters)
2734 : {
2735 : // If this node exists we have added Histogram data -- complete with
2736 : // some defaults.
2737 6 : poEntry->SetIntField("SkipFactorX", 1);
2738 6 : poEntry->SetIntField("SkipFactorY", 1);
2739 :
2740 6 : const int nNumBins = poEntry->GetIntField("BinFunction.numBins");
2741 : const double dMinLimit =
2742 6 : poEntry->GetDoubleField("BinFunction.minLimit");
2743 : const double dMaxLimit =
2744 6 : poEntry->GetDoubleField("BinFunction.maxLimit");
2745 :
2746 : // Fill the descriptor table - check it isn't there already.
2747 6 : poEntry = poNode->GetNamedChild("Descriptor_Table");
2748 6 : if (poEntry == nullptr || !EQUAL(poEntry->GetType(), "Edsc_Table"))
2749 1 : poEntry = HFAEntry::New(hHFA, "Descriptor_Table", "Edsc_Table",
2750 : poNode);
2751 :
2752 6 : poEntry->SetIntField("numRows", nNumBins);
2753 :
2754 : // Bin function.
2755 6 : HFAEntry *poBinFunc = poEntry->GetNamedChild("#Bin_Function#");
2756 10 : if (poBinFunc == nullptr ||
2757 4 : !EQUAL(poBinFunc->GetType(), "Edsc_BinFunction"))
2758 2 : poBinFunc = HFAEntry::New(hHFA, "#Bin_Function#",
2759 : "Edsc_BinFunction", poEntry);
2760 :
2761 6 : poBinFunc->MakeData(30);
2762 6 : poBinFunc->SetIntField("numBins", nNumBins);
2763 6 : poBinFunc->SetDoubleField("minLimit", dMinLimit);
2764 6 : poBinFunc->SetDoubleField("maxLimit", dMaxLimit);
2765 : // Direct for thematic layers, linear otherwise.
2766 6 : if (STARTS_WITH_CI(poNode->GetStringField("layerType"), "thematic"))
2767 1 : poBinFunc->SetStringField("binFunctionType", "direct");
2768 : else
2769 5 : poBinFunc->SetStringField("binFunctionType", "linear");
2770 :
2771 : // We need a child named histogram.
2772 6 : HFAEntry *poHisto = poEntry->GetNamedChild("Histogram");
2773 6 : if (poHisto == nullptr || !EQUAL(poHisto->GetType(), "Edsc_Column"))
2774 : poHisto =
2775 1 : HFAEntry::New(hHFA, "Histogram", "Edsc_Column", poEntry);
2776 :
2777 6 : poHisto->SetIntField("numRows", nNumBins);
2778 : // Allocate space for the bin values.
2779 6 : const auto nOffset64 = HFAAllocateSpace(hHFA, nNumBins * 8);
2780 6 : if (nOffset64 >= static_cast<unsigned>(INT_MAX))
2781 : {
2782 0 : return CE_Failure;
2783 : }
2784 6 : const int nOffset = static_cast<int>(nOffset64);
2785 6 : poHisto->SetIntField("columnDataPtr", nOffset);
2786 6 : poHisto->SetStringField("dataType", "real");
2787 6 : poHisto->SetIntField("maxNumChars", 0);
2788 : // Write out histogram data.
2789 6 : char *pszWork = pszBinValues;
2790 1531 : for (int nBin = 0; nBin < nNumBins; ++nBin)
2791 : {
2792 1525 : char *pszEnd = strchr(pszWork, '|');
2793 1525 : if (pszEnd != nullptr)
2794 : {
2795 1525 : *pszEnd = 0;
2796 1525 : bRet &=
2797 3050 : VSIFSeekL(hHFA->fp,
2798 1525 : nOffset + 8 * static_cast<vsi_l_offset>(nBin),
2799 1525 : SEEK_SET) >= 0;
2800 1525 : double nValue = CPLAtof(pszWork);
2801 : HFAStandard(8, &nValue);
2802 :
2803 1525 : bRet &= VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
2804 1525 : pszWork = pszEnd + 1;
2805 : }
2806 6 : }
2807 : }
2808 1 : else if (poEntry != nullptr)
2809 : {
2810 : // In this case, there are HistogramParameters present, but we did
2811 : // not create them. However, we might be modifying them, in the case
2812 : // where the data has changed and the histogram counts need to be
2813 : // updated. It could be worse than that, but that is all we are
2814 : // going to cope with for now. We are assuming that we did not
2815 : // change any of the other stuff, like skip factors and so
2816 : // forth. The main need for this case is for programs (such as
2817 : // Imagine itself) which will happily modify the pixel values
2818 : // without re-calculating the histogram counts.
2819 1 : int nNumBins = poEntry->GetIntField("BinFunction.numBins");
2820 : HFAEntry *poEntryDescrTbl =
2821 1 : poNode->GetNamedChild("Descriptor_Table");
2822 1 : HFAEntry *poHisto = nullptr;
2823 1 : if (poEntryDescrTbl != nullptr)
2824 : {
2825 1 : poHisto = poEntryDescrTbl->GetNamedChild("Histogram");
2826 : }
2827 1 : if (poHisto != nullptr)
2828 : {
2829 1 : int nOffset = poHisto->GetIntField("columnDataPtr");
2830 : // Write out histogram data.
2831 1 : char *pszWork = pszBinValues;
2832 :
2833 : // Check whether histogram counts were written as int or double
2834 1 : bool bCountIsInt = true;
2835 1 : const char *pszDataType = poHisto->GetStringField("dataType");
2836 1 : if (STARTS_WITH_CI(pszDataType, "real"))
2837 : {
2838 1 : bCountIsInt = false;
2839 : }
2840 257 : for (int nBin = 0; nBin < nNumBins; ++nBin)
2841 : {
2842 256 : char *pszEnd = strchr(pszWork, '|');
2843 256 : if (pszEnd != nullptr)
2844 : {
2845 256 : *pszEnd = 0;
2846 256 : if (bCountIsInt)
2847 : {
2848 : // Histogram counts were written as ints, so
2849 : // re-write them the same way.
2850 0 : bRet &= VSIFSeekL(hHFA->fp, nOffset + 4 * nBin,
2851 0 : SEEK_SET) >= 0;
2852 0 : int nValue = atoi(pszWork);
2853 : HFAStandard(4, &nValue);
2854 0 : bRet &=
2855 0 : VSIFWriteL((void *)&nValue, 4, 1, hHFA->fp) > 0;
2856 : }
2857 : else
2858 : {
2859 : // Histogram were written as doubles, as is now the
2860 : // default behavior.
2861 256 : bRet &= VSIFSeekL(hHFA->fp, nOffset + 8 * nBin,
2862 256 : SEEK_SET) >= 0;
2863 256 : double nValue = CPLAtof(pszWork);
2864 : HFAStandard(8, &nValue);
2865 256 : bRet &=
2866 256 : VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
2867 : }
2868 256 : pszWork = pszEnd + 1;
2869 : }
2870 : }
2871 : }
2872 : }
2873 7 : CPLFree(pszBinValues);
2874 : }
2875 :
2876 : // If we created a statistics node then try to create a
2877 : // StatisticsParameters node too.
2878 65 : if (bCreatedStatistics)
2879 : {
2880 : HFAEntry *poEntry =
2881 8 : HFAEntry::New(hHFA, "StatisticsParameters",
2882 : "Eimg_StatisticsParameters830", poNode);
2883 :
2884 8 : poEntry->MakeData(70);
2885 : // poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
2886 :
2887 8 : poEntry->SetIntField("SkipFactorX", 1);
2888 8 : poEntry->SetIntField("SkipFactorY", 1);
2889 : }
2890 :
2891 : // Write out metadata items without a special place.
2892 65 : if (bRet && CSLCount(papszGDALMD) != 0)
2893 : {
2894 42 : CPLErr eErr = HFASetGDALMetadata(hHFA, nBand, papszGDALMD);
2895 :
2896 42 : CSLDestroy(papszGDALMD);
2897 42 : return eErr;
2898 : }
2899 : else
2900 : {
2901 23 : CSLDestroy(papszGDALMD);
2902 23 : return CE_Failure;
2903 : }
2904 : }
2905 :
2906 : /************************************************************************/
2907 : /* HFAGetIGEFilename() */
2908 : /* */
2909 : /* Returns the .ige filename if one is associated with this */
2910 : /* object. For files not newly created we need to scan the */
2911 : /* bands for spill files. Presumably there will only be one. */
2912 : /* */
2913 : /* NOTE: Returns full path, not just the filename portion. */
2914 : /************************************************************************/
2915 :
2916 115 : std::string HFAGetIGEFilename(HFAHandle hHFA)
2917 :
2918 : {
2919 115 : if (hHFA->pszIGEFilename == nullptr)
2920 : {
2921 : std::vector<HFAEntry *> apoDMSList =
2922 112 : hHFA->poRoot->FindChildren(nullptr, "ImgExternalRaster");
2923 :
2924 112 : HFAEntry *poDMS = apoDMSList.empty() ? nullptr : apoDMSList[0];
2925 :
2926 : // Get the IGE filename from if we have an ExternalRasterDMS.
2927 112 : if (poDMS)
2928 : {
2929 : const char *pszRawFilename =
2930 32 : poDMS->GetStringField("fileName.string");
2931 :
2932 32 : if (pszRawFilename != nullptr)
2933 : {
2934 32 : if (CPLHasPathTraversal(pszRawFilename))
2935 : {
2936 0 : CPLError(CE_Failure, CPLE_AppDefined,
2937 : "Path traversal detected in %s", pszRawFilename);
2938 0 : return std::string();
2939 : }
2940 :
2941 : VSIStatBufL sStatBuf;
2942 : std::string osFullFilename =
2943 64 : CPLFormFilenameSafe(hHFA->pszPath, pszRawFilename, nullptr);
2944 :
2945 32 : if (VSIStatL(osFullFilename.c_str(), &sStatBuf) != 0)
2946 : {
2947 : const CPLString osExtension =
2948 0 : CPLGetExtensionSafe(pszRawFilename);
2949 : const CPLString osBasename =
2950 0 : CPLGetBasenameSafe(hHFA->pszFilename);
2951 0 : osFullFilename = CPLFormFilenameSafe(
2952 0 : hHFA->pszPath, osBasename, osExtension);
2953 :
2954 0 : if (VSIStatL(osFullFilename.c_str(), &sStatBuf) == 0)
2955 0 : hHFA->pszIGEFilename =
2956 0 : CPLStrdup(CPLFormFilenameSafe(nullptr, osBasename,
2957 : osExtension)
2958 : .c_str());
2959 : else
2960 0 : hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
2961 : }
2962 : else
2963 : {
2964 32 : hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
2965 : }
2966 : }
2967 : }
2968 : }
2969 :
2970 : // Return the full filename.
2971 115 : if (hHFA->pszIGEFilename)
2972 35 : return CPLFormFilenameSafe(hHFA->pszPath, hHFA->pszIGEFilename,
2973 35 : nullptr);
2974 :
2975 80 : return std::string();
2976 : }
2977 :
2978 : /************************************************************************/
2979 : /* HFACreateSpillStack() */
2980 : /* */
2981 : /* Create a new stack of raster layers in the spill (.ige) */
2982 : /* file. Create the spill file if it didn't exist before. */
2983 : /************************************************************************/
2984 :
2985 7 : bool HFACreateSpillStack(HFAInfo_t *psInfo, int nXSize, int nYSize, int nLayers,
2986 : int nBlockSize, EPTType eDataType,
2987 : GIntBig *pnValidFlagsOffset, GIntBig *pnDataOffset)
2988 :
2989 : {
2990 : // Form .ige filename.
2991 7 : if (nBlockSize <= 0)
2992 : {
2993 0 : CPLError(CE_Failure, CPLE_IllegalArg,
2994 : "HFACreateSpillStack: nBlockXSize < 0");
2995 0 : return false;
2996 : }
2997 :
2998 7 : if (psInfo->pszIGEFilename == nullptr)
2999 : {
3000 14 : const auto osExt = CPLGetExtensionSafe(psInfo->pszFilename);
3001 7 : if (EQUAL(osExt.c_str(), "rrd"))
3002 0 : psInfo->pszIGEFilename = CPLStrdup(
3003 0 : CPLResetExtensionSafe(psInfo->pszFilename, "rde").c_str());
3004 7 : else if (EQUAL(osExt.c_str(), "aux"))
3005 0 : psInfo->pszIGEFilename = CPLStrdup(
3006 0 : CPLResetExtensionSafe(psInfo->pszFilename, "axe").c_str());
3007 : else
3008 7 : psInfo->pszIGEFilename = CPLStrdup(
3009 14 : CPLResetExtensionSafe(psInfo->pszFilename, "ige").c_str());
3010 : }
3011 :
3012 7 : char *pszFullFilename = CPLStrdup(
3013 14 : CPLFormFilenameSafe(psInfo->pszPath, psInfo->pszIGEFilename, nullptr)
3014 : .c_str());
3015 :
3016 : // Try and open it. If we fail, create it and write the magic header.
3017 : static const char *const pszMagick = "ERDAS_IMG_EXTERNAL_RASTER";
3018 :
3019 7 : bool bRet = true;
3020 7 : VSILFILE *fpVSIL = VSIFOpenL(pszFullFilename, "r+b");
3021 7 : if (fpVSIL == nullptr)
3022 : {
3023 7 : fpVSIL = VSIFOpenL(pszFullFilename, "w+");
3024 7 : if (fpVSIL == nullptr)
3025 : {
3026 0 : CPLError(CE_Failure, CPLE_OpenFailed,
3027 : "Failed to create spill file %s.\n%s",
3028 0 : psInfo->pszIGEFilename, VSIStrerror(errno));
3029 0 : return false;
3030 : }
3031 :
3032 7 : bRet &=
3033 7 : VSIFWriteL((void *)pszMagick, strlen(pszMagick) + 1, 1, fpVSIL) > 0;
3034 : }
3035 :
3036 7 : CPLFree(pszFullFilename);
3037 :
3038 : // Work out some details about the tiling scheme.
3039 7 : const int nBlocksPerRow = DIV_ROUND_UP(nXSize, nBlockSize);
3040 7 : const int nBlocksPerColumn = DIV_ROUND_UP(nYSize, nBlockSize);
3041 : // const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
3042 : const int nBytesPerBlock =
3043 7 : (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
3044 :
3045 7 : const int nBytesPerRow = (nBlocksPerRow + 7) / 8;
3046 7 : const int nBlockMapSize = nBytesPerRow * nBlocksPerColumn;
3047 : // const int iFlagsSize = nBlockMapSize + 20;
3048 :
3049 : // Write stack prefix information.
3050 7 : bRet &= VSIFSeekL(fpVSIL, 0, SEEK_END) >= 0;
3051 :
3052 7 : GByte bUnknown = 1;
3053 7 : bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
3054 :
3055 7 : GInt32 nValue32 = nLayers;
3056 : HFAStandard(4, &nValue32);
3057 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3058 7 : nValue32 = nXSize;
3059 : HFAStandard(4, &nValue32);
3060 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3061 7 : nValue32 = nYSize;
3062 : HFAStandard(4, &nValue32);
3063 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3064 7 : nValue32 = nBlockSize;
3065 : HFAStandard(4, &nValue32);
3066 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3067 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3068 7 : bUnknown = 3;
3069 7 : bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
3070 7 : bUnknown = 0;
3071 7 : bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
3072 :
3073 : // Write out ValidFlags section(s).
3074 7 : *pnValidFlagsOffset = VSIFTellL(fpVSIL);
3075 :
3076 : unsigned char *pabyBlockMap =
3077 7 : static_cast<unsigned char *>(VSI_MALLOC_VERBOSE(nBlockMapSize));
3078 7 : if (pabyBlockMap == nullptr)
3079 : {
3080 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpVSIL));
3081 0 : return false;
3082 : }
3083 :
3084 7 : memset(pabyBlockMap, 0xff, nBlockMapSize);
3085 14 : for (int iBand = 0; iBand < nLayers; iBand++)
3086 : {
3087 7 : nValue32 = 1; // Unknown
3088 : HFAStandard(4, &nValue32);
3089 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3090 7 : nValue32 = 0; // Unknown
3091 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3092 7 : nValue32 = nBlocksPerColumn;
3093 : HFAStandard(4, &nValue32);
3094 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3095 7 : nValue32 = nBlocksPerRow;
3096 : HFAStandard(4, &nValue32);
3097 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3098 7 : nValue32 = 0x30000; // Unknown
3099 : HFAStandard(4, &nValue32);
3100 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3101 :
3102 7 : const int iRemainder = nBlocksPerRow % 8;
3103 7 : CPLDebug("HFACreate",
3104 : "Block map size %d, bytes per row %d, remainder %d.",
3105 : nBlockMapSize, nBytesPerRow, iRemainder);
3106 7 : if (iRemainder)
3107 : {
3108 14 : for (int i = nBytesPerRow - 1; i < nBlockMapSize; i += nBytesPerRow)
3109 7 : pabyBlockMap[i] = static_cast<GByte>((1 << iRemainder) - 1);
3110 : }
3111 :
3112 7 : bRet &= VSIFWriteL(pabyBlockMap, nBlockMapSize, 1, fpVSIL) > 0;
3113 : }
3114 7 : CPLFree(pabyBlockMap);
3115 7 : pabyBlockMap = nullptr;
3116 :
3117 : // Extend the file to account for all the imagery space.
3118 7 : const GIntBig nTileDataSize = static_cast<GIntBig>(nBytesPerBlock) *
3119 7 : nBlocksPerRow * nBlocksPerColumn * nLayers;
3120 :
3121 7 : *pnDataOffset = VSIFTellL(fpVSIL);
3122 :
3123 7 : if (!bRet || VSIFTruncateL(fpVSIL, nTileDataSize + *pnDataOffset) != 0)
3124 : {
3125 0 : CPLError(CE_Failure, CPLE_FileIO,
3126 : "Failed to extend %s to full size (" CPL_FRMT_GIB " bytes), "
3127 : "likely out of disk space.\n%s",
3128 0 : psInfo->pszIGEFilename, nTileDataSize + *pnDataOffset,
3129 0 : VSIStrerror(errno));
3130 :
3131 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpVSIL));
3132 0 : return false;
3133 : }
3134 :
3135 7 : if (VSIFCloseL(fpVSIL) != 0)
3136 0 : return false;
3137 :
3138 7 : return true;
3139 : }
3140 :
3141 : /************************************************************************/
3142 : /* HFAReadAndValidatePoly() */
3143 : /************************************************************************/
3144 :
3145 3 : static bool HFAReadAndValidatePoly(HFAEntry *poTarget, const char *pszName,
3146 : Efga_Polynomial *psRetPoly)
3147 :
3148 : {
3149 3 : memset(psRetPoly, 0, sizeof(Efga_Polynomial));
3150 :
3151 6 : CPLString osFldName;
3152 3 : osFldName.Printf("%sorder", pszName);
3153 3 : psRetPoly->order = poTarget->GetIntField(osFldName);
3154 :
3155 3 : if (psRetPoly->order < 1 || psRetPoly->order > 3)
3156 0 : return false;
3157 :
3158 : // Validate that things are in a "well known" form.
3159 3 : osFldName.Printf("%snumdimtransform", pszName);
3160 3 : const int numdimtransform = poTarget->GetIntField(osFldName);
3161 :
3162 3 : osFldName.Printf("%snumdimpolynomial", pszName);
3163 3 : const int numdimpolynomial = poTarget->GetIntField(osFldName);
3164 :
3165 3 : osFldName.Printf("%stermcount", pszName);
3166 3 : const int termcount = poTarget->GetIntField(osFldName);
3167 :
3168 3 : if (numdimtransform != 2 || numdimpolynomial != 2)
3169 0 : return false;
3170 :
3171 3 : if ((psRetPoly->order == 1 && termcount != 3) ||
3172 3 : (psRetPoly->order == 2 && termcount != 6) ||
3173 3 : (psRetPoly->order == 3 && termcount != 10))
3174 0 : return false;
3175 :
3176 : // We don't check the exponent organization for now. Hopefully
3177 : // it is always standard.
3178 :
3179 : // Get coefficients.
3180 43 : for (int i = 0; i < termcount * 2 - 2; i++)
3181 : {
3182 40 : osFldName.Printf("%spolycoefmtx[%d]", pszName, i);
3183 40 : psRetPoly->polycoefmtx[i] = poTarget->GetDoubleField(osFldName);
3184 : }
3185 :
3186 9 : for (int i = 0; i < 2; i++)
3187 : {
3188 6 : osFldName.Printf("%spolycoefvector[%d]", pszName, i);
3189 6 : psRetPoly->polycoefvector[i] = poTarget->GetDoubleField(osFldName);
3190 : }
3191 :
3192 3 : return true;
3193 : }
3194 :
3195 : /************************************************************************/
3196 : /* HFAReadXFormStack() */
3197 : /************************************************************************/
3198 :
3199 293 : int HFAReadXFormStack(HFAHandle hHFA, Efga_Polynomial **ppasPolyListForward,
3200 : Efga_Polynomial **ppasPolyListReverse)
3201 :
3202 : {
3203 293 : if (hHFA->nBands == 0)
3204 0 : return 0;
3205 :
3206 : // Get the HFA node.
3207 : HFAEntry *poXFormHeader =
3208 293 : hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm");
3209 293 : if (poXFormHeader == nullptr)
3210 291 : return 0;
3211 :
3212 : // Loop over children, collecting XForms.
3213 2 : int nStepCount = 0;
3214 2 : *ppasPolyListForward = nullptr;
3215 2 : *ppasPolyListReverse = nullptr;
3216 :
3217 5 : for (HFAEntry *poXForm = poXFormHeader->GetChild(); poXForm != nullptr;
3218 3 : poXForm = poXForm->GetNext())
3219 : {
3220 3 : bool bSuccess = false;
3221 : Efga_Polynomial sForward;
3222 : Efga_Polynomial sReverse;
3223 3 : memset(&sForward, 0, sizeof(sForward));
3224 3 : memset(&sReverse, 0, sizeof(sReverse));
3225 :
3226 3 : if (EQUAL(poXForm->GetType(), "Efga_Polynomial"))
3227 : {
3228 1 : bSuccess = HFAReadAndValidatePoly(poXForm, "", &sForward);
3229 :
3230 1 : if (bSuccess)
3231 : {
3232 : double adfGT[6] = {
3233 1 : sForward.polycoefvector[0], sForward.polycoefmtx[0],
3234 1 : sForward.polycoefmtx[2], sForward.polycoefvector[1],
3235 1 : sForward.polycoefmtx[1], sForward.polycoefmtx[3]};
3236 :
3237 1 : double adfInvGT[6] = {};
3238 1 : bSuccess = HFAInvGeoTransform(adfGT, adfInvGT);
3239 1 : if (!bSuccess)
3240 0 : memset(adfInvGT, 0, sizeof(adfInvGT));
3241 :
3242 1 : sReverse.order = sForward.order;
3243 1 : sReverse.polycoefvector[0] = adfInvGT[0];
3244 1 : sReverse.polycoefmtx[0] = adfInvGT[1];
3245 1 : sReverse.polycoefmtx[2] = adfInvGT[2];
3246 1 : sReverse.polycoefvector[1] = adfInvGT[3];
3247 1 : sReverse.polycoefmtx[1] = adfInvGT[4];
3248 1 : sReverse.polycoefmtx[3] = adfInvGT[5];
3249 : }
3250 : }
3251 2 : else if (EQUAL(poXForm->GetType(), "GM_PolyPair"))
3252 : {
3253 2 : bSuccess = HFAReadAndValidatePoly(poXForm, "forward.", &sForward) &&
3254 1 : HFAReadAndValidatePoly(poXForm, "reverse.", &sReverse);
3255 : }
3256 :
3257 3 : if (bSuccess)
3258 : {
3259 2 : nStepCount++;
3260 4 : *ppasPolyListForward = static_cast<Efga_Polynomial *>(CPLRealloc(
3261 2 : *ppasPolyListForward, sizeof(Efga_Polynomial) * nStepCount));
3262 2 : memcpy(*ppasPolyListForward + nStepCount - 1, &sForward,
3263 : sizeof(sForward));
3264 :
3265 4 : *ppasPolyListReverse = static_cast<Efga_Polynomial *>(CPLRealloc(
3266 2 : *ppasPolyListReverse, sizeof(Efga_Polynomial) * nStepCount));
3267 2 : memcpy(*ppasPolyListReverse + nStepCount - 1, &sReverse,
3268 : sizeof(sReverse));
3269 : }
3270 : }
3271 :
3272 2 : return nStepCount;
3273 : }
3274 :
3275 : /************************************************************************/
3276 : /* HFAEvaluateXFormStack() */
3277 : /************************************************************************/
3278 :
3279 36 : int HFAEvaluateXFormStack(int nStepCount, int bForward,
3280 : Efga_Polynomial *pasPolyList, double *pdfX,
3281 : double *pdfY)
3282 :
3283 : {
3284 108 : for (int iStep = 0; iStep < nStepCount; iStep++)
3285 : {
3286 72 : const Efga_Polynomial *psStep =
3287 0 : bForward ? pasPolyList + iStep
3288 72 : : pasPolyList + nStepCount - iStep - 1;
3289 :
3290 72 : if (psStep->order == 1)
3291 : {
3292 36 : const double dfXOut = psStep->polycoefvector[0] +
3293 36 : psStep->polycoefmtx[0] * *pdfX +
3294 36 : psStep->polycoefmtx[2] * *pdfY;
3295 :
3296 36 : const double dfYOut = psStep->polycoefvector[1] +
3297 36 : psStep->polycoefmtx[1] * *pdfX +
3298 36 : psStep->polycoefmtx[3] * *pdfY;
3299 :
3300 36 : *pdfX = dfXOut;
3301 36 : *pdfY = dfYOut;
3302 : }
3303 36 : else if (psStep->order == 2)
3304 : {
3305 0 : const double dfXOut = psStep->polycoefvector[0] +
3306 0 : psStep->polycoefmtx[0] * *pdfX +
3307 0 : psStep->polycoefmtx[2] * *pdfY +
3308 0 : psStep->polycoefmtx[4] * *pdfX * *pdfX +
3309 0 : psStep->polycoefmtx[6] * *pdfX * *pdfY +
3310 0 : psStep->polycoefmtx[8] * *pdfY * *pdfY;
3311 0 : const double dfYOut = psStep->polycoefvector[1] +
3312 0 : psStep->polycoefmtx[1] * *pdfX +
3313 0 : psStep->polycoefmtx[3] * *pdfY +
3314 0 : psStep->polycoefmtx[5] * *pdfX * *pdfX +
3315 0 : psStep->polycoefmtx[7] * *pdfX * *pdfY +
3316 0 : psStep->polycoefmtx[9] * *pdfY * *pdfY;
3317 :
3318 0 : *pdfX = dfXOut;
3319 0 : *pdfY = dfYOut;
3320 : }
3321 36 : else if (psStep->order == 3)
3322 : {
3323 36 : const double dfXOut =
3324 36 : psStep->polycoefvector[0] + psStep->polycoefmtx[0] * *pdfX +
3325 36 : psStep->polycoefmtx[2] * *pdfY +
3326 36 : psStep->polycoefmtx[4] * *pdfX * *pdfX +
3327 36 : psStep->polycoefmtx[6] * *pdfX * *pdfY +
3328 36 : psStep->polycoefmtx[8] * *pdfY * *pdfY +
3329 36 : psStep->polycoefmtx[10] * *pdfX * *pdfX * *pdfX +
3330 36 : psStep->polycoefmtx[12] * *pdfX * *pdfX * *pdfY +
3331 36 : psStep->polycoefmtx[14] * *pdfX * *pdfY * *pdfY +
3332 36 : psStep->polycoefmtx[16] * *pdfY * *pdfY * *pdfY;
3333 36 : const double dfYOut =
3334 36 : psStep->polycoefvector[1] + psStep->polycoefmtx[1] * *pdfX +
3335 36 : psStep->polycoefmtx[3] * *pdfY +
3336 36 : psStep->polycoefmtx[5] * *pdfX * *pdfX +
3337 36 : psStep->polycoefmtx[7] * *pdfX * *pdfY +
3338 36 : psStep->polycoefmtx[9] * *pdfY * *pdfY +
3339 36 : psStep->polycoefmtx[11] * *pdfX * *pdfX * *pdfX +
3340 36 : psStep->polycoefmtx[13] * *pdfX * *pdfX * *pdfY +
3341 36 : psStep->polycoefmtx[15] * *pdfX * *pdfY * *pdfY +
3342 36 : psStep->polycoefmtx[17] * *pdfY * *pdfY * *pdfY;
3343 :
3344 36 : *pdfX = dfXOut;
3345 36 : *pdfY = dfYOut;
3346 : }
3347 : else
3348 0 : return FALSE;
3349 : }
3350 :
3351 36 : return TRUE;
3352 : }
3353 :
3354 : /************************************************************************/
3355 : /* HFAWriteXFormStack() */
3356 : /************************************************************************/
3357 :
3358 2 : CPLErr HFAWriteXFormStack(HFAHandle hHFA, int nBand, int nXFormCount,
3359 : Efga_Polynomial **ppasPolyListForward,
3360 : Efga_Polynomial **ppasPolyListReverse)
3361 :
3362 : {
3363 2 : if (nXFormCount == 0)
3364 0 : return CE_None;
3365 :
3366 2 : if (ppasPolyListForward[0]->order != 1)
3367 : {
3368 0 : CPLError(
3369 : CE_Failure, CPLE_AppDefined,
3370 : "For now HFAWriteXFormStack() only supports order 1 polynomials");
3371 0 : return CE_Failure;
3372 : }
3373 :
3374 2 : if (nBand < 0 || nBand > hHFA->nBands)
3375 0 : return CE_Failure;
3376 :
3377 : // If no band number is provided, operate on all bands.
3378 2 : if (nBand == 0)
3379 : {
3380 2 : for (nBand = 1; nBand <= hHFA->nBands; nBand++)
3381 : {
3382 : CPLErr eErr =
3383 1 : HFAWriteXFormStack(hHFA, nBand, nXFormCount,
3384 : ppasPolyListForward, ppasPolyListReverse);
3385 1 : if (eErr != CE_None)
3386 0 : return eErr;
3387 : }
3388 :
3389 1 : return CE_None;
3390 : }
3391 :
3392 : // Fetch our band node.
3393 1 : HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
3394 1 : HFAEntry *poXFormHeader = poBandNode->GetNamedChild("MapToPixelXForm");
3395 1 : if (poXFormHeader == nullptr)
3396 : {
3397 1 : poXFormHeader = HFAEntry::New(hHFA, "MapToPixelXForm",
3398 : "Exfr_GenericXFormHeader", poBandNode);
3399 1 : poXFormHeader->MakeData(23);
3400 1 : poXFormHeader->SetPosition();
3401 1 : poXFormHeader->SetStringField("titleList.string", "Affine");
3402 : }
3403 :
3404 : // Loop over XForms.
3405 2 : for (int iXForm = 0; iXForm < nXFormCount; iXForm++)
3406 : {
3407 1 : Efga_Polynomial *psForward = *ppasPolyListForward + iXForm;
3408 2 : CPLString osXFormName;
3409 1 : osXFormName.Printf("XForm%d", iXForm);
3410 :
3411 1 : HFAEntry *poXForm = poXFormHeader->GetNamedChild(osXFormName);
3412 :
3413 1 : if (poXForm == nullptr)
3414 : {
3415 1 : poXForm = HFAEntry::New(hHFA, osXFormName, "Efga_Polynomial",
3416 : poXFormHeader);
3417 1 : poXForm->MakeData(136);
3418 1 : poXForm->SetPosition();
3419 : }
3420 :
3421 1 : poXForm->SetIntField("order", 1);
3422 1 : poXForm->SetIntField("numdimtransform", 2);
3423 1 : poXForm->SetIntField("numdimpolynomial", 2);
3424 1 : poXForm->SetIntField("termcount", 3);
3425 1 : poXForm->SetIntField("exponentlist[0]", 0);
3426 1 : poXForm->SetIntField("exponentlist[1]", 0);
3427 1 : poXForm->SetIntField("exponentlist[2]", 1);
3428 1 : poXForm->SetIntField("exponentlist[3]", 0);
3429 1 : poXForm->SetIntField("exponentlist[4]", 0);
3430 1 : poXForm->SetIntField("exponentlist[5]", 1);
3431 :
3432 1 : poXForm->SetIntField("polycoefmtx[-3]", EPT_f64);
3433 1 : poXForm->SetIntField("polycoefmtx[-2]", 2);
3434 1 : poXForm->SetIntField("polycoefmtx[-1]", 2);
3435 1 : poXForm->SetDoubleField("polycoefmtx[0]", psForward->polycoefmtx[0]);
3436 1 : poXForm->SetDoubleField("polycoefmtx[1]", psForward->polycoefmtx[1]);
3437 1 : poXForm->SetDoubleField("polycoefmtx[2]", psForward->polycoefmtx[2]);
3438 1 : poXForm->SetDoubleField("polycoefmtx[3]", psForward->polycoefmtx[3]);
3439 :
3440 1 : poXForm->SetIntField("polycoefvector[-3]", EPT_f64);
3441 1 : poXForm->SetIntField("polycoefvector[-2]", 1);
3442 1 : poXForm->SetIntField("polycoefvector[-1]", 2);
3443 1 : poXForm->SetDoubleField("polycoefvector[0]",
3444 : psForward->polycoefvector[0]);
3445 1 : poXForm->SetDoubleField("polycoefvector[1]",
3446 : psForward->polycoefvector[1]);
3447 : }
3448 :
3449 1 : return CE_None;
3450 : }
3451 :
3452 : /************************************************************************/
3453 : /* HFAReadCameraModel() */
3454 : /************************************************************************/
3455 :
3456 543 : char **HFAReadCameraModel(HFAHandle hHFA)
3457 :
3458 : {
3459 543 : if (hHFA->nBands == 0)
3460 0 : return nullptr;
3461 :
3462 : // Get the camera model node, and confirm its type.
3463 : HFAEntry *poXForm =
3464 543 : hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
3465 543 : if (poXForm == nullptr)
3466 537 : return nullptr;
3467 :
3468 6 : if (!EQUAL(poXForm->GetType(), "Camera_ModelX"))
3469 5 : return nullptr;
3470 :
3471 : // Convert the values to metadata.
3472 1 : char **papszMD = nullptr;
3473 : static const char *const apszFields[] = {"direction",
3474 : "refType",
3475 : "demsource",
3476 : "PhotoDirection",
3477 : "RotationSystem",
3478 : "demfilename",
3479 : "demzunits",
3480 : "forSrcAffine[0]",
3481 : "forSrcAffine[1]",
3482 : "forSrcAffine[2]",
3483 : "forSrcAffine[3]",
3484 : "forSrcAffine[4]",
3485 : "forSrcAffine[5]",
3486 : "forDstAffine[0]",
3487 : "forDstAffine[1]",
3488 : "forDstAffine[2]",
3489 : "forDstAffine[3]",
3490 : "forDstAffine[4]",
3491 : "forDstAffine[5]",
3492 : "invSrcAffine[0]",
3493 : "invSrcAffine[1]",
3494 : "invSrcAffine[2]",
3495 : "invSrcAffine[3]",
3496 : "invSrcAffine[4]",
3497 : "invSrcAffine[5]",
3498 : "invDstAffine[0]",
3499 : "invDstAffine[1]",
3500 : "invDstAffine[2]",
3501 : "invDstAffine[3]",
3502 : "invDstAffine[4]",
3503 : "invDstAffine[5]",
3504 : "z_mean",
3505 : "lat0",
3506 : "lon0",
3507 : "coeffs[0]",
3508 : "coeffs[1]",
3509 : "coeffs[2]",
3510 : "coeffs[3]",
3511 : "coeffs[4]",
3512 : "coeffs[5]",
3513 : "coeffs[6]",
3514 : "coeffs[7]",
3515 : "coeffs[8]",
3516 : "LensDistortion[0]",
3517 : "LensDistortion[1]",
3518 : "LensDistortion[2]",
3519 : nullptr};
3520 :
3521 1 : const char *pszValue = nullptr;
3522 47 : for (int i = 0; apszFields[i] != nullptr; i++)
3523 : {
3524 46 : pszValue = poXForm->GetStringField(apszFields[i]);
3525 46 : if (pszValue == nullptr)
3526 1 : pszValue = "";
3527 :
3528 46 : papszMD = CSLSetNameValue(papszMD, apszFields[i], pszValue);
3529 : }
3530 :
3531 : // Create a pseudo-entry for the MIFObject with the outputProjection.
3532 : HFAEntry *poProjInfo =
3533 1 : HFAEntry::BuildEntryFromMIFObject(poXForm, "outputProjection");
3534 1 : if (poProjInfo)
3535 : {
3536 : // Fetch the datum.
3537 : Eprj_Datum sDatum;
3538 :
3539 1 : memset(&sDatum, 0, sizeof(sDatum));
3540 :
3541 1 : sDatum.datumname =
3542 1 : (char *)poProjInfo->GetStringField("earthModel.datum.datumname");
3543 :
3544 1 : const int nDatumType = poProjInfo->GetIntField("earthModel.datum.type");
3545 1 : if (nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE)
3546 : {
3547 0 : CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
3548 0 : sDatum.type = EPRJ_DATUM_NONE;
3549 : }
3550 : else
3551 : {
3552 1 : sDatum.type = static_cast<Eprj_DatumType>(nDatumType);
3553 : }
3554 :
3555 8 : for (int i = 0; i < 7; i++)
3556 : {
3557 7 : char szFieldName[60] = {};
3558 :
3559 7 : snprintf(szFieldName, sizeof(szFieldName),
3560 : "earthModel.datum.params[%d]", i);
3561 7 : sDatum.params[i] = poProjInfo->GetDoubleField(szFieldName);
3562 : }
3563 :
3564 1 : sDatum.gridname =
3565 1 : (char *)poProjInfo->GetStringField("earthModel.datum.gridname");
3566 :
3567 : // Fetch the projection parameters.
3568 : Eprj_ProParameters sPro;
3569 :
3570 1 : memset(&sPro, 0, sizeof(sPro));
3571 :
3572 1 : sPro.proType =
3573 1 : (Eprj_ProType)poProjInfo->GetIntField("projectionObject.proType");
3574 1 : sPro.proNumber = poProjInfo->GetIntField("projectionObject.proNumber");
3575 1 : sPro.proExeName =
3576 1 : (char *)poProjInfo->GetStringField("projectionObject.proExeName");
3577 1 : sPro.proName =
3578 1 : (char *)poProjInfo->GetStringField("projectionObject.proName");
3579 1 : sPro.proZone = poProjInfo->GetIntField("projectionObject.proZone");
3580 :
3581 16 : for (int i = 0; i < 15; i++)
3582 : {
3583 15 : char szFieldName[40] = {};
3584 :
3585 15 : snprintf(szFieldName, sizeof(szFieldName),
3586 : "projectionObject.proParams[%d]", i);
3587 15 : sPro.proParams[i] = poProjInfo->GetDoubleField(szFieldName);
3588 : }
3589 :
3590 : // Fetch the spheroid.
3591 1 : sPro.proSpheroid.sphereName = (char *)poProjInfo->GetStringField(
3592 : "earthModel.proSpheroid.sphereName");
3593 1 : sPro.proSpheroid.a =
3594 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.a");
3595 1 : sPro.proSpheroid.b =
3596 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.b");
3597 1 : sPro.proSpheroid.eSquared =
3598 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.eSquared");
3599 1 : sPro.proSpheroid.radius =
3600 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.radius");
3601 :
3602 : // Fetch the projection info.
3603 : // poProjInfo->DumpFieldValues( stdout, "" );
3604 :
3605 2 : auto poSRS = HFAPCSStructToOSR(&sDatum, &sPro, nullptr, nullptr);
3606 :
3607 1 : if (poSRS)
3608 : {
3609 1 : char *pszProjection = nullptr;
3610 1 : if (poSRS->exportToWkt(&pszProjection) == OGRERR_NONE)
3611 : {
3612 : papszMD =
3613 1 : CSLSetNameValue(papszMD, "outputProjection", pszProjection);
3614 : }
3615 1 : CPLFree(pszProjection);
3616 : }
3617 :
3618 1 : delete poProjInfo;
3619 : }
3620 :
3621 : // Fetch the horizontal units.
3622 1 : pszValue = poXForm->GetStringField("outputHorizontalUnits.string");
3623 1 : if (pszValue == nullptr)
3624 0 : pszValue = "";
3625 :
3626 1 : papszMD = CSLSetNameValue(papszMD, "outputHorizontalUnits", pszValue);
3627 :
3628 : // Fetch the elevationinfo.
3629 : HFAEntry *poElevInfo =
3630 1 : HFAEntry::BuildEntryFromMIFObject(poXForm, "outputElevationInfo");
3631 1 : if (poElevInfo)
3632 : {
3633 : // poElevInfo->DumpFieldValues( stdout, "" );
3634 :
3635 1 : if (poElevInfo->GetDataSize() != 0)
3636 : {
3637 : static const char *const apszEFields[] = {
3638 : "verticalDatum.datumname", "verticalDatum.type",
3639 : "elevationUnit", "elevationType", nullptr};
3640 :
3641 5 : for (int i = 0; apszEFields[i] != nullptr; i++)
3642 : {
3643 4 : pszValue = poElevInfo->GetStringField(apszEFields[i]);
3644 4 : if (pszValue == nullptr)
3645 0 : pszValue = "";
3646 :
3647 4 : papszMD = CSLSetNameValue(papszMD, apszEFields[i], pszValue);
3648 : }
3649 : }
3650 :
3651 1 : delete poElevInfo;
3652 : }
3653 :
3654 1 : return papszMD;
3655 : }
3656 :
3657 : /************************************************************************/
3658 : /* HFAReadElevationUnit() */
3659 : /************************************************************************/
3660 :
3661 623 : const char *HFAReadElevationUnit(HFAHandle hHFA, int iBand)
3662 : {
3663 623 : if (hHFA->nBands <= iBand)
3664 0 : return nullptr;
3665 :
3666 623 : HFABand *poBand(hHFA->papoBand[iBand]);
3667 623 : if (poBand == nullptr || poBand->poNode == nullptr)
3668 : {
3669 0 : return nullptr;
3670 : }
3671 623 : HFAEntry *poElevInfo = poBand->poNode->GetNamedChild("Elevation_Info");
3672 623 : if (poElevInfo == nullptr)
3673 : {
3674 620 : return nullptr;
3675 : }
3676 3 : return poElevInfo->GetStringField("elevationUnit");
3677 : }
3678 :
3679 : /************************************************************************/
3680 : /* HFASetGeoTransform() */
3681 : /* */
3682 : /* Set a MapInformation and XForm block. Allows for rotated */
3683 : /* and shared geotransforms. */
3684 : /************************************************************************/
3685 :
3686 1 : CPLErr HFASetGeoTransform(HFAHandle hHFA, const char *pszProName,
3687 : const char *pszUnits, double *padfGeoTransform)
3688 :
3689 : {
3690 : // Write MapInformation.
3691 2 : for (int nBand = 1; nBand <= hHFA->nBands; nBand++)
3692 : {
3693 1 : HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
3694 :
3695 1 : HFAEntry *poMI = poBandNode->GetNamedChild("MapInformation");
3696 1 : if (poMI == nullptr)
3697 : {
3698 1 : poMI = HFAEntry::New(hHFA, "MapInformation", "Eimg_MapInformation",
3699 : poBandNode);
3700 1 : poMI->MakeData(
3701 1 : static_cast<int>(18 + strlen(pszProName) + strlen(pszUnits)));
3702 1 : poMI->SetPosition();
3703 : }
3704 :
3705 1 : poMI->SetStringField("projection.string", pszProName);
3706 1 : poMI->SetStringField("units.string", pszUnits);
3707 : }
3708 :
3709 : // Write XForm.
3710 1 : double adfAdjTransform[6] = {};
3711 :
3712 : // Offset by half pixel.
3713 :
3714 1 : memcpy(adfAdjTransform, padfGeoTransform, sizeof(double) * 6);
3715 1 : adfAdjTransform[0] += adfAdjTransform[1] * 0.5;
3716 1 : adfAdjTransform[0] += adfAdjTransform[2] * 0.5;
3717 1 : adfAdjTransform[3] += adfAdjTransform[4] * 0.5;
3718 1 : adfAdjTransform[3] += adfAdjTransform[5] * 0.5;
3719 :
3720 : // Invert.
3721 1 : double adfRevTransform[6] = {};
3722 1 : if (!HFAInvGeoTransform(adfAdjTransform, adfRevTransform))
3723 0 : memset(adfRevTransform, 0, sizeof(adfRevTransform));
3724 :
3725 : // Assign to polynomial object.
3726 :
3727 : Efga_Polynomial sForward;
3728 1 : memset(&sForward, 0, sizeof(sForward));
3729 1 : Efga_Polynomial *psForward = &sForward;
3730 1 : sForward.order = 1;
3731 1 : sForward.polycoefvector[0] = adfRevTransform[0];
3732 1 : sForward.polycoefmtx[0] = adfRevTransform[1];
3733 1 : sForward.polycoefmtx[1] = adfRevTransform[4];
3734 1 : sForward.polycoefvector[1] = adfRevTransform[3];
3735 1 : sForward.polycoefmtx[2] = adfRevTransform[2];
3736 1 : sForward.polycoefmtx[3] = adfRevTransform[5];
3737 :
3738 1 : Efga_Polynomial sReverse = sForward;
3739 1 : Efga_Polynomial *psReverse = &sReverse;
3740 :
3741 2 : return HFAWriteXFormStack(hHFA, 0, 1, &psForward, &psReverse);
3742 : }
3743 :
3744 : /************************************************************************/
3745 : /* HFARenameReferences() */
3746 : /* */
3747 : /* Rename references in this .img file from the old basename to */
3748 : /* a new basename. This should be passed on to .aux and .rrd */
3749 : /* files and should include references to .aux, .rrd and .ige. */
3750 : /************************************************************************/
3751 :
3752 4 : CPLErr HFARenameReferences(HFAHandle hHFA, const char *pszNewBase,
3753 : const char *pszOldBase)
3754 :
3755 : {
3756 : // Handle RRDNamesList updates.
3757 : std::vector<HFAEntry *> apoNodeList =
3758 4 : hHFA->poRoot->FindChildren("RRDNamesList", nullptr);
3759 :
3760 6 : for (size_t iNode = 0; iNode < apoNodeList.size(); iNode++)
3761 : {
3762 2 : HFAEntry *poRRDNL = apoNodeList[iNode];
3763 4 : std::vector<CPLString> aosNL;
3764 :
3765 : // Collect all the existing names.
3766 2 : const int nNameCount = poRRDNL->GetFieldCount("nameList");
3767 2 : if (nNameCount >= 0)
3768 : {
3769 4 : CPLString osAlgorithm = poRRDNL->GetStringField("algorithm.string");
3770 4 : for (int i = 0; i < nNameCount; i++)
3771 : {
3772 2 : CPLString osFN;
3773 2 : osFN.Printf("nameList[%d].string", i);
3774 2 : aosNL.push_back(poRRDNL->GetStringField(osFN));
3775 : }
3776 :
3777 : // Adjust the names to the new form.
3778 4 : for (int i = 0; i < nNameCount; i++)
3779 : {
3780 2 : if (strncmp(aosNL[i], pszOldBase, strlen(pszOldBase)) == 0)
3781 : {
3782 2 : std::string osNew = pszNewBase;
3783 2 : osNew += aosNL[i].c_str() + strlen(pszOldBase);
3784 2 : aosNL[i] = std::move(osNew);
3785 : }
3786 : }
3787 :
3788 : // Try to make sure the RRDNamesList is big enough to hold the
3789 : // adjusted name list.
3790 2 : if (strlen(pszNewBase) > strlen(pszOldBase))
3791 : {
3792 1 : CPLDebug("HFA", "Growing RRDNamesList to hold new names");
3793 1 : poRRDNL->MakeData(static_cast<int>(
3794 1 : poRRDNL->GetDataSize() +
3795 1 : nNameCount * (strlen(pszNewBase) - strlen(pszOldBase))));
3796 : }
3797 :
3798 : // Initialize the whole thing to zeros for a clean start.
3799 2 : memset(poRRDNL->GetData(), 0, poRRDNL->GetDataSize());
3800 :
3801 : // Write the updates back to the file.
3802 2 : poRRDNL->SetStringField("algorithm.string", osAlgorithm);
3803 4 : for (int i = 0; i < nNameCount; i++)
3804 : {
3805 4 : CPLString osFN;
3806 2 : osFN.Printf("nameList[%d].string", i);
3807 2 : poRRDNL->SetStringField(osFN, aosNL[i]);
3808 : }
3809 : }
3810 : }
3811 :
3812 : // Spill file references.
3813 : apoNodeList =
3814 4 : hHFA->poRoot->FindChildren("ExternalRasterDMS", "ImgExternalRaster");
3815 :
3816 8 : for (size_t iNode = 0; iNode < apoNodeList.size(); iNode++)
3817 : {
3818 4 : HFAEntry *poERDMS = apoNodeList[iNode];
3819 :
3820 4 : if (poERDMS == nullptr)
3821 0 : continue;
3822 :
3823 : // Fetch all existing values.
3824 8 : CPLString osFileName = poERDMS->GetStringField("fileName.string");
3825 :
3826 : GInt32 anValidFlagsOffset[2] = {
3827 4 : poERDMS->GetIntField("layerStackValidFlagsOffset[0]"),
3828 4 : poERDMS->GetIntField("layerStackValidFlagsOffset[1]")};
3829 :
3830 : GInt32 anStackDataOffset[2] = {
3831 4 : poERDMS->GetIntField("layerStackDataOffset[0]"),
3832 4 : poERDMS->GetIntField("layerStackDataOffset[1]")};
3833 :
3834 4 : const GInt32 nStackCount = poERDMS->GetIntField("layerStackCount");
3835 4 : const GInt32 nStackIndex = poERDMS->GetIntField("layerStackIndex");
3836 :
3837 : // Update the filename.
3838 4 : if (strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0)
3839 : {
3840 4 : std::string osNew = pszNewBase;
3841 4 : osNew += osFileName.c_str() + strlen(pszOldBase);
3842 4 : osFileName = std::move(osNew);
3843 : }
3844 :
3845 : // Grow the node if needed.
3846 4 : if (strlen(pszNewBase) > strlen(pszOldBase))
3847 : {
3848 2 : CPLDebug("HFA", "Growing ExternalRasterDMS to hold new names");
3849 2 : poERDMS->MakeData(
3850 2 : static_cast<int>(poERDMS->GetDataSize() +
3851 2 : (strlen(pszNewBase) - strlen(pszOldBase))));
3852 : }
3853 :
3854 : // Initialize the whole thing to zeros for a clean start.
3855 4 : memset(poERDMS->GetData(), 0, poERDMS->GetDataSize());
3856 :
3857 : // Write it all out again, this may change the size of the node.
3858 4 : poERDMS->SetStringField("fileName.string", osFileName);
3859 4 : poERDMS->SetIntField("layerStackValidFlagsOffset[0]",
3860 : anValidFlagsOffset[0]);
3861 4 : poERDMS->SetIntField("layerStackValidFlagsOffset[1]",
3862 : anValidFlagsOffset[1]);
3863 :
3864 4 : poERDMS->SetIntField("layerStackDataOffset[0]", anStackDataOffset[0]);
3865 4 : poERDMS->SetIntField("layerStackDataOffset[1]", anStackDataOffset[1]);
3866 :
3867 4 : poERDMS->SetIntField("layerStackCount", nStackCount);
3868 4 : poERDMS->SetIntField("layerStackIndex", nStackIndex);
3869 : }
3870 :
3871 : // DependentFile.
3872 : apoNodeList =
3873 4 : hHFA->poRoot->FindChildren("DependentFile", "Eimg_DependentFile");
3874 :
3875 6 : for (size_t iNode = 0; iNode < apoNodeList.size(); iNode++)
3876 : {
3877 : CPLString osFileName =
3878 4 : apoNodeList[iNode]->GetStringField("dependent.string");
3879 :
3880 : // Grow the node if needed.
3881 2 : if (strlen(pszNewBase) > strlen(pszOldBase))
3882 : {
3883 1 : CPLDebug("HFA", "Growing DependentFile to hold new names");
3884 1 : apoNodeList[iNode]->MakeData(
3885 1 : static_cast<int>(apoNodeList[iNode]->GetDataSize() +
3886 1 : (strlen(pszNewBase) - strlen(pszOldBase))));
3887 : }
3888 :
3889 : // Update the filename.
3890 2 : if (strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0)
3891 : {
3892 2 : std::string osNew = pszNewBase;
3893 2 : osNew += (osFileName.c_str() + strlen(pszOldBase));
3894 2 : osFileName = std::move(osNew);
3895 : }
3896 :
3897 2 : apoNodeList[iNode]->SetStringField("dependent.string", osFileName);
3898 : }
3899 :
3900 8 : return CE_None;
3901 : }
3902 :
3903 : /* ==================================================================== */
3904 : /* Table relating USGS and ESRI state plane zones. */
3905 : /* ==================================================================== */
3906 : constexpr int anUsgsEsriZones[] = {
3907 : 101, 3101, 102, 3126, 201, 3151, 202, 3176, 203, 3201, 301, 3226,
3908 : 302, 3251, 401, 3276, 402, 3301, 403, 3326, 404, 3351, 405, 3376,
3909 : 406, 3401, 407, 3426, 501, 3451, 502, 3476, 503, 3501, 600, 3526,
3910 : 700, 3551, 901, 3601, 902, 3626, 903, 3576, 1001, 3651, 1002, 3676,
3911 : 1101, 3701, 1102, 3726, 1103, 3751, 1201, 3776, 1202, 3801, 1301, 3826,
3912 : 1302, 3851, 1401, 3876, 1402, 3901, 1501, 3926, 1502, 3951, 1601, 3976,
3913 : 1602, 4001, 1701, 4026, 1702, 4051, 1703, 6426, 1801, 4076, 1802, 4101,
3914 : 1900, 4126, 2001, 4151, 2002, 4176, 2101, 4201, 2102, 4226, 2103, 4251,
3915 : 2111, 6351, 2112, 6376, 2113, 6401, 2201, 4276, 2202, 4301, 2203, 4326,
3916 : 2301, 4351, 2302, 4376, 2401, 4401, 2402, 4426, 2403, 4451, 2500, 0,
3917 : 2501, 4476, 2502, 4501, 2503, 4526, 2600, 0, 2601, 4551, 2602, 4576,
3918 : 2701, 4601, 2702, 4626, 2703, 4651, 2800, 4676, 2900, 4701, 3001, 4726,
3919 : 3002, 4751, 3003, 4776, 3101, 4801, 3102, 4826, 3103, 4851, 3104, 4876,
3920 : 3200, 4901, 3301, 4926, 3302, 4951, 3401, 4976, 3402, 5001, 3501, 5026,
3921 : 3502, 5051, 3601, 5076, 3602, 5101, 3701, 5126, 3702, 5151, 3800, 5176,
3922 : 3900, 0, 3901, 5201, 3902, 5226, 4001, 5251, 4002, 5276, 4100, 5301,
3923 : 4201, 5326, 4202, 5351, 4203, 5376, 4204, 5401, 4205, 5426, 4301, 5451,
3924 : 4302, 5476, 4303, 5501, 4400, 5526, 4501, 5551, 4502, 5576, 4601, 5601,
3925 : 4602, 5626, 4701, 5651, 4702, 5676, 4801, 5701, 4802, 5726, 4803, 5751,
3926 : 4901, 5776, 4902, 5801, 4903, 5826, 4904, 5851, 5001, 6101, 5002, 6126,
3927 : 5003, 6151, 5004, 6176, 5005, 6201, 5006, 6226, 5007, 6251, 5008, 6276,
3928 : 5009, 6301, 5010, 6326, 5101, 5876, 5102, 5901, 5103, 5926, 5104, 5951,
3929 : 5105, 5976, 5201, 6001, 5200, 6026, 5200, 6076, 5201, 6051, 5202, 6051,
3930 : 5300, 0, 5400, 0};
3931 :
3932 : /************************************************************************/
3933 : /* ESRIToUSGSZone() */
3934 : /* */
3935 : /* Convert ESRI style state plane zones to USGS style state */
3936 : /* plane zones. */
3937 : /************************************************************************/
3938 :
3939 3 : static int ESRIToUSGSZone(int nESRIZone)
3940 :
3941 : {
3942 3 : if (nESRIZone == INT_MIN)
3943 0 : return 0;
3944 3 : if (nESRIZone < 0)
3945 1 : return std::abs(nESRIZone);
3946 :
3947 2 : const int nPairs = sizeof(anUsgsEsriZones) / (2 * sizeof(int));
3948 212 : for (int i = 0; i < nPairs; i++)
3949 : {
3950 212 : if (anUsgsEsriZones[i * 2 + 1] == nESRIZone)
3951 2 : return anUsgsEsriZones[i * 2];
3952 : }
3953 :
3954 0 : return 0;
3955 : }
3956 :
3957 : static const char *const apszDatumMap[] = {
3958 : // Imagine name, WKT name.
3959 : "NAD27",
3960 : "North_American_Datum_1927",
3961 : "NAD83",
3962 : "North_American_Datum_1983",
3963 : "WGS 84",
3964 : "WGS_1984",
3965 : "WGS 1972",
3966 : "WGS_1972",
3967 : "GDA94",
3968 : "Geocentric_Datum_of_Australia_1994",
3969 : "Pulkovo 1942",
3970 : "Pulkovo_1942",
3971 : "Geodetic Datum 1949",
3972 : "New_Zealand_Geodetic_Datum_1949",
3973 : nullptr,
3974 : nullptr};
3975 :
3976 337 : const char *const *HFAGetDatumMap()
3977 : {
3978 337 : return apszDatumMap;
3979 : }
3980 :
3981 : static const char *const apszUnitMap[] = {"meters",
3982 : "1.0",
3983 : "meter",
3984 : "1.0",
3985 : "m",
3986 : "1.0",
3987 : "centimeters",
3988 : "0.01",
3989 : "centimeter",
3990 : "0.01",
3991 : "cm",
3992 : "0.01",
3993 : "millimeters",
3994 : "0.001",
3995 : "millimeter",
3996 : "0.001",
3997 : "mm",
3998 : "0.001",
3999 : "kilometers",
4000 : "1000.0",
4001 : "kilometer",
4002 : "1000.0",
4003 : "km",
4004 : "1000.0",
4005 : "us_survey_feet",
4006 : "0.3048006096012192",
4007 : "us_survey_foot",
4008 : "0.3048006096012192",
4009 : "feet",
4010 : "0.3048006096012192",
4011 : "foot",
4012 : "0.3048006096012192",
4013 : "ft",
4014 : "0.3048006096012192",
4015 : "international_feet",
4016 : "0.3048",
4017 : "international_foot",
4018 : "0.3048",
4019 : "inches",
4020 : "0.0254000508001",
4021 : "inch",
4022 : "0.0254000508001",
4023 : "in",
4024 : "0.0254000508001",
4025 : "yards",
4026 : "0.9144",
4027 : "yard",
4028 : "0.9144",
4029 : "yd",
4030 : "0.9144",
4031 : "clarke_yard",
4032 : "0.9143917962",
4033 : "miles",
4034 : "1304.544",
4035 : "mile",
4036 : "1304.544",
4037 : "mi",
4038 : "1304.544",
4039 : "modified_american_feet",
4040 : "0.3048122530",
4041 : "modified_american_foot",
4042 : "0.3048122530",
4043 : "clarke_feet",
4044 : "0.3047972651",
4045 : "clarke_foot",
4046 : "0.3047972651",
4047 : "indian_feet",
4048 : "0.3047995142",
4049 : "indian_foot",
4050 : "0.3047995142",
4051 : nullptr,
4052 : nullptr};
4053 :
4054 197 : const char *const *HFAGetUnitMap()
4055 : {
4056 197 : return apszUnitMap;
4057 : }
4058 :
4059 : /************************************************************************/
4060 : /* HFAPCSStructToOSR() */
4061 : /* */
4062 : /* Convert the datum, proparameters and mapinfo structures into */
4063 : /* WKT format. */
4064 : /************************************************************************/
4065 :
4066 : std::unique_ptr<OGRSpatialReference>
4067 210 : HFAPCSStructToOSR(const Eprj_Datum *psDatum, const Eprj_ProParameters *psPro,
4068 : const Eprj_MapInfo *psMapInfo, HFAEntry *poMapInformation)
4069 :
4070 : {
4071 : // General case for Erdas style projections.
4072 :
4073 : // We make a particular effort to adapt the mapinfo->proname as
4074 : // the PROJCS[] name per #2422.
4075 420 : auto poSRS = std::make_unique<OGRSpatialReference>();
4076 210 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4077 :
4078 210 : if (psPro == nullptr && psMapInfo != nullptr)
4079 : {
4080 6 : poSRS->SetLocalCS(psMapInfo->proName);
4081 : }
4082 204 : else if (psPro == nullptr)
4083 : {
4084 0 : return nullptr;
4085 : }
4086 204 : else if (psPro->proType == EPRJ_EXTERNAL)
4087 : {
4088 4 : if (EQUALN(psPro->proExeName, EPRJ_EXTERNAL_NZMG, 4))
4089 : {
4090 : // Handle New Zealand Map Grid (NZMG) external projection. See:
4091 : // http://www.linz.govt.nz/
4092 : //
4093 : // Is there a better way that doesn't require hardcoding
4094 : // of these numbers?
4095 4 : poSRS->SetNZMG(-41.0, 173.0, 2510000, 6023150);
4096 : }
4097 : else
4098 : {
4099 0 : poSRS->SetLocalCS(psPro->proName);
4100 : }
4101 : }
4102 200 : else if (psPro->proNumber != EPRJ_LATLONG && psMapInfo != nullptr)
4103 : {
4104 167 : poSRS->SetProjCS(psMapInfo->proName);
4105 : }
4106 33 : else if (psPro->proNumber != EPRJ_LATLONG)
4107 : {
4108 7 : poSRS->SetProjCS(psPro->proName);
4109 : }
4110 :
4111 : // Handle units. It is important to deal with this first so
4112 : // that the projection Set methods will automatically do
4113 : // translation of linear values (like false easting) to PROJCS
4114 : // units from meters. Erdas linear projection values are
4115 : // always in meters.
4116 210 : if (poSRS->IsProjected() || poSRS->IsLocal())
4117 : {
4118 184 : const char *pszUnits = nullptr;
4119 :
4120 184 : if (psMapInfo)
4121 177 : pszUnits = psMapInfo->units;
4122 7 : else if (poMapInformation != nullptr)
4123 6 : pszUnits = poMapInformation->GetStringField("units.string");
4124 :
4125 184 : if (pszUnits != nullptr)
4126 : {
4127 183 : const char *const *papszUnitMap = HFAGetUnitMap();
4128 183 : int iUnitIndex = 0; // Used after for.
4129 304 : for (; papszUnitMap[iUnitIndex] != nullptr; iUnitIndex += 2)
4130 : {
4131 304 : if (EQUAL(papszUnitMap[iUnitIndex], pszUnits))
4132 183 : break;
4133 : }
4134 :
4135 183 : if (papszUnitMap[iUnitIndex] == nullptr)
4136 0 : iUnitIndex = 0;
4137 :
4138 366 : poSRS->SetLinearUnits(pszUnits,
4139 183 : CPLAtof(papszUnitMap[iUnitIndex + 1]));
4140 : }
4141 : else
4142 : {
4143 1 : poSRS->SetLinearUnits(SRS_UL_METER, 1.0);
4144 : }
4145 : }
4146 :
4147 210 : if (psPro == nullptr)
4148 : {
4149 6 : if (poSRS->IsLocal())
4150 : {
4151 6 : return poSRS;
4152 : }
4153 : else
4154 0 : return nullptr;
4155 : }
4156 :
4157 : // Try to work out ellipsoid and datum information.
4158 204 : const char *pszDatumName = psPro->proSpheroid.sphereName;
4159 204 : const char *pszEllipsoidName = psPro->proSpheroid.sphereName;
4160 :
4161 204 : if (psDatum != nullptr)
4162 : {
4163 204 : pszDatumName = psDatum->datumname;
4164 :
4165 : // Imagine to WKT translation.
4166 204 : const char *const *papszDatumMap = HFAGetDatumMap();
4167 741 : for (int i = 0; papszDatumMap[i] != nullptr; i += 2)
4168 : {
4169 689 : if (EQUAL(pszDatumName, papszDatumMap[i]))
4170 : {
4171 152 : pszDatumName = papszDatumMap[i + 1];
4172 152 : break;
4173 : }
4174 : }
4175 : }
4176 :
4177 204 : if (psPro->proSpheroid.a == 0.0)
4178 0 : ((Eprj_ProParameters *)psPro)->proSpheroid.a = 6378137.0;
4179 204 : if (psPro->proSpheroid.b == 0.0)
4180 0 : ((Eprj_ProParameters *)psPro)->proSpheroid.b = 6356752.3;
4181 :
4182 : const double dfInvFlattening =
4183 204 : OSRCalcInvFlattening(psPro->proSpheroid.a, psPro->proSpheroid.b);
4184 :
4185 : // Handle different projection methods.
4186 204 : switch (psPro->proNumber)
4187 : {
4188 30 : case EPRJ_LATLONG:
4189 30 : break;
4190 :
4191 90 : case EPRJ_UTM:
4192 : // We change this to unnamed so that SetUTM will set the long
4193 : // UTM description.
4194 90 : poSRS->SetProjCS("unnamed");
4195 90 : poSRS->SetUTM(psPro->proZone, psPro->proParams[3] >= 0.0);
4196 :
4197 : // The PCS name from the above function may be different with the
4198 : // input name. If there is a PCS name in psMapInfo that is
4199 : // different with the one in psPro, just use it as the PCS name.
4200 : // This case happens if the dataset's SR was written by the new
4201 : // GDAL.
4202 90 : if (psMapInfo && strlen(psMapInfo->proName) > 0 &&
4203 88 : strlen(psPro->proName) > 0 &&
4204 88 : !EQUAL(psMapInfo->proName, psPro->proName))
4205 63 : poSRS->SetProjCS(psMapInfo->proName);
4206 90 : break;
4207 :
4208 3 : case EPRJ_STATE_PLANE:
4209 : {
4210 3 : CPLString osUnitsName;
4211 : double dfLinearUnits;
4212 : {
4213 3 : const char *pszUnitsName = nullptr;
4214 3 : dfLinearUnits = poSRS->GetLinearUnits(&pszUnitsName);
4215 3 : if (pszUnitsName)
4216 3 : osUnitsName = pszUnitsName;
4217 : }
4218 :
4219 : // Historically, hfa used esri state plane zone code. Try esri pe
4220 : // string first.
4221 3 : const int zoneCode = ESRIToUSGSZone(psPro->proZone);
4222 : const char *pszDatum;
4223 3 : if (psDatum)
4224 3 : pszDatum = psDatum->datumname;
4225 : else
4226 0 : pszDatum = "HARN";
4227 : const char *pszUnits;
4228 3 : if (psMapInfo)
4229 0 : pszUnits = psMapInfo->units;
4230 3 : else if (!osUnitsName.empty())
4231 3 : pszUnits = osUnitsName;
4232 : else
4233 0 : pszUnits = "meters";
4234 3 : const int proNu = psPro->proNumber;
4235 3 : if (poSRS->ImportFromESRIStatePlaneWKT(zoneCode, pszDatum, pszUnits,
4236 3 : proNu) == OGRERR_NONE)
4237 : {
4238 3 : poSRS->AutoIdentifyEPSG();
4239 :
4240 3 : return poSRS;
4241 : }
4242 :
4243 : // Set state plane zone. Set NAD83/27 on basis of spheroid.
4244 0 : poSRS->SetStatePlane(ESRIToUSGSZone(psPro->proZone),
4245 0 : fabs(psPro->proSpheroid.a - 6378137.0) < 1.0,
4246 0 : osUnitsName.empty() ? nullptr
4247 0 : : osUnitsName.c_str(),
4248 : dfLinearUnits);
4249 :
4250 : // Same as the UTM, The following is needed.
4251 0 : if (psMapInfo && strlen(psMapInfo->proName) > 0 &&
4252 0 : strlen(psPro->proName) > 0 &&
4253 0 : !EQUAL(psMapInfo->proName, psPro->proName))
4254 0 : poSRS->SetProjCS(psMapInfo->proName);
4255 : }
4256 0 : break;
4257 :
4258 1 : case EPRJ_ALBERS_CONIC_EQUAL_AREA:
4259 1 : poSRS->SetACEA(psPro->proParams[2] * R2D, psPro->proParams[3] * R2D,
4260 1 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4261 1 : psPro->proParams[6], psPro->proParams[7]);
4262 1 : break;
4263 :
4264 6 : case EPRJ_LAMBERT_CONFORMAL_CONIC:
4265 : // Check the possible Wisconsin first.
4266 6 : if (psDatum && psMapInfo && EQUAL(psDatum->datumname, "HARN"))
4267 : {
4268 : // ERO: I doubt this works. Wisconsin LCC is LCC_1SP whereas
4269 : // we are here in the LCC_2SP case...
4270 0 : if (poSRS->ImportFromESRIWisconsinWKT(
4271 0 : "Lambert_Conformal_Conic", psPro->proParams[4] * R2D,
4272 0 : psPro->proParams[5] * R2D,
4273 0 : psMapInfo->units) == OGRERR_NONE)
4274 : {
4275 0 : return poSRS;
4276 : }
4277 : }
4278 6 : poSRS->SetLCC(psPro->proParams[2] * R2D, psPro->proParams[3] * R2D,
4279 6 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4280 6 : psPro->proParams[6], psPro->proParams[7]);
4281 6 : break;
4282 :
4283 1 : case EPRJ_MERCATOR:
4284 1 : poSRS->SetMercator(psPro->proParams[5] * R2D,
4285 1 : psPro->proParams[4] * R2D, 1.0,
4286 1 : psPro->proParams[6], psPro->proParams[7]);
4287 1 : break;
4288 :
4289 2 : case EPRJ_POLAR_STEREOGRAPHIC:
4290 2 : poSRS->SetPS(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4291 2 : 1.0, psPro->proParams[6], psPro->proParams[7]);
4292 2 : break;
4293 :
4294 2 : case EPRJ_POLYCONIC:
4295 2 : poSRS->SetPolyconic(psPro->proParams[5] * R2D,
4296 2 : psPro->proParams[4] * R2D, psPro->proParams[6],
4297 2 : psPro->proParams[7]);
4298 2 : break;
4299 :
4300 1 : case EPRJ_EQUIDISTANT_CONIC:
4301 : {
4302 2 : const double dfStdParallel2 = psPro->proParams[8] != 0.0
4303 1 : ? psPro->proParams[3] * R2D
4304 0 : : psPro->proParams[2] * R2D;
4305 1 : poSRS->SetEC(psPro->proParams[2] * R2D, dfStdParallel2,
4306 1 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4307 1 : psPro->proParams[6], psPro->proParams[7]);
4308 1 : break;
4309 : }
4310 18 : case EPRJ_TRANSVERSE_MERCATOR:
4311 : case EPRJ_GAUSS_KRUGER:
4312 : // Check the possible Wisconsin first.
4313 18 : if (psDatum && psMapInfo && EQUAL(psDatum->datumname, "HARN"))
4314 : {
4315 1 : if (poSRS->ImportFromESRIWisconsinWKT(
4316 1 : "Transverse_Mercator", psPro->proParams[4] * R2D,
4317 1 : psPro->proParams[5] * R2D,
4318 1 : psMapInfo->units) == OGRERR_NONE)
4319 : {
4320 1 : return poSRS;
4321 : }
4322 : }
4323 17 : poSRS->SetTM(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4324 17 : psPro->proParams[2], psPro->proParams[6],
4325 17 : psPro->proParams[7]);
4326 17 : break;
4327 :
4328 0 : case EPRJ_STEREOGRAPHIC:
4329 0 : poSRS->SetStereographic(psPro->proParams[5] * R2D,
4330 0 : psPro->proParams[4] * R2D, 1.0,
4331 0 : psPro->proParams[6], psPro->proParams[7]);
4332 0 : break;
4333 :
4334 1 : case EPRJ_LAMBERT_AZIMUTHAL_EQUAL_AREA:
4335 1 : poSRS->SetLAEA(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4336 1 : psPro->proParams[6], psPro->proParams[7]);
4337 1 : break;
4338 :
4339 1 : case EPRJ_AZIMUTHAL_EQUIDISTANT:
4340 1 : poSRS->SetAE(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4341 1 : psPro->proParams[6], psPro->proParams[7]);
4342 1 : break;
4343 :
4344 1 : case EPRJ_GNOMONIC:
4345 1 : poSRS->SetGnomonic(psPro->proParams[5] * R2D,
4346 1 : psPro->proParams[4] * R2D, psPro->proParams[6],
4347 1 : psPro->proParams[7]);
4348 1 : break;
4349 :
4350 2 : case EPRJ_ORTHOGRAPHIC:
4351 2 : poSRS->SetOrthographic(psPro->proParams[5] * R2D,
4352 2 : psPro->proParams[4] * R2D,
4353 2 : psPro->proParams[6], psPro->proParams[7]);
4354 2 : break;
4355 :
4356 1 : case EPRJ_SINUSOIDAL:
4357 1 : poSRS->SetSinusoidal(psPro->proParams[4] * R2D, psPro->proParams[6],
4358 1 : psPro->proParams[7]);
4359 1 : break;
4360 :
4361 3 : case EPRJ_PLATE_CARREE:
4362 : case EPRJ_EQUIRECTANGULAR:
4363 3 : poSRS->SetEquirectangular2(
4364 3 : 0.0, psPro->proParams[4] * R2D, psPro->proParams[5] * R2D,
4365 3 : psPro->proParams[6], psPro->proParams[7]);
4366 3 : break;
4367 :
4368 0 : case EPRJ_EQUIDISTANT_CYLINDRICAL:
4369 0 : poSRS->SetEquirectangular2(
4370 0 : 0.0, psPro->proParams[4] * R2D, psPro->proParams[2] * R2D,
4371 0 : psPro->proParams[6], psPro->proParams[7]);
4372 0 : break;
4373 :
4374 1 : case EPRJ_MILLER_CYLINDRICAL:
4375 1 : poSRS->SetMC(0.0, psPro->proParams[4] * R2D, psPro->proParams[6],
4376 1 : psPro->proParams[7]);
4377 1 : break;
4378 :
4379 1 : case EPRJ_VANDERGRINTEN:
4380 1 : poSRS->SetVDG(psPro->proParams[4] * R2D, psPro->proParams[6],
4381 1 : psPro->proParams[7]);
4382 1 : break;
4383 :
4384 0 : case EPRJ_HOTINE_OBLIQUE_MERCATOR:
4385 0 : if (psPro->proParams[12] > 0.0)
4386 0 : poSRS->SetHOM(
4387 0 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4388 0 : psPro->proParams[3] * R2D, 0.0, psPro->proParams[2],
4389 0 : psPro->proParams[6], psPro->proParams[7]);
4390 0 : break;
4391 :
4392 2 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER:
4393 2 : poSRS->SetHOMAC(
4394 2 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4395 2 : psPro->proParams[3] * R2D,
4396 2 : psPro->proParams[3] *
4397 : R2D, // We reuse azimuth as rectified_grid_angle
4398 2 : psPro->proParams[2], psPro->proParams[6], psPro->proParams[7]);
4399 2 : break;
4400 :
4401 1 : case EPRJ_ROBINSON:
4402 1 : poSRS->SetRobinson(psPro->proParams[4] * R2D, psPro->proParams[6],
4403 1 : psPro->proParams[7]);
4404 1 : break;
4405 :
4406 1 : case EPRJ_MOLLWEIDE:
4407 1 : poSRS->SetMollweide(psPro->proParams[4] * R2D, psPro->proParams[6],
4408 1 : psPro->proParams[7]);
4409 1 : break;
4410 :
4411 1 : case EPRJ_GALL_STEREOGRAPHIC:
4412 1 : poSRS->SetGS(psPro->proParams[4] * R2D, psPro->proParams[6],
4413 1 : psPro->proParams[7]);
4414 1 : break;
4415 :
4416 1 : case EPRJ_ECKERT_I:
4417 1 : poSRS->SetEckert(1, psPro->proParams[4] * R2D, psPro->proParams[6],
4418 1 : psPro->proParams[7]);
4419 1 : break;
4420 :
4421 1 : case EPRJ_ECKERT_II:
4422 1 : poSRS->SetEckert(2, psPro->proParams[4] * R2D, psPro->proParams[6],
4423 1 : psPro->proParams[7]);
4424 1 : break;
4425 :
4426 1 : case EPRJ_ECKERT_III:
4427 1 : poSRS->SetEckert(3, psPro->proParams[4] * R2D, psPro->proParams[6],
4428 1 : psPro->proParams[7]);
4429 1 : break;
4430 :
4431 1 : case EPRJ_ECKERT_IV:
4432 1 : poSRS->SetEckert(4, psPro->proParams[4] * R2D, psPro->proParams[6],
4433 1 : psPro->proParams[7]);
4434 1 : break;
4435 :
4436 1 : case EPRJ_ECKERT_V:
4437 1 : poSRS->SetEckert(5, psPro->proParams[4] * R2D, psPro->proParams[6],
4438 1 : psPro->proParams[7]);
4439 1 : break;
4440 :
4441 1 : case EPRJ_ECKERT_VI:
4442 1 : poSRS->SetEckert(6, psPro->proParams[4] * R2D, psPro->proParams[6],
4443 1 : psPro->proParams[7]);
4444 1 : break;
4445 :
4446 2 : case EPRJ_CASSINI:
4447 2 : poSRS->SetCS(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4448 2 : psPro->proParams[6], psPro->proParams[7]);
4449 2 : break;
4450 :
4451 1 : case EPRJ_TWO_POINT_EQUIDISTANT:
4452 1 : poSRS->SetTPED(psPro->proParams[9] * R2D, psPro->proParams[8] * R2D,
4453 1 : psPro->proParams[11] * R2D,
4454 1 : psPro->proParams[10] * R2D, psPro->proParams[6],
4455 1 : psPro->proParams[7]);
4456 1 : break;
4457 :
4458 0 : case EPRJ_STEREOGRAPHIC_EXTENDED:
4459 0 : poSRS->SetStereographic(
4460 0 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4461 0 : psPro->proParams[2], psPro->proParams[6], psPro->proParams[7]);
4462 0 : break;
4463 :
4464 1 : case EPRJ_BONNE:
4465 1 : poSRS->SetBonne(psPro->proParams[2] * R2D,
4466 1 : psPro->proParams[4] * R2D, psPro->proParams[6],
4467 1 : psPro->proParams[7]);
4468 1 : break;
4469 :
4470 1 : case EPRJ_LOXIMUTHAL:
4471 : {
4472 1 : poSRS->SetProjection("Loximuthal");
4473 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4474 1 : psPro->proParams[4] * R2D);
4475 1 : poSRS->SetNormProjParm("central_parallel",
4476 1 : psPro->proParams[5] * R2D);
4477 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4478 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4479 : }
4480 1 : break;
4481 :
4482 1 : case EPRJ_QUARTIC_AUTHALIC:
4483 : {
4484 1 : poSRS->SetProjection("Quartic_Authalic");
4485 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4486 1 : psPro->proParams[4] * R2D);
4487 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4488 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4489 : }
4490 1 : break;
4491 :
4492 1 : case EPRJ_WINKEL_I:
4493 : {
4494 1 : poSRS->SetProjection("Winkel_I");
4495 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4496 1 : psPro->proParams[4] * R2D);
4497 1 : poSRS->SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,
4498 1 : psPro->proParams[2] * R2D);
4499 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4500 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4501 : }
4502 1 : break;
4503 :
4504 1 : case EPRJ_WINKEL_II:
4505 : {
4506 1 : poSRS->SetProjection("Winkel_II");
4507 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4508 1 : psPro->proParams[4] * R2D);
4509 1 : poSRS->SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,
4510 1 : psPro->proParams[2] * R2D);
4511 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4512 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4513 : }
4514 1 : break;
4515 :
4516 0 : case EPRJ_BEHRMANN:
4517 : {
4518 0 : poSRS->SetProjection("Behrmann");
4519 0 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4520 0 : psPro->proParams[4] * R2D);
4521 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4522 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4523 : }
4524 0 : break;
4525 :
4526 1 : case EPRJ_KROVAK:
4527 1 : poSRS->SetKrovak(
4528 1 : psPro->proParams[4] * R2D, psPro->proParams[5] * R2D,
4529 1 : psPro->proParams[3] * R2D, psPro->proParams[9] * R2D,
4530 1 : psPro->proParams[2], psPro->proParams[6], psPro->proParams[7]);
4531 1 : break;
4532 :
4533 0 : case EPRJ_DOUBLE_STEREOGRAPHIC:
4534 : {
4535 0 : poSRS->SetProjection("Double_Stereographic");
4536 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,
4537 0 : psPro->proParams[5] * R2D);
4538 0 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4539 0 : psPro->proParams[4] * R2D);
4540 0 : poSRS->SetNormProjParm(SRS_PP_SCALE_FACTOR, psPro->proParams[2]);
4541 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4542 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4543 : }
4544 0 : break;
4545 :
4546 1 : case EPRJ_AITOFF:
4547 : {
4548 1 : poSRS->SetProjection("Aitoff");
4549 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4550 1 : psPro->proParams[4] * R2D);
4551 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4552 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4553 : }
4554 1 : break;
4555 :
4556 1 : case EPRJ_CRASTER_PARABOLIC:
4557 : {
4558 1 : poSRS->SetProjection("Craster_Parabolic");
4559 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4560 1 : psPro->proParams[4] * R2D);
4561 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4562 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4563 : }
4564 1 : break;
4565 :
4566 1 : case EPRJ_CYLINDRICAL_EQUAL_AREA:
4567 1 : poSRS->SetCEA(psPro->proParams[2] * R2D, psPro->proParams[4] * R2D,
4568 1 : psPro->proParams[6], psPro->proParams[7]);
4569 1 : break;
4570 :
4571 1 : case EPRJ_FLAT_POLAR_QUARTIC:
4572 : {
4573 1 : poSRS->SetProjection("Flat_Polar_Quartic");
4574 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4575 1 : psPro->proParams[4] * R2D);
4576 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4577 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4578 : }
4579 1 : break;
4580 :
4581 1 : case EPRJ_TIMES:
4582 : {
4583 1 : poSRS->SetProjection("Times");
4584 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4585 1 : psPro->proParams[4] * R2D);
4586 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4587 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4588 : }
4589 1 : break;
4590 :
4591 1 : case EPRJ_WINKEL_TRIPEL:
4592 : {
4593 1 : poSRS->SetProjection("Winkel_Tripel");
4594 1 : poSRS->SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,
4595 1 : psPro->proParams[2] * R2D);
4596 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4597 1 : psPro->proParams[4] * R2D);
4598 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4599 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4600 : }
4601 1 : break;
4602 :
4603 1 : case EPRJ_HAMMER_AITOFF:
4604 : {
4605 1 : poSRS->SetProjection("Hammer_Aitoff");
4606 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4607 1 : psPro->proParams[4] * R2D);
4608 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4609 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4610 : }
4611 1 : break;
4612 :
4613 1 : case EPRJ_VERTICAL_NEAR_SIDE_PERSPECTIVE:
4614 : {
4615 1 : poSRS->SetVerticalPerspective(
4616 1 : psPro->proParams[5] * R2D, // dfTopoOriginLat
4617 1 : psPro->proParams[4] * R2D, // dfTopoOriginLon
4618 : 0, // dfTopoOriginHeight
4619 1 : psPro->proParams[2], // dfViewPointHeight
4620 1 : psPro->proParams[6], // dfFalseEasting
4621 1 : psPro->proParams[7]); // dfFalseNorthing
4622 : }
4623 1 : break;
4624 :
4625 0 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_CENTER:
4626 : {
4627 0 : poSRS->SetProjection("Hotine_Oblique_Mercator_Twp_Point_Center");
4628 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_CENTER,
4629 0 : psPro->proParams[5] * R2D);
4630 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_1ST_POINT,
4631 0 : psPro->proParams[9] * R2D);
4632 0 : poSRS->SetNormProjParm(SRS_PP_LONGITUDE_OF_1ST_POINT,
4633 0 : psPro->proParams[8] * R2D);
4634 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_2ND_POINT,
4635 0 : psPro->proParams[11] * R2D);
4636 0 : poSRS->SetNormProjParm(SRS_PP_LONGITUDE_OF_2ND_POINT,
4637 0 : psPro->proParams[10] * R2D);
4638 0 : poSRS->SetNormProjParm(SRS_PP_SCALE_FACTOR, psPro->proParams[2]);
4639 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4640 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4641 : }
4642 0 : break;
4643 :
4644 1 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN:
4645 1 : poSRS->SetHOM2PNO(
4646 1 : psPro->proParams[5] * R2D, psPro->proParams[8] * R2D,
4647 1 : psPro->proParams[9] * R2D, psPro->proParams[10] * R2D,
4648 1 : psPro->proParams[11] * R2D, psPro->proParams[2],
4649 1 : psPro->proParams[6], psPro->proParams[7]);
4650 1 : break;
4651 :
4652 0 : case EPRJ_LAMBERT_CONFORMAL_CONIC_1SP:
4653 0 : poSRS->SetLCC1SP(psPro->proParams[3] * R2D,
4654 0 : psPro->proParams[2] * R2D, psPro->proParams[4],
4655 0 : psPro->proParams[5], psPro->proParams[6]);
4656 0 : break;
4657 :
4658 2 : case EPRJ_MERCATOR_VARIANT_A:
4659 2 : poSRS->SetMercator(psPro->proParams[5] * R2D,
4660 2 : psPro->proParams[4] * R2D, psPro->proParams[2],
4661 2 : psPro->proParams[6], psPro->proParams[7]);
4662 2 : break;
4663 :
4664 0 : case EPRJ_PSEUDO_MERCATOR: // Likely this is google mercator?
4665 0 : poSRS->SetMercator(psPro->proParams[5] * R2D,
4666 0 : psPro->proParams[4] * R2D, 1.0,
4667 0 : psPro->proParams[6], psPro->proParams[7]);
4668 0 : break;
4669 :
4670 5 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_VARIANT_A:
4671 5 : poSRS->SetHOM(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4672 5 : psPro->proParams[3] * R2D, psPro->proParams[8] * R2D,
4673 5 : psPro->proParams[2], psPro->proParams[6],
4674 5 : psPro->proParams[7]);
4675 5 : break;
4676 :
4677 3 : case EPRJ_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED:
4678 3 : poSRS->SetTMSO(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4679 3 : psPro->proParams[2], psPro->proParams[6],
4680 3 : psPro->proParams[7]);
4681 3 : break;
4682 :
4683 0 : default:
4684 0 : if (poSRS->IsProjected())
4685 0 : poSRS->GetRoot()->SetValue("LOCAL_CS");
4686 : else
4687 0 : poSRS->SetLocalCS(psPro->proName);
4688 0 : break;
4689 : }
4690 :
4691 : // Try and set the GeogCS information.
4692 200 : if (!poSRS->IsLocal())
4693 : {
4694 200 : bool bWellKnownDatum = false;
4695 200 : if (pszDatumName == nullptr)
4696 0 : poSRS->SetGeogCS(pszDatumName, pszDatumName, pszEllipsoidName,
4697 0 : psPro->proSpheroid.a, dfInvFlattening);
4698 200 : else if (EQUAL(pszDatumName, "WGS 84") ||
4699 200 : EQUAL(pszDatumName, "WGS_1984"))
4700 : {
4701 36 : bWellKnownDatum = true;
4702 36 : poSRS->SetWellKnownGeogCS("WGS84");
4703 : }
4704 164 : else if (strstr(pszDatumName, "NAD27") != nullptr ||
4705 162 : EQUAL(pszDatumName, "North_American_Datum_1927"))
4706 : {
4707 88 : bWellKnownDatum = true;
4708 88 : poSRS->SetWellKnownGeogCS("NAD27");
4709 : }
4710 76 : else if (EQUAL(pszDatumName, "NAD83") ||
4711 76 : EQUAL(pszDatumName, "North_American_Datum_1983"))
4712 : {
4713 7 : bWellKnownDatum = true;
4714 7 : poSRS->SetWellKnownGeogCS("NAD83");
4715 : }
4716 : else
4717 : {
4718 138 : CPLString osGeogCRSName(pszDatumName);
4719 :
4720 69 : if (poSRS->IsProjected())
4721 : {
4722 68 : PJ_CONTEXT *ctxt = OSRGetProjTLSContext();
4723 68 : const PJ_TYPE type = PJ_TYPE_PROJECTED_CRS;
4724 : PJ_OBJ_LIST *list =
4725 68 : proj_create_from_name(ctxt, nullptr, poSRS->GetName(),
4726 : &type, 1, false, 1, nullptr);
4727 68 : if (list)
4728 : {
4729 68 : const auto listSize = proj_list_get_count(list);
4730 68 : if (listSize == 1)
4731 : {
4732 43 : auto crs = proj_list_get(ctxt, list, 0);
4733 43 : if (crs)
4734 : {
4735 43 : auto geogCRS = proj_crs_get_geodetic_crs(ctxt, crs);
4736 43 : if (geogCRS)
4737 : {
4738 43 : const char *pszName = proj_get_name(geogCRS);
4739 43 : if (pszName)
4740 43 : osGeogCRSName = pszName;
4741 43 : proj_destroy(geogCRS);
4742 : }
4743 43 : proj_destroy(crs);
4744 : }
4745 : }
4746 68 : proj_list_destroy(list);
4747 : }
4748 : }
4749 :
4750 138 : poSRS->SetGeogCS(osGeogCRSName, pszDatumName, pszEllipsoidName,
4751 69 : psPro->proSpheroid.a, dfInvFlattening);
4752 : }
4753 :
4754 200 : if (psDatum != nullptr && psDatum->type == EPRJ_DATUM_PARAMETRIC)
4755 : {
4756 159 : if (bWellKnownDatum &&
4757 45 : CPLTestBool(CPLGetConfigOption("OSR_STRIP_TOWGS84", "YES")))
4758 : {
4759 45 : CPLDebug("OSR",
4760 : "TOWGS84 information has been removed. "
4761 : "It can be kept by setting the OSR_STRIP_TOWGS84 "
4762 : "configuration option to NO");
4763 : }
4764 : else
4765 : {
4766 69 : poSRS->SetTOWGS84(psDatum->params[0], psDatum->params[1],
4767 69 : psDatum->params[2],
4768 69 : -psDatum->params[3] * RAD2ARCSEC,
4769 69 : -psDatum->params[4] * RAD2ARCSEC,
4770 69 : -psDatum->params[5] * RAD2ARCSEC,
4771 69 : psDatum->params[6] * 1e+6);
4772 69 : poSRS->StripTOWGS84IfKnownDatumAndAllowed();
4773 : }
4774 : }
4775 : }
4776 :
4777 : // Try to insert authority information if possible.
4778 200 : poSRS->AutoIdentifyEPSG();
4779 :
4780 200 : auto poSRSBestMatch = poSRS->FindBestMatch(90, nullptr, nullptr);
4781 200 : if (poSRSBestMatch)
4782 : {
4783 137 : poSRS.reset(poSRSBestMatch);
4784 : }
4785 :
4786 200 : return poSRS;
4787 : }
|
