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