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