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