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 680 : const char *const *GetHFAAuxMetaDataList()
115 : {
116 680 : return apszAuxMetadataItems;
117 : }
118 :
119 : /************************************************************************/
120 : /* HFAGetDictionary() */
121 : /************************************************************************/
122 :
123 554 : static char *HFAGetDictionary(HFAHandle hHFA)
124 :
125 : {
126 554 : int nDictMax = 100;
127 554 : char *pszDictionary = static_cast<char *>(CPLMalloc(nDictMax));
128 554 : int nDictSize = 0;
129 :
130 554 : 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 2045100 : if (nDictSize >= nDictMax - 1)
139 : {
140 2751 : nDictMax = nDictSize * 2 + 100;
141 : pszDictionary =
142 2751 : static_cast<char *>(CPLRealloc(pszDictionary, nDictMax));
143 : }
144 :
145 2045100 : if (VSIFReadL(pszDictionary + nDictSize, 1, 1, hHFA->fp) < 1 ||
146 4088000 : pszDictionary[nDictSize] == '\0' ||
147 2042890 : (nDictSize > 2 && pszDictionary[nDictSize - 2] == ',' &&
148 159645 : pszDictionary[nDictSize - 1] == '.'))
149 554 : break;
150 :
151 2044550 : nDictSize++;
152 : }
153 :
154 554 : pszDictionary[nDictSize] = '\0';
155 :
156 554 : return pszDictionary;
157 : }
158 :
159 : /************************************************************************/
160 : /* HFAOpen() */
161 : /************************************************************************/
162 :
163 554 : HFAHandle HFAOpen(const char *pszFilename, const char *pszAccess)
164 :
165 : {
166 554 : VSILFILE *fp = VSIFOpenL(
167 : pszFilename,
168 554 : (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 554 : 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 554 : char szHeader[16] = {};
182 554 : 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 554 : 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 554 : static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
202 :
203 554 : psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
204 554 : psInfo->pszPath = CPLStrdup(CPLGetPathSafe(pszFilename).c_str());
205 554 : psInfo->fp = fp;
206 554 : if (EQUAL(pszAccess, "r") || EQUAL(pszAccess, "rb"))
207 341 : psInfo->eAccess = HFA_ReadOnly;
208 : else
209 213 : psInfo->eAccess = HFA_Update;
210 554 : psInfo->bTreeDirty = false;
211 :
212 : // Where is the header?
213 554 : GUInt32 nHeaderPos = 0;
214 554 : bool bRet = VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, fp) > 0;
215 : HFAStandard(4, &nHeaderPos);
216 :
217 : // Read the header.
218 554 : bRet &= VSIFSeekL(fp, nHeaderPos, SEEK_SET) >= 0;
219 :
220 554 : bRet &= VSIFReadL(&(psInfo->nVersion), sizeof(GInt32), 1, fp) > 0;
221 : HFAStandard(4, &(psInfo->nVersion));
222 :
223 554 : bRet &= VSIFReadL(szHeader, 4, 1, fp) > 0; // Skip freeList.
224 :
225 554 : bRet &= VSIFReadL(&(psInfo->nRootPos), sizeof(GInt32), 1, fp) > 0;
226 : HFAStandard(4, &(psInfo->nRootPos));
227 :
228 554 : bRet &= VSIFReadL(&(psInfo->nEntryHeaderLength), sizeof(GInt16), 1, fp) > 0;
229 : HFAStandard(2, &(psInfo->nEntryHeaderLength));
230 :
231 554 : bRet &= VSIFReadL(&(psInfo->nDictionaryPos), sizeof(GInt32), 1, fp) > 0;
232 : HFAStandard(4, &(psInfo->nDictionaryPos));
233 :
234 : // Collect file size.
235 554 : bRet &= VSIFSeekL(fp, 0, SEEK_END) >= 0;
236 554 : 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 554 : psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
245 :
246 : // Instantiate the root entry.
247 554 : psInfo->poRoot = HFAEntry::New(psInfo, psInfo->nRootPos, nullptr, nullptr);
248 554 : 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 554 : psInfo->pszDictionary = HFAGetDictionary(psInfo);
259 554 : psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
260 :
261 : // Collect band definitions.
262 554 : HFAParseBandInfo(psInfo);
263 :
264 554 : 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 49 : HFAInfo_t *HFAGetDependent(HFAInfo_t *psBase, const char *pszFilename)
327 :
328 : {
329 49 : if (EQUAL(pszFilename, psBase->pszFilename))
330 30 : return psBase;
331 :
332 19 : if (psBase->psDependent != nullptr)
333 : {
334 2 : if (EQUAL(pszFilename, psBase->psDependent->pszFilename))
335 2 : return psBase->psDependent;
336 : else
337 0 : return nullptr;
338 : }
339 :
340 : // Try to open the dependent file.
341 17 : const char *pszMode = psBase->eAccess == HFA_Update ? "r+b" : "rb";
342 :
343 17 : char *pszDependent = CPLStrdup(
344 34 : CPLFormFilenameSafe(psBase->pszPath, pszFilename, nullptr).c_str());
345 :
346 17 : VSILFILE *fp = VSIFOpenL(pszDependent, pszMode);
347 17 : if (fp != nullptr)
348 : {
349 13 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
350 13 : psBase->psDependent = HFAOpen(pszDependent, pszMode);
351 : }
352 :
353 17 : CPLFree(pszDependent);
354 :
355 17 : 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 743 : CPLErr HFAParseBandInfo(HFAInfo_t *psInfo)
366 :
367 : {
368 : // Find the first band node.
369 743 : psInfo->nBands = 0;
370 743 : HFAEntry *poNode = psInfo->poRoot->GetChild();
371 2430 : while (poNode != nullptr)
372 : {
373 1687 : if (EQUAL(poNode->GetType(), "Eimg_Layer") &&
374 2530 : poNode->GetIntField("width") > 0 &&
375 843 : poNode->GetIntField("height") > 0)
376 : {
377 843 : if (psInfo->nBands == 0)
378 : {
379 725 : psInfo->nXSize = poNode->GetIntField("width");
380 725 : 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 1686 : psInfo->papoBand = static_cast<HFABand **>(CPLRealloc(
389 843 : psInfo->papoBand, sizeof(HFABand *) * (psInfo->nBands + 1)));
390 843 : psInfo->papoBand[psInfo->nBands] = new HFABand(psInfo, poNode);
391 843 : if (psInfo->papoBand[psInfo->nBands]->nWidth == 0)
392 : {
393 0 : delete psInfo->papoBand[psInfo->nBands];
394 0 : return CE_Failure;
395 : }
396 843 : psInfo->nBands++;
397 : }
398 :
399 1687 : poNode = poNode->GetNext();
400 : }
401 :
402 743 : return CE_None;
403 : }
404 :
405 : /************************************************************************/
406 : /* HFAClose() */
407 : /************************************************************************/
408 :
409 752 : int HFAClose(HFAHandle hHFA)
410 :
411 : {
412 752 : if (hHFA->eAccess == HFA_Update &&
413 411 : (hHFA->bTreeDirty || (hHFA->poDictionary != nullptr &&
414 14 : hHFA->poDictionary->bDictionaryTextDirty)))
415 393 : HFAFlush(hHFA);
416 :
417 752 : int nRet = 0;
418 752 : if (hHFA->psDependent != nullptr)
419 : {
420 13 : if (HFAClose(hHFA->psDependent) != 0)
421 0 : nRet = -1;
422 : }
423 :
424 752 : delete hHFA->poRoot;
425 :
426 752 : if (VSIFCloseL(hHFA->fp) != 0)
427 0 : nRet = -1;
428 :
429 752 : if (hHFA->poDictionary != nullptr)
430 748 : delete hHFA->poDictionary;
431 :
432 752 : CPLFree(hHFA->pszDictionary);
433 752 : CPLFree(hHFA->pszFilename);
434 752 : CPLFree(hHFA->pszIGEFilename);
435 752 : CPLFree(hHFA->pszPath);
436 :
437 1595 : for (int i = 0; i < hHFA->nBands; i++)
438 : {
439 843 : delete hHFA->papoBand[i];
440 : }
441 :
442 752 : CPLFree(hHFA->papoBand);
443 :
444 752 : 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 752 : 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 752 : 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 752 : CPLFree(hHFA);
471 752 : 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 539 : CPLErr HFAGetRasterInfo(HFAHandle hHFA, int *pnXSize, int *pnYSize,
549 : int *pnBands)
550 :
551 : {
552 539 : if (pnXSize != nullptr)
553 539 : *pnXSize = hHFA->nXSize;
554 539 : if (pnYSize != nullptr)
555 539 : *pnYSize = hHFA->nYSize;
556 539 : if (pnBands != nullptr)
557 539 : *pnBands = hHFA->nBands;
558 539 : return CE_None;
559 : }
560 :
561 : /************************************************************************/
562 : /* HFAGetBandInfo() */
563 : /************************************************************************/
564 :
565 669 : CPLErr HFAGetBandInfo(HFAHandle hHFA, int nBand, EPTType *peDataType,
566 : int *pnBlockXSize, int *pnBlockYSize,
567 : int *pnCompressionType)
568 :
569 : {
570 669 : if (nBand < 0 || nBand > hHFA->nBands)
571 : {
572 0 : CPLAssert(false);
573 : return CE_Failure;
574 : }
575 :
576 669 : HFABand *poBand = hHFA->papoBand[nBand - 1];
577 :
578 669 : if (peDataType != nullptr)
579 669 : *peDataType = poBand->eDataType;
580 :
581 669 : if (pnBlockXSize != nullptr)
582 669 : *pnBlockXSize = poBand->nBlockXSize;
583 :
584 669 : if (pnBlockYSize != nullptr)
585 669 : *pnBlockYSize = poBand->nBlockYSize;
586 :
587 : // Get compression code from RasterDMS.
588 669 : if (pnCompressionType != nullptr)
589 : {
590 669 : *pnCompressionType = 0;
591 :
592 669 : HFAEntry *poDMS = poBand->poNode->GetNamedChild("RasterDMS");
593 :
594 669 : if (poDMS != nullptr)
595 579 : *pnCompressionType = poDMS->GetIntField("compressionType");
596 : }
597 :
598 669 : 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 9 : CPLErr HFASetBandNoData(HFAHandle hHFA, int nBand, double dfValue)
636 :
637 : {
638 9 : if (nBand < 0 || nBand > hHFA->nBands)
639 : {
640 0 : CPLAssert(false);
641 : return CE_Failure;
642 : }
643 :
644 9 : HFABand *poBand = hHFA->papoBand[nBand - 1];
645 :
646 9 : return poBand->SetNoDataValue(dfValue);
647 : }
648 :
649 : /************************************************************************/
650 : /* HFAGetOverviewCount() */
651 : /************************************************************************/
652 :
653 123 : int HFAGetOverviewCount(HFAHandle hHFA, int nBand)
654 :
655 : {
656 123 : if (nBand < 0 || nBand > hHFA->nBands)
657 : {
658 0 : CPLAssert(false);
659 : return CE_Failure;
660 : }
661 :
662 123 : HFABand *poBand = hHFA->papoBand[nBand - 1];
663 123 : poBand->LoadOverviews();
664 :
665 123 : return poBand->nOverviews;
666 : }
667 :
668 : /************************************************************************/
669 : /* HFAGetOverviewInfo() */
670 : /************************************************************************/
671 :
672 54 : CPLErr HFAGetOverviewInfo(HFAHandle hHFA, int nBand, int iOverview,
673 : int *pnXSize, int *pnYSize, int *pnBlockXSize,
674 : int *pnBlockYSize, EPTType *peHFADataType)
675 :
676 : {
677 54 : if (nBand < 0 || nBand > hHFA->nBands)
678 : {
679 0 : CPLAssert(false);
680 : return CE_Failure;
681 : }
682 :
683 54 : HFABand *poBand = hHFA->papoBand[nBand - 1];
684 54 : poBand->LoadOverviews();
685 :
686 54 : if (iOverview < 0 || iOverview >= poBand->nOverviews)
687 : {
688 0 : CPLAssert(false);
689 : return CE_Failure;
690 : }
691 54 : poBand = poBand->papoOverviews[iOverview];
692 54 : if (poBand == nullptr)
693 : {
694 0 : return CE_Failure;
695 : }
696 :
697 54 : if (pnXSize != nullptr)
698 54 : *pnXSize = poBand->nWidth;
699 :
700 54 : if (pnYSize != nullptr)
701 54 : *pnYSize = poBand->nHeight;
702 :
703 54 : if (pnBlockXSize != nullptr)
704 54 : *pnBlockXSize = poBand->nBlockXSize;
705 :
706 54 : if (pnBlockYSize != nullptr)
707 54 : *pnBlockYSize = poBand->nBlockYSize;
708 :
709 54 : if (peHFADataType != nullptr)
710 54 : *peHFADataType = poBand->eDataType;
711 :
712 54 : 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 1330 : CPLErr HFAGetRasterBlockEx(HFAHandle hHFA, int nBand, int nXBlock, int nYBlock,
731 : void *pData, int nDataSize)
732 :
733 : {
734 1330 : if (nBand < 1 || nBand > hHFA->nBands)
735 0 : return CE_Failure;
736 :
737 1330 : return hHFA->papoBand[nBand - 1]->GetRasterBlock(nXBlock, nYBlock, pData,
738 1330 : 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 105 : CPLErr HFASetRasterBlock(HFAHandle hHFA, int nBand, int nXBlock, int nYBlock,
777 : void *pData)
778 :
779 : {
780 105 : if (nBand < 1 || nBand > hHFA->nBands)
781 0 : return CE_Failure;
782 :
783 105 : 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 212 : const char *HFAGetBandName(HFAHandle hHFA, int nBand)
809 : {
810 212 : if (nBand < 1 || nBand > hHFA->nBands)
811 0 : return "";
812 :
813 212 : 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 4712 : int HFAGetDataTypeBits(EPTType eDataType)
833 :
834 : {
835 4712 : switch (eDataType)
836 : {
837 76 : case EPT_u1:
838 76 : return 1;
839 :
840 42 : case EPT_u2:
841 42 : return 2;
842 :
843 18 : case EPT_u4:
844 18 : return 4;
845 :
846 1328 : case EPT_u8:
847 : case EPT_s8:
848 1328 : return 8;
849 :
850 1148 : case EPT_u16:
851 : case EPT_s16:
852 1148 : return 16;
853 :
854 268 : case EPT_u32:
855 : case EPT_s32:
856 : case EPT_f32:
857 268 : return 32;
858 :
859 1793 : case EPT_f64:
860 : case EPT_c64:
861 1793 : 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 1070 : const Eprj_MapInfo *HFAGetMapInfo(HFAHandle hHFA)
930 :
931 : {
932 1070 : if (hHFA->nBands < 1)
933 0 : return nullptr;
934 :
935 : // Do we already have it?
936 1070 : 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 825 : HFAEntry *poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild("Map_Info");
942 825 : if (poMIEntry == nullptr)
943 : {
944 580 : for (HFAEntry *poChild = hHFA->papoBand[0]->poNode->GetChild();
945 2008 : 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 825 : if (poMIEntry == nullptr)
954 : {
955 580 : 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 535 : int HFAGetGeoTransform(HFAHandle hHFA, double *padfGeoTransform)
1030 :
1031 : {
1032 535 : const Eprj_MapInfo *psMapInfo = HFAGetMapInfo(hHFA);
1033 :
1034 535 : padfGeoTransform[0] = 0.0;
1035 535 : padfGeoTransform[1] = 1.0;
1036 535 : padfGeoTransform[2] = 0.0;
1037 535 : padfGeoTransform[3] = 0.0;
1038 535 : padfGeoTransform[4] = 0.0;
1039 535 : padfGeoTransform[5] = 1.0;
1040 :
1041 : // Simple (north up) MapInfo approach.
1042 535 : 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 290 : if (hHFA->nBands == 0)
1078 0 : return FALSE;
1079 :
1080 : HFAEntry *poXForm0 =
1081 290 : hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
1082 :
1083 290 : if (poXForm0 == nullptr)
1084 284 : 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 145 : CPLErr HFASetMapInfo(HFAHandle hHFA, const Eprj_MapInfo *poMapInfo)
1128 :
1129 : {
1130 : // Loop over bands, setting information on each one.
1131 328 : 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 183 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Map_Info");
1136 183 : if (poMIEntry == nullptr)
1137 : {
1138 182 : poMIEntry = HFAEntry::New(hHFA, "Map_Info", "Eprj_MapInfo",
1139 182 : hHFA->papoBand[iBand]->poNode);
1140 : }
1141 :
1142 183 : poMIEntry->MarkDirty();
1143 :
1144 : // Ensure we have enough space for all the data.
1145 : // TODO(schwehr): Explain 48 and 40 constants.
1146 183 : const int nSize =
1147 183 : static_cast<int>(48 + 40 + strlen(poMapInfo->proName) + 1 +
1148 183 : strlen(poMapInfo->units) + 1);
1149 :
1150 183 : GByte *pabyData = poMIEntry->MakeData(nSize);
1151 183 : memset(pabyData, 0, nSize);
1152 :
1153 183 : poMIEntry->SetPosition();
1154 :
1155 : // Write the various fields.
1156 183 : poMIEntry->SetStringField("proName", poMapInfo->proName);
1157 :
1158 183 : poMIEntry->SetDoubleField("upperLeftCenter.x",
1159 183 : poMapInfo->upperLeftCenter.x);
1160 183 : poMIEntry->SetDoubleField("upperLeftCenter.y",
1161 183 : poMapInfo->upperLeftCenter.y);
1162 :
1163 183 : poMIEntry->SetDoubleField("lowerRightCenter.x",
1164 183 : poMapInfo->lowerRightCenter.x);
1165 183 : poMIEntry->SetDoubleField("lowerRightCenter.y",
1166 183 : poMapInfo->lowerRightCenter.y);
1167 :
1168 183 : poMIEntry->SetDoubleField("pixelSize.width",
1169 183 : poMapInfo->pixelSize.width);
1170 183 : poMIEntry->SetDoubleField("pixelSize.height",
1171 183 : poMapInfo->pixelSize.height);
1172 :
1173 183 : poMIEntry->SetStringField("units", poMapInfo->units);
1174 : }
1175 :
1176 145 : 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 133 : CPLErr HFASetPEString(HFAHandle hHFA, const char *pszPEString)
1229 :
1230 : {
1231 133 : if (!CPLTestBool(CPLGetConfigOption("HFA_WRITE_PE_STRING", "YES")))
1232 0 : return CE_None;
1233 :
1234 : // Loop over bands, setting information on each one.
1235 304 : 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 171 : hHFA->papoBand[iBand]->poNode->GetNamedChild("ProjectionX");
1241 :
1242 : // If we are setting an empty string then a missing entry is equivalent.
1243 171 : if (strlen(pszPEString) == 0 && poProX == nullptr)
1244 1 : continue;
1245 :
1246 : // Create the node.
1247 170 : if (poProX == nullptr)
1248 : {
1249 338 : poProX = HFAEntry::New(hHFA, "ProjectionX", "Eprj_MapProjection842",
1250 169 : hHFA->papoBand[iBand]->poNode);
1251 169 : 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 170 : poProX->MakeData(static_cast<int>(700 + strlen(pszPEString)));
1258 170 : if (!pabyData)
1259 0 : return CE_Failure;
1260 :
1261 170 : memset(pabyData, 0, 250 + strlen(pszPEString));
1262 :
1263 170 : poProX->SetPosition();
1264 :
1265 170 : poProX->SetStringField("projection.type.string", "PE_COORDSYS");
1266 170 : 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 170 : pabyData = poProX->GetData();
1274 170 : int nDataSize = poProX->GetDataSize();
1275 170 : GUInt32 iOffset = poProX->GetDataPos();
1276 :
1277 15810 : while (nDataSize > 10 &&
1278 15810 : !STARTS_WITH_CI((const char *)pabyData, "PE_COORDSYS,."))
1279 : {
1280 15640 : pabyData++;
1281 15640 : nDataSize--;
1282 15640 : iOffset++;
1283 : }
1284 :
1285 170 : CPLAssert(nDataSize > static_cast<int>(strlen(pszPEString)) + 10);
1286 :
1287 170 : pabyData += 14;
1288 170 : iOffset += 14;
1289 :
1290 : // Set the size and offset of the mifobject.
1291 170 : iOffset += 8;
1292 :
1293 170 : GUInt32 nSize = static_cast<GUInt32>(strlen(pszPEString) + 9);
1294 :
1295 : HFAStandard(4, &nSize);
1296 170 : memcpy(pabyData, &nSize, 4);
1297 170 : pabyData += 4;
1298 :
1299 : HFAStandard(4, &iOffset);
1300 170 : memcpy(pabyData, &iOffset, 4);
1301 170 : pabyData += 4;
1302 :
1303 : // Set the size and offset of the string value.
1304 170 : nSize = static_cast<GUInt32>(strlen(pszPEString) + 1);
1305 :
1306 : HFAStandard(4, &nSize);
1307 170 : memcpy(pabyData, &nSize, 4);
1308 170 : pabyData += 4;
1309 :
1310 170 : iOffset = 8;
1311 : HFAStandard(4, &iOffset);
1312 170 : memcpy(pabyData, &iOffset, 4);
1313 170 : pabyData += 4;
1314 :
1315 : // Place the string itself.
1316 170 : memcpy(pabyData, pszPEString, strlen(pszPEString) + 1);
1317 :
1318 170 : poProX->SetStringField("title.string", "PE");
1319 : }
1320 :
1321 133 : return CE_None;
1322 : }
1323 :
1324 : /************************************************************************/
1325 : /* HFAGetProParameters() */
1326 : /************************************************************************/
1327 :
1328 535 : const Eprj_ProParameters *HFAGetProParameters(HFAHandle hHFA)
1329 :
1330 : {
1331 535 : if (hHFA->nBands < 1)
1332 0 : return nullptr;
1333 :
1334 : // Do we already have it?
1335 535 : if (hHFA->pProParameters != nullptr)
1336 0 : return (Eprj_ProParameters *)hHFA->pProParameters;
1337 :
1338 : // Get the HFA node.
1339 : HFAEntry *poMIEntry =
1340 535 : hHFA->papoBand[0]->poNode->GetNamedChild("Projection");
1341 535 : if (poMIEntry == nullptr)
1342 333 : 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 131 : CPLErr HFASetProParameters(HFAHandle hHFA, const Eprj_ProParameters *poPro)
1390 :
1391 : {
1392 : // Loop over bands, setting information on each one.
1393 300 : 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 169 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1398 169 : if (poMIEntry == nullptr)
1399 : {
1400 168 : poMIEntry = HFAEntry::New(hHFA, "Projection", "Eprj_ProParameters",
1401 168 : hHFA->papoBand[iBand]->poNode);
1402 : }
1403 :
1404 169 : poMIEntry->MarkDirty();
1405 :
1406 : // Ensure we have enough space for all the data.
1407 : // TODO(schwehr): Explain all these constants.
1408 169 : int nSize =
1409 169 : static_cast<int>(34 + 15 * 8 + 8 + strlen(poPro->proName) + 1 + 32 +
1410 169 : 8 + strlen(poPro->proSpheroid.sphereName) + 1);
1411 :
1412 169 : if (poPro->proExeName != nullptr)
1413 1 : nSize += static_cast<int>(strlen(poPro->proExeName) + 1);
1414 :
1415 169 : GByte *pabyData = poMIEntry->MakeData(nSize);
1416 169 : if (!pabyData)
1417 0 : return CE_Failure;
1418 :
1419 169 : poMIEntry->SetPosition();
1420 :
1421 : // Initialize the whole thing to zeros for a clean start.
1422 169 : memset(poMIEntry->GetData(), 0, poMIEntry->GetDataSize());
1423 :
1424 : // Write the various fields.
1425 169 : poMIEntry->SetIntField("proType", poPro->proType);
1426 :
1427 169 : poMIEntry->SetIntField("proNumber", poPro->proNumber);
1428 :
1429 169 : poMIEntry->SetStringField("proExeName", poPro->proExeName);
1430 169 : poMIEntry->SetStringField("proName", poPro->proName);
1431 169 : poMIEntry->SetIntField("proZone", poPro->proZone);
1432 169 : poMIEntry->SetDoubleField("proParams[0]", poPro->proParams[0]);
1433 169 : poMIEntry->SetDoubleField("proParams[1]", poPro->proParams[1]);
1434 169 : poMIEntry->SetDoubleField("proParams[2]", poPro->proParams[2]);
1435 169 : poMIEntry->SetDoubleField("proParams[3]", poPro->proParams[3]);
1436 169 : poMIEntry->SetDoubleField("proParams[4]", poPro->proParams[4]);
1437 169 : poMIEntry->SetDoubleField("proParams[5]", poPro->proParams[5]);
1438 169 : poMIEntry->SetDoubleField("proParams[6]", poPro->proParams[6]);
1439 169 : poMIEntry->SetDoubleField("proParams[7]", poPro->proParams[7]);
1440 169 : poMIEntry->SetDoubleField("proParams[8]", poPro->proParams[8]);
1441 169 : poMIEntry->SetDoubleField("proParams[9]", poPro->proParams[9]);
1442 169 : poMIEntry->SetDoubleField("proParams[10]", poPro->proParams[10]);
1443 169 : poMIEntry->SetDoubleField("proParams[11]", poPro->proParams[11]);
1444 169 : poMIEntry->SetDoubleField("proParams[12]", poPro->proParams[12]);
1445 169 : poMIEntry->SetDoubleField("proParams[13]", poPro->proParams[13]);
1446 169 : poMIEntry->SetDoubleField("proParams[14]", poPro->proParams[14]);
1447 169 : poMIEntry->SetStringField("proSpheroid.sphereName",
1448 169 : poPro->proSpheroid.sphereName);
1449 169 : poMIEntry->SetDoubleField("proSpheroid.a", poPro->proSpheroid.a);
1450 169 : poMIEntry->SetDoubleField("proSpheroid.b", poPro->proSpheroid.b);
1451 169 : poMIEntry->SetDoubleField("proSpheroid.eSquared",
1452 169 : poPro->proSpheroid.eSquared);
1453 169 : poMIEntry->SetDoubleField("proSpheroid.radius",
1454 169 : poPro->proSpheroid.radius);
1455 : }
1456 :
1457 131 : return CE_None;
1458 : }
1459 :
1460 : /************************************************************************/
1461 : /* HFAGetDatum() */
1462 : /************************************************************************/
1463 :
1464 535 : const Eprj_Datum *HFAGetDatum(HFAHandle hHFA)
1465 :
1466 : {
1467 535 : if (hHFA->nBands < 1)
1468 0 : return nullptr;
1469 :
1470 : // Do we already have it?
1471 535 : if (hHFA->pDatum != nullptr)
1472 0 : return (Eprj_Datum *)hHFA->pDatum;
1473 :
1474 : // Get the HFA node.
1475 : HFAEntry *poMIEntry =
1476 535 : hHFA->papoBand[0]->poNode->GetNamedChild("Projection.Datum");
1477 535 : if (poMIEntry == nullptr)
1478 333 : 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 131 : CPLErr HFASetDatum(HFAHandle hHFA, const Eprj_Datum *poDatum)
1514 :
1515 : {
1516 : // Loop over bands, setting information on each one.
1517 300 : 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 169 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1522 169 : 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 169 : HFAEntry *poDatumEntry = poProParams->GetNamedChild("Datum");
1530 169 : if (poDatumEntry == nullptr)
1531 : {
1532 : poDatumEntry =
1533 168 : HFAEntry::New(hHFA, "Datum", "Eprj_Datum", poProParams);
1534 : }
1535 :
1536 169 : poDatumEntry->MarkDirty();
1537 :
1538 : // Ensure we have enough space for all the data.
1539 : // TODO(schwehr): Explain constants.
1540 169 : int nSize =
1541 169 : static_cast<int>(26 + strlen(poDatum->datumname) + 1 + 7 * 8);
1542 :
1543 169 : if (poDatum->gridname != nullptr)
1544 30 : nSize += static_cast<int>(strlen(poDatum->gridname) + 1);
1545 :
1546 169 : GByte *pabyData = poDatumEntry->MakeData(nSize);
1547 169 : if (!pabyData)
1548 0 : return CE_Failure;
1549 :
1550 169 : poDatumEntry->SetPosition();
1551 :
1552 : // Initialize the whole thing to zeros for a clean start.
1553 169 : memset(poDatumEntry->GetData(), 0, poDatumEntry->GetDataSize());
1554 :
1555 : // Write the various fields.
1556 169 : poDatumEntry->SetStringField("datumname", poDatum->datumname);
1557 169 : poDatumEntry->SetIntField("type", poDatum->type);
1558 :
1559 169 : poDatumEntry->SetDoubleField("params[0]", poDatum->params[0]);
1560 169 : poDatumEntry->SetDoubleField("params[1]", poDatum->params[1]);
1561 169 : poDatumEntry->SetDoubleField("params[2]", poDatum->params[2]);
1562 169 : poDatumEntry->SetDoubleField("params[3]", poDatum->params[3]);
1563 169 : poDatumEntry->SetDoubleField("params[4]", poDatum->params[4]);
1564 169 : poDatumEntry->SetDoubleField("params[5]", poDatum->params[5]);
1565 169 : poDatumEntry->SetDoubleField("params[6]", poDatum->params[6]);
1566 :
1567 169 : poDatumEntry->SetStringField("gridname", poDatum->gridname);
1568 : }
1569 :
1570 131 : return CE_None;
1571 : }
1572 :
1573 : /************************************************************************/
1574 : /* HFAGetPCT() */
1575 : /* */
1576 : /* Read the PCT from a band, if it has one. */
1577 : /************************************************************************/
1578 :
1579 615 : CPLErr HFAGetPCT(HFAHandle hHFA, int nBand, int *pnColors, double **ppadfRed,
1580 : double **ppadfGreen, double **ppadfBlue, double **ppadfAlpha,
1581 : double **ppadfBins)
1582 :
1583 : {
1584 615 : if (nBand < 1 || nBand > hHFA->nBands)
1585 0 : return CE_Failure;
1586 :
1587 615 : return hHFA->papoBand[nBand - 1]->GetPCT(pnColors, ppadfRed, ppadfGreen,
1588 615 : 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 201 : HFAHandle HFACreateLL(const char *pszFilename)
1786 :
1787 : {
1788 : // Create the file in the file system.
1789 201 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
1790 201 : 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 198 : static_cast<HFAInfo_t *>(CPLCalloc(sizeof(HFAInfo_t), 1));
1800 :
1801 198 : psInfo->fp = fp;
1802 198 : psInfo->eAccess = HFA_Update;
1803 198 : psInfo->nXSize = 0;
1804 198 : psInfo->nYSize = 0;
1805 198 : psInfo->nBands = 0;
1806 198 : psInfo->papoBand = nullptr;
1807 198 : psInfo->pMapInfo = nullptr;
1808 198 : psInfo->pDatum = nullptr;
1809 198 : psInfo->pProParameters = nullptr;
1810 198 : psInfo->bTreeDirty = false;
1811 198 : psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
1812 198 : psInfo->pszPath = CPLStrdup(CPLGetPathSafe(pszFilename).c_str());
1813 :
1814 : // Write out the Ehfa_HeaderTag.
1815 198 : bool bRet = VSIFWriteL((void *)"EHFA_HEADER_TAG", 1, 16, fp) > 0;
1816 :
1817 198 : GInt32 nHeaderPos = 20;
1818 : HFAStandard(4, &nHeaderPos);
1819 198 : bRet &= VSIFWriteL(&nHeaderPos, 4, 1, fp) > 0;
1820 :
1821 : // Write the Ehfa_File node, locked in at offset 20.
1822 198 : GInt32 nVersion = 1;
1823 198 : GInt32 nFreeList = 0;
1824 198 : GInt32 nRootEntry = 0;
1825 198 : GInt16 nEntryHeaderLength = 128;
1826 198 : GInt32 nDictionaryPtr = 38;
1827 :
1828 198 : psInfo->nEntryHeaderLength = nEntryHeaderLength;
1829 198 : psInfo->nRootPos = 0;
1830 198 : psInfo->nDictionaryPos = nDictionaryPtr;
1831 198 : 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 198 : bRet &= VSIFWriteL(&nVersion, 4, 1, fp) > 0;
1840 198 : bRet &= VSIFWriteL(&nFreeList, 4, 1, fp) > 0;
1841 198 : bRet &= VSIFWriteL(&nRootEntry, 4, 1, fp) > 0;
1842 198 : bRet &= VSIFWriteL(&nEntryHeaderLength, 2, 1, fp) > 0;
1843 198 : 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 198 : int nDictLen = 0;
1850 :
1851 2178 : for (int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++)
1852 1980 : nDictLen += static_cast<int>(strlen(aszDefaultDD[iChunk]));
1853 :
1854 198 : psInfo->pszDictionary = static_cast<char *>(CPLMalloc(nDictLen + 1));
1855 198 : psInfo->pszDictionary[0] = '\0';
1856 :
1857 2178 : for (int iChunk = 0; aszDefaultDD[iChunk] != nullptr; iChunk++)
1858 1980 : strcat(psInfo->pszDictionary, aszDefaultDD[iChunk]);
1859 :
1860 396 : bRet &= VSIFWriteL((void *)psInfo->pszDictionary,
1861 198 : strlen(psInfo->pszDictionary) + 1, 1, fp) > 0;
1862 198 : if (!bRet)
1863 : {
1864 4 : CPL_IGNORE_RET_VAL(HFAClose(psInfo));
1865 4 : return nullptr;
1866 : }
1867 :
1868 194 : psInfo->poDictionary = new HFADictionary(psInfo->pszDictionary);
1869 :
1870 194 : psInfo->nEndOfFile = static_cast<GUInt32>(VSIFTellL(fp));
1871 :
1872 : // Create a root entry.
1873 194 : 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 194 : CPLString osExtension = CPLGetExtensionSafe(pszFilename);
1878 194 : if (!EQUAL(osExtension, "rrd") && !EQUAL(osExtension, "aux"))
1879 : {
1880 378 : CPLString osPath = CPLGetPathSafe(pszFilename);
1881 378 : CPLString osBasename = CPLGetBasenameSafe(pszFilename);
1882 : VSIStatBufL sStatBuf;
1883 378 : CPLString osSupFile = CPLFormCIFilenameSafe(osPath, osBasename, "rrd");
1884 :
1885 189 : if (VSIStatL(osSupFile, &sStatBuf) == 0)
1886 0 : VSIUnlink(osSupFile);
1887 :
1888 189 : osSupFile = CPLFormCIFilenameSafe(osPath, osBasename, "ige");
1889 :
1890 189 : if (VSIStatL(osSupFile, &sStatBuf) == 0)
1891 0 : VSIUnlink(osSupFile);
1892 : }
1893 :
1894 194 : 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 2843 : GUInt32 HFAAllocateSpace(HFAInfo_t *psInfo, GUInt32 nBytes)
1911 :
1912 : {
1913 : // TODO(schwehr): Check if this will wrap over 2GB limit.
1914 :
1915 2843 : psInfo->nEndOfFile += nBytes;
1916 2843 : 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 393 : CPLErr HFAFlush(HFAHandle hHFA)
1927 :
1928 : {
1929 393 : if (!hHFA->bTreeDirty && !hHFA->poDictionary->bDictionaryTextDirty)
1930 0 : return CE_None;
1931 :
1932 393 : CPLAssert(hHFA->poRoot != nullptr);
1933 :
1934 : // Flush HFAEntry tree to disk.
1935 393 : if (hHFA->bTreeDirty)
1936 : {
1937 393 : const CPLErr eErr = hHFA->poRoot->FlushToDisk();
1938 393 : if (eErr != CE_None)
1939 6 : return eErr;
1940 :
1941 387 : hHFA->bTreeDirty = false;
1942 : }
1943 :
1944 : // Flush Dictionary to disk.
1945 387 : GUInt32 nNewDictionaryPos = hHFA->nDictionaryPos;
1946 387 : bool bRet = true;
1947 387 : 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 585 : if (hHFA->nRootPos != hHFA->poRoot->GetFilePos() ||
1960 198 : nNewDictionaryPos != hHFA->nDictionaryPos)
1961 : {
1962 189 : GUInt32 nHeaderPos = 0;
1963 :
1964 189 : bRet &= VSIFSeekL(hHFA->fp, 16, SEEK_SET) >= 0;
1965 189 : bRet &= VSIFReadL(&nHeaderPos, sizeof(GInt32), 1, hHFA->fp) > 0;
1966 : HFAStandard(4, &nHeaderPos);
1967 :
1968 189 : GUInt32 nOffset = hHFA->poRoot->GetFilePos();
1969 189 : hHFA->nRootPos = nOffset;
1970 : HFAStandard(4, &nOffset);
1971 189 : bRet &= VSIFSeekL(hHFA->fp, nHeaderPos + 8, SEEK_SET) >= 0;
1972 189 : bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
1973 :
1974 189 : nOffset = nNewDictionaryPos;
1975 189 : hHFA->nDictionaryPos = nNewDictionaryPos;
1976 : HFAStandard(4, &nOffset);
1977 189 : bRet &= VSIFSeekL(hHFA->fp, nHeaderPos + 14, SEEK_SET) >= 0;
1978 189 : bRet &= VSIFWriteL(&nOffset, 4, 1, hHFA->fp) > 0;
1979 : }
1980 :
1981 387 : 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 241 : 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 241 : const char *pszLayerType =
2002 241 : bOverview ? "Eimg_Layer_SubSample" : "Eimg_Layer";
2003 :
2004 241 : 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 241 : const int nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2013 241 : const int nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2014 241 : const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
2015 : const int nBytesPerBlock =
2016 241 : (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
2017 :
2018 : // Create the Eimg_Layer for the band.
2019 : HFAEntry *poEimg_Layer =
2020 241 : HFAEntry::New(psInfo, pszLayerName, pszLayerType, poParent);
2021 :
2022 241 : poEimg_Layer->SetIntField("width", nXSize);
2023 241 : poEimg_Layer->SetIntField("height", nYSize);
2024 241 : poEimg_Layer->SetStringField("layerType", "athematic");
2025 241 : poEimg_Layer->SetIntField("pixelType", eDataType);
2026 241 : poEimg_Layer->SetIntField("blockWidth", nBlockSize);
2027 241 : 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 241 : if (!bCreateLargeRaster && !bDependentLayer)
2034 : {
2035 : HFAEntry *poEdms_State =
2036 230 : HFAEntry::New(psInfo, "RasterDMS", "Edms_State", poEimg_Layer);
2037 :
2038 : // TODO(schwehr): Explain constants.
2039 230 : const int nDmsSize = 14 * nBlocks + 38;
2040 230 : GByte *pabyData = poEdms_State->MakeData(nDmsSize);
2041 :
2042 : // Set some simple values.
2043 230 : poEdms_State->SetIntField("numvirtualblocks", nBlocks);
2044 230 : poEdms_State->SetIntField("numobjectsperblock",
2045 : nBlockSize * nBlockSize);
2046 230 : poEdms_State->SetIntField("nextobjectnum",
2047 230 : nBlockSize * nBlockSize * nBlocks);
2048 :
2049 : // Is file compressed or not?
2050 230 : if (bCreateCompressed)
2051 : {
2052 10 : poEdms_State->SetStringField("compressionType", "RLC compression");
2053 : }
2054 : else
2055 : {
2056 220 : poEdms_State->SetStringField("compressionType", "no compression");
2057 : }
2058 :
2059 : // We need to hardcode file offset into the data, so locate it now.
2060 230 : poEdms_State->SetPosition();
2061 :
2062 : // Set block info headers.
2063 :
2064 : // Blockinfo count.
2065 230 : GUInt32 nValue = nBlocks;
2066 : HFAStandard(4, &nValue);
2067 230 : memcpy(pabyData + 14, &nValue, 4);
2068 :
2069 : // Blockinfo position.
2070 230 : nValue = poEdms_State->GetDataPos() + 22;
2071 : HFAStandard(4, &nValue);
2072 230 : memcpy(pabyData + 18, &nValue, 4);
2073 :
2074 : // Set each blockinfo.
2075 728 : for (int iBlock = 0; iBlock < nBlocks; iBlock++)
2076 : {
2077 498 : int nOffset = 22 + 14 * iBlock;
2078 :
2079 : // fileCode.
2080 498 : GInt16 nValue16 = 0;
2081 : HFAStandard(2, &nValue16);
2082 498 : memcpy(pabyData + nOffset, &nValue16, 2);
2083 :
2084 : // Offset.
2085 498 : if (bCreateCompressed)
2086 : {
2087 : // Flag it with zero offset. Allocate space when we compress it.
2088 18 : nValue = 0;
2089 : }
2090 : else
2091 : {
2092 480 : nValue = HFAAllocateSpace(psInfo, nBytesPerBlock);
2093 : }
2094 : HFAStandard(4, &nValue);
2095 498 : memcpy(pabyData + nOffset + 2, &nValue, 4);
2096 :
2097 : // Size.
2098 498 : if (bCreateCompressed)
2099 : {
2100 : // Flag with zero size. Don't know until we compress it.
2101 18 : nValue = 0;
2102 : }
2103 : else
2104 : {
2105 480 : nValue = nBytesPerBlock;
2106 : }
2107 : HFAStandard(4, &nValue);
2108 498 : memcpy(pabyData + nOffset + 6, &nValue, 4);
2109 :
2110 : // logValid (false).
2111 498 : nValue16 = 0;
2112 : HFAStandard(2, &nValue16);
2113 498 : memcpy(pabyData + nOffset + 10, &nValue16, 2);
2114 :
2115 : // compressionType.
2116 498 : if (bCreateCompressed)
2117 18 : nValue16 = 1;
2118 : else
2119 480 : nValue16 = 0;
2120 :
2121 : HFAStandard(2, &nValue16);
2122 498 : memcpy(pabyData + nOffset + 12, &nValue16, 2);
2123 230 : }
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 241 : char chBandType = '\0';
2165 :
2166 241 : if (eDataType == EPT_u1)
2167 2 : chBandType = '1';
2168 239 : else if (eDataType == EPT_u2)
2169 2 : chBandType = '2';
2170 237 : else if (eDataType == EPT_u4)
2171 2 : chBandType = '4';
2172 235 : else if (eDataType == EPT_u8)
2173 147 : chBandType = 'c';
2174 88 : else if (eDataType == EPT_s8)
2175 4 : chBandType = 'C';
2176 84 : else if (eDataType == EPT_u16)
2177 13 : chBandType = 's';
2178 71 : else if (eDataType == EPT_s16)
2179 9 : chBandType = 'S';
2180 62 : 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 52 : else if (eDataType == EPT_s32)
2185 11 : chBandType = 'L';
2186 41 : else if (eDataType == EPT_f32)
2187 11 : 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 241 : char szLDict[128] = {};
2203 241 : snprintf(szLDict, sizeof(szLDict), "{%d:%cdata,}RasterDMS,.",
2204 : nBlockSize * nBlockSize, chBandType);
2205 :
2206 : HFAEntry *poEhfa_Layer =
2207 241 : HFAEntry::New(psInfo, "Ehfa_Layer", "Ehfa_Layer", poEimg_Layer);
2208 241 : poEhfa_Layer->MakeData();
2209 241 : poEhfa_Layer->SetPosition();
2210 : const GUInt32 nLDict =
2211 241 : HFAAllocateSpace(psInfo, static_cast<GUInt32>(strlen(szLDict) + 1));
2212 :
2213 241 : poEhfa_Layer->SetStringField("type", "raster");
2214 241 : poEhfa_Layer->SetIntField("dictionaryPtr", nLDict);
2215 :
2216 241 : bool bRet = VSIFSeekL(psInfo->fp, nLDict, SEEK_SET) >= 0;
2217 241 : bRet &= VSIFWriteL((void *)szLDict, strlen(szLDict) + 1, 1, psInfo->fp) > 0;
2218 :
2219 241 : return bRet;
2220 : }
2221 :
2222 : /************************************************************************/
2223 : /* HFACreate() */
2224 : /************************************************************************/
2225 :
2226 201 : HFAHandle HFACreate(const char *pszFilename, int nXSize, int nYSize, int nBands,
2227 : EPTType eDataType, char **papszOptions)
2228 :
2229 : {
2230 201 : if (nXSize == 0 || nYSize == 0)
2231 : {
2232 1 : CPLError(CE_Failure, CPLE_NotSupported,
2233 : "nXSize == 0 || nYSize == 0 not supported");
2234 1 : return nullptr;
2235 : }
2236 200 : int nBlockSize = 64;
2237 200 : const char *pszValue = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
2238 :
2239 200 : if (pszValue != nullptr)
2240 : {
2241 1 : nBlockSize = atoi(pszValue);
2242 : // Check for sane values.
2243 2 : if (nBlockSize == 0 ||
2244 1 : (((nBlockSize < 32) || (nBlockSize > 2048)) &&
2245 0 : !CPLTestBool(CPLGetConfigOption("FORCE_BLOCKSIZE", "NO"))))
2246 : {
2247 0 : if (nBlockSize != 0)
2248 0 : CPLError(CE_Warning, CPLE_AppDefined, "Forcing BLOCKSIZE to %d",
2249 : 64);
2250 0 : nBlockSize = 64;
2251 : }
2252 : }
2253 200 : bool bCreateLargeRaster = CPLFetchBool(papszOptions, "USE_SPILL", false);
2254 396 : bool bCreateCompressed = CPLFetchBool(papszOptions, "COMPRESS", false) ||
2255 196 : CPLFetchBool(papszOptions, "COMPRESSED", false);
2256 200 : const bool bCreateAux = CPLFetchBool(papszOptions, "AUX", false);
2257 :
2258 200 : char *pszFullFilename = nullptr;
2259 200 : char *pszRawFilename = nullptr;
2260 :
2261 : // Work out some details about the tiling scheme.
2262 200 : const int nBlocksPerRow = DIV_ROUND_UP(nXSize, nBlockSize);
2263 200 : const int nBlocksPerColumn = DIV_ROUND_UP(nYSize, nBlockSize);
2264 200 : if (nBlocksPerRow > INT_MAX / nBlocksPerColumn)
2265 : {
2266 0 : CPLError(CE_Failure, CPLE_NotSupported, "Too many blocks");
2267 0 : return nullptr;
2268 : }
2269 200 : const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
2270 : const GInt64 nBytesPerBlock64 =
2271 400 : (static_cast<GInt64>(nBlockSize) * nBlockSize *
2272 200 : HFAGetDataTypeBits(eDataType) +
2273 : 7) /
2274 200 : 8;
2275 200 : if (nBytesPerBlock64 > INT_MAX)
2276 : {
2277 0 : CPLError(CE_Failure, CPLE_NotSupported, "Too large block");
2278 0 : return nullptr;
2279 : }
2280 200 : const int nBytesPerBlock = static_cast<int>(nBytesPerBlock64);
2281 :
2282 : // Create the low level structure.
2283 200 : HFAHandle psInfo = HFACreateLL(pszFilename);
2284 200 : if (psInfo == nullptr)
2285 7 : return nullptr;
2286 :
2287 : // Create the DependentFile node if requested.
2288 : const char *pszDependentFile =
2289 193 : CSLFetchNameValue(papszOptions, "DEPENDENT_FILE");
2290 :
2291 193 : if (pszDependentFile != nullptr)
2292 : {
2293 4 : HFAEntry *poDF = HFAEntry::New(psInfo, "DependentFile",
2294 : "Eimg_DependentFile", psInfo->poRoot);
2295 :
2296 4 : poDF->MakeData(static_cast<int>(strlen(pszDependentFile) + 50));
2297 4 : poDF->SetPosition();
2298 4 : poDF->SetStringField("dependent.string", pszDependentFile);
2299 : }
2300 :
2301 193 : CPLDebug("HFACreate",
2302 : "Blocks per row %d, blocks per column %d, "
2303 : "total number of blocks %d, bytes per block %d.",
2304 : nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock);
2305 :
2306 : // Check whether we should create external large file with
2307 : // image. We create a spill file if the amount of imagery is
2308 : // close to 2GB. We don't check the amount of auxiliary
2309 : // information, so in theory if there were an awful lot of
2310 : // non-imagery data our approximate size could be smaller than
2311 : // the file will actually we be. We leave room for 10MB of
2312 : // auxiliary data.
2313 : // We can also force spill file creation using option
2314 : // SPILL_FILE=YES.
2315 193 : const double dfApproxSize = static_cast<double>(nBytesPerBlock) *
2316 193 : static_cast<double>(nBlocks) *
2317 193 : static_cast<double>(nBands) +
2318 : 10000000.0;
2319 :
2320 193 : if (dfApproxSize > 2147483648.0 && !bCreateAux)
2321 0 : bCreateLargeRaster = true;
2322 :
2323 : // Erdas Imagine creates this entry even if an external spill file is used.
2324 193 : if (!bCreateAux)
2325 : {
2326 189 : HFAEntry *poImgFormat = HFAEntry::New(
2327 : psInfo, "IMGFormatInfo", "ImgFormatInfo831", psInfo->poRoot);
2328 189 : poImgFormat->MakeData();
2329 189 : if (bCreateLargeRaster)
2330 : {
2331 7 : poImgFormat->SetIntField("spaceUsedForRasterData", 0);
2332 : // Can't be compressed if we are creating a spillfile.
2333 7 : bCreateCompressed = false;
2334 : }
2335 : else
2336 : {
2337 182 : poImgFormat->SetIntField("spaceUsedForRasterData",
2338 182 : nBytesPerBlock * nBlocks * nBands);
2339 : }
2340 : }
2341 :
2342 : // Create external file and write its header.
2343 193 : GIntBig nValidFlagsOffset = 0;
2344 193 : GIntBig nDataOffset = 0;
2345 :
2346 193 : if (bCreateLargeRaster)
2347 : {
2348 7 : if (!HFACreateSpillStack(psInfo, nXSize, nYSize, nBands, nBlockSize,
2349 : eDataType, &nValidFlagsOffset, &nDataOffset))
2350 : {
2351 0 : CPLFree(pszRawFilename);
2352 0 : CPLFree(pszFullFilename);
2353 0 : return nullptr;
2354 : }
2355 : }
2356 :
2357 : // Create each band (layer).
2358 419 : for (int iBand = 0; iBand < nBands; iBand++)
2359 : {
2360 230 : char szName[128] = {};
2361 :
2362 230 : snprintf(szName, sizeof(szName), "Layer_%d", iBand + 1);
2363 :
2364 230 : if (!HFACreateLayer(psInfo, psInfo->poRoot, szName, FALSE, nBlockSize,
2365 : bCreateCompressed, bCreateLargeRaster, bCreateAux,
2366 : nXSize, nYSize, eDataType, papszOptions,
2367 : nValidFlagsOffset, nDataOffset, nBands, iBand))
2368 : {
2369 4 : CPL_IGNORE_RET_VAL(HFAClose(psInfo));
2370 4 : return nullptr;
2371 : }
2372 : }
2373 :
2374 : // Initialize the band information.
2375 189 : HFAParseBandInfo(psInfo);
2376 :
2377 189 : return psInfo;
2378 : }
2379 :
2380 : /************************************************************************/
2381 : /* HFACreateOverview() */
2382 : /* */
2383 : /* Create an overview layer object for a band. */
2384 : /************************************************************************/
2385 :
2386 11 : int HFACreateOverview(HFAHandle hHFA, int nBand, int nOverviewLevel,
2387 : const char *pszResampling)
2388 :
2389 : {
2390 11 : if (nBand < 1 || nBand > hHFA->nBands)
2391 0 : return -1;
2392 :
2393 11 : HFABand *poBand = hHFA->papoBand[nBand - 1];
2394 11 : return poBand->CreateOverview(nOverviewLevel, pszResampling);
2395 : }
2396 :
2397 : /************************************************************************/
2398 : /* HFAGetMetadata() */
2399 : /* */
2400 : /* Read metadata structured in a table called GDAL_MetaData. */
2401 : /************************************************************************/
2402 :
2403 1150 : char **HFAGetMetadata(HFAHandle hHFA, int nBand)
2404 :
2405 : {
2406 1150 : HFAEntry *poTable = nullptr;
2407 :
2408 1150 : if (nBand > 0 && nBand <= hHFA->nBands)
2409 615 : poTable = hHFA->papoBand[nBand - 1]->poNode->GetChild();
2410 535 : else if (nBand == 0)
2411 535 : poTable = hHFA->poRoot->GetChild();
2412 : else
2413 0 : return nullptr;
2414 :
2415 4721 : for (; poTable != nullptr && !EQUAL(poTable->GetName(), "GDAL_MetaData");
2416 3571 : poTable = poTable->GetNext())
2417 : {
2418 : }
2419 :
2420 1150 : if (poTable == nullptr || !EQUAL(poTable->GetType(), "Edsc_Table"))
2421 1055 : return nullptr;
2422 :
2423 95 : if (poTable->GetIntField("numRows") != 1)
2424 : {
2425 0 : CPLDebug("HFADataset", "GDAL_MetaData.numRows = %d, expected 1!",
2426 : poTable->GetIntField("numRows"));
2427 0 : return nullptr;
2428 : }
2429 :
2430 : // Loop over each column. Each column will be one metadata
2431 : // entry, with the title being the key, and the row value being
2432 : // the value. There is only ever one row in GDAL_MetaData tables.
2433 95 : char **papszMD = nullptr;
2434 :
2435 286 : for (HFAEntry *poColumn = poTable->GetChild(); poColumn != nullptr;
2436 191 : poColumn = poColumn->GetNext())
2437 : {
2438 : // Skip the #Bin_Function# entry.
2439 191 : if (STARTS_WITH_CI(poColumn->GetName(), "#"))
2440 95 : continue;
2441 :
2442 96 : const char *pszValue = poColumn->GetStringField("dataType");
2443 96 : if (pszValue == nullptr || !EQUAL(pszValue, "string"))
2444 0 : continue;
2445 :
2446 96 : const int columnDataPtr = poColumn->GetIntField("columnDataPtr");
2447 96 : if (columnDataPtr <= 0)
2448 0 : continue;
2449 :
2450 : // Read up to nMaxNumChars bytes from the indicated location.
2451 : // allocate required space temporarily nMaxNumChars should have been
2452 : // set by GDAL originally so we should trust it, but who knows.
2453 96 : const int nMaxNumChars = poColumn->GetIntField("maxNumChars");
2454 :
2455 96 : if (nMaxNumChars <= 0)
2456 : {
2457 0 : papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), "");
2458 : }
2459 : else
2460 : {
2461 : char *pszMDValue =
2462 96 : static_cast<char *>(VSI_MALLOC_VERBOSE(nMaxNumChars));
2463 96 : if (pszMDValue == nullptr)
2464 : {
2465 0 : continue;
2466 : }
2467 :
2468 96 : if (VSIFSeekL(hHFA->fp, columnDataPtr, SEEK_SET) != 0)
2469 : {
2470 0 : CPLFree(pszMDValue);
2471 0 : continue;
2472 : }
2473 :
2474 : const int nMDBytes = static_cast<int>(
2475 96 : VSIFReadL(pszMDValue, 1, nMaxNumChars, hHFA->fp));
2476 96 : if (nMDBytes == 0)
2477 : {
2478 0 : CPLFree(pszMDValue);
2479 0 : continue;
2480 : }
2481 :
2482 96 : pszMDValue[nMaxNumChars - 1] = '\0';
2483 :
2484 96 : papszMD = CSLSetNameValue(papszMD, poColumn->GetName(), pszMDValue);
2485 96 : CPLFree(pszMDValue);
2486 : }
2487 : }
2488 :
2489 95 : return papszMD;
2490 : }
2491 :
2492 : /************************************************************************/
2493 : /* HFASetGDALMetadata() */
2494 : /* */
2495 : /* This function is used to set metadata in a table called */
2496 : /* GDAL_MetaData. It is called by HFASetMetadata() for all */
2497 : /* metadata items that aren't some specific supported */
2498 : /* information (like histogram or stats info). */
2499 : /************************************************************************/
2500 :
2501 42 : static CPLErr HFASetGDALMetadata(HFAHandle hHFA, int nBand, char **papszMD)
2502 :
2503 : {
2504 42 : if (papszMD == nullptr)
2505 0 : return CE_None;
2506 :
2507 42 : HFAEntry *poNode = nullptr;
2508 :
2509 42 : if (nBand > 0 && nBand <= hHFA->nBands)
2510 6 : poNode = hHFA->papoBand[nBand - 1]->poNode;
2511 36 : else if (nBand == 0)
2512 36 : poNode = hHFA->poRoot;
2513 : else
2514 0 : return CE_Failure;
2515 :
2516 : // Create the Descriptor table.
2517 : // Check we have no table with this name already.
2518 42 : HFAEntry *poEdsc_Table = poNode->GetNamedChild("GDAL_MetaData");
2519 :
2520 43 : if (poEdsc_Table == nullptr ||
2521 1 : !EQUAL(poEdsc_Table->GetType(), "Edsc_Table"))
2522 : poEdsc_Table =
2523 41 : HFAEntry::New(hHFA, "GDAL_MetaData", "Edsc_Table", poNode);
2524 :
2525 42 : poEdsc_Table->SetIntField("numrows", 1);
2526 :
2527 : // Create the Binning function node. Do we really need this though?
2528 : // Check it doesn't exist already.
2529 : HFAEntry *poEdsc_BinFunction =
2530 42 : poEdsc_Table->GetNamedChild("#Bin_Function#");
2531 :
2532 43 : if (poEdsc_BinFunction == nullptr ||
2533 1 : !EQUAL(poEdsc_BinFunction->GetType(), "Edsc_BinFunction"))
2534 41 : poEdsc_BinFunction = HFAEntry::New(hHFA, "#Bin_Function#",
2535 : "Edsc_BinFunction", poEdsc_Table);
2536 :
2537 : // Because of the BaseData we have to hardcode the size.
2538 42 : poEdsc_BinFunction->MakeData(30);
2539 :
2540 42 : poEdsc_BinFunction->SetIntField("numBins", 1);
2541 42 : poEdsc_BinFunction->SetStringField("binFunction", "direct");
2542 42 : poEdsc_BinFunction->SetDoubleField("minLimit", 0.0);
2543 42 : poEdsc_BinFunction->SetDoubleField("maxLimit", 0.0);
2544 :
2545 : // Process each metadata item as a separate column.
2546 42 : bool bRet = true;
2547 85 : for (int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++)
2548 : {
2549 43 : char *pszKey = nullptr;
2550 43 : const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
2551 43 : if (pszValue == nullptr)
2552 0 : continue;
2553 :
2554 : // Create the Edsc_Column.
2555 : // Check it doesn't exist already.
2556 43 : HFAEntry *poEdsc_Column = poEdsc_Table->GetNamedChild(pszKey);
2557 :
2558 44 : if (poEdsc_Column == nullptr ||
2559 1 : !EQUAL(poEdsc_Column->GetType(), "Edsc_Column"))
2560 : poEdsc_Column =
2561 42 : HFAEntry::New(hHFA, pszKey, "Edsc_Column", poEdsc_Table);
2562 :
2563 43 : poEdsc_Column->SetIntField("numRows", 1);
2564 43 : poEdsc_Column->SetStringField("dataType", "string");
2565 43 : poEdsc_Column->SetIntField("maxNumChars",
2566 43 : static_cast<GUInt32>(strlen(pszValue) + 1));
2567 :
2568 : // Write the data out.
2569 : const int nOffset =
2570 43 : HFAAllocateSpace(hHFA, static_cast<GUInt32>(strlen(pszValue) + 1));
2571 :
2572 43 : poEdsc_Column->SetIntField("columnDataPtr", nOffset);
2573 :
2574 43 : bRet &= VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) >= 0;
2575 43 : bRet &=
2576 43 : VSIFWriteL((void *)pszValue, strlen(pszValue) + 1, 1, hHFA->fp) > 0;
2577 :
2578 43 : CPLFree(pszKey);
2579 : }
2580 :
2581 42 : return bRet ? CE_None : CE_Failure;
2582 : }
2583 :
2584 : /************************************************************************/
2585 : /* HFASetMetadata() */
2586 : /************************************************************************/
2587 :
2588 65 : CPLErr HFASetMetadata(HFAHandle hHFA, int nBand, char **papszMD)
2589 :
2590 : {
2591 65 : char **papszGDALMD = nullptr;
2592 :
2593 65 : if (CSLCount(papszMD) == 0)
2594 0 : return CE_None;
2595 :
2596 65 : HFAEntry *poNode = nullptr;
2597 :
2598 65 : if (nBand > 0 && nBand <= hHFA->nBands)
2599 15 : poNode = hHFA->papoBand[nBand - 1]->poNode;
2600 50 : else if (nBand == 0)
2601 50 : poNode = hHFA->poRoot;
2602 : else
2603 0 : return CE_Failure;
2604 : #ifdef DEBUG
2605 : // To please Clang Static Analyzer (CSA).
2606 65 : if (poNode == nullptr)
2607 : {
2608 0 : CPLAssert(false);
2609 : return CE_Failure;
2610 : }
2611 : #endif
2612 :
2613 : // Check if the Metadata is an "known" entity which should be
2614 : // stored in a better place.
2615 65 : char *pszBinValues = nullptr;
2616 65 : bool bCreatedHistogramParameters = false;
2617 65 : bool bCreatedStatistics = false;
2618 65 : const char *const *pszAuxMetaData = GetHFAAuxMetaDataList();
2619 : // Check each metadata item.
2620 229 : for (int iColumn = 0; papszMD[iColumn] != nullptr; iColumn++)
2621 : {
2622 164 : char *pszKey = nullptr;
2623 164 : const char *pszValue = CPLParseNameValue(papszMD[iColumn], &pszKey);
2624 164 : if (pszValue == nullptr)
2625 28 : continue;
2626 :
2627 : // Know look if its known.
2628 150 : int i = 0; // Used after for.
2629 1408 : for (; pszAuxMetaData[i] != nullptr; i += 4)
2630 : {
2631 1358 : if (EQUALN(pszAuxMetaData[i + 2], pszKey, strlen(pszKey)))
2632 100 : break;
2633 : }
2634 150 : if (pszAuxMetaData[i] != nullptr)
2635 : {
2636 : // Found one, get the right entry.
2637 100 : HFAEntry *poEntry = nullptr;
2638 :
2639 100 : if (strlen(pszAuxMetaData[i]) > 0)
2640 90 : poEntry = poNode->GetNamedChild(pszAuxMetaData[i]);
2641 : else
2642 10 : poEntry = poNode;
2643 :
2644 100 : if (poEntry == nullptr && strlen(pszAuxMetaData[i + 3]) > 0)
2645 : {
2646 : // Child does not yet exist --> create it,
2647 28 : poEntry = HFAEntry::New(hHFA, pszAuxMetaData[i],
2648 14 : pszAuxMetaData[i + 3], poNode);
2649 :
2650 14 : if (STARTS_WITH_CI(pszAuxMetaData[i], "Statistics"))
2651 8 : bCreatedStatistics = true;
2652 :
2653 14 : if (STARTS_WITH_CI(pszAuxMetaData[i], "HistogramParameters"))
2654 : {
2655 : // A bit nasty. Need to set the string field for the object
2656 : // first because the SetStringField sets the count for the
2657 : // object BinFunction to the length of the string.
2658 6 : poEntry->MakeData(70);
2659 6 : poEntry->SetStringField("BinFunction.binFunctionType",
2660 : "direct");
2661 :
2662 6 : bCreatedHistogramParameters = true;
2663 : }
2664 : }
2665 100 : if (poEntry == nullptr)
2666 : {
2667 14 : CPLFree(pszKey);
2668 14 : continue;
2669 : }
2670 :
2671 86 : const char *pszFieldName = pszAuxMetaData[i + 1] + 1;
2672 86 : switch (pszAuxMetaData[i + 1][0])
2673 : {
2674 67 : case 'd':
2675 : {
2676 67 : double dfValue = CPLAtof(pszValue);
2677 67 : poEntry->SetDoubleField(pszFieldName, dfValue);
2678 : }
2679 67 : break;
2680 9 : case 'i':
2681 : case 'l':
2682 : {
2683 9 : int nValue = atoi(pszValue);
2684 9 : poEntry->SetIntField(pszFieldName, nValue);
2685 : }
2686 9 : break;
2687 10 : case 's':
2688 : case 'e':
2689 : {
2690 10 : poEntry->SetStringField(pszFieldName, pszValue);
2691 : }
2692 10 : break;
2693 0 : default:
2694 0 : CPLAssert(false);
2695 : }
2696 : }
2697 50 : else if (STARTS_WITH_CI(pszKey, "STATISTICS_HISTOBINVALUES"))
2698 : {
2699 7 : CPLFree(pszBinValues);
2700 7 : pszBinValues = CPLStrdup(pszValue);
2701 : }
2702 : else
2703 : {
2704 43 : papszGDALMD = CSLAddString(papszGDALMD, papszMD[iColumn]);
2705 : }
2706 :
2707 136 : CPLFree(pszKey);
2708 : }
2709 :
2710 : // Special case to write out the histogram.
2711 65 : bool bRet = true;
2712 65 : if (pszBinValues != nullptr)
2713 : {
2714 7 : HFAEntry *poEntry = poNode->GetNamedChild("HistogramParameters");
2715 7 : if (poEntry != nullptr && bCreatedHistogramParameters)
2716 : {
2717 : // If this node exists we have added Histogram data -- complete with
2718 : // some defaults.
2719 6 : poEntry->SetIntField("SkipFactorX", 1);
2720 6 : poEntry->SetIntField("SkipFactorY", 1);
2721 :
2722 6 : const int nNumBins = poEntry->GetIntField("BinFunction.numBins");
2723 : const double dMinLimit =
2724 6 : poEntry->GetDoubleField("BinFunction.minLimit");
2725 : const double dMaxLimit =
2726 6 : poEntry->GetDoubleField("BinFunction.maxLimit");
2727 :
2728 : // Fill the descriptor table - check it isn't there already.
2729 6 : poEntry = poNode->GetNamedChild("Descriptor_Table");
2730 6 : if (poEntry == nullptr || !EQUAL(poEntry->GetType(), "Edsc_Table"))
2731 1 : poEntry = HFAEntry::New(hHFA, "Descriptor_Table", "Edsc_Table",
2732 : poNode);
2733 :
2734 6 : poEntry->SetIntField("numRows", nNumBins);
2735 :
2736 : // Bin function.
2737 6 : HFAEntry *poBinFunc = poEntry->GetNamedChild("#Bin_Function#");
2738 10 : if (poBinFunc == nullptr ||
2739 4 : !EQUAL(poBinFunc->GetType(), "Edsc_BinFunction"))
2740 2 : poBinFunc = HFAEntry::New(hHFA, "#Bin_Function#",
2741 : "Edsc_BinFunction", poEntry);
2742 :
2743 6 : poBinFunc->MakeData(30);
2744 6 : poBinFunc->SetIntField("numBins", nNumBins);
2745 6 : poBinFunc->SetDoubleField("minLimit", dMinLimit);
2746 6 : poBinFunc->SetDoubleField("maxLimit", dMaxLimit);
2747 : // Direct for thematic layers, linear otherwise.
2748 6 : if (STARTS_WITH_CI(poNode->GetStringField("layerType"), "thematic"))
2749 1 : poBinFunc->SetStringField("binFunctionType", "direct");
2750 : else
2751 5 : poBinFunc->SetStringField("binFunctionType", "linear");
2752 :
2753 : // We need a child named histogram.
2754 6 : HFAEntry *poHisto = poEntry->GetNamedChild("Histogram");
2755 6 : if (poHisto == nullptr || !EQUAL(poHisto->GetType(), "Edsc_Column"))
2756 : poHisto =
2757 1 : HFAEntry::New(hHFA, "Histogram", "Edsc_Column", poEntry);
2758 :
2759 6 : poHisto->SetIntField("numRows", nNumBins);
2760 : // Allocate space for the bin values.
2761 6 : GUInt32 nOffset = HFAAllocateSpace(hHFA, nNumBins * 8);
2762 6 : poHisto->SetIntField("columnDataPtr", nOffset);
2763 6 : poHisto->SetStringField("dataType", "real");
2764 6 : poHisto->SetIntField("maxNumChars", 0);
2765 : // Write out histogram data.
2766 6 : char *pszWork = pszBinValues;
2767 1531 : for (int nBin = 0; nBin < nNumBins; ++nBin)
2768 : {
2769 1525 : char *pszEnd = strchr(pszWork, '|');
2770 1525 : if (pszEnd != nullptr)
2771 : {
2772 1525 : *pszEnd = 0;
2773 1525 : bRet &=
2774 1525 : VSIFSeekL(hHFA->fp, nOffset + 8 * nBin, SEEK_SET) >= 0;
2775 1525 : double nValue = CPLAtof(pszWork);
2776 : HFAStandard(8, &nValue);
2777 :
2778 1525 : bRet &= VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
2779 1525 : pszWork = pszEnd + 1;
2780 : }
2781 6 : }
2782 : }
2783 1 : else if (poEntry != nullptr)
2784 : {
2785 : // In this case, there are HistogramParameters present, but we did
2786 : // not create them. However, we might be modifying them, in the case
2787 : // where the data has changed and the histogram counts need to be
2788 : // updated. It could be worse than that, but that is all we are
2789 : // going to cope with for now. We are assuming that we did not
2790 : // change any of the other stuff, like skip factors and so
2791 : // forth. The main need for this case is for programs (such as
2792 : // Imagine itself) which will happily modify the pixel values
2793 : // without re-calculating the histogram counts.
2794 1 : int nNumBins = poEntry->GetIntField("BinFunction.numBins");
2795 : HFAEntry *poEntryDescrTbl =
2796 1 : poNode->GetNamedChild("Descriptor_Table");
2797 1 : HFAEntry *poHisto = nullptr;
2798 1 : if (poEntryDescrTbl != nullptr)
2799 : {
2800 1 : poHisto = poEntryDescrTbl->GetNamedChild("Histogram");
2801 : }
2802 1 : if (poHisto != nullptr)
2803 : {
2804 1 : int nOffset = poHisto->GetIntField("columnDataPtr");
2805 : // Write out histogram data.
2806 1 : char *pszWork = pszBinValues;
2807 :
2808 : // Check whether histogram counts were written as int or double
2809 1 : bool bCountIsInt = true;
2810 1 : const char *pszDataType = poHisto->GetStringField("dataType");
2811 1 : if (STARTS_WITH_CI(pszDataType, "real"))
2812 : {
2813 1 : bCountIsInt = false;
2814 : }
2815 257 : for (int nBin = 0; nBin < nNumBins; ++nBin)
2816 : {
2817 256 : char *pszEnd = strchr(pszWork, '|');
2818 256 : if (pszEnd != nullptr)
2819 : {
2820 256 : *pszEnd = 0;
2821 256 : if (bCountIsInt)
2822 : {
2823 : // Histogram counts were written as ints, so
2824 : // re-write them the same way.
2825 0 : bRet &= VSIFSeekL(hHFA->fp, nOffset + 4 * nBin,
2826 0 : SEEK_SET) >= 0;
2827 0 : int nValue = atoi(pszWork);
2828 : HFAStandard(4, &nValue);
2829 0 : bRet &=
2830 0 : VSIFWriteL((void *)&nValue, 4, 1, hHFA->fp) > 0;
2831 : }
2832 : else
2833 : {
2834 : // Histogram were written as doubles, as is now the
2835 : // default behavior.
2836 256 : bRet &= VSIFSeekL(hHFA->fp, nOffset + 8 * nBin,
2837 256 : SEEK_SET) >= 0;
2838 256 : double nValue = CPLAtof(pszWork);
2839 : HFAStandard(8, &nValue);
2840 256 : bRet &=
2841 256 : VSIFWriteL((void *)&nValue, 8, 1, hHFA->fp) > 0;
2842 : }
2843 256 : pszWork = pszEnd + 1;
2844 : }
2845 : }
2846 : }
2847 : }
2848 7 : CPLFree(pszBinValues);
2849 : }
2850 :
2851 : // If we created a statistics node then try to create a
2852 : // StatisticsParameters node too.
2853 65 : if (bCreatedStatistics)
2854 : {
2855 : HFAEntry *poEntry =
2856 8 : HFAEntry::New(hHFA, "StatisticsParameters",
2857 : "Eimg_StatisticsParameters830", poNode);
2858 :
2859 8 : poEntry->MakeData(70);
2860 : // poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
2861 :
2862 8 : poEntry->SetIntField("SkipFactorX", 1);
2863 8 : poEntry->SetIntField("SkipFactorY", 1);
2864 : }
2865 :
2866 : // Write out metadata items without a special place.
2867 65 : if (bRet && CSLCount(papszGDALMD) != 0)
2868 : {
2869 42 : CPLErr eErr = HFASetGDALMetadata(hHFA, nBand, papszGDALMD);
2870 :
2871 42 : CSLDestroy(papszGDALMD);
2872 42 : return eErr;
2873 : }
2874 : else
2875 : {
2876 23 : CSLDestroy(papszGDALMD);
2877 23 : return CE_Failure;
2878 : }
2879 : }
2880 :
2881 : /************************************************************************/
2882 : /* HFAGetIGEFilename() */
2883 : /* */
2884 : /* Returns the .ige filename if one is associated with this */
2885 : /* object. For files not newly created we need to scan the */
2886 : /* bands for spill files. Presumably there will only be one. */
2887 : /* */
2888 : /* NOTE: Returns full path, not just the filename portion. */
2889 : /************************************************************************/
2890 :
2891 113 : std::string HFAGetIGEFilename(HFAHandle hHFA)
2892 :
2893 : {
2894 113 : if (hHFA->pszIGEFilename == nullptr)
2895 : {
2896 : std::vector<HFAEntry *> apoDMSList =
2897 220 : hHFA->poRoot->FindChildren(nullptr, "ImgExternalRaster");
2898 :
2899 110 : HFAEntry *poDMS = apoDMSList.empty() ? nullptr : apoDMSList[0];
2900 :
2901 : // Get the IGE filename from if we have an ExternalRasterDMS.
2902 110 : if (poDMS)
2903 : {
2904 : const char *pszRawFilename =
2905 32 : poDMS->GetStringField("fileName.string");
2906 :
2907 32 : if (pszRawFilename != nullptr)
2908 : {
2909 : VSIStatBufL sStatBuf;
2910 : std::string osFullFilename =
2911 64 : CPLFormFilenameSafe(hHFA->pszPath, pszRawFilename, nullptr);
2912 :
2913 32 : if (VSIStatL(osFullFilename.c_str(), &sStatBuf) != 0)
2914 : {
2915 : const CPLString osExtension =
2916 0 : CPLGetExtensionSafe(pszRawFilename);
2917 : const CPLString osBasename =
2918 0 : CPLGetBasenameSafe(hHFA->pszFilename);
2919 0 : osFullFilename = CPLFormFilenameSafe(
2920 0 : hHFA->pszPath, osBasename, osExtension);
2921 :
2922 0 : if (VSIStatL(osFullFilename.c_str(), &sStatBuf) == 0)
2923 0 : hHFA->pszIGEFilename =
2924 0 : CPLStrdup(CPLFormFilenameSafe(nullptr, osBasename,
2925 : osExtension)
2926 : .c_str());
2927 : else
2928 0 : hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
2929 : }
2930 : else
2931 : {
2932 32 : hHFA->pszIGEFilename = CPLStrdup(pszRawFilename);
2933 : }
2934 : }
2935 : }
2936 : }
2937 :
2938 : // Return the full filename.
2939 113 : if (hHFA->pszIGEFilename)
2940 35 : return CPLFormFilenameSafe(hHFA->pszPath, hHFA->pszIGEFilename,
2941 35 : nullptr);
2942 :
2943 78 : return std::string();
2944 : }
2945 :
2946 : /************************************************************************/
2947 : /* HFACreateSpillStack() */
2948 : /* */
2949 : /* Create a new stack of raster layers in the spill (.ige) */
2950 : /* file. Create the spill file if it didn't exist before. */
2951 : /************************************************************************/
2952 :
2953 7 : bool HFACreateSpillStack(HFAInfo_t *psInfo, int nXSize, int nYSize, int nLayers,
2954 : int nBlockSize, EPTType eDataType,
2955 : GIntBig *pnValidFlagsOffset, GIntBig *pnDataOffset)
2956 :
2957 : {
2958 : // Form .ige filename.
2959 7 : if (nBlockSize <= 0)
2960 : {
2961 0 : CPLError(CE_Failure, CPLE_IllegalArg,
2962 : "HFACreateSpillStack: nBlockXSize < 0");
2963 0 : return false;
2964 : }
2965 :
2966 7 : if (psInfo->pszIGEFilename == nullptr)
2967 : {
2968 14 : const auto osExt = CPLGetExtensionSafe(psInfo->pszFilename);
2969 7 : if (EQUAL(osExt.c_str(), "rrd"))
2970 0 : psInfo->pszIGEFilename = CPLStrdup(
2971 0 : CPLResetExtensionSafe(psInfo->pszFilename, "rde").c_str());
2972 7 : else if (EQUAL(osExt.c_str(), "aux"))
2973 0 : psInfo->pszIGEFilename = CPLStrdup(
2974 0 : CPLResetExtensionSafe(psInfo->pszFilename, "axe").c_str());
2975 : else
2976 7 : psInfo->pszIGEFilename = CPLStrdup(
2977 14 : CPLResetExtensionSafe(psInfo->pszFilename, "ige").c_str());
2978 : }
2979 :
2980 7 : char *pszFullFilename = CPLStrdup(
2981 14 : CPLFormFilenameSafe(psInfo->pszPath, psInfo->pszIGEFilename, nullptr)
2982 : .c_str());
2983 :
2984 : // Try and open it. If we fail, create it and write the magic header.
2985 : static const char *const pszMagick = "ERDAS_IMG_EXTERNAL_RASTER";
2986 :
2987 7 : bool bRet = true;
2988 7 : VSILFILE *fpVSIL = VSIFOpenL(pszFullFilename, "r+b");
2989 7 : if (fpVSIL == nullptr)
2990 : {
2991 7 : fpVSIL = VSIFOpenL(pszFullFilename, "w+");
2992 7 : if (fpVSIL == nullptr)
2993 : {
2994 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2995 : "Failed to create spill file %s.\n%s",
2996 0 : psInfo->pszIGEFilename, VSIStrerror(errno));
2997 0 : return false;
2998 : }
2999 :
3000 7 : bRet &=
3001 7 : VSIFWriteL((void *)pszMagick, strlen(pszMagick) + 1, 1, fpVSIL) > 0;
3002 : }
3003 :
3004 7 : CPLFree(pszFullFilename);
3005 :
3006 : // Work out some details about the tiling scheme.
3007 7 : const int nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
3008 7 : const int nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
3009 : // const int nBlocks = nBlocksPerRow * nBlocksPerColumn;
3010 : const int nBytesPerBlock =
3011 7 : (nBlockSize * nBlockSize * HFAGetDataTypeBits(eDataType) + 7) / 8;
3012 :
3013 7 : const int nBytesPerRow = (nBlocksPerRow + 7) / 8;
3014 7 : const int nBlockMapSize = nBytesPerRow * nBlocksPerColumn;
3015 : // const int iFlagsSize = nBlockMapSize + 20;
3016 :
3017 : // Write stack prefix information.
3018 7 : bRet &= VSIFSeekL(fpVSIL, 0, SEEK_END) >= 0;
3019 :
3020 7 : GByte bUnknown = 1;
3021 7 : bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
3022 :
3023 7 : GInt32 nValue32 = nLayers;
3024 : HFAStandard(4, &nValue32);
3025 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3026 7 : nValue32 = nXSize;
3027 : HFAStandard(4, &nValue32);
3028 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3029 7 : nValue32 = nYSize;
3030 : HFAStandard(4, &nValue32);
3031 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3032 7 : nValue32 = nBlockSize;
3033 : HFAStandard(4, &nValue32);
3034 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3035 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3036 7 : bUnknown = 3;
3037 7 : bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
3038 7 : bUnknown = 0;
3039 7 : bRet &= VSIFWriteL(&bUnknown, 1, 1, fpVSIL) > 0;
3040 :
3041 : // Write out ValidFlags section(s).
3042 7 : *pnValidFlagsOffset = VSIFTellL(fpVSIL);
3043 :
3044 : unsigned char *pabyBlockMap =
3045 7 : static_cast<unsigned char *>(VSI_MALLOC_VERBOSE(nBlockMapSize));
3046 7 : if (pabyBlockMap == nullptr)
3047 : {
3048 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpVSIL));
3049 0 : return false;
3050 : }
3051 :
3052 7 : memset(pabyBlockMap, 0xff, nBlockMapSize);
3053 14 : for (int iBand = 0; iBand < nLayers; iBand++)
3054 : {
3055 7 : nValue32 = 1; // Unknown
3056 : HFAStandard(4, &nValue32);
3057 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3058 7 : nValue32 = 0; // Unknown
3059 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3060 7 : nValue32 = nBlocksPerColumn;
3061 : HFAStandard(4, &nValue32);
3062 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3063 7 : nValue32 = nBlocksPerRow;
3064 : HFAStandard(4, &nValue32);
3065 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3066 7 : nValue32 = 0x30000; // Unknown
3067 : HFAStandard(4, &nValue32);
3068 7 : bRet &= VSIFWriteL(&nValue32, 4, 1, fpVSIL) > 0;
3069 :
3070 7 : const int iRemainder = nBlocksPerRow % 8;
3071 7 : CPLDebug("HFACreate",
3072 : "Block map size %d, bytes per row %d, remainder %d.",
3073 : nBlockMapSize, nBytesPerRow, iRemainder);
3074 7 : if (iRemainder)
3075 : {
3076 14 : for (int i = nBytesPerRow - 1; i < nBlockMapSize; i += nBytesPerRow)
3077 7 : pabyBlockMap[i] = static_cast<GByte>((1 << iRemainder) - 1);
3078 : }
3079 :
3080 7 : bRet &= VSIFWriteL(pabyBlockMap, nBlockMapSize, 1, fpVSIL) > 0;
3081 : }
3082 7 : CPLFree(pabyBlockMap);
3083 7 : pabyBlockMap = nullptr;
3084 :
3085 : // Extend the file to account for all the imagery space.
3086 7 : const GIntBig nTileDataSize = static_cast<GIntBig>(nBytesPerBlock) *
3087 7 : nBlocksPerRow * nBlocksPerColumn * nLayers;
3088 :
3089 7 : *pnDataOffset = VSIFTellL(fpVSIL);
3090 :
3091 7 : if (!bRet || VSIFTruncateL(fpVSIL, nTileDataSize + *pnDataOffset) != 0)
3092 : {
3093 0 : CPLError(CE_Failure, CPLE_FileIO,
3094 : "Failed to extend %s to full size (" CPL_FRMT_GIB " bytes), "
3095 : "likely out of disk space.\n%s",
3096 0 : psInfo->pszIGEFilename, nTileDataSize + *pnDataOffset,
3097 0 : VSIStrerror(errno));
3098 :
3099 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpVSIL));
3100 0 : return false;
3101 : }
3102 :
3103 7 : if (VSIFCloseL(fpVSIL) != 0)
3104 0 : return false;
3105 :
3106 7 : return true;
3107 : }
3108 :
3109 : /************************************************************************/
3110 : /* HFAReadAndValidatePoly() */
3111 : /************************************************************************/
3112 :
3113 3 : static bool HFAReadAndValidatePoly(HFAEntry *poTarget, const char *pszName,
3114 : Efga_Polynomial *psRetPoly)
3115 :
3116 : {
3117 3 : memset(psRetPoly, 0, sizeof(Efga_Polynomial));
3118 :
3119 6 : CPLString osFldName;
3120 3 : osFldName.Printf("%sorder", pszName);
3121 3 : psRetPoly->order = poTarget->GetIntField(osFldName);
3122 :
3123 3 : if (psRetPoly->order < 1 || psRetPoly->order > 3)
3124 0 : return false;
3125 :
3126 : // Validate that things are in a "well known" form.
3127 3 : osFldName.Printf("%snumdimtransform", pszName);
3128 3 : const int numdimtransform = poTarget->GetIntField(osFldName);
3129 :
3130 3 : osFldName.Printf("%snumdimpolynomial", pszName);
3131 3 : const int numdimpolynomial = poTarget->GetIntField(osFldName);
3132 :
3133 3 : osFldName.Printf("%stermcount", pszName);
3134 3 : const int termcount = poTarget->GetIntField(osFldName);
3135 :
3136 3 : if (numdimtransform != 2 || numdimpolynomial != 2)
3137 0 : return false;
3138 :
3139 3 : if ((psRetPoly->order == 1 && termcount != 3) ||
3140 3 : (psRetPoly->order == 2 && termcount != 6) ||
3141 3 : (psRetPoly->order == 3 && termcount != 10))
3142 0 : return false;
3143 :
3144 : // We don't check the exponent organization for now. Hopefully
3145 : // it is always standard.
3146 :
3147 : // Get coefficients.
3148 43 : for (int i = 0; i < termcount * 2 - 2; i++)
3149 : {
3150 40 : osFldName.Printf("%spolycoefmtx[%d]", pszName, i);
3151 40 : psRetPoly->polycoefmtx[i] = poTarget->GetDoubleField(osFldName);
3152 : }
3153 :
3154 9 : for (int i = 0; i < 2; i++)
3155 : {
3156 6 : osFldName.Printf("%spolycoefvector[%d]", pszName, i);
3157 6 : psRetPoly->polycoefvector[i] = poTarget->GetDoubleField(osFldName);
3158 : }
3159 :
3160 3 : return true;
3161 : }
3162 :
3163 : /************************************************************************/
3164 : /* HFAReadXFormStack() */
3165 : /************************************************************************/
3166 :
3167 286 : int HFAReadXFormStack(HFAHandle hHFA, Efga_Polynomial **ppasPolyListForward,
3168 : Efga_Polynomial **ppasPolyListReverse)
3169 :
3170 : {
3171 286 : if (hHFA->nBands == 0)
3172 0 : return 0;
3173 :
3174 : // Get the HFA node.
3175 : HFAEntry *poXFormHeader =
3176 286 : hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm");
3177 286 : if (poXFormHeader == nullptr)
3178 284 : return 0;
3179 :
3180 : // Loop over children, collecting XForms.
3181 2 : int nStepCount = 0;
3182 2 : *ppasPolyListForward = nullptr;
3183 2 : *ppasPolyListReverse = nullptr;
3184 :
3185 5 : for (HFAEntry *poXForm = poXFormHeader->GetChild(); poXForm != nullptr;
3186 3 : poXForm = poXForm->GetNext())
3187 : {
3188 3 : bool bSuccess = false;
3189 : Efga_Polynomial sForward;
3190 : Efga_Polynomial sReverse;
3191 3 : memset(&sForward, 0, sizeof(sForward));
3192 3 : memset(&sReverse, 0, sizeof(sReverse));
3193 :
3194 3 : if (EQUAL(poXForm->GetType(), "Efga_Polynomial"))
3195 : {
3196 1 : bSuccess = HFAReadAndValidatePoly(poXForm, "", &sForward);
3197 :
3198 1 : if (bSuccess)
3199 : {
3200 : double adfGT[6] = {
3201 1 : sForward.polycoefvector[0], sForward.polycoefmtx[0],
3202 1 : sForward.polycoefmtx[2], sForward.polycoefvector[1],
3203 1 : sForward.polycoefmtx[1], sForward.polycoefmtx[3]};
3204 :
3205 1 : double adfInvGT[6] = {};
3206 1 : bSuccess = HFAInvGeoTransform(adfGT, adfInvGT);
3207 1 : if (!bSuccess)
3208 0 : memset(adfInvGT, 0, sizeof(adfInvGT));
3209 :
3210 1 : sReverse.order = sForward.order;
3211 1 : sReverse.polycoefvector[0] = adfInvGT[0];
3212 1 : sReverse.polycoefmtx[0] = adfInvGT[1];
3213 1 : sReverse.polycoefmtx[2] = adfInvGT[2];
3214 1 : sReverse.polycoefvector[1] = adfInvGT[3];
3215 1 : sReverse.polycoefmtx[1] = adfInvGT[4];
3216 1 : sReverse.polycoefmtx[3] = adfInvGT[5];
3217 : }
3218 : }
3219 2 : else if (EQUAL(poXForm->GetType(), "GM_PolyPair"))
3220 : {
3221 2 : bSuccess = HFAReadAndValidatePoly(poXForm, "forward.", &sForward) &&
3222 1 : HFAReadAndValidatePoly(poXForm, "reverse.", &sReverse);
3223 : }
3224 :
3225 3 : if (bSuccess)
3226 : {
3227 2 : nStepCount++;
3228 4 : *ppasPolyListForward = static_cast<Efga_Polynomial *>(CPLRealloc(
3229 2 : *ppasPolyListForward, sizeof(Efga_Polynomial) * nStepCount));
3230 2 : memcpy(*ppasPolyListForward + nStepCount - 1, &sForward,
3231 : sizeof(sForward));
3232 :
3233 4 : *ppasPolyListReverse = static_cast<Efga_Polynomial *>(CPLRealloc(
3234 2 : *ppasPolyListReverse, sizeof(Efga_Polynomial) * nStepCount));
3235 2 : memcpy(*ppasPolyListReverse + nStepCount - 1, &sReverse,
3236 : sizeof(sReverse));
3237 : }
3238 : }
3239 :
3240 2 : return nStepCount;
3241 : }
3242 :
3243 : /************************************************************************/
3244 : /* HFAEvaluateXFormStack() */
3245 : /************************************************************************/
3246 :
3247 36 : int HFAEvaluateXFormStack(int nStepCount, int bForward,
3248 : Efga_Polynomial *pasPolyList, double *pdfX,
3249 : double *pdfY)
3250 :
3251 : {
3252 108 : for (int iStep = 0; iStep < nStepCount; iStep++)
3253 : {
3254 72 : const Efga_Polynomial *psStep =
3255 0 : bForward ? pasPolyList + iStep
3256 72 : : pasPolyList + nStepCount - iStep - 1;
3257 :
3258 72 : if (psStep->order == 1)
3259 : {
3260 36 : const double dfXOut = psStep->polycoefvector[0] +
3261 36 : psStep->polycoefmtx[0] * *pdfX +
3262 36 : psStep->polycoefmtx[2] * *pdfY;
3263 :
3264 36 : const double dfYOut = psStep->polycoefvector[1] +
3265 36 : psStep->polycoefmtx[1] * *pdfX +
3266 36 : psStep->polycoefmtx[3] * *pdfY;
3267 :
3268 36 : *pdfX = dfXOut;
3269 36 : *pdfY = dfYOut;
3270 : }
3271 36 : else if (psStep->order == 2)
3272 : {
3273 0 : const double dfXOut = psStep->polycoefvector[0] +
3274 0 : psStep->polycoefmtx[0] * *pdfX +
3275 0 : psStep->polycoefmtx[2] * *pdfY +
3276 0 : psStep->polycoefmtx[4] * *pdfX * *pdfX +
3277 0 : psStep->polycoefmtx[6] * *pdfX * *pdfY +
3278 0 : psStep->polycoefmtx[8] * *pdfY * *pdfY;
3279 0 : const double dfYOut = psStep->polycoefvector[1] +
3280 0 : psStep->polycoefmtx[1] * *pdfX +
3281 0 : psStep->polycoefmtx[3] * *pdfY +
3282 0 : psStep->polycoefmtx[5] * *pdfX * *pdfX +
3283 0 : psStep->polycoefmtx[7] * *pdfX * *pdfY +
3284 0 : psStep->polycoefmtx[9] * *pdfY * *pdfY;
3285 :
3286 0 : *pdfX = dfXOut;
3287 0 : *pdfY = dfYOut;
3288 : }
3289 36 : else if (psStep->order == 3)
3290 : {
3291 36 : const double dfXOut =
3292 36 : psStep->polycoefvector[0] + psStep->polycoefmtx[0] * *pdfX +
3293 36 : psStep->polycoefmtx[2] * *pdfY +
3294 36 : psStep->polycoefmtx[4] * *pdfX * *pdfX +
3295 36 : psStep->polycoefmtx[6] * *pdfX * *pdfY +
3296 36 : psStep->polycoefmtx[8] * *pdfY * *pdfY +
3297 36 : psStep->polycoefmtx[10] * *pdfX * *pdfX * *pdfX +
3298 36 : psStep->polycoefmtx[12] * *pdfX * *pdfX * *pdfY +
3299 36 : psStep->polycoefmtx[14] * *pdfX * *pdfY * *pdfY +
3300 36 : psStep->polycoefmtx[16] * *pdfY * *pdfY * *pdfY;
3301 36 : const double dfYOut =
3302 36 : psStep->polycoefvector[1] + psStep->polycoefmtx[1] * *pdfX +
3303 36 : psStep->polycoefmtx[3] * *pdfY +
3304 36 : psStep->polycoefmtx[5] * *pdfX * *pdfX +
3305 36 : psStep->polycoefmtx[7] * *pdfX * *pdfY +
3306 36 : psStep->polycoefmtx[9] * *pdfY * *pdfY +
3307 36 : psStep->polycoefmtx[11] * *pdfX * *pdfX * *pdfX +
3308 36 : psStep->polycoefmtx[13] * *pdfX * *pdfX * *pdfY +
3309 36 : psStep->polycoefmtx[15] * *pdfX * *pdfY * *pdfY +
3310 36 : psStep->polycoefmtx[17] * *pdfY * *pdfY * *pdfY;
3311 :
3312 36 : *pdfX = dfXOut;
3313 36 : *pdfY = dfYOut;
3314 : }
3315 : else
3316 0 : return FALSE;
3317 : }
3318 :
3319 36 : return TRUE;
3320 : }
3321 :
3322 : /************************************************************************/
3323 : /* HFAWriteXFormStack() */
3324 : /************************************************************************/
3325 :
3326 2 : CPLErr HFAWriteXFormStack(HFAHandle hHFA, int nBand, int nXFormCount,
3327 : Efga_Polynomial **ppasPolyListForward,
3328 : Efga_Polynomial **ppasPolyListReverse)
3329 :
3330 : {
3331 2 : if (nXFormCount == 0)
3332 0 : return CE_None;
3333 :
3334 2 : if (ppasPolyListForward[0]->order != 1)
3335 : {
3336 0 : CPLError(
3337 : CE_Failure, CPLE_AppDefined,
3338 : "For now HFAWriteXFormStack() only supports order 1 polynomials");
3339 0 : return CE_Failure;
3340 : }
3341 :
3342 2 : if (nBand < 0 || nBand > hHFA->nBands)
3343 0 : return CE_Failure;
3344 :
3345 : // If no band number is provided, operate on all bands.
3346 2 : if (nBand == 0)
3347 : {
3348 2 : for (nBand = 1; nBand <= hHFA->nBands; nBand++)
3349 : {
3350 : CPLErr eErr =
3351 1 : HFAWriteXFormStack(hHFA, nBand, nXFormCount,
3352 : ppasPolyListForward, ppasPolyListReverse);
3353 1 : if (eErr != CE_None)
3354 0 : return eErr;
3355 : }
3356 :
3357 1 : return CE_None;
3358 : }
3359 :
3360 : // Fetch our band node.
3361 1 : HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
3362 1 : HFAEntry *poXFormHeader = poBandNode->GetNamedChild("MapToPixelXForm");
3363 1 : if (poXFormHeader == nullptr)
3364 : {
3365 1 : poXFormHeader = HFAEntry::New(hHFA, "MapToPixelXForm",
3366 : "Exfr_GenericXFormHeader", poBandNode);
3367 1 : poXFormHeader->MakeData(23);
3368 1 : poXFormHeader->SetPosition();
3369 1 : poXFormHeader->SetStringField("titleList.string", "Affine");
3370 : }
3371 :
3372 : // Loop over XForms.
3373 2 : for (int iXForm = 0; iXForm < nXFormCount; iXForm++)
3374 : {
3375 1 : Efga_Polynomial *psForward = *ppasPolyListForward + iXForm;
3376 2 : CPLString osXFormName;
3377 1 : osXFormName.Printf("XForm%d", iXForm);
3378 :
3379 1 : HFAEntry *poXForm = poXFormHeader->GetNamedChild(osXFormName);
3380 :
3381 1 : if (poXForm == nullptr)
3382 : {
3383 1 : poXForm = HFAEntry::New(hHFA, osXFormName, "Efga_Polynomial",
3384 : poXFormHeader);
3385 1 : poXForm->MakeData(136);
3386 1 : poXForm->SetPosition();
3387 : }
3388 :
3389 1 : poXForm->SetIntField("order", 1);
3390 1 : poXForm->SetIntField("numdimtransform", 2);
3391 1 : poXForm->SetIntField("numdimpolynomial", 2);
3392 1 : poXForm->SetIntField("termcount", 3);
3393 1 : poXForm->SetIntField("exponentlist[0]", 0);
3394 1 : poXForm->SetIntField("exponentlist[1]", 0);
3395 1 : poXForm->SetIntField("exponentlist[2]", 1);
3396 1 : poXForm->SetIntField("exponentlist[3]", 0);
3397 1 : poXForm->SetIntField("exponentlist[4]", 0);
3398 1 : poXForm->SetIntField("exponentlist[5]", 1);
3399 :
3400 1 : poXForm->SetIntField("polycoefmtx[-3]", EPT_f64);
3401 1 : poXForm->SetIntField("polycoefmtx[-2]", 2);
3402 1 : poXForm->SetIntField("polycoefmtx[-1]", 2);
3403 1 : poXForm->SetDoubleField("polycoefmtx[0]", psForward->polycoefmtx[0]);
3404 1 : poXForm->SetDoubleField("polycoefmtx[1]", psForward->polycoefmtx[1]);
3405 1 : poXForm->SetDoubleField("polycoefmtx[2]", psForward->polycoefmtx[2]);
3406 1 : poXForm->SetDoubleField("polycoefmtx[3]", psForward->polycoefmtx[3]);
3407 :
3408 1 : poXForm->SetIntField("polycoefvector[-3]", EPT_f64);
3409 1 : poXForm->SetIntField("polycoefvector[-2]", 1);
3410 1 : poXForm->SetIntField("polycoefvector[-1]", 2);
3411 1 : poXForm->SetDoubleField("polycoefvector[0]",
3412 : psForward->polycoefvector[0]);
3413 1 : poXForm->SetDoubleField("polycoefvector[1]",
3414 : psForward->polycoefvector[1]);
3415 : }
3416 :
3417 1 : return CE_None;
3418 : }
3419 :
3420 : /************************************************************************/
3421 : /* HFAReadCameraModel() */
3422 : /************************************************************************/
3423 :
3424 535 : char **HFAReadCameraModel(HFAHandle hHFA)
3425 :
3426 : {
3427 535 : if (hHFA->nBands == 0)
3428 0 : return nullptr;
3429 :
3430 : // Get the camera model node, and confirm its type.
3431 : HFAEntry *poXForm =
3432 535 : hHFA->papoBand[0]->poNode->GetNamedChild("MapToPixelXForm.XForm0");
3433 535 : if (poXForm == nullptr)
3434 529 : return nullptr;
3435 :
3436 6 : if (!EQUAL(poXForm->GetType(), "Camera_ModelX"))
3437 5 : return nullptr;
3438 :
3439 : // Convert the values to metadata.
3440 1 : char **papszMD = nullptr;
3441 : static const char *const apszFields[] = {"direction",
3442 : "refType",
3443 : "demsource",
3444 : "PhotoDirection",
3445 : "RotationSystem",
3446 : "demfilename",
3447 : "demzunits",
3448 : "forSrcAffine[0]",
3449 : "forSrcAffine[1]",
3450 : "forSrcAffine[2]",
3451 : "forSrcAffine[3]",
3452 : "forSrcAffine[4]",
3453 : "forSrcAffine[5]",
3454 : "forDstAffine[0]",
3455 : "forDstAffine[1]",
3456 : "forDstAffine[2]",
3457 : "forDstAffine[3]",
3458 : "forDstAffine[4]",
3459 : "forDstAffine[5]",
3460 : "invSrcAffine[0]",
3461 : "invSrcAffine[1]",
3462 : "invSrcAffine[2]",
3463 : "invSrcAffine[3]",
3464 : "invSrcAffine[4]",
3465 : "invSrcAffine[5]",
3466 : "invDstAffine[0]",
3467 : "invDstAffine[1]",
3468 : "invDstAffine[2]",
3469 : "invDstAffine[3]",
3470 : "invDstAffine[4]",
3471 : "invDstAffine[5]",
3472 : "z_mean",
3473 : "lat0",
3474 : "lon0",
3475 : "coeffs[0]",
3476 : "coeffs[1]",
3477 : "coeffs[2]",
3478 : "coeffs[3]",
3479 : "coeffs[4]",
3480 : "coeffs[5]",
3481 : "coeffs[6]",
3482 : "coeffs[7]",
3483 : "coeffs[8]",
3484 : "LensDistortion[0]",
3485 : "LensDistortion[1]",
3486 : "LensDistortion[2]",
3487 : nullptr};
3488 :
3489 1 : const char *pszValue = nullptr;
3490 47 : for (int i = 0; apszFields[i] != nullptr; i++)
3491 : {
3492 46 : pszValue = poXForm->GetStringField(apszFields[i]);
3493 46 : if (pszValue == nullptr)
3494 1 : pszValue = "";
3495 :
3496 46 : papszMD = CSLSetNameValue(papszMD, apszFields[i], pszValue);
3497 : }
3498 :
3499 : // Create a pseudo-entry for the MIFObject with the outputProjection.
3500 : HFAEntry *poProjInfo =
3501 1 : HFAEntry::BuildEntryFromMIFObject(poXForm, "outputProjection");
3502 1 : if (poProjInfo)
3503 : {
3504 : // Fetch the datum.
3505 : Eprj_Datum sDatum;
3506 :
3507 1 : memset(&sDatum, 0, sizeof(sDatum));
3508 :
3509 1 : sDatum.datumname =
3510 1 : (char *)poProjInfo->GetStringField("earthModel.datum.datumname");
3511 :
3512 1 : const int nDatumType = poProjInfo->GetIntField("earthModel.datum.type");
3513 1 : if (nDatumType < 0 || nDatumType > EPRJ_DATUM_NONE)
3514 : {
3515 0 : CPLDebug("HFA", "Invalid value for datum type: %d", nDatumType);
3516 0 : sDatum.type = EPRJ_DATUM_NONE;
3517 : }
3518 : else
3519 : {
3520 1 : sDatum.type = static_cast<Eprj_DatumType>(nDatumType);
3521 : }
3522 :
3523 8 : for (int i = 0; i < 7; i++)
3524 : {
3525 7 : char szFieldName[60] = {};
3526 :
3527 7 : snprintf(szFieldName, sizeof(szFieldName),
3528 : "earthModel.datum.params[%d]", i);
3529 7 : sDatum.params[i] = poProjInfo->GetDoubleField(szFieldName);
3530 : }
3531 :
3532 1 : sDatum.gridname =
3533 1 : (char *)poProjInfo->GetStringField("earthModel.datum.gridname");
3534 :
3535 : // Fetch the projection parameters.
3536 : Eprj_ProParameters sPro;
3537 :
3538 1 : memset(&sPro, 0, sizeof(sPro));
3539 :
3540 1 : sPro.proType =
3541 1 : (Eprj_ProType)poProjInfo->GetIntField("projectionObject.proType");
3542 1 : sPro.proNumber = poProjInfo->GetIntField("projectionObject.proNumber");
3543 1 : sPro.proExeName =
3544 1 : (char *)poProjInfo->GetStringField("projectionObject.proExeName");
3545 1 : sPro.proName =
3546 1 : (char *)poProjInfo->GetStringField("projectionObject.proName");
3547 1 : sPro.proZone = poProjInfo->GetIntField("projectionObject.proZone");
3548 :
3549 16 : for (int i = 0; i < 15; i++)
3550 : {
3551 15 : char szFieldName[40] = {};
3552 :
3553 15 : snprintf(szFieldName, sizeof(szFieldName),
3554 : "projectionObject.proParams[%d]", i);
3555 15 : sPro.proParams[i] = poProjInfo->GetDoubleField(szFieldName);
3556 : }
3557 :
3558 : // Fetch the spheroid.
3559 1 : sPro.proSpheroid.sphereName = (char *)poProjInfo->GetStringField(
3560 : "earthModel.proSpheroid.sphereName");
3561 1 : sPro.proSpheroid.a =
3562 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.a");
3563 1 : sPro.proSpheroid.b =
3564 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.b");
3565 1 : sPro.proSpheroid.eSquared =
3566 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.eSquared");
3567 1 : sPro.proSpheroid.radius =
3568 1 : poProjInfo->GetDoubleField("earthModel.proSpheroid.radius");
3569 :
3570 : // Fetch the projection info.
3571 : // poProjInfo->DumpFieldValues( stdout, "" );
3572 :
3573 2 : auto poSRS = HFAPCSStructToOSR(&sDatum, &sPro, nullptr, nullptr);
3574 :
3575 1 : if (poSRS)
3576 : {
3577 1 : char *pszProjection = nullptr;
3578 1 : if (poSRS->exportToWkt(&pszProjection) == OGRERR_NONE)
3579 : {
3580 : papszMD =
3581 1 : CSLSetNameValue(papszMD, "outputProjection", pszProjection);
3582 : }
3583 1 : CPLFree(pszProjection);
3584 : }
3585 :
3586 1 : delete poProjInfo;
3587 : }
3588 :
3589 : // Fetch the horizontal units.
3590 1 : pszValue = poXForm->GetStringField("outputHorizontalUnits.string");
3591 1 : if (pszValue == nullptr)
3592 0 : pszValue = "";
3593 :
3594 1 : papszMD = CSLSetNameValue(papszMD, "outputHorizontalUnits", pszValue);
3595 :
3596 : // Fetch the elevationinfo.
3597 : HFAEntry *poElevInfo =
3598 1 : HFAEntry::BuildEntryFromMIFObject(poXForm, "outputElevationInfo");
3599 1 : if (poElevInfo)
3600 : {
3601 : // poElevInfo->DumpFieldValues( stdout, "" );
3602 :
3603 1 : if (poElevInfo->GetDataSize() != 0)
3604 : {
3605 : static const char *const apszEFields[] = {
3606 : "verticalDatum.datumname", "verticalDatum.type",
3607 : "elevationUnit", "elevationType", nullptr};
3608 :
3609 5 : for (int i = 0; apszEFields[i] != nullptr; i++)
3610 : {
3611 4 : pszValue = poElevInfo->GetStringField(apszEFields[i]);
3612 4 : if (pszValue == nullptr)
3613 0 : pszValue = "";
3614 :
3615 4 : papszMD = CSLSetNameValue(papszMD, apszEFields[i], pszValue);
3616 : }
3617 : }
3618 :
3619 1 : delete poElevInfo;
3620 : }
3621 :
3622 1 : return papszMD;
3623 : }
3624 :
3625 : /************************************************************************/
3626 : /* HFAReadElevationUnit() */
3627 : /************************************************************************/
3628 :
3629 615 : const char *HFAReadElevationUnit(HFAHandle hHFA, int iBand)
3630 : {
3631 615 : if (hHFA->nBands <= iBand)
3632 0 : return nullptr;
3633 :
3634 615 : HFABand *poBand(hHFA->papoBand[iBand]);
3635 615 : if (poBand == nullptr || poBand->poNode == nullptr)
3636 : {
3637 0 : return nullptr;
3638 : }
3639 615 : HFAEntry *poElevInfo = poBand->poNode->GetNamedChild("Elevation_Info");
3640 615 : if (poElevInfo == nullptr)
3641 : {
3642 612 : return nullptr;
3643 : }
3644 3 : return poElevInfo->GetStringField("elevationUnit");
3645 : }
3646 :
3647 : /************************************************************************/
3648 : /* HFASetGeoTransform() */
3649 : /* */
3650 : /* Set a MapInformation and XForm block. Allows for rotated */
3651 : /* and shared geotransforms. */
3652 : /************************************************************************/
3653 :
3654 1 : CPLErr HFASetGeoTransform(HFAHandle hHFA, const char *pszProName,
3655 : const char *pszUnits, double *padfGeoTransform)
3656 :
3657 : {
3658 : // Write MapInformation.
3659 2 : for (int nBand = 1; nBand <= hHFA->nBands; nBand++)
3660 : {
3661 1 : HFAEntry *poBandNode = hHFA->papoBand[nBand - 1]->poNode;
3662 :
3663 1 : HFAEntry *poMI = poBandNode->GetNamedChild("MapInformation");
3664 1 : if (poMI == nullptr)
3665 : {
3666 1 : poMI = HFAEntry::New(hHFA, "MapInformation", "Eimg_MapInformation",
3667 : poBandNode);
3668 1 : poMI->MakeData(
3669 1 : static_cast<int>(18 + strlen(pszProName) + strlen(pszUnits)));
3670 1 : poMI->SetPosition();
3671 : }
3672 :
3673 1 : poMI->SetStringField("projection.string", pszProName);
3674 1 : poMI->SetStringField("units.string", pszUnits);
3675 : }
3676 :
3677 : // Write XForm.
3678 1 : double adfAdjTransform[6] = {};
3679 :
3680 : // Offset by half pixel.
3681 :
3682 1 : memcpy(adfAdjTransform, padfGeoTransform, sizeof(double) * 6);
3683 1 : adfAdjTransform[0] += adfAdjTransform[1] * 0.5;
3684 1 : adfAdjTransform[0] += adfAdjTransform[2] * 0.5;
3685 1 : adfAdjTransform[3] += adfAdjTransform[4] * 0.5;
3686 1 : adfAdjTransform[3] += adfAdjTransform[5] * 0.5;
3687 :
3688 : // Invert.
3689 1 : double adfRevTransform[6] = {};
3690 1 : if (!HFAInvGeoTransform(adfAdjTransform, adfRevTransform))
3691 0 : memset(adfRevTransform, 0, sizeof(adfRevTransform));
3692 :
3693 : // Assign to polynomial object.
3694 :
3695 : Efga_Polynomial sForward;
3696 1 : memset(&sForward, 0, sizeof(sForward));
3697 1 : Efga_Polynomial *psForward = &sForward;
3698 1 : sForward.order = 1;
3699 1 : sForward.polycoefvector[0] = adfRevTransform[0];
3700 1 : sForward.polycoefmtx[0] = adfRevTransform[1];
3701 1 : sForward.polycoefmtx[1] = adfRevTransform[4];
3702 1 : sForward.polycoefvector[1] = adfRevTransform[3];
3703 1 : sForward.polycoefmtx[2] = adfRevTransform[2];
3704 1 : sForward.polycoefmtx[3] = adfRevTransform[5];
3705 :
3706 1 : Efga_Polynomial sReverse = sForward;
3707 1 : Efga_Polynomial *psReverse = &sReverse;
3708 :
3709 2 : return HFAWriteXFormStack(hHFA, 0, 1, &psForward, &psReverse);
3710 : }
3711 :
3712 : /************************************************************************/
3713 : /* HFARenameReferences() */
3714 : /* */
3715 : /* Rename references in this .img file from the old basename to */
3716 : /* a new basename. This should be passed on to .aux and .rrd */
3717 : /* files and should include references to .aux, .rrd and .ige. */
3718 : /************************************************************************/
3719 :
3720 4 : CPLErr HFARenameReferences(HFAHandle hHFA, const char *pszNewBase,
3721 : const char *pszOldBase)
3722 :
3723 : {
3724 : // Handle RRDNamesList updates.
3725 : std::vector<HFAEntry *> apoNodeList =
3726 4 : hHFA->poRoot->FindChildren("RRDNamesList", nullptr);
3727 :
3728 6 : for (size_t iNode = 0; iNode < apoNodeList.size(); iNode++)
3729 : {
3730 2 : HFAEntry *poRRDNL = apoNodeList[iNode];
3731 4 : std::vector<CPLString> aosNL;
3732 :
3733 : // Collect all the existing names.
3734 2 : const int nNameCount = poRRDNL->GetFieldCount("nameList");
3735 :
3736 4 : CPLString osAlgorithm = poRRDNL->GetStringField("algorithm.string");
3737 4 : for (int i = 0; i < nNameCount; i++)
3738 : {
3739 2 : CPLString osFN;
3740 2 : osFN.Printf("nameList[%d].string", i);
3741 2 : aosNL.push_back(poRRDNL->GetStringField(osFN));
3742 : }
3743 :
3744 : // Adjust the names to the new form.
3745 4 : for (int i = 0; i < nNameCount; i++)
3746 : {
3747 2 : if (strncmp(aosNL[i], pszOldBase, strlen(pszOldBase)) == 0)
3748 : {
3749 2 : std::string osNew = pszNewBase;
3750 2 : osNew += aosNL[i].c_str() + strlen(pszOldBase);
3751 2 : aosNL[i] = std::move(osNew);
3752 : }
3753 : }
3754 :
3755 : // Try to make sure the RRDNamesList is big enough to hold the
3756 : // adjusted name list.
3757 2 : if (strlen(pszNewBase) > strlen(pszOldBase))
3758 : {
3759 1 : CPLDebug("HFA", "Growing RRDNamesList to hold new names");
3760 1 : poRRDNL->MakeData(static_cast<int>(
3761 1 : poRRDNL->GetDataSize() +
3762 1 : nNameCount * (strlen(pszNewBase) - strlen(pszOldBase))));
3763 : }
3764 :
3765 : // Initialize the whole thing to zeros for a clean start.
3766 2 : memset(poRRDNL->GetData(), 0, poRRDNL->GetDataSize());
3767 :
3768 : // Write the updates back to the file.
3769 2 : poRRDNL->SetStringField("algorithm.string", osAlgorithm);
3770 4 : for (int i = 0; i < nNameCount; i++)
3771 : {
3772 4 : CPLString osFN;
3773 2 : osFN.Printf("nameList[%d].string", i);
3774 2 : poRRDNL->SetStringField(osFN, aosNL[i]);
3775 : }
3776 : }
3777 :
3778 : // Spill file references.
3779 : apoNodeList =
3780 4 : hHFA->poRoot->FindChildren("ExternalRasterDMS", "ImgExternalRaster");
3781 :
3782 8 : for (size_t iNode = 0; iNode < apoNodeList.size(); iNode++)
3783 : {
3784 4 : HFAEntry *poERDMS = apoNodeList[iNode];
3785 :
3786 4 : if (poERDMS == nullptr)
3787 0 : continue;
3788 :
3789 : // Fetch all existing values.
3790 8 : CPLString osFileName = poERDMS->GetStringField("fileName.string");
3791 :
3792 : GInt32 anValidFlagsOffset[2] = {
3793 4 : poERDMS->GetIntField("layerStackValidFlagsOffset[0]"),
3794 4 : poERDMS->GetIntField("layerStackValidFlagsOffset[1]")};
3795 :
3796 : GInt32 anStackDataOffset[2] = {
3797 4 : poERDMS->GetIntField("layerStackDataOffset[0]"),
3798 4 : poERDMS->GetIntField("layerStackDataOffset[1]")};
3799 :
3800 4 : const GInt32 nStackCount = poERDMS->GetIntField("layerStackCount");
3801 4 : const GInt32 nStackIndex = poERDMS->GetIntField("layerStackIndex");
3802 :
3803 : // Update the filename.
3804 4 : if (strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0)
3805 : {
3806 4 : std::string osNew = pszNewBase;
3807 4 : osNew += osFileName.c_str() + strlen(pszOldBase);
3808 4 : osFileName = std::move(osNew);
3809 : }
3810 :
3811 : // Grow the node if needed.
3812 4 : if (strlen(pszNewBase) > strlen(pszOldBase))
3813 : {
3814 2 : CPLDebug("HFA", "Growing ExternalRasterDMS to hold new names");
3815 2 : poERDMS->MakeData(
3816 2 : static_cast<int>(poERDMS->GetDataSize() +
3817 2 : (strlen(pszNewBase) - strlen(pszOldBase))));
3818 : }
3819 :
3820 : // Initialize the whole thing to zeros for a clean start.
3821 4 : memset(poERDMS->GetData(), 0, poERDMS->GetDataSize());
3822 :
3823 : // Write it all out again, this may change the size of the node.
3824 4 : poERDMS->SetStringField("fileName.string", osFileName);
3825 4 : poERDMS->SetIntField("layerStackValidFlagsOffset[0]",
3826 : anValidFlagsOffset[0]);
3827 4 : poERDMS->SetIntField("layerStackValidFlagsOffset[1]",
3828 : anValidFlagsOffset[1]);
3829 :
3830 4 : poERDMS->SetIntField("layerStackDataOffset[0]", anStackDataOffset[0]);
3831 4 : poERDMS->SetIntField("layerStackDataOffset[1]", anStackDataOffset[1]);
3832 :
3833 4 : poERDMS->SetIntField("layerStackCount", nStackCount);
3834 4 : poERDMS->SetIntField("layerStackIndex", nStackIndex);
3835 : }
3836 :
3837 : // DependentFile.
3838 : apoNodeList =
3839 4 : hHFA->poRoot->FindChildren("DependentFile", "Eimg_DependentFile");
3840 :
3841 6 : for (size_t iNode = 0; iNode < apoNodeList.size(); iNode++)
3842 : {
3843 : CPLString osFileName =
3844 4 : apoNodeList[iNode]->GetStringField("dependent.string");
3845 :
3846 : // Grow the node if needed.
3847 2 : if (strlen(pszNewBase) > strlen(pszOldBase))
3848 : {
3849 1 : CPLDebug("HFA", "Growing DependentFile to hold new names");
3850 1 : apoNodeList[iNode]->MakeData(
3851 1 : static_cast<int>(apoNodeList[iNode]->GetDataSize() +
3852 1 : (strlen(pszNewBase) - strlen(pszOldBase))));
3853 : }
3854 :
3855 : // Update the filename.
3856 2 : if (strncmp(osFileName, pszOldBase, strlen(pszOldBase)) == 0)
3857 : {
3858 2 : std::string osNew = pszNewBase;
3859 2 : osNew += (osFileName.c_str() + strlen(pszOldBase));
3860 2 : osFileName = std::move(osNew);
3861 : }
3862 :
3863 2 : apoNodeList[iNode]->SetStringField("dependent.string", osFileName);
3864 : }
3865 :
3866 8 : return CE_None;
3867 : }
3868 :
3869 : /* ==================================================================== */
3870 : /* Table relating USGS and ESRI state plane zones. */
3871 : /* ==================================================================== */
3872 : constexpr int anUsgsEsriZones[] = {
3873 : 101, 3101, 102, 3126, 201, 3151, 202, 3176, 203, 3201, 301, 3226,
3874 : 302, 3251, 401, 3276, 402, 3301, 403, 3326, 404, 3351, 405, 3376,
3875 : 406, 3401, 407, 3426, 501, 3451, 502, 3476, 503, 3501, 600, 3526,
3876 : 700, 3551, 901, 3601, 902, 3626, 903, 3576, 1001, 3651, 1002, 3676,
3877 : 1101, 3701, 1102, 3726, 1103, 3751, 1201, 3776, 1202, 3801, 1301, 3826,
3878 : 1302, 3851, 1401, 3876, 1402, 3901, 1501, 3926, 1502, 3951, 1601, 3976,
3879 : 1602, 4001, 1701, 4026, 1702, 4051, 1703, 6426, 1801, 4076, 1802, 4101,
3880 : 1900, 4126, 2001, 4151, 2002, 4176, 2101, 4201, 2102, 4226, 2103, 4251,
3881 : 2111, 6351, 2112, 6376, 2113, 6401, 2201, 4276, 2202, 4301, 2203, 4326,
3882 : 2301, 4351, 2302, 4376, 2401, 4401, 2402, 4426, 2403, 4451, 2500, 0,
3883 : 2501, 4476, 2502, 4501, 2503, 4526, 2600, 0, 2601, 4551, 2602, 4576,
3884 : 2701, 4601, 2702, 4626, 2703, 4651, 2800, 4676, 2900, 4701, 3001, 4726,
3885 : 3002, 4751, 3003, 4776, 3101, 4801, 3102, 4826, 3103, 4851, 3104, 4876,
3886 : 3200, 4901, 3301, 4926, 3302, 4951, 3401, 4976, 3402, 5001, 3501, 5026,
3887 : 3502, 5051, 3601, 5076, 3602, 5101, 3701, 5126, 3702, 5151, 3800, 5176,
3888 : 3900, 0, 3901, 5201, 3902, 5226, 4001, 5251, 4002, 5276, 4100, 5301,
3889 : 4201, 5326, 4202, 5351, 4203, 5376, 4204, 5401, 4205, 5426, 4301, 5451,
3890 : 4302, 5476, 4303, 5501, 4400, 5526, 4501, 5551, 4502, 5576, 4601, 5601,
3891 : 4602, 5626, 4701, 5651, 4702, 5676, 4801, 5701, 4802, 5726, 4803, 5751,
3892 : 4901, 5776, 4902, 5801, 4903, 5826, 4904, 5851, 5001, 6101, 5002, 6126,
3893 : 5003, 6151, 5004, 6176, 5005, 6201, 5006, 6226, 5007, 6251, 5008, 6276,
3894 : 5009, 6301, 5010, 6326, 5101, 5876, 5102, 5901, 5103, 5926, 5104, 5951,
3895 : 5105, 5976, 5201, 6001, 5200, 6026, 5200, 6076, 5201, 6051, 5202, 6051,
3896 : 5300, 0, 5400, 0};
3897 :
3898 : /************************************************************************/
3899 : /* ESRIToUSGSZone() */
3900 : /* */
3901 : /* Convert ESRI style state plane zones to USGS style state */
3902 : /* plane zones. */
3903 : /************************************************************************/
3904 :
3905 3 : static int ESRIToUSGSZone(int nESRIZone)
3906 :
3907 : {
3908 3 : if (nESRIZone == INT_MIN)
3909 0 : return 0;
3910 3 : if (nESRIZone < 0)
3911 1 : return std::abs(nESRIZone);
3912 :
3913 2 : const int nPairs = sizeof(anUsgsEsriZones) / (2 * sizeof(int));
3914 212 : for (int i = 0; i < nPairs; i++)
3915 : {
3916 212 : if (anUsgsEsriZones[i * 2 + 1] == nESRIZone)
3917 2 : return anUsgsEsriZones[i * 2];
3918 : }
3919 :
3920 0 : return 0;
3921 : }
3922 :
3923 : static const char *const apszDatumMap[] = {
3924 : // Imagine name, WKT name.
3925 : "NAD27",
3926 : "North_American_Datum_1927",
3927 : "NAD83",
3928 : "North_American_Datum_1983",
3929 : "WGS 84",
3930 : "WGS_1984",
3931 : "WGS 1972",
3932 : "WGS_1972",
3933 : "GDA94",
3934 : "Geocentric_Datum_of_Australia_1994",
3935 : "Pulkovo 1942",
3936 : "Pulkovo_1942",
3937 : "Geodetic Datum 1949",
3938 : "New_Zealand_Geodetic_Datum_1949",
3939 : nullptr,
3940 : nullptr};
3941 :
3942 336 : const char *const *HFAGetDatumMap()
3943 : {
3944 336 : return apszDatumMap;
3945 : }
3946 :
3947 : static const char *const apszUnitMap[] = {"meters",
3948 : "1.0",
3949 : "meter",
3950 : "1.0",
3951 : "m",
3952 : "1.0",
3953 : "centimeters",
3954 : "0.01",
3955 : "centimeter",
3956 : "0.01",
3957 : "cm",
3958 : "0.01",
3959 : "millimeters",
3960 : "0.001",
3961 : "millimeter",
3962 : "0.001",
3963 : "mm",
3964 : "0.001",
3965 : "kilometers",
3966 : "1000.0",
3967 : "kilometer",
3968 : "1000.0",
3969 : "km",
3970 : "1000.0",
3971 : "us_survey_feet",
3972 : "0.3048006096012192",
3973 : "us_survey_foot",
3974 : "0.3048006096012192",
3975 : "feet",
3976 : "0.3048006096012192",
3977 : "foot",
3978 : "0.3048006096012192",
3979 : "ft",
3980 : "0.3048006096012192",
3981 : "international_feet",
3982 : "0.3048",
3983 : "international_foot",
3984 : "0.3048",
3985 : "inches",
3986 : "0.0254000508001",
3987 : "inch",
3988 : "0.0254000508001",
3989 : "in",
3990 : "0.0254000508001",
3991 : "yards",
3992 : "0.9144",
3993 : "yard",
3994 : "0.9144",
3995 : "yd",
3996 : "0.9144",
3997 : "clarke_yard",
3998 : "0.9143917962",
3999 : "miles",
4000 : "1304.544",
4001 : "mile",
4002 : "1304.544",
4003 : "mi",
4004 : "1304.544",
4005 : "modified_american_feet",
4006 : "0.3048122530",
4007 : "modified_american_foot",
4008 : "0.3048122530",
4009 : "clarke_feet",
4010 : "0.3047972651",
4011 : "clarke_foot",
4012 : "0.3047972651",
4013 : "indian_feet",
4014 : "0.3047995142",
4015 : "indian_foot",
4016 : "0.3047995142",
4017 : nullptr,
4018 : nullptr};
4019 :
4020 196 : const char *const *HFAGetUnitMap()
4021 : {
4022 196 : return apszUnitMap;
4023 : }
4024 :
4025 : /************************************************************************/
4026 : /* HFAPCSStructToOSR() */
4027 : /* */
4028 : /* Convert the datum, proparameters and mapinfo structures into */
4029 : /* WKT format. */
4030 : /************************************************************************/
4031 :
4032 : std::unique_ptr<OGRSpatialReference>
4033 209 : HFAPCSStructToOSR(const Eprj_Datum *psDatum, const Eprj_ProParameters *psPro,
4034 : const Eprj_MapInfo *psMapInfo, HFAEntry *poMapInformation)
4035 :
4036 : {
4037 : // General case for Erdas style projections.
4038 :
4039 : // We make a particular effort to adapt the mapinfo->proname as
4040 : // the PROJCS[] name per #2422.
4041 418 : auto poSRS = std::make_unique<OGRSpatialReference>();
4042 209 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4043 :
4044 209 : if (psPro == nullptr && psMapInfo != nullptr)
4045 : {
4046 6 : poSRS->SetLocalCS(psMapInfo->proName);
4047 : }
4048 203 : else if (psPro == nullptr)
4049 : {
4050 0 : return nullptr;
4051 : }
4052 203 : else if (psPro->proType == EPRJ_EXTERNAL)
4053 : {
4054 3 : if (EQUALN(psPro->proExeName, EPRJ_EXTERNAL_NZMG, 4))
4055 : {
4056 : // Handle New Zealand Map Grid (NZMG) external projection. See:
4057 : // http://www.linz.govt.nz/
4058 : //
4059 : // Is there a better way that doesn't require hardcoding
4060 : // of these numbers?
4061 3 : poSRS->SetNZMG(-41.0, 173.0, 2510000, 6023150);
4062 : }
4063 : else
4064 : {
4065 0 : poSRS->SetLocalCS(psPro->proName);
4066 : }
4067 : }
4068 200 : else if (psPro->proNumber != EPRJ_LATLONG && psMapInfo != nullptr)
4069 : {
4070 167 : poSRS->SetProjCS(psMapInfo->proName);
4071 : }
4072 33 : else if (psPro->proNumber != EPRJ_LATLONG)
4073 : {
4074 7 : poSRS->SetProjCS(psPro->proName);
4075 : }
4076 :
4077 : // Handle units. It is important to deal with this first so
4078 : // that the projection Set methods will automatically do
4079 : // translation of linear values (like false easting) to PROJCS
4080 : // units from meters. Erdas linear projection values are
4081 : // always in meters.
4082 209 : if (poSRS->IsProjected() || poSRS->IsLocal())
4083 : {
4084 183 : const char *pszUnits = nullptr;
4085 :
4086 183 : if (psMapInfo)
4087 176 : pszUnits = psMapInfo->units;
4088 7 : else if (poMapInformation != nullptr)
4089 6 : pszUnits = poMapInformation->GetStringField("units.string");
4090 :
4091 183 : if (pszUnits != nullptr)
4092 : {
4093 182 : const char *const *papszUnitMap = HFAGetUnitMap();
4094 182 : int iUnitIndex = 0; // Used after for.
4095 303 : for (; papszUnitMap[iUnitIndex] != nullptr; iUnitIndex += 2)
4096 : {
4097 303 : if (EQUAL(papszUnitMap[iUnitIndex], pszUnits))
4098 182 : break;
4099 : }
4100 :
4101 182 : if (papszUnitMap[iUnitIndex] == nullptr)
4102 0 : iUnitIndex = 0;
4103 :
4104 364 : poSRS->SetLinearUnits(pszUnits,
4105 182 : CPLAtof(papszUnitMap[iUnitIndex + 1]));
4106 : }
4107 : else
4108 : {
4109 1 : poSRS->SetLinearUnits(SRS_UL_METER, 1.0);
4110 : }
4111 : }
4112 :
4113 209 : if (psPro == nullptr)
4114 : {
4115 6 : if (poSRS->IsLocal())
4116 : {
4117 6 : return poSRS;
4118 : }
4119 : else
4120 0 : return nullptr;
4121 : }
4122 :
4123 : // Try to work out ellipsoid and datum information.
4124 203 : const char *pszDatumName = psPro->proSpheroid.sphereName;
4125 203 : const char *pszEllipsoidName = psPro->proSpheroid.sphereName;
4126 :
4127 203 : if (psDatum != nullptr)
4128 : {
4129 203 : pszDatumName = psDatum->datumname;
4130 :
4131 : // Imagine to WKT translation.
4132 203 : const char *const *papszDatumMap = HFAGetDatumMap();
4133 734 : for (int i = 0; papszDatumMap[i] != nullptr; i += 2)
4134 : {
4135 682 : if (EQUAL(pszDatumName, papszDatumMap[i]))
4136 : {
4137 151 : pszDatumName = papszDatumMap[i + 1];
4138 151 : break;
4139 : }
4140 : }
4141 : }
4142 :
4143 203 : if (psPro->proSpheroid.a == 0.0)
4144 0 : ((Eprj_ProParameters *)psPro)->proSpheroid.a = 6378137.0;
4145 203 : if (psPro->proSpheroid.b == 0.0)
4146 0 : ((Eprj_ProParameters *)psPro)->proSpheroid.b = 6356752.3;
4147 :
4148 : const double dfInvFlattening =
4149 203 : OSRCalcInvFlattening(psPro->proSpheroid.a, psPro->proSpheroid.b);
4150 :
4151 : // Handle different projection methods.
4152 203 : switch (psPro->proNumber)
4153 : {
4154 29 : case EPRJ_LATLONG:
4155 29 : break;
4156 :
4157 90 : case EPRJ_UTM:
4158 : // We change this to unnamed so that SetUTM will set the long
4159 : // UTM description.
4160 90 : poSRS->SetProjCS("unnamed");
4161 90 : poSRS->SetUTM(psPro->proZone, psPro->proParams[3] >= 0.0);
4162 :
4163 : // The PCS name from the above function may be different with the
4164 : // input name. If there is a PCS name in psMapInfo that is
4165 : // different with the one in psPro, just use it as the PCS name.
4166 : // This case happens if the dataset's SR was written by the new
4167 : // GDAL.
4168 90 : if (psMapInfo && strlen(psMapInfo->proName) > 0 &&
4169 88 : strlen(psPro->proName) > 0 &&
4170 88 : !EQUAL(psMapInfo->proName, psPro->proName))
4171 63 : poSRS->SetProjCS(psMapInfo->proName);
4172 90 : break;
4173 :
4174 3 : case EPRJ_STATE_PLANE:
4175 : {
4176 3 : CPLString osUnitsName;
4177 : double dfLinearUnits;
4178 : {
4179 3 : const char *pszUnitsName = nullptr;
4180 3 : dfLinearUnits = poSRS->GetLinearUnits(&pszUnitsName);
4181 3 : if (pszUnitsName)
4182 3 : osUnitsName = pszUnitsName;
4183 : }
4184 :
4185 : // Historically, hfa used esri state plane zone code. Try esri pe
4186 : // string first.
4187 3 : const int zoneCode = ESRIToUSGSZone(psPro->proZone);
4188 : const char *pszDatum;
4189 3 : if (psDatum)
4190 3 : pszDatum = psDatum->datumname;
4191 : else
4192 0 : pszDatum = "HARN";
4193 : const char *pszUnits;
4194 3 : if (psMapInfo)
4195 0 : pszUnits = psMapInfo->units;
4196 3 : else if (!osUnitsName.empty())
4197 3 : pszUnits = osUnitsName;
4198 : else
4199 0 : pszUnits = "meters";
4200 3 : const int proNu = psPro->proNumber;
4201 3 : if (poSRS->ImportFromESRIStatePlaneWKT(zoneCode, pszDatum, pszUnits,
4202 3 : proNu) == OGRERR_NONE)
4203 : {
4204 3 : poSRS->AutoIdentifyEPSG();
4205 :
4206 3 : return poSRS;
4207 : }
4208 :
4209 : // Set state plane zone. Set NAD83/27 on basis of spheroid.
4210 0 : poSRS->SetStatePlane(ESRIToUSGSZone(psPro->proZone),
4211 0 : fabs(psPro->proSpheroid.a - 6378137.0) < 1.0,
4212 0 : osUnitsName.empty() ? nullptr
4213 0 : : osUnitsName.c_str(),
4214 : dfLinearUnits);
4215 :
4216 : // Same as the UTM, The following is needed.
4217 0 : if (psMapInfo && strlen(psMapInfo->proName) > 0 &&
4218 0 : strlen(psPro->proName) > 0 &&
4219 0 : !EQUAL(psMapInfo->proName, psPro->proName))
4220 0 : poSRS->SetProjCS(psMapInfo->proName);
4221 : }
4222 0 : break;
4223 :
4224 1 : case EPRJ_ALBERS_CONIC_EQUAL_AREA:
4225 1 : poSRS->SetACEA(psPro->proParams[2] * R2D, psPro->proParams[3] * R2D,
4226 1 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4227 1 : psPro->proParams[6], psPro->proParams[7]);
4228 1 : break;
4229 :
4230 6 : case EPRJ_LAMBERT_CONFORMAL_CONIC:
4231 : // Check the possible Wisconsin first.
4232 6 : if (psDatum && psMapInfo && EQUAL(psDatum->datumname, "HARN"))
4233 : {
4234 : // ERO: I doubt this works. Wisconsin LCC is LCC_1SP whereas
4235 : // we are here in the LCC_2SP case...
4236 0 : if (poSRS->ImportFromESRIWisconsinWKT(
4237 0 : "Lambert_Conformal_Conic", psPro->proParams[4] * R2D,
4238 0 : psPro->proParams[5] * R2D,
4239 0 : psMapInfo->units) == OGRERR_NONE)
4240 : {
4241 0 : return poSRS;
4242 : }
4243 : }
4244 6 : poSRS->SetLCC(psPro->proParams[2] * R2D, psPro->proParams[3] * R2D,
4245 6 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4246 6 : psPro->proParams[6], psPro->proParams[7]);
4247 6 : break;
4248 :
4249 1 : case EPRJ_MERCATOR:
4250 1 : poSRS->SetMercator(psPro->proParams[5] * R2D,
4251 1 : psPro->proParams[4] * R2D, 1.0,
4252 1 : psPro->proParams[6], psPro->proParams[7]);
4253 1 : break;
4254 :
4255 2 : case EPRJ_POLAR_STEREOGRAPHIC:
4256 2 : poSRS->SetPS(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4257 2 : 1.0, psPro->proParams[6], psPro->proParams[7]);
4258 2 : break;
4259 :
4260 2 : case EPRJ_POLYCONIC:
4261 2 : poSRS->SetPolyconic(psPro->proParams[5] * R2D,
4262 2 : psPro->proParams[4] * R2D, psPro->proParams[6],
4263 2 : psPro->proParams[7]);
4264 2 : break;
4265 :
4266 1 : case EPRJ_EQUIDISTANT_CONIC:
4267 : {
4268 2 : const double dfStdParallel2 = psPro->proParams[8] != 0.0
4269 1 : ? psPro->proParams[3] * R2D
4270 0 : : psPro->proParams[2] * R2D;
4271 1 : poSRS->SetEC(psPro->proParams[2] * R2D, dfStdParallel2,
4272 1 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4273 1 : psPro->proParams[6], psPro->proParams[7]);
4274 1 : break;
4275 : }
4276 18 : case EPRJ_TRANSVERSE_MERCATOR:
4277 : case EPRJ_GAUSS_KRUGER:
4278 : // Check the possible Wisconsin first.
4279 18 : if (psDatum && psMapInfo && EQUAL(psDatum->datumname, "HARN"))
4280 : {
4281 1 : if (poSRS->ImportFromESRIWisconsinWKT(
4282 1 : "Transverse_Mercator", psPro->proParams[4] * R2D,
4283 1 : psPro->proParams[5] * R2D,
4284 1 : psMapInfo->units) == OGRERR_NONE)
4285 : {
4286 1 : return poSRS;
4287 : }
4288 : }
4289 17 : poSRS->SetTM(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4290 17 : psPro->proParams[2], psPro->proParams[6],
4291 17 : psPro->proParams[7]);
4292 17 : break;
4293 :
4294 0 : case EPRJ_STEREOGRAPHIC:
4295 0 : poSRS->SetStereographic(psPro->proParams[5] * R2D,
4296 0 : psPro->proParams[4] * R2D, 1.0,
4297 0 : psPro->proParams[6], psPro->proParams[7]);
4298 0 : break;
4299 :
4300 1 : case EPRJ_LAMBERT_AZIMUTHAL_EQUAL_AREA:
4301 1 : poSRS->SetLAEA(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4302 1 : psPro->proParams[6], psPro->proParams[7]);
4303 1 : break;
4304 :
4305 1 : case EPRJ_AZIMUTHAL_EQUIDISTANT:
4306 1 : poSRS->SetAE(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4307 1 : psPro->proParams[6], psPro->proParams[7]);
4308 1 : break;
4309 :
4310 1 : case EPRJ_GNOMONIC:
4311 1 : poSRS->SetGnomonic(psPro->proParams[5] * R2D,
4312 1 : psPro->proParams[4] * R2D, psPro->proParams[6],
4313 1 : psPro->proParams[7]);
4314 1 : break;
4315 :
4316 2 : case EPRJ_ORTHOGRAPHIC:
4317 2 : poSRS->SetOrthographic(psPro->proParams[5] * R2D,
4318 2 : psPro->proParams[4] * R2D,
4319 2 : psPro->proParams[6], psPro->proParams[7]);
4320 2 : break;
4321 :
4322 1 : case EPRJ_SINUSOIDAL:
4323 1 : poSRS->SetSinusoidal(psPro->proParams[4] * R2D, psPro->proParams[6],
4324 1 : psPro->proParams[7]);
4325 1 : break;
4326 :
4327 3 : case EPRJ_PLATE_CARREE:
4328 : case EPRJ_EQUIRECTANGULAR:
4329 3 : poSRS->SetEquirectangular2(
4330 3 : 0.0, psPro->proParams[4] * R2D, psPro->proParams[5] * R2D,
4331 3 : psPro->proParams[6], psPro->proParams[7]);
4332 3 : break;
4333 :
4334 0 : case EPRJ_EQUIDISTANT_CYLINDRICAL:
4335 0 : poSRS->SetEquirectangular2(
4336 0 : 0.0, psPro->proParams[4] * R2D, psPro->proParams[2] * R2D,
4337 0 : psPro->proParams[6], psPro->proParams[7]);
4338 0 : break;
4339 :
4340 1 : case EPRJ_MILLER_CYLINDRICAL:
4341 1 : poSRS->SetMC(0.0, psPro->proParams[4] * R2D, psPro->proParams[6],
4342 1 : psPro->proParams[7]);
4343 1 : break;
4344 :
4345 1 : case EPRJ_VANDERGRINTEN:
4346 1 : poSRS->SetVDG(psPro->proParams[4] * R2D, psPro->proParams[6],
4347 1 : psPro->proParams[7]);
4348 1 : break;
4349 :
4350 0 : case EPRJ_HOTINE_OBLIQUE_MERCATOR:
4351 0 : if (psPro->proParams[12] > 0.0)
4352 0 : poSRS->SetHOM(
4353 0 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4354 0 : psPro->proParams[3] * R2D, 0.0, psPro->proParams[2],
4355 0 : psPro->proParams[6], psPro->proParams[7]);
4356 0 : break;
4357 :
4358 2 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER:
4359 2 : poSRS->SetHOMAC(
4360 2 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4361 2 : psPro->proParams[3] * R2D,
4362 2 : psPro->proParams[3] *
4363 : R2D, // We reuse azimuth as rectified_grid_angle
4364 2 : psPro->proParams[2], psPro->proParams[6], psPro->proParams[7]);
4365 2 : break;
4366 :
4367 1 : case EPRJ_ROBINSON:
4368 1 : poSRS->SetRobinson(psPro->proParams[4] * R2D, psPro->proParams[6],
4369 1 : psPro->proParams[7]);
4370 1 : break;
4371 :
4372 1 : case EPRJ_MOLLWEIDE:
4373 1 : poSRS->SetMollweide(psPro->proParams[4] * R2D, psPro->proParams[6],
4374 1 : psPro->proParams[7]);
4375 1 : break;
4376 :
4377 1 : case EPRJ_GALL_STEREOGRAPHIC:
4378 1 : poSRS->SetGS(psPro->proParams[4] * R2D, psPro->proParams[6],
4379 1 : psPro->proParams[7]);
4380 1 : break;
4381 :
4382 1 : case EPRJ_ECKERT_I:
4383 1 : poSRS->SetEckert(1, psPro->proParams[4] * R2D, psPro->proParams[6],
4384 1 : psPro->proParams[7]);
4385 1 : break;
4386 :
4387 1 : case EPRJ_ECKERT_II:
4388 1 : poSRS->SetEckert(2, psPro->proParams[4] * R2D, psPro->proParams[6],
4389 1 : psPro->proParams[7]);
4390 1 : break;
4391 :
4392 1 : case EPRJ_ECKERT_III:
4393 1 : poSRS->SetEckert(3, psPro->proParams[4] * R2D, psPro->proParams[6],
4394 1 : psPro->proParams[7]);
4395 1 : break;
4396 :
4397 1 : case EPRJ_ECKERT_IV:
4398 1 : poSRS->SetEckert(4, psPro->proParams[4] * R2D, psPro->proParams[6],
4399 1 : psPro->proParams[7]);
4400 1 : break;
4401 :
4402 1 : case EPRJ_ECKERT_V:
4403 1 : poSRS->SetEckert(5, psPro->proParams[4] * R2D, psPro->proParams[6],
4404 1 : psPro->proParams[7]);
4405 1 : break;
4406 :
4407 1 : case EPRJ_ECKERT_VI:
4408 1 : poSRS->SetEckert(6, psPro->proParams[4] * R2D, psPro->proParams[6],
4409 1 : psPro->proParams[7]);
4410 1 : break;
4411 :
4412 2 : case EPRJ_CASSINI:
4413 2 : poSRS->SetCS(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4414 2 : psPro->proParams[6], psPro->proParams[7]);
4415 2 : break;
4416 :
4417 1 : case EPRJ_TWO_POINT_EQUIDISTANT:
4418 1 : poSRS->SetTPED(psPro->proParams[9] * R2D, psPro->proParams[8] * R2D,
4419 1 : psPro->proParams[11] * R2D,
4420 1 : psPro->proParams[10] * R2D, psPro->proParams[6],
4421 1 : psPro->proParams[7]);
4422 1 : break;
4423 :
4424 0 : case EPRJ_STEREOGRAPHIC_EXTENDED:
4425 0 : poSRS->SetStereographic(
4426 0 : psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4427 0 : psPro->proParams[2], psPro->proParams[6], psPro->proParams[7]);
4428 0 : break;
4429 :
4430 1 : case EPRJ_BONNE:
4431 1 : poSRS->SetBonne(psPro->proParams[2] * R2D,
4432 1 : psPro->proParams[4] * R2D, psPro->proParams[6],
4433 1 : psPro->proParams[7]);
4434 1 : break;
4435 :
4436 1 : case EPRJ_LOXIMUTHAL:
4437 : {
4438 1 : poSRS->SetProjection("Loximuthal");
4439 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4440 1 : psPro->proParams[4] * R2D);
4441 1 : poSRS->SetNormProjParm("central_parallel",
4442 1 : psPro->proParams[5] * R2D);
4443 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4444 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4445 : }
4446 1 : break;
4447 :
4448 1 : case EPRJ_QUARTIC_AUTHALIC:
4449 : {
4450 1 : poSRS->SetProjection("Quartic_Authalic");
4451 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4452 1 : psPro->proParams[4] * R2D);
4453 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4454 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4455 : }
4456 1 : break;
4457 :
4458 1 : case EPRJ_WINKEL_I:
4459 : {
4460 1 : poSRS->SetProjection("Winkel_I");
4461 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4462 1 : psPro->proParams[4] * R2D);
4463 1 : poSRS->SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,
4464 1 : psPro->proParams[2] * R2D);
4465 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4466 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4467 : }
4468 1 : break;
4469 :
4470 1 : case EPRJ_WINKEL_II:
4471 : {
4472 1 : poSRS->SetProjection("Winkel_II");
4473 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4474 1 : psPro->proParams[4] * R2D);
4475 1 : poSRS->SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,
4476 1 : psPro->proParams[2] * R2D);
4477 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4478 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4479 : }
4480 1 : break;
4481 :
4482 0 : case EPRJ_BEHRMANN:
4483 : {
4484 0 : poSRS->SetProjection("Behrmann");
4485 0 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4486 0 : psPro->proParams[4] * R2D);
4487 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4488 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4489 : }
4490 0 : break;
4491 :
4492 1 : case EPRJ_KROVAK:
4493 1 : poSRS->SetKrovak(
4494 1 : psPro->proParams[4] * R2D, psPro->proParams[5] * R2D,
4495 1 : psPro->proParams[3] * R2D, psPro->proParams[9] * R2D,
4496 1 : psPro->proParams[2], psPro->proParams[6], psPro->proParams[7]);
4497 1 : break;
4498 :
4499 0 : case EPRJ_DOUBLE_STEREOGRAPHIC:
4500 : {
4501 0 : poSRS->SetProjection("Double_Stereographic");
4502 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,
4503 0 : psPro->proParams[5] * R2D);
4504 0 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4505 0 : psPro->proParams[4] * R2D);
4506 0 : poSRS->SetNormProjParm(SRS_PP_SCALE_FACTOR, psPro->proParams[2]);
4507 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4508 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4509 : }
4510 0 : break;
4511 :
4512 1 : case EPRJ_AITOFF:
4513 : {
4514 1 : poSRS->SetProjection("Aitoff");
4515 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4516 1 : psPro->proParams[4] * R2D);
4517 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4518 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4519 : }
4520 1 : break;
4521 :
4522 1 : case EPRJ_CRASTER_PARABOLIC:
4523 : {
4524 1 : poSRS->SetProjection("Craster_Parabolic");
4525 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4526 1 : psPro->proParams[4] * R2D);
4527 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4528 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4529 : }
4530 1 : break;
4531 :
4532 1 : case EPRJ_CYLINDRICAL_EQUAL_AREA:
4533 1 : poSRS->SetCEA(psPro->proParams[2] * R2D, psPro->proParams[4] * R2D,
4534 1 : psPro->proParams[6], psPro->proParams[7]);
4535 1 : break;
4536 :
4537 1 : case EPRJ_FLAT_POLAR_QUARTIC:
4538 : {
4539 1 : poSRS->SetProjection("Flat_Polar_Quartic");
4540 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4541 1 : psPro->proParams[4] * R2D);
4542 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4543 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4544 : }
4545 1 : break;
4546 :
4547 1 : case EPRJ_TIMES:
4548 : {
4549 1 : poSRS->SetProjection("Times");
4550 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4551 1 : psPro->proParams[4] * R2D);
4552 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4553 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4554 : }
4555 1 : break;
4556 :
4557 1 : case EPRJ_WINKEL_TRIPEL:
4558 : {
4559 1 : poSRS->SetProjection("Winkel_Tripel");
4560 1 : poSRS->SetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,
4561 1 : psPro->proParams[2] * R2D);
4562 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4563 1 : psPro->proParams[4] * R2D);
4564 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4565 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4566 : }
4567 1 : break;
4568 :
4569 1 : case EPRJ_HAMMER_AITOFF:
4570 : {
4571 1 : poSRS->SetProjection("Hammer_Aitoff");
4572 1 : poSRS->SetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,
4573 1 : psPro->proParams[4] * R2D);
4574 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4575 1 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4576 : }
4577 1 : break;
4578 :
4579 1 : case EPRJ_VERTICAL_NEAR_SIDE_PERSPECTIVE:
4580 : {
4581 1 : poSRS->SetVerticalPerspective(
4582 1 : psPro->proParams[5] * R2D, // dfTopoOriginLat
4583 1 : psPro->proParams[4] * R2D, // dfTopoOriginLon
4584 : 0, // dfTopoOriginHeight
4585 1 : psPro->proParams[2], // dfViewPointHeight
4586 1 : psPro->proParams[6], // dfFalseEasting
4587 1 : psPro->proParams[7]); // dfFalseNorthing
4588 : }
4589 1 : break;
4590 :
4591 0 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_CENTER:
4592 : {
4593 0 : poSRS->SetProjection("Hotine_Oblique_Mercator_Twp_Point_Center");
4594 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_CENTER,
4595 0 : psPro->proParams[5] * R2D);
4596 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_1ST_POINT,
4597 0 : psPro->proParams[9] * R2D);
4598 0 : poSRS->SetNormProjParm(SRS_PP_LONGITUDE_OF_1ST_POINT,
4599 0 : psPro->proParams[8] * R2D);
4600 0 : poSRS->SetNormProjParm(SRS_PP_LATITUDE_OF_2ND_POINT,
4601 0 : psPro->proParams[11] * R2D);
4602 0 : poSRS->SetNormProjParm(SRS_PP_LONGITUDE_OF_2ND_POINT,
4603 0 : psPro->proParams[10] * R2D);
4604 0 : poSRS->SetNormProjParm(SRS_PP_SCALE_FACTOR, psPro->proParams[2]);
4605 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_EASTING, psPro->proParams[6]);
4606 0 : poSRS->SetNormProjParm(SRS_PP_FALSE_NORTHING, psPro->proParams[7]);
4607 : }
4608 0 : break;
4609 :
4610 1 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN:
4611 1 : poSRS->SetHOM2PNO(
4612 1 : psPro->proParams[5] * R2D, psPro->proParams[8] * R2D,
4613 1 : psPro->proParams[9] * R2D, psPro->proParams[10] * R2D,
4614 1 : psPro->proParams[11] * R2D, psPro->proParams[2],
4615 1 : psPro->proParams[6], psPro->proParams[7]);
4616 1 : break;
4617 :
4618 0 : case EPRJ_LAMBERT_CONFORMAL_CONIC_1SP:
4619 0 : poSRS->SetLCC1SP(psPro->proParams[3] * R2D,
4620 0 : psPro->proParams[2] * R2D, psPro->proParams[4],
4621 0 : psPro->proParams[5], psPro->proParams[6]);
4622 0 : break;
4623 :
4624 2 : case EPRJ_MERCATOR_VARIANT_A:
4625 2 : poSRS->SetMercator(psPro->proParams[5] * R2D,
4626 2 : psPro->proParams[4] * R2D, psPro->proParams[2],
4627 2 : psPro->proParams[6], psPro->proParams[7]);
4628 2 : break;
4629 :
4630 0 : case EPRJ_PSEUDO_MERCATOR: // Likely this is google mercator?
4631 0 : poSRS->SetMercator(psPro->proParams[5] * R2D,
4632 0 : psPro->proParams[4] * R2D, 1.0,
4633 0 : psPro->proParams[6], psPro->proParams[7]);
4634 0 : break;
4635 :
4636 5 : case EPRJ_HOTINE_OBLIQUE_MERCATOR_VARIANT_A:
4637 5 : poSRS->SetHOM(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4638 5 : psPro->proParams[3] * R2D, psPro->proParams[8] * R2D,
4639 5 : psPro->proParams[2], psPro->proParams[6],
4640 5 : psPro->proParams[7]);
4641 5 : break;
4642 :
4643 3 : case EPRJ_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED:
4644 3 : poSRS->SetTMSO(psPro->proParams[5] * R2D, psPro->proParams[4] * R2D,
4645 3 : psPro->proParams[2], psPro->proParams[6],
4646 3 : psPro->proParams[7]);
4647 3 : break;
4648 :
4649 0 : default:
4650 0 : if (poSRS->IsProjected())
4651 0 : poSRS->GetRoot()->SetValue("LOCAL_CS");
4652 : else
4653 0 : poSRS->SetLocalCS(psPro->proName);
4654 0 : break;
4655 : }
4656 :
4657 : // Try and set the GeogCS information.
4658 199 : if (!poSRS->IsLocal())
4659 : {
4660 199 : bool bWellKnownDatum = false;
4661 199 : if (pszDatumName == nullptr)
4662 0 : poSRS->SetGeogCS(pszDatumName, pszDatumName, pszEllipsoidName,
4663 0 : psPro->proSpheroid.a, dfInvFlattening);
4664 199 : else if (EQUAL(pszDatumName, "WGS 84") ||
4665 199 : EQUAL(pszDatumName, "WGS_1984"))
4666 : {
4667 36 : bWellKnownDatum = true;
4668 36 : poSRS->SetWellKnownGeogCS("WGS84");
4669 : }
4670 163 : else if (strstr(pszDatumName, "NAD27") != nullptr ||
4671 161 : EQUAL(pszDatumName, "North_American_Datum_1927"))
4672 : {
4673 88 : bWellKnownDatum = true;
4674 88 : poSRS->SetWellKnownGeogCS("NAD27");
4675 : }
4676 75 : else if (EQUAL(pszDatumName, "NAD83") ||
4677 75 : EQUAL(pszDatumName, "North_American_Datum_1983"))
4678 : {
4679 7 : bWellKnownDatum = true;
4680 7 : poSRS->SetWellKnownGeogCS("NAD83");
4681 : }
4682 : else
4683 : {
4684 136 : CPLString osGeogCRSName(pszDatumName);
4685 :
4686 68 : if (poSRS->IsProjected())
4687 : {
4688 67 : PJ_CONTEXT *ctxt = OSRGetProjTLSContext();
4689 67 : const PJ_TYPE type = PJ_TYPE_PROJECTED_CRS;
4690 : PJ_OBJ_LIST *list =
4691 67 : proj_create_from_name(ctxt, nullptr, poSRS->GetName(),
4692 : &type, 1, false, 1, nullptr);
4693 67 : if (list)
4694 : {
4695 67 : const auto listSize = proj_list_get_count(list);
4696 67 : if (listSize == 1)
4697 : {
4698 43 : auto crs = proj_list_get(ctxt, list, 0);
4699 43 : if (crs)
4700 : {
4701 43 : auto geogCRS = proj_crs_get_geodetic_crs(ctxt, crs);
4702 43 : if (geogCRS)
4703 : {
4704 43 : const char *pszName = proj_get_name(geogCRS);
4705 43 : if (pszName)
4706 43 : osGeogCRSName = pszName;
4707 43 : proj_destroy(geogCRS);
4708 : }
4709 43 : proj_destroy(crs);
4710 : }
4711 : }
4712 67 : proj_list_destroy(list);
4713 : }
4714 : }
4715 :
4716 136 : poSRS->SetGeogCS(osGeogCRSName, pszDatumName, pszEllipsoidName,
4717 68 : psPro->proSpheroid.a, dfInvFlattening);
4718 : }
4719 :
4720 199 : if (psDatum != nullptr && psDatum->type == EPRJ_DATUM_PARAMETRIC)
4721 : {
4722 158 : if (bWellKnownDatum &&
4723 45 : CPLTestBool(CPLGetConfigOption("OSR_STRIP_TOWGS84", "YES")))
4724 : {
4725 45 : CPLDebug("OSR",
4726 : "TOWGS84 information has been removed. "
4727 : "It can be kept by setting the OSR_STRIP_TOWGS84 "
4728 : "configuration option to NO");
4729 : }
4730 : else
4731 : {
4732 68 : poSRS->SetTOWGS84(psDatum->params[0], psDatum->params[1],
4733 68 : psDatum->params[2],
4734 68 : -psDatum->params[3] * RAD2ARCSEC,
4735 68 : -psDatum->params[4] * RAD2ARCSEC,
4736 68 : -psDatum->params[5] * RAD2ARCSEC,
4737 68 : psDatum->params[6] * 1e+6);
4738 68 : poSRS->StripTOWGS84IfKnownDatumAndAllowed();
4739 : }
4740 : }
4741 : }
4742 :
4743 : // Try to insert authority information if possible.
4744 199 : poSRS->AutoIdentifyEPSG();
4745 :
4746 199 : auto poSRSBestMatch = poSRS->FindBestMatch(90, nullptr, nullptr);
4747 199 : if (poSRSBestMatch)
4748 : {
4749 136 : poSRS.reset(poSRSBestMatch);
4750 : }
4751 :
4752 199 : return poSRS;
4753 : }
|
