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