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