Line data Source code
1 : /******************************************************************************
2 : *
3 : * Name: hfadataset.cpp
4 : * Project: Erdas Imagine Driver
5 : * Purpose: Main driver for Erdas Imagine format.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam
10 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "cpl_port.h"
16 : #include "hfadataset.h"
17 : #include "hfa_p.h"
18 :
19 : #include <cassert>
20 : #include <climits>
21 : #include <cmath>
22 : #include <cstddef>
23 : #include <cstdio>
24 : #include <cstdlib>
25 : #include <cstring>
26 : #if HAVE_FCNTL_H
27 : #include <fcntl.h>
28 : #endif
29 : #include <algorithm>
30 : #include <limits>
31 : #include <memory>
32 : #include <string>
33 : #include <vector>
34 :
35 : #include "cpl_conv.h"
36 : #include "cpl_error.h"
37 : #include "cpl_minixml.h"
38 : #include "cpl_progress.h"
39 : #include "cpl_string.h"
40 : #include "cpl_vsi.h"
41 : #include "gdal.h"
42 : #include "gdal_frmts.h"
43 : #include "gdal_pam.h"
44 : #include "gdal_priv.h"
45 : #include "gdal_rat.h"
46 : #include "hfa.h"
47 : #include "ogr_core.h"
48 : #include "ogr_spatialref.h"
49 : #include "ogr_srs_api.h"
50 :
51 : constexpr double D2R = M_PI / 180.0;
52 :
53 : constexpr double ARCSEC2RAD = M_PI / 648000.0;
54 :
55 : int WritePeStringIfNeeded(const OGRSpatialReference *poSRS, HFAHandle hHFA);
56 : void ClearSR(HFAHandle hHFA);
57 :
58 : /************************************************************************/
59 : /* HFARasterAttributeTable() */
60 : /************************************************************************/
61 :
62 622 : HFARasterAttributeTable::HFARasterAttributeTable(HFARasterBand *poBand,
63 622 : const char *pszName)
64 622 : : hHFA(poBand->hHFA),
65 1244 : poDT(poBand->hHFA->papoBand[poBand->nBand - 1]->poNode->GetNamedChild(
66 : pszName)),
67 1244 : osName(pszName), nBand(poBand->nBand), eAccess(poBand->GetAccess()),
68 : nRows(0), bLinearBinning(false), dfRow0Min(0.0), dfBinSize(0.0),
69 622 : eTableType(GRTT_THEMATIC)
70 : {
71 622 : if (poDT != nullptr)
72 : {
73 105 : nRows = std::max(0, poDT->GetIntField("numRows"));
74 :
75 : // Scan under table for columns.
76 447 : for (HFAEntry *poDTChild = poDT->GetChild(); poDTChild != nullptr;
77 342 : poDTChild = poDTChild->GetNext())
78 : {
79 342 : if (EQUAL(poDTChild->GetType(), "Edsc_BinFunction"))
80 : {
81 73 : const double dfMax = poDTChild->GetDoubleField("maxLimit");
82 73 : const double dfMin = poDTChild->GetDoubleField("minLimit");
83 73 : const int nBinCount = poDTChild->GetIntField("numBins");
84 :
85 73 : if (nBinCount == nRows && dfMax != dfMin && nBinCount > 1)
86 : {
87 : // Can't call SetLinearBinning since it will re-write
88 : // which we might not have permission to do.
89 69 : bLinearBinning = true;
90 69 : dfRow0Min = dfMin;
91 69 : dfBinSize = (dfMax - dfMin) / (nBinCount - 1);
92 : }
93 : }
94 :
95 342 : if (EQUAL(poDTChild->GetType(), "Edsc_BinFunction840"))
96 : {
97 : const char *pszValue =
98 13 : poDTChild->GetStringField("binFunction.type.string");
99 13 : if (pszValue && EQUAL(pszValue, "BFUnique"))
100 : {
101 13 : AddColumn("BinValues", GFT_Real, GFU_MinMax, 0, 0,
102 : poDTChild, true);
103 : }
104 : }
105 :
106 342 : if (!EQUAL(poDTChild->GetType(), "Edsc_Column"))
107 86 : continue;
108 :
109 256 : const int nOffset = poDTChild->GetIntField("columnDataPtr");
110 256 : const char *pszType = poDTChild->GetStringField("dataType");
111 256 : GDALRATFieldUsage eUsage = GFU_Generic;
112 256 : bool bConvertColors = false;
113 :
114 256 : if (pszType == nullptr || nOffset == 0)
115 0 : continue;
116 :
117 : GDALRATFieldType eType;
118 256 : if (EQUAL(pszType, "real"))
119 170 : eType = GFT_Real;
120 86 : else if (EQUAL(pszType, "string"))
121 42 : eType = GFT_String;
122 44 : else if (STARTS_WITH_CI(pszType, "int"))
123 44 : eType = GFT_Integer;
124 : else
125 0 : continue;
126 :
127 256 : if (EQUAL(poDTChild->GetName(), "Histogram"))
128 83 : eUsage = GFU_PixelCount;
129 173 : else if (EQUAL(poDTChild->GetName(), "Red"))
130 : {
131 11 : eUsage = GFU_Red;
132 : // Treat color columns as ints regardless
133 : // of how they are stored.
134 11 : bConvertColors = eType == GFT_Real;
135 11 : eType = GFT_Integer;
136 : }
137 162 : else if (EQUAL(poDTChild->GetName(), "Green"))
138 : {
139 11 : eUsage = GFU_Green;
140 11 : bConvertColors = eType == GFT_Real;
141 11 : eType = GFT_Integer;
142 : }
143 151 : else if (EQUAL(poDTChild->GetName(), "Blue"))
144 : {
145 11 : eUsage = GFU_Blue;
146 11 : bConvertColors = eType == GFT_Real;
147 11 : eType = GFT_Integer;
148 : }
149 140 : else if (EQUAL(poDTChild->GetName(), "Opacity"))
150 : {
151 11 : eUsage = GFU_Alpha;
152 11 : bConvertColors = eType == GFT_Real;
153 11 : eType = GFT_Integer;
154 : }
155 129 : else if (EQUAL(poDTChild->GetName(), "Class_Names"))
156 0 : eUsage = GFU_Name;
157 :
158 256 : if (eType == GFT_Real)
159 : {
160 126 : AddColumn(poDTChild->GetName(), GFT_Real, eUsage, nOffset,
161 : sizeof(double), poDTChild);
162 : }
163 130 : else if (eType == GFT_String)
164 : {
165 42 : int nMaxNumChars = poDTChild->GetIntField("maxNumChars");
166 42 : if (nMaxNumChars <= 0)
167 : {
168 0 : CPLError(CE_Failure, CPLE_AppDefined,
169 : "Invalid nMaxNumChars = %d for column %s",
170 : nMaxNumChars, poDTChild->GetName());
171 0 : nMaxNumChars = 1;
172 : }
173 42 : AddColumn(poDTChild->GetName(), GFT_String, eUsage, nOffset,
174 : nMaxNumChars, poDTChild);
175 : }
176 88 : else if (eType == GFT_Integer)
177 : {
178 88 : int nSize = sizeof(GInt32);
179 88 : if (bConvertColors)
180 44 : nSize = sizeof(double);
181 88 : AddColumn(poDTChild->GetName(), GFT_Integer, eUsage, nOffset,
182 : nSize, poDTChild, false, bConvertColors);
183 : }
184 : }
185 : }
186 622 : }
187 :
188 : /************************************************************************/
189 : /* ~HFARasterAttributeTable() */
190 : /************************************************************************/
191 :
192 1244 : HFARasterAttributeTable::~HFARasterAttributeTable()
193 : {
194 1244 : }
195 :
196 : /************************************************************************/
197 : /* Clone() */
198 : /************************************************************************/
199 :
200 17 : GDALRasterAttributeTable *HFARasterAttributeTable::Clone() const
201 : {
202 17 : if (nRows > 0 && GetColumnCount() > RAT_MAX_ELEM_FOR_CLONE / nRows)
203 0 : return nullptr;
204 :
205 : GDALDefaultRasterAttributeTable *poRAT =
206 17 : new GDALDefaultRasterAttributeTable();
207 :
208 44 : for (int iCol = 0; iCol < static_cast<int>(aoFields.size()); iCol++)
209 : {
210 27 : poRAT->CreateColumn(aoFields[iCol].sName, aoFields[iCol].eType,
211 27 : aoFields[iCol].eUsage);
212 27 : poRAT->SetRowCount(nRows);
213 :
214 27 : if (aoFields[iCol].eType == GFT_Integer)
215 : {
216 : int *panColData =
217 5 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nRows));
218 5 : if (panColData == nullptr)
219 : {
220 0 : delete poRAT;
221 0 : return nullptr;
222 : }
223 :
224 10 : if (const_cast<HFARasterAttributeTable *>(this)->ValuesIO(
225 5 : GF_Read, iCol, 0, nRows, panColData) != CE_None)
226 : {
227 0 : CPLFree(panColData);
228 0 : delete poRAT;
229 0 : return nullptr;
230 : }
231 :
232 501 : for (int iRow = 0; iRow < nRows; iRow++)
233 : {
234 496 : poRAT->SetValue(iRow, iCol, panColData[iRow]);
235 : }
236 5 : CPLFree(panColData);
237 : }
238 27 : if (aoFields[iCol].eType == GFT_Real)
239 : {
240 : double *padfColData = static_cast<double *>(
241 17 : VSI_MALLOC2_VERBOSE(sizeof(double), nRows));
242 17 : if (padfColData == nullptr)
243 : {
244 0 : delete poRAT;
245 0 : return nullptr;
246 : }
247 :
248 34 : if (const_cast<HFARasterAttributeTable *>(this)->ValuesIO(
249 17 : GF_Read, iCol, 0, nRows, padfColData) != CE_None)
250 : {
251 0 : CPLFree(padfColData);
252 0 : delete poRAT;
253 0 : return nullptr;
254 : }
255 :
256 3509 : for (int iRow = 0; iRow < nRows; iRow++)
257 : {
258 3492 : poRAT->SetValue(iRow, iCol, padfColData[iRow]);
259 : }
260 17 : CPLFree(padfColData);
261 : }
262 27 : if (aoFields[iCol].eType == GFT_String)
263 : {
264 : char **papszColData = static_cast<char **>(
265 5 : VSI_MALLOC2_VERBOSE(sizeof(char *), nRows));
266 5 : if (papszColData == nullptr)
267 : {
268 0 : delete poRAT;
269 0 : return nullptr;
270 : }
271 :
272 10 : if (const_cast<HFARasterAttributeTable *>(this)->ValuesIO(
273 5 : GF_Read, iCol, 0, nRows, papszColData) != CE_None)
274 : {
275 0 : CPLFree(papszColData);
276 0 : delete poRAT;
277 0 : return nullptr;
278 : }
279 :
280 501 : for (int iRow = 0; iRow < nRows; iRow++)
281 : {
282 496 : poRAT->SetValue(iRow, iCol, papszColData[iRow]);
283 496 : CPLFree(papszColData[iRow]);
284 : }
285 5 : CPLFree(papszColData);
286 : }
287 : }
288 :
289 17 : if (bLinearBinning)
290 11 : poRAT->SetLinearBinning(dfRow0Min, dfBinSize);
291 :
292 17 : poRAT->SetTableType(this->GetTableType());
293 :
294 17 : return poRAT;
295 : }
296 :
297 : /************************************************************************/
298 : /* GetColumnCount() */
299 : /************************************************************************/
300 :
301 39 : int HFARasterAttributeTable::GetColumnCount() const
302 : {
303 39 : return static_cast<int>(aoFields.size());
304 : }
305 :
306 : /************************************************************************/
307 : /* GetNameOfCol() */
308 : /************************************************************************/
309 :
310 9 : const char *HFARasterAttributeTable::GetNameOfCol(int nCol) const
311 : {
312 9 : if (nCol < 0 || nCol >= static_cast<int>(aoFields.size()))
313 0 : return nullptr;
314 :
315 9 : return aoFields[nCol].sName;
316 : }
317 :
318 : /************************************************************************/
319 : /* GetUsageOfCol() */
320 : /************************************************************************/
321 :
322 95 : GDALRATFieldUsage HFARasterAttributeTable::GetUsageOfCol(int nCol) const
323 : {
324 95 : if (nCol < 0 || nCol >= static_cast<int>(aoFields.size()))
325 0 : return GFU_Generic;
326 :
327 95 : return aoFields[nCol].eUsage;
328 : }
329 :
330 : /************************************************************************/
331 : /* GetTypeOfCol() */
332 : /************************************************************************/
333 :
334 3741 : GDALRATFieldType HFARasterAttributeTable::GetTypeOfCol(int nCol) const
335 : {
336 3741 : if (nCol < 0 || nCol >= static_cast<int>(aoFields.size()))
337 0 : return GFT_Integer;
338 :
339 3741 : return aoFields[nCol].eType;
340 : }
341 :
342 : /************************************************************************/
343 : /* GetColOfUsage() */
344 : /************************************************************************/
345 :
346 1 : int HFARasterAttributeTable::GetColOfUsage(GDALRATFieldUsage eUsage) const
347 : {
348 1 : for (unsigned int i = 0; i < aoFields.size(); i++)
349 : {
350 1 : if (aoFields[i].eUsage == eUsage)
351 1 : return i;
352 : }
353 :
354 0 : return -1;
355 : }
356 :
357 : /************************************************************************/
358 : /* GetRowCount() */
359 : /************************************************************************/
360 :
361 39 : int HFARasterAttributeTable::GetRowCount() const
362 : {
363 39 : return nRows;
364 : }
365 :
366 : /************************************************************************/
367 : /* GetValueAsString() */
368 : /************************************************************************/
369 :
370 10 : const char *HFARasterAttributeTable::GetValueAsString(int iRow,
371 : int iField) const
372 : {
373 : // Let ValuesIO do the work.
374 10 : char *apszStrList[1] = {nullptr};
375 10 : if (const_cast<HFARasterAttributeTable *>(this)->ValuesIO(
376 10 : GF_Read, iField, iRow, 1, apszStrList) != CE_None)
377 : {
378 0 : return "";
379 : }
380 :
381 : const_cast<HFARasterAttributeTable *>(this)->osWorkingResult =
382 10 : apszStrList[0];
383 10 : CPLFree(apszStrList[0]);
384 :
385 10 : return osWorkingResult;
386 : }
387 :
388 : /************************************************************************/
389 : /* GetValueAsInt() */
390 : /************************************************************************/
391 :
392 3125 : int HFARasterAttributeTable::GetValueAsInt(int iRow, int iField) const
393 : {
394 : // Let ValuesIO do the work.
395 3125 : int nValue = 0;
396 3125 : if (const_cast<HFARasterAttributeTable *>(this)->ValuesIO(
397 3125 : GF_Read, iField, iRow, 1, &nValue) != CE_None)
398 : {
399 0 : return 0;
400 : }
401 :
402 3125 : return nValue;
403 : }
404 :
405 : /************************************************************************/
406 : /* GetValueAsDouble() */
407 : /************************************************************************/
408 :
409 1592 : double HFARasterAttributeTable::GetValueAsDouble(int iRow, int iField) const
410 : {
411 : // Let ValuesIO do the work.
412 1592 : double dfValue = 0.0;
413 1592 : if (const_cast<HFARasterAttributeTable *>(this)->ValuesIO(
414 1592 : GF_Read, iField, iRow, 1, &dfValue) != CE_None)
415 : {
416 0 : return 0.0;
417 : }
418 :
419 1592 : return dfValue;
420 : }
421 :
422 : /************************************************************************/
423 : /* SetValue() */
424 : /************************************************************************/
425 :
426 20 : CPLErr HFARasterAttributeTable::SetValue(int iRow, int iField,
427 : const char *pszValue)
428 : {
429 : // Let ValuesIO do the work.
430 20 : char *apszValues[1] = {const_cast<char *>(pszValue)};
431 40 : return ValuesIO(GF_Write, iField, iRow, 1, apszValues);
432 : }
433 :
434 : /************************************************************************/
435 : /* SetValue() */
436 : /************************************************************************/
437 :
438 20 : CPLErr HFARasterAttributeTable::SetValue(int iRow, int iField, double dfValue)
439 : {
440 : // Let ValuesIO do the work.
441 20 : return ValuesIO(GF_Write, iField, iRow, 1, &dfValue);
442 : }
443 :
444 : /************************************************************************/
445 : /* SetValue() */
446 : /************************************************************************/
447 :
448 20 : CPLErr HFARasterAttributeTable::SetValue(int iRow, int iField, int nValue)
449 : {
450 : // Let ValuesIO do the work.
451 20 : return ValuesIO(GF_Write, iField, iRow, 1, &nValue);
452 : }
453 :
454 : /************************************************************************/
455 : /* ValuesIO() */
456 : /************************************************************************/
457 :
458 1651 : CPLErr HFARasterAttributeTable::ValuesIO(GDALRWFlag eRWFlag, int iField,
459 : int iStartRow, int iLength,
460 : double *pdfData)
461 : {
462 1651 : if (eRWFlag == GF_Write && eAccess == GA_ReadOnly)
463 : {
464 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
465 : "Dataset not open in update mode");
466 0 : return CE_Failure;
467 : }
468 :
469 1651 : if (iField < 0 || iField >= static_cast<int>(aoFields.size()))
470 : {
471 0 : CPLError(CE_Failure, CPLE_AppDefined, "iField (%d) out of range.",
472 : iField);
473 :
474 0 : return CE_Failure;
475 : }
476 :
477 1651 : if (iStartRow < 0 || iLength >= INT_MAX - iStartRow ||
478 1651 : (iStartRow + iLength) > nRows)
479 : {
480 0 : CPLError(CE_Failure, CPLE_AppDefined,
481 : "iStartRow (%d) + iLength(%d) out of range.", iStartRow,
482 : iLength);
483 :
484 0 : return CE_Failure;
485 : }
486 :
487 1651 : if (aoFields[iField].bConvertColors)
488 : {
489 : // Convert to/from float color field.
490 : int *panColData =
491 0 : static_cast<int *>(VSI_MALLOC2_VERBOSE(iLength, sizeof(int)));
492 0 : if (panColData == nullptr)
493 : {
494 0 : CPLFree(panColData);
495 0 : return CE_Failure;
496 : }
497 :
498 0 : if (eRWFlag == GF_Write)
499 : {
500 0 : for (int i = 0; i < iLength; i++)
501 0 : panColData[i] = static_cast<int>(pdfData[i]);
502 : }
503 :
504 : const CPLErr ret =
505 0 : ColorsIO(eRWFlag, iField, iStartRow, iLength, panColData);
506 :
507 0 : if (eRWFlag == GF_Read)
508 : {
509 : // Copy them back to doubles.
510 0 : for (int i = 0; i < iLength; i++)
511 0 : pdfData[i] = panColData[i];
512 : }
513 :
514 0 : CPLFree(panColData);
515 0 : return ret;
516 : }
517 :
518 1651 : switch (aoFields[iField].eType)
519 : {
520 1 : case GFT_Integer:
521 : {
522 : // Allocate space for ints.
523 : int *panColData =
524 1 : static_cast<int *>(VSI_MALLOC2_VERBOSE(iLength, sizeof(int)));
525 1 : if (panColData == nullptr)
526 : {
527 0 : CPLFree(panColData);
528 0 : return CE_Failure;
529 : }
530 :
531 1 : if (eRWFlag == GF_Write)
532 : {
533 : // Copy the application supplied doubles to ints.
534 11 : for (int i = 0; i < iLength; i++)
535 10 : panColData[i] = static_cast<int>(pdfData[i]);
536 : }
537 :
538 : // Do the ValuesIO as ints.
539 : const CPLErr eVal =
540 1 : ValuesIO(eRWFlag, iField, iStartRow, iLength, panColData);
541 1 : if (eVal != CE_None)
542 : {
543 0 : CPLFree(panColData);
544 0 : return eVal;
545 : }
546 :
547 1 : if (eRWFlag == GF_Read)
548 : {
549 : // Copy them back to doubles.
550 0 : for (int i = 0; i < iLength; i++)
551 0 : pdfData[i] = panColData[i];
552 : }
553 :
554 1 : CPLFree(panColData);
555 : }
556 1 : break;
557 1649 : case GFT_Real:
558 : {
559 1649 : if ((eRWFlag == GF_Read) && aoFields[iField].bIsBinValues)
560 : {
561 : // Probably could change HFAReadBFUniqueBins to only read needed
562 : // rows.
563 193 : double *padfBinValues = HFAReadBFUniqueBins(
564 193 : aoFields[iField].poColumn, iStartRow + iLength);
565 193 : if (padfBinValues == nullptr)
566 0 : return CE_Failure;
567 193 : memcpy(pdfData, &padfBinValues[iStartRow],
568 193 : sizeof(double) * iLength);
569 193 : CPLFree(padfBinValues);
570 : }
571 : else
572 : {
573 2912 : if (VSIFSeekL(hHFA->fp,
574 1456 : aoFields[iField].nDataOffset +
575 2912 : (static_cast<vsi_l_offset>(iStartRow) *
576 1456 : aoFields[iField].nElementSize),
577 1456 : SEEK_SET) != 0)
578 : {
579 0 : return CE_Failure;
580 : }
581 :
582 1456 : if (eRWFlag == GF_Read)
583 : {
584 2864 : if (static_cast<int>(VSIFReadL(pdfData, sizeof(double),
585 1432 : iLength, hHFA->fp)) !=
586 : iLength)
587 : {
588 0 : CPLError(CE_Failure, CPLE_AppDefined,
589 : "HFARasterAttributeTable::ValuesIO: "
590 : "Cannot read values");
591 0 : return CE_Failure;
592 : }
593 : #ifdef CPL_MSB
594 : GDALSwapWords(pdfData, 8, iLength, 8);
595 : #endif
596 : }
597 : else
598 : {
599 : #ifdef CPL_MSB
600 : GDALSwapWords(pdfData, 8, iLength, 8);
601 : #endif
602 : // Note: HFAAllocateSpace now called by CreateColumn so
603 : // space should exist.
604 48 : if (static_cast<int>(VSIFWriteL(pdfData, sizeof(double),
605 24 : iLength, hHFA->fp)) !=
606 : iLength)
607 : {
608 0 : CPLError(CE_Failure, CPLE_AppDefined,
609 : "HFARasterAttributeTable::ValuesIO: "
610 : "Cannot write values");
611 0 : return CE_Failure;
612 : }
613 : #ifdef CPL_MSB
614 : // Swap back.
615 : GDALSwapWords(pdfData, 8, iLength, 8);
616 : #endif
617 : }
618 : }
619 : }
620 1649 : break;
621 1 : case GFT_String:
622 : {
623 : // Allocate space for string pointers.
624 : char **papszColData = static_cast<char **>(
625 1 : VSI_MALLOC2_VERBOSE(iLength, sizeof(char *)));
626 1 : if (papszColData == nullptr)
627 : {
628 0 : return CE_Failure;
629 : }
630 :
631 1 : if (eRWFlag == GF_Write)
632 : {
633 : // Copy the application supplied doubles to strings.
634 11 : for (int i = 0; i < iLength; i++)
635 : {
636 10 : osWorkingResult.Printf("%.16g", pdfData[i]);
637 10 : papszColData[i] = CPLStrdup(osWorkingResult);
638 : }
639 : }
640 :
641 : // Do the ValuesIO as strings.
642 : const CPLErr eVal =
643 1 : ValuesIO(eRWFlag, iField, iStartRow, iLength, papszColData);
644 1 : if (eVal != CE_None)
645 : {
646 0 : if (eRWFlag == GF_Write)
647 : {
648 0 : for (int i = 0; i < iLength; i++)
649 0 : CPLFree(papszColData[i]);
650 : }
651 0 : CPLFree(papszColData);
652 0 : return eVal;
653 : }
654 :
655 1 : if (eRWFlag == GF_Read)
656 : {
657 : // Copy them back to doubles.
658 0 : for (int i = 0; i < iLength; i++)
659 0 : pdfData[i] = CPLAtof(papszColData[i]);
660 : }
661 :
662 : // Either we allocated them for write, or they were allocated
663 : // by ValuesIO on read.
664 11 : for (int i = 0; i < iLength; i++)
665 10 : CPLFree(papszColData[i]);
666 :
667 1 : CPLFree(papszColData);
668 : }
669 1 : break;
670 : }
671 :
672 1651 : return CE_None;
673 : }
674 :
675 : /************************************************************************/
676 : /* ValuesIO() */
677 : /************************************************************************/
678 :
679 3169 : CPLErr HFARasterAttributeTable::ValuesIO(GDALRWFlag eRWFlag, int iField,
680 : int iStartRow, int iLength,
681 : int *pnData)
682 : {
683 3169 : if (eRWFlag == GF_Write && eAccess == GA_ReadOnly)
684 : {
685 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
686 : "Dataset not open in update mode");
687 0 : return CE_Failure;
688 : }
689 :
690 3169 : if (iField < 0 || iField >= static_cast<int>(aoFields.size()))
691 : {
692 0 : CPLError(CE_Failure, CPLE_AppDefined, "iField (%d) out of range.",
693 : iField);
694 :
695 0 : return CE_Failure;
696 : }
697 :
698 3169 : if (iStartRow < 0 || iLength >= INT_MAX - iStartRow ||
699 3169 : (iStartRow + iLength) > nRows)
700 : {
701 0 : CPLError(CE_Failure, CPLE_AppDefined,
702 : "iStartRow (%d) + iLength(%d) out of range.", iStartRow,
703 : iLength);
704 :
705 0 : return CE_Failure;
706 : }
707 :
708 3169 : if (aoFields[iField].bConvertColors)
709 : {
710 : // Convert to/from float color field.
711 3112 : return ColorsIO(eRWFlag, iField, iStartRow, iLength, pnData);
712 : }
713 :
714 57 : switch (aoFields[iField].eType)
715 : {
716 52 : case GFT_Integer:
717 : {
718 104 : if (VSIFSeekL(hHFA->fp,
719 52 : aoFields[iField].nDataOffset +
720 104 : (static_cast<vsi_l_offset>(iStartRow) *
721 52 : aoFields[iField].nElementSize),
722 52 : SEEK_SET) != 0)
723 : {
724 0 : return CE_Failure;
725 : }
726 : GInt32 *panColData = static_cast<GInt32 *>(
727 52 : VSI_MALLOC2_VERBOSE(iLength, sizeof(GInt32)));
728 52 : if (panColData == nullptr)
729 : {
730 0 : return CE_Failure;
731 : }
732 :
733 52 : if (eRWFlag == GF_Read)
734 : {
735 56 : if (static_cast<int>(VSIFReadL(panColData, sizeof(GInt32),
736 28 : iLength, hHFA->fp)) != iLength)
737 : {
738 0 : CPLError(CE_Failure, CPLE_AppDefined,
739 : "HFARasterAttributeTable::ValuesIO: "
740 : "Cannot read values");
741 0 : CPLFree(panColData);
742 0 : return CE_Failure;
743 : }
744 : #ifdef CPL_MSB
745 : GDALSwapWords(panColData, 4, iLength, 4);
746 : #endif
747 : // Now copy into application buffer. This extra step
748 : // may not be necessary if sizeof(int) == sizeof(GInt32).
749 664 : for (int i = 0; i < iLength; i++)
750 636 : pnData[i] = panColData[i];
751 : }
752 : else
753 : {
754 : // Copy from application buffer.
755 84 : for (int i = 0; i < iLength; i++)
756 60 : panColData[i] = pnData[i];
757 :
758 : #ifdef CPL_MSB
759 : GDALSwapWords(panColData, 4, iLength, 4);
760 : #endif
761 : // Note: HFAAllocateSpace now called by CreateColumn so space
762 : // should exist.
763 48 : if (static_cast<int>(VSIFWriteL(panColData, sizeof(GInt32),
764 24 : iLength, hHFA->fp)) != iLength)
765 : {
766 0 : CPLError(CE_Failure, CPLE_AppDefined,
767 : "HFARasterAttributeTable::ValuesIO: "
768 : "Cannot write values");
769 0 : CPLFree(panColData);
770 0 : return CE_Failure;
771 : }
772 : }
773 52 : CPLFree(panColData);
774 : }
775 52 : break;
776 4 : case GFT_Real:
777 : {
778 : // Allocate space for doubles.
779 : double *padfColData = static_cast<double *>(
780 4 : VSI_MALLOC2_VERBOSE(iLength, sizeof(double)));
781 4 : if (padfColData == nullptr)
782 : {
783 0 : return CE_Failure;
784 : }
785 :
786 4 : if (eRWFlag == GF_Write)
787 : {
788 : // Copy the application supplied ints to doubles.
789 11 : for (int i = 0; i < iLength; i++)
790 10 : padfColData[i] = pnData[i];
791 : }
792 :
793 : // Do the ValuesIO as doubles.
794 : const CPLErr eVal =
795 4 : ValuesIO(eRWFlag, iField, iStartRow, iLength, padfColData);
796 4 : if (eVal != CE_None)
797 : {
798 0 : CPLFree(padfColData);
799 0 : return eVal;
800 : }
801 :
802 4 : if (eRWFlag == GF_Read)
803 : {
804 : // Copy them back to ints.
805 6 : for (int i = 0; i < iLength; i++)
806 3 : pnData[i] = static_cast<int>(padfColData[i]);
807 : }
808 :
809 4 : CPLFree(padfColData);
810 : }
811 4 : break;
812 1 : case GFT_String:
813 : {
814 : // Allocate space for string pointers.
815 : char **papszColData = static_cast<char **>(
816 1 : VSI_MALLOC2_VERBOSE(iLength, sizeof(char *)));
817 1 : if (papszColData == nullptr)
818 : {
819 0 : return CE_Failure;
820 : }
821 :
822 1 : if (eRWFlag == GF_Write)
823 : {
824 : // Copy the application supplied ints to strings.
825 11 : for (int i = 0; i < iLength; i++)
826 : {
827 10 : osWorkingResult.Printf("%d", pnData[i]);
828 10 : papszColData[i] = CPLStrdup(osWorkingResult);
829 : }
830 : }
831 :
832 : // Do the ValuesIO as strings.
833 : const CPLErr eVal =
834 1 : ValuesIO(eRWFlag, iField, iStartRow, iLength, papszColData);
835 1 : if (eVal != CE_None)
836 : {
837 0 : if (eRWFlag == GF_Write)
838 : {
839 0 : for (int i = 0; i < iLength; i++)
840 0 : CPLFree(papszColData[i]);
841 : }
842 0 : CPLFree(papszColData);
843 0 : return eVal;
844 : }
845 :
846 1 : if (eRWFlag == GF_Read)
847 : {
848 : // Copy them back to ints.
849 0 : for (int i = 0; i < iLength; i++)
850 0 : pnData[i] = atoi(papszColData[i]);
851 : }
852 :
853 : // Either we allocated them for write, or they were allocated
854 : // by ValuesIO on read.
855 11 : for (int i = 0; i < iLength; i++)
856 10 : CPLFree(papszColData[i]);
857 :
858 1 : CPLFree(papszColData);
859 : }
860 1 : break;
861 : }
862 :
863 57 : return CE_None;
864 : }
865 :
866 : /************************************************************************/
867 : /* ValuesIO() */
868 : /************************************************************************/
869 :
870 56 : CPLErr HFARasterAttributeTable::ValuesIO(GDALRWFlag eRWFlag, int iField,
871 : int iStartRow, int iLength,
872 : char **papszStrList)
873 : {
874 56 : if (eRWFlag == GF_Write && eAccess == GA_ReadOnly)
875 : {
876 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
877 : "Dataset not open in update mode");
878 0 : return CE_Failure;
879 : }
880 :
881 56 : if (iField < 0 || iField >= static_cast<int>(aoFields.size()))
882 : {
883 0 : CPLError(CE_Failure, CPLE_AppDefined, "iField (%d) out of range.",
884 : iField);
885 :
886 0 : return CE_Failure;
887 : }
888 :
889 56 : if (iStartRow < 0 || iLength >= INT_MAX - iStartRow ||
890 56 : (iStartRow + iLength) > nRows)
891 : {
892 0 : CPLError(CE_Failure, CPLE_AppDefined,
893 : "iStartRow (%d) + iLength(%d) out of range.", iStartRow,
894 : iLength);
895 :
896 0 : return CE_Failure;
897 : }
898 :
899 56 : if (aoFields[iField].bConvertColors)
900 : {
901 : // Convert to/from float color field.
902 : int *panColData =
903 0 : static_cast<int *>(VSI_MALLOC2_VERBOSE(iLength, sizeof(int)));
904 0 : if (panColData == nullptr)
905 : {
906 0 : CPLFree(panColData);
907 0 : return CE_Failure;
908 : }
909 :
910 0 : if (eRWFlag == GF_Write)
911 : {
912 0 : for (int i = 0; i < iLength; i++)
913 0 : panColData[i] = atoi(papszStrList[i]);
914 : }
915 :
916 : const CPLErr ret =
917 0 : ColorsIO(eRWFlag, iField, iStartRow, iLength, panColData);
918 :
919 0 : if (eRWFlag == GF_Read)
920 : {
921 : // Copy them back to strings.
922 0 : for (int i = 0; i < iLength; i++)
923 : {
924 0 : osWorkingResult.Printf("%d", panColData[i]);
925 0 : papszStrList[i] = CPLStrdup(osWorkingResult);
926 : }
927 : }
928 :
929 0 : CPLFree(panColData);
930 0 : return ret;
931 : }
932 :
933 56 : switch (aoFields[iField].eType)
934 : {
935 1 : case GFT_Integer:
936 : {
937 : // Allocate space for ints.
938 : int *panColData =
939 1 : static_cast<int *>(VSI_MALLOC2_VERBOSE(iLength, sizeof(int)));
940 1 : if (panColData == nullptr)
941 : {
942 0 : return CE_Failure;
943 : }
944 :
945 1 : if (eRWFlag == GF_Write)
946 : {
947 : // Convert user supplied strings to ints.
948 11 : for (int i = 0; i < iLength; i++)
949 10 : panColData[i] = atoi(papszStrList[i]);
950 : }
951 :
952 : // Call values IO to read/write ints.
953 : const CPLErr eVal =
954 1 : ValuesIO(eRWFlag, iField, iStartRow, iLength, panColData);
955 1 : if (eVal != CE_None)
956 : {
957 0 : CPLFree(panColData);
958 0 : return eVal;
959 : }
960 :
961 1 : if (eRWFlag == GF_Read)
962 : {
963 : // Convert ints back to strings.
964 0 : for (int i = 0; i < iLength; i++)
965 : {
966 0 : osWorkingResult.Printf("%d", panColData[i]);
967 0 : papszStrList[i] = CPLStrdup(osWorkingResult);
968 : }
969 : }
970 1 : CPLFree(panColData);
971 : }
972 1 : break;
973 1 : case GFT_Real:
974 : {
975 : // Allocate space for doubles.
976 : double *padfColData = static_cast<double *>(
977 1 : VSI_MALLOC2_VERBOSE(iLength, sizeof(double)));
978 1 : if (padfColData == nullptr)
979 : {
980 0 : return CE_Failure;
981 : }
982 :
983 1 : if (eRWFlag == GF_Write)
984 : {
985 : // Convert user supplied strings to doubles.
986 11 : for (int i = 0; i < iLength; i++)
987 10 : padfColData[i] = CPLAtof(papszStrList[i]);
988 : }
989 :
990 : // Call value IO to read/write doubles.
991 : const CPLErr eVal =
992 1 : ValuesIO(eRWFlag, iField, iStartRow, iLength, padfColData);
993 1 : if (eVal != CE_None)
994 : {
995 0 : CPLFree(padfColData);
996 0 : return eVal;
997 : }
998 :
999 1 : if (eRWFlag == GF_Read)
1000 : {
1001 : // Convert doubles back to strings.
1002 0 : for (int i = 0; i < iLength; i++)
1003 : {
1004 0 : osWorkingResult.Printf("%.16g", padfColData[i]);
1005 0 : papszStrList[i] = CPLStrdup(osWorkingResult);
1006 : }
1007 : }
1008 1 : CPLFree(padfColData);
1009 : }
1010 1 : break;
1011 54 : case GFT_String:
1012 : {
1013 108 : if (VSIFSeekL(hHFA->fp,
1014 54 : aoFields[iField].nDataOffset +
1015 108 : (static_cast<vsi_l_offset>(iStartRow) *
1016 54 : aoFields[iField].nElementSize),
1017 54 : SEEK_SET) != 0)
1018 : {
1019 0 : return CE_Failure;
1020 : }
1021 : char *pachColData = static_cast<char *>(
1022 54 : VSI_MALLOC2_VERBOSE(iLength, aoFields[iField].nElementSize));
1023 54 : if (pachColData == nullptr)
1024 : {
1025 0 : return CE_Failure;
1026 : }
1027 :
1028 54 : if (eRWFlag == GF_Read)
1029 : {
1030 87 : if (static_cast<int>(VSIFReadL(pachColData,
1031 29 : aoFields[iField].nElementSize,
1032 58 : iLength, hHFA->fp)) != iLength)
1033 : {
1034 0 : CPLError(CE_Failure, CPLE_AppDefined,
1035 : "HFARasterAttributeTable::ValuesIO: "
1036 : "Cannot read values");
1037 0 : CPLFree(pachColData);
1038 0 : return CE_Failure;
1039 : }
1040 :
1041 : // Now copy into application buffer.
1042 675 : for (int i = 0; i < iLength; i++)
1043 : {
1044 : osWorkingResult.assign(
1045 1292 : pachColData + aoFields[iField].nElementSize * i,
1046 646 : aoFields[iField].nElementSize);
1047 646 : papszStrList[i] = CPLStrdup(osWorkingResult);
1048 : }
1049 : }
1050 : else
1051 : {
1052 : // We need to check that these strings will fit in the allocated
1053 : // space.
1054 25 : int nNewMaxChars = aoFields[iField].nElementSize;
1055 95 : for (int i = 0; i < iLength; i++)
1056 : {
1057 70 : const int nStringSize =
1058 70 : static_cast<int>(strlen(papszStrList[i])) + 1;
1059 70 : if (nStringSize > nNewMaxChars)
1060 3 : nNewMaxChars = nStringSize;
1061 : }
1062 :
1063 25 : if (nNewMaxChars > aoFields[iField].nElementSize)
1064 : {
1065 : // OK we have a problem: The allocated space is not big
1066 : // enough we need to re-allocate the space and update the
1067 : // pointers and copy across the old data.
1068 : const int nNewOffset =
1069 4 : HFAAllocateSpace(hHFA->papoBand[nBand - 1]->psInfo,
1070 2 : nRows * nNewMaxChars);
1071 4 : char *pszBuffer = static_cast<char *>(VSIMalloc2(
1072 2 : aoFields[iField].nElementSize, sizeof(char)));
1073 32 : for (int i = 0; i < nRows; i++)
1074 : {
1075 : // Seek to the old place.
1076 30 : CPL_IGNORE_RET_VAL(
1077 30 : VSIFSeekL(hHFA->fp,
1078 30 : aoFields[iField].nDataOffset +
1079 60 : (static_cast<vsi_l_offset>(i) *
1080 30 : aoFields[iField].nElementSize),
1081 : SEEK_SET));
1082 : // Read in old data.
1083 30 : CPL_IGNORE_RET_VAL(
1084 30 : VSIFReadL(pszBuffer, aoFields[iField].nElementSize,
1085 30 : 1, hHFA->fp));
1086 : // Seek to new place.
1087 60 : bool bOK = VSIFSeekL(hHFA->fp,
1088 30 : nNewOffset +
1089 30 : (static_cast<vsi_l_offset>(i) *
1090 30 : nNewMaxChars),
1091 30 : SEEK_SET) == 0;
1092 : // Write data to new place.
1093 30 : bOK &=
1094 30 : VSIFWriteL(pszBuffer, aoFields[iField].nElementSize,
1095 30 : 1, hHFA->fp) == 1;
1096 : // Make sure there is a terminating null byte just to be
1097 : // safe.
1098 30 : const char cNullByte = '\0';
1099 60 : bOK &= VSIFWriteL(&cNullByte, sizeof(char), 1,
1100 30 : hHFA->fp) == 1;
1101 30 : if (!bOK)
1102 : {
1103 0 : CPLFree(pszBuffer);
1104 0 : CPLFree(pachColData);
1105 0 : CPLError(CE_Failure, CPLE_AppDefined,
1106 : "HFARasterAttributeTable::ValuesIO: "
1107 : "Cannot write values");
1108 0 : return CE_Failure;
1109 : }
1110 : }
1111 : // Update our data structures.
1112 2 : aoFields[iField].nElementSize = nNewMaxChars;
1113 2 : aoFields[iField].nDataOffset = nNewOffset;
1114 : // Update file.
1115 2 : aoFields[iField].poColumn->SetIntField("columnDataPtr",
1116 : nNewOffset);
1117 2 : aoFields[iField].poColumn->SetIntField("maxNumChars",
1118 : nNewMaxChars);
1119 :
1120 : // Note: There isn't an HFAFreeSpace so we can't un-allocate
1121 : // the old space in the file.
1122 2 : CPLFree(pszBuffer);
1123 :
1124 : // Re-allocate our buffer.
1125 2 : CPLFree(pachColData);
1126 : pachColData = static_cast<char *>(
1127 2 : VSI_MALLOC2_VERBOSE(iLength, nNewMaxChars));
1128 2 : if (pachColData == nullptr)
1129 : {
1130 0 : return CE_Failure;
1131 : }
1132 :
1133 : // Lastly seek to the right place in the new space ready to
1134 : // write.
1135 4 : if (VSIFSeekL(hHFA->fp,
1136 2 : nNewOffset +
1137 2 : (static_cast<vsi_l_offset>(iStartRow) *
1138 2 : nNewMaxChars),
1139 2 : SEEK_SET) != 0)
1140 : {
1141 0 : VSIFree(pachColData);
1142 0 : return CE_Failure;
1143 : }
1144 : }
1145 :
1146 : // Copy from application buffer.
1147 95 : for (int i = 0; i < iLength; i++)
1148 70 : strcpy(&pachColData[nNewMaxChars * i], papszStrList[i]);
1149 :
1150 : // Note: HFAAllocateSpace now called by CreateColumn so space
1151 : // should exist.
1152 75 : if (static_cast<int>(VSIFWriteL(pachColData,
1153 25 : aoFields[iField].nElementSize,
1154 50 : iLength, hHFA->fp)) != iLength)
1155 : {
1156 0 : CPLError(CE_Failure, CPLE_AppDefined,
1157 : "HFARasterAttributeTable::ValuesIO: "
1158 : "Cannot write values");
1159 0 : CPLFree(pachColData);
1160 0 : return CE_Failure;
1161 : }
1162 : }
1163 54 : CPLFree(pachColData);
1164 : }
1165 54 : break;
1166 : }
1167 :
1168 56 : return CE_None;
1169 : }
1170 :
1171 : /************************************************************************/
1172 : /* ColorsIO() */
1173 : /************************************************************************/
1174 :
1175 : // Handle the fact that HFA stores colours as floats, but we need to
1176 : // read them in as ints 0...255.
1177 3112 : CPLErr HFARasterAttributeTable::ColorsIO(GDALRWFlag eRWFlag, int iField,
1178 : int iStartRow, int iLength,
1179 : int *pnData)
1180 : {
1181 : // Allocate space for doubles.
1182 : double *padfData =
1183 3112 : static_cast<double *>(VSI_MALLOC2_VERBOSE(iLength, sizeof(double)));
1184 3112 : if (padfData == nullptr)
1185 : {
1186 0 : return CE_Failure;
1187 : }
1188 :
1189 3112 : if (eRWFlag == GF_Write)
1190 : {
1191 : // Copy the application supplied ints to doubles
1192 : // and convert 0..255 to 0..1 in the same manner
1193 : // as the color table.
1194 0 : for (int i = 0; i < iLength; i++)
1195 0 : padfData[i] = pnData[i] / 255.0;
1196 : }
1197 :
1198 6224 : if (VSIFSeekL(hHFA->fp,
1199 3112 : aoFields[iField].nDataOffset +
1200 6224 : (static_cast<vsi_l_offset>(iStartRow) *
1201 3112 : aoFields[iField].nElementSize),
1202 3112 : SEEK_SET) != 0)
1203 : {
1204 0 : CPLFree(padfData);
1205 0 : return CE_Failure;
1206 : }
1207 :
1208 3112 : if (eRWFlag == GF_Read)
1209 : {
1210 6224 : if (static_cast<int>(VSIFReadL(padfData, sizeof(double), iLength,
1211 3112 : hHFA->fp)) != iLength)
1212 : {
1213 0 : CPLError(CE_Failure, CPLE_AppDefined,
1214 : "HFARasterAttributeTable::ColorsIO: Cannot read values");
1215 0 : CPLFree(padfData);
1216 0 : return CE_Failure;
1217 : }
1218 : #ifdef CPL_MSB
1219 : GDALSwapWords(padfData, 8, iLength, 8);
1220 : #endif
1221 : }
1222 : else
1223 : {
1224 : #ifdef CPL_MSB
1225 : GDALSwapWords(padfData, 8, iLength, 8);
1226 : #endif
1227 : // Note: HFAAllocateSpace now called by CreateColumn so space should
1228 : // exist.
1229 0 : if (static_cast<int>(VSIFWriteL(padfData, sizeof(double), iLength,
1230 0 : hHFA->fp)) != iLength)
1231 : {
1232 0 : CPLError(CE_Failure, CPLE_AppDefined,
1233 : "HFARasterAttributeTable::ColorsIO: Cannot write values");
1234 0 : CPLFree(padfData);
1235 0 : return CE_Failure;
1236 : }
1237 : }
1238 :
1239 3112 : if (eRWFlag == GF_Read)
1240 : {
1241 : // Copy them back to ints converting 0..1 to 0..255 in
1242 : // the same manner as the color table.
1243 : // TODO(schwehr): Symbolic constants for 255 and 256.
1244 6224 : for (int i = 0; i < iLength; i++)
1245 3112 : pnData[i] = std::min(255, static_cast<int>(padfData[i] * 256));
1246 : }
1247 :
1248 3112 : CPLFree(padfData);
1249 :
1250 3112 : return CE_None;
1251 : }
1252 :
1253 : /************************************************************************/
1254 : /* ChangesAreWrittenToFile() */
1255 : /************************************************************************/
1256 :
1257 1 : int HFARasterAttributeTable::ChangesAreWrittenToFile()
1258 : {
1259 1 : return TRUE;
1260 : }
1261 :
1262 : /************************************************************************/
1263 : /* SetRowCount() */
1264 : /************************************************************************/
1265 :
1266 2 : void HFARasterAttributeTable::SetRowCount(int iCount)
1267 : {
1268 2 : if (eAccess == GA_ReadOnly)
1269 : {
1270 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
1271 : "Dataset not open in update mode");
1272 0 : return;
1273 : }
1274 :
1275 2 : if (iCount > nRows)
1276 : {
1277 : // Making the RAT larger - a bit hard.
1278 : // We need to re-allocate space on disc.
1279 20 : for (int iCol = 0; iCol < static_cast<int>(aoFields.size()); iCol++)
1280 : {
1281 : // New space.
1282 : const int nNewOffset =
1283 36 : HFAAllocateSpace(hHFA->papoBand[nBand - 1]->psInfo,
1284 18 : iCount * aoFields[iCol].nElementSize);
1285 :
1286 : // Only need to bother if there are actually rows.
1287 18 : if (nRows > 0)
1288 : {
1289 : // Temp buffer for this column.
1290 : void *pData =
1291 9 : VSI_MALLOC2_VERBOSE(nRows, aoFields[iCol].nElementSize);
1292 9 : if (pData == nullptr)
1293 : {
1294 0 : return;
1295 : }
1296 : // Read old data.
1297 9 : if (VSIFSeekL(hHFA->fp, aoFields[iCol].nDataOffset, SEEK_SET) !=
1298 18 : 0 ||
1299 27 : static_cast<int>(VSIFReadL(pData,
1300 9 : aoFields[iCol].nElementSize,
1301 18 : nRows, hHFA->fp)) != nRows)
1302 : {
1303 0 : CPLError(CE_Failure, CPLE_AppDefined,
1304 : "HFARasterAttributeTable::SetRowCount: "
1305 : "Cannot read values");
1306 0 : CPLFree(pData);
1307 0 : return;
1308 : }
1309 :
1310 : // Write data - new space will be uninitialised.
1311 18 : if (VSIFSeekL(hHFA->fp, nNewOffset, SEEK_SET) != 0 ||
1312 27 : static_cast<int>(VSIFWriteL(pData,
1313 9 : aoFields[iCol].nElementSize,
1314 18 : nRows, hHFA->fp)) != nRows)
1315 : {
1316 0 : CPLError(CE_Failure, CPLE_AppDefined,
1317 : "HFARasterAttributeTable::SetRowCount: "
1318 : "Cannot write values");
1319 0 : CPLFree(pData);
1320 0 : return;
1321 : }
1322 9 : CPLFree(pData);
1323 : }
1324 :
1325 : // Update our data structures.
1326 18 : aoFields[iCol].nDataOffset = nNewOffset;
1327 : // Update file.
1328 18 : aoFields[iCol].poColumn->SetIntField("columnDataPtr", nNewOffset);
1329 18 : aoFields[iCol].poColumn->SetIntField("numRows", iCount);
1330 : }
1331 : }
1332 0 : else if (iCount < nRows)
1333 : {
1334 : // Update the numRows.
1335 0 : for (int iCol = 0; iCol < static_cast<int>(aoFields.size()); iCol++)
1336 : {
1337 0 : aoFields[iCol].poColumn->SetIntField("numRows", iCount);
1338 : }
1339 : }
1340 :
1341 2 : nRows = iCount;
1342 :
1343 2 : if (poDT != nullptr && EQUAL(poDT->GetType(), "Edsc_Table"))
1344 : {
1345 2 : poDT->SetIntField("numrows", iCount);
1346 : }
1347 : }
1348 :
1349 : /************************************************************************/
1350 : /* GetRowOfValue() */
1351 : /************************************************************************/
1352 1 : int HFARasterAttributeTable::GetRowOfValue(double dfValue) const
1353 : {
1354 : // Handle case of regular binning.
1355 1 : if (bLinearBinning)
1356 : {
1357 1 : const int iBin =
1358 1 : static_cast<int>(floor((dfValue - dfRow0Min) / dfBinSize));
1359 1 : if (iBin < 0 || iBin >= nRows)
1360 0 : return -1;
1361 1 : return iBin;
1362 : }
1363 : // Do we have any information?
1364 0 : int nMinCol = GetColOfUsage(GFU_Min);
1365 0 : if (nMinCol == -1)
1366 0 : nMinCol = GetColOfUsage(GFU_MinMax);
1367 0 : int nMaxCol = GetColOfUsage(GFU_Max);
1368 0 : if (nMaxCol == -1)
1369 0 : nMaxCol = GetColOfUsage(GFU_MinMax);
1370 0 : if (nMinCol == -1 && nMaxCol == -1)
1371 0 : return -1;
1372 : // Search through rows for match.
1373 0 : for (int iRow = 0; iRow < nRows; iRow++)
1374 : {
1375 0 : if (nMinCol != -1)
1376 : {
1377 0 : while (iRow < nRows && dfValue < GetValueAsDouble(iRow, nMinCol))
1378 0 : iRow++;
1379 0 : if (iRow == nRows)
1380 0 : break;
1381 : }
1382 0 : if (nMaxCol != -1)
1383 : {
1384 0 : if (dfValue > GetValueAsDouble(iRow, nMaxCol))
1385 0 : continue;
1386 : }
1387 0 : return iRow;
1388 : }
1389 0 : return -1;
1390 : }
1391 :
1392 : /************************************************************************/
1393 : /* GetRowOfValue() */
1394 : /* */
1395 : /* Int arg for now just converted to double. Perhaps we will */
1396 : /* handle this in a special way some day? */
1397 : /************************************************************************/
1398 0 : int HFARasterAttributeTable::GetRowOfValue(int nValue) const
1399 : {
1400 0 : return GetRowOfValue(static_cast<double>(nValue));
1401 : }
1402 :
1403 : /************************************************************************/
1404 : /* CreateColumn() */
1405 : /************************************************************************/
1406 :
1407 9 : CPLErr HFARasterAttributeTable::CreateColumn(const char *pszFieldName,
1408 : GDALRATFieldType eFieldType,
1409 : GDALRATFieldUsage eFieldUsage)
1410 : {
1411 9 : if (eAccess == GA_ReadOnly)
1412 : {
1413 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
1414 : "Dataset not open in update mode");
1415 0 : return CE_Failure;
1416 : }
1417 :
1418 : // Do we have a descriptor table already?
1419 9 : if (poDT == nullptr || !EQUAL(poDT->GetType(), "Edsc_Table"))
1420 1 : CreateDT();
1421 :
1422 9 : bool bConvertColors = false;
1423 :
1424 : // Imagine doesn't have a concept of usage - works of the names instead.
1425 : // Must make sure name matches use.
1426 9 : if (eFieldUsage == GFU_Red)
1427 : {
1428 0 : pszFieldName = "Red";
1429 : // Create a real column in the file, but make it
1430 : // available as int to GDAL.
1431 0 : bConvertColors = true;
1432 0 : eFieldType = GFT_Real;
1433 : }
1434 9 : else if (eFieldUsage == GFU_Green)
1435 : {
1436 0 : pszFieldName = "Green";
1437 0 : bConvertColors = true;
1438 0 : eFieldType = GFT_Real;
1439 : }
1440 9 : else if (eFieldUsage == GFU_Blue)
1441 : {
1442 0 : pszFieldName = "Blue";
1443 0 : bConvertColors = true;
1444 0 : eFieldType = GFT_Real;
1445 : }
1446 9 : else if (eFieldUsage == GFU_Alpha)
1447 : {
1448 0 : pszFieldName = "Opacity";
1449 0 : bConvertColors = true;
1450 0 : eFieldType = GFT_Real;
1451 : }
1452 9 : else if (eFieldUsage == GFU_PixelCount)
1453 : {
1454 0 : pszFieldName = "Histogram";
1455 : // Histogram is always float in HFA.
1456 0 : eFieldType = GFT_Real;
1457 : }
1458 9 : else if (eFieldUsage == GFU_Name)
1459 : {
1460 0 : pszFieldName = "Class_Names";
1461 : }
1462 :
1463 : // Check to see if a column with pszFieldName exists and create it
1464 : // if necessary.
1465 9 : HFAEntry *poColumn = poDT->GetNamedChild(pszFieldName);
1466 :
1467 9 : if (poColumn == nullptr || !EQUAL(poColumn->GetType(), "Edsc_Column"))
1468 9 : poColumn = HFAEntry::New(hHFA->papoBand[nBand - 1]->psInfo,
1469 : pszFieldName, "Edsc_Column", poDT);
1470 :
1471 9 : poColumn->SetIntField("numRows", nRows);
1472 9 : int nElementSize = 0;
1473 :
1474 9 : if (eFieldType == GFT_Integer)
1475 : {
1476 3 : nElementSize = sizeof(GInt32);
1477 3 : poColumn->SetStringField("dataType", "integer");
1478 : }
1479 6 : else if (eFieldType == GFT_Real)
1480 : {
1481 3 : nElementSize = sizeof(double);
1482 3 : poColumn->SetStringField("dataType", "real");
1483 : }
1484 3 : else if (eFieldType == GFT_String)
1485 : {
1486 : // Just have to guess here since we don't have any strings to check.
1487 3 : nElementSize = 10;
1488 3 : poColumn->SetStringField("dataType", "string");
1489 3 : poColumn->SetIntField("maxNumChars", nElementSize);
1490 : }
1491 : else
1492 : {
1493 : // Cannot deal with any of the others yet.
1494 0 : CPLError(CE_Failure, CPLE_NotSupported,
1495 : "Writing this data type in a column is not supported "
1496 : "for this Raster Attribute Table.");
1497 0 : return CE_Failure;
1498 : }
1499 :
1500 18 : const int nOffset = HFAAllocateSpace(hHFA->papoBand[nBand - 1]->psInfo,
1501 9 : nRows * nElementSize);
1502 9 : poColumn->SetIntField("columnDataPtr", nOffset);
1503 :
1504 9 : if (bConvertColors)
1505 : {
1506 : // GDAL Int column
1507 0 : eFieldType = GFT_Integer;
1508 : }
1509 :
1510 9 : AddColumn(pszFieldName, eFieldType, eFieldUsage, nOffset, nElementSize,
1511 : poColumn, false, bConvertColors);
1512 :
1513 9 : return CE_None;
1514 : }
1515 :
1516 : /************************************************************************/
1517 : /* SetLinearBinning() */
1518 : /************************************************************************/
1519 :
1520 1 : CPLErr HFARasterAttributeTable::SetLinearBinning(double dfRow0MinIn,
1521 : double dfBinSizeIn)
1522 : {
1523 1 : if (eAccess == GA_ReadOnly)
1524 : {
1525 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
1526 : "Dataset not open in update mode");
1527 0 : return CE_Failure;
1528 : }
1529 :
1530 1 : bLinearBinning = true;
1531 1 : dfRow0Min = dfRow0MinIn;
1532 1 : dfBinSize = dfBinSizeIn;
1533 :
1534 : // Do we have a descriptor table already?
1535 1 : if (poDT == nullptr || !EQUAL(poDT->GetType(), "Edsc_Table"))
1536 0 : CreateDT();
1537 :
1538 : // We should have an Edsc_BinFunction.
1539 1 : HFAEntry *poBinFunction = poDT->GetNamedChild("#Bin_Function#");
1540 1 : if (poBinFunction == nullptr ||
1541 0 : !EQUAL(poBinFunction->GetType(), "Edsc_BinFunction"))
1542 : {
1543 : poBinFunction =
1544 1 : HFAEntry::New(hHFA->papoBand[nBand - 1]->psInfo, "#Bin_Function#",
1545 : "Edsc_BinFunction", poDT);
1546 : }
1547 :
1548 : // Because of the BaseData we have to hardcode the size.
1549 1 : poBinFunction->MakeData(30);
1550 :
1551 1 : poBinFunction->SetStringField("binFunction", "direct");
1552 1 : poBinFunction->SetDoubleField("minLimit", dfRow0Min);
1553 1 : poBinFunction->SetDoubleField("maxLimit",
1554 1 : (nRows - 1) * dfBinSize + dfRow0Min);
1555 1 : poBinFunction->SetIntField("numBins", nRows);
1556 :
1557 1 : return CE_None;
1558 : }
1559 :
1560 : /************************************************************************/
1561 : /* GetLinearBinning() */
1562 : /************************************************************************/
1563 :
1564 9 : int HFARasterAttributeTable::GetLinearBinning(double *pdfRow0Min,
1565 : double *pdfBinSize) const
1566 : {
1567 9 : if (!bLinearBinning)
1568 2 : return FALSE;
1569 :
1570 7 : *pdfRow0Min = dfRow0Min;
1571 7 : *pdfBinSize = dfBinSize;
1572 :
1573 7 : return TRUE;
1574 : }
1575 :
1576 : /************************************************************************/
1577 : /* Serialize() */
1578 : /************************************************************************/
1579 :
1580 1 : CPLXMLNode *HFARasterAttributeTable::Serialize() const
1581 : {
1582 2 : if (GetRowCount() != 0 &&
1583 1 : GetColumnCount() > RAT_MAX_ELEM_FOR_CLONE / GetRowCount())
1584 0 : return nullptr;
1585 :
1586 1 : return GDALRasterAttributeTable::Serialize();
1587 : }
1588 :
1589 : /************************************************************************/
1590 : /* SetTableType() */
1591 : /************************************************************************/
1592 :
1593 : CPLErr
1594 617 : HFARasterAttributeTable::SetTableType(const GDALRATTableType eInTableType)
1595 : {
1596 617 : eTableType = eInTableType;
1597 617 : return CE_None;
1598 : }
1599 :
1600 : /************************************************************************/
1601 : /* GetTableType() */
1602 : /************************************************************************/
1603 :
1604 25 : GDALRATTableType HFARasterAttributeTable::GetTableType() const
1605 : {
1606 25 : return eTableType;
1607 : }
1608 :
1609 0 : void HFARasterAttributeTable::RemoveStatistics()
1610 : {
1611 : // since we are storing the fields in a vector it will generally
1612 : // be faster to create a new vector and replace the old one
1613 : // rather than actually erasing columns.
1614 0 : std::vector<HFAAttributeField> aoNewFields;
1615 0 : for (const auto &field : aoFields)
1616 : {
1617 0 : switch (field.eUsage)
1618 : {
1619 0 : case GFU_PixelCount:
1620 : case GFU_Min:
1621 : case GFU_Max:
1622 : case GFU_RedMin:
1623 : case GFU_GreenMin:
1624 : case GFU_BlueMin:
1625 : case GFU_AlphaMin:
1626 : case GFU_RedMax:
1627 : case GFU_GreenMax:
1628 : case GFU_BlueMax:
1629 : case GFU_AlphaMax:
1630 : {
1631 0 : break;
1632 : }
1633 :
1634 0 : default:
1635 0 : if (field.sName != "Histogram")
1636 : {
1637 0 : aoNewFields.push_back(field);
1638 : }
1639 : }
1640 : }
1641 0 : aoFields = std::move(aoNewFields);
1642 0 : }
1643 :
1644 : /************************************************************************/
1645 : /* HFARasterBand() */
1646 : /************************************************************************/
1647 :
1648 : namespace
1649 : {
1650 :
1651 : // Convert 0..1 input color range to 0..255.
1652 : // Clamp overflow and underflow.
1653 11260 : short ColorToShort(double val)
1654 : {
1655 11260 : const double dfScaled = val * 256.0;
1656 : // Clamp to [0..255].
1657 11260 : const double dfClamped = std::max(0.0, std::min(255.0, dfScaled));
1658 11260 : return static_cast<short>(dfClamped);
1659 : }
1660 :
1661 : } // namespace
1662 :
1663 670 : HFARasterBand::HFARasterBand(HFADataset *poDSIn, int nBandIn, int iOverview)
1664 : : poCT(nullptr),
1665 : // eHFADataType
1666 : nOverviews(-1), nThisOverview(iOverview), papoOverviewBands(nullptr),
1667 670 : hHFA(poDSIn->hHFA), bMetadataDirty(false), poDefaultRAT(nullptr)
1668 : {
1669 670 : if (iOverview == -1)
1670 616 : poDS = poDSIn;
1671 : else
1672 54 : poDS = nullptr;
1673 :
1674 670 : nBand = nBandIn;
1675 670 : eAccess = poDSIn->GetAccess();
1676 :
1677 670 : int nCompression = 0;
1678 670 : HFAGetBandInfo(hHFA, nBand, &eHFADataType, &nBlockXSize, &nBlockYSize,
1679 : &nCompression);
1680 :
1681 : // If this is an overview, we need to fetch the actual size,
1682 : // and block size.
1683 670 : if (iOverview > -1)
1684 : {
1685 : EPTType eHFADataTypeO;
1686 :
1687 54 : nOverviews = 0;
1688 54 : if (HFAGetOverviewInfo(hHFA, nBand, iOverview, &nRasterXSize,
1689 : &nRasterYSize, &nBlockXSize, &nBlockYSize,
1690 54 : &eHFADataTypeO) != CE_None)
1691 : {
1692 0 : nRasterXSize = 0;
1693 0 : nRasterYSize = 0;
1694 0 : return;
1695 : }
1696 :
1697 : // If we are an 8bit overview of a 1bit layer, we need to mark
1698 : // ourselves as being "resample: average_bit2grayscale".
1699 54 : if (eHFADataType == EPT_u1 && eHFADataTypeO == EPT_u8)
1700 : {
1701 5 : GDALMajorObject::SetMetadataItem("RESAMPLING",
1702 : "AVERAGE_BIT2GRAYSCALE");
1703 5 : GDALMajorObject::SetMetadataItem("NBITS", "8");
1704 : }
1705 54 : eHFADataType = eHFADataTypeO;
1706 : }
1707 :
1708 : // Set some other information.
1709 670 : if (nCompression != 0)
1710 73 : GDALMajorObject::SetMetadataItem("COMPRESSION", "RLE",
1711 : "IMAGE_STRUCTURE");
1712 :
1713 670 : switch (eHFADataType)
1714 : {
1715 428 : case EPT_u1:
1716 : case EPT_u2:
1717 : case EPT_u4:
1718 : case EPT_u8:
1719 428 : eDataType = GDT_Byte;
1720 428 : break;
1721 :
1722 11 : case EPT_s8:
1723 11 : eDataType = GDT_Int8;
1724 11 : break;
1725 :
1726 36 : case EPT_u16:
1727 36 : eDataType = GDT_UInt16;
1728 36 : break;
1729 :
1730 23 : case EPT_s16:
1731 23 : eDataType = GDT_Int16;
1732 23 : break;
1733 :
1734 25 : case EPT_u32:
1735 25 : eDataType = GDT_UInt32;
1736 25 : break;
1737 :
1738 40 : case EPT_s32:
1739 40 : eDataType = GDT_Int32;
1740 40 : break;
1741 :
1742 33 : case EPT_f32:
1743 33 : eDataType = GDT_Float32;
1744 33 : break;
1745 :
1746 32 : case EPT_f64:
1747 32 : eDataType = GDT_Float64;
1748 32 : break;
1749 :
1750 21 : case EPT_c64:
1751 21 : eDataType = GDT_CFloat32;
1752 21 : break;
1753 :
1754 21 : case EPT_c128:
1755 21 : eDataType = GDT_CFloat64;
1756 21 : break;
1757 :
1758 0 : default:
1759 0 : eDataType = GDT_Byte;
1760 : // This should really report an error, but this isn't
1761 : // so easy from within constructors.
1762 0 : CPLDebug("GDAL", "Unsupported pixel type in HFARasterBand: %d.",
1763 0 : eHFADataType);
1764 0 : break;
1765 : }
1766 :
1767 670 : if (HFAGetDataTypeBits(eHFADataType) < 8)
1768 : {
1769 18 : GDALMajorObject::SetMetadataItem(
1770 36 : "NBITS", CPLString().Printf("%d", HFAGetDataTypeBits(eHFADataType)),
1771 : "IMAGE_STRUCTURE");
1772 : }
1773 :
1774 : // Collect color table if present.
1775 670 : double *padfRed = nullptr;
1776 670 : double *padfGreen = nullptr;
1777 670 : double *padfBlue = nullptr;
1778 670 : double *padfAlpha = nullptr;
1779 670 : double *padfBins = nullptr;
1780 670 : int nColors = 0;
1781 :
1782 1286 : if (iOverview == -1 &&
1783 616 : HFAGetPCT(hHFA, nBand, &nColors, &padfRed, &padfGreen, &padfBlue,
1784 1286 : &padfAlpha, &padfBins) == CE_None &&
1785 10 : nColors > 0)
1786 : {
1787 10 : poCT = new GDALColorTable();
1788 2825 : for (int iColor = 0; iColor < nColors; iColor++)
1789 : {
1790 : // The following mapping assigns "equal sized" section of
1791 : // the [0...1] range to each possible output value and avoid
1792 : // rounding issues for the "normal" values generated using n/255.
1793 : // See bug #1732 for some discussion.
1794 2815 : GDALColorEntry sEntry = {ColorToShort(padfRed[iColor]),
1795 5630 : ColorToShort(padfGreen[iColor]),
1796 5630 : ColorToShort(padfBlue[iColor]),
1797 2815 : ColorToShort(padfAlpha[iColor])};
1798 :
1799 2815 : if (padfBins != nullptr)
1800 : {
1801 300 : const double dfIdx = padfBins[iColor];
1802 300 : if (!(dfIdx >= 0.0 && dfIdx <= 65535.0))
1803 : {
1804 0 : CPLError(CE_Failure, CPLE_NotSupported,
1805 : "Invalid index padfBins[%d] = %g", iColor, dfIdx);
1806 0 : break;
1807 : }
1808 : else
1809 : {
1810 300 : poCT->SetColorEntry(static_cast<int>(dfIdx), &sEntry);
1811 : }
1812 : }
1813 : else
1814 : {
1815 2515 : poCT->SetColorEntry(iColor, &sEntry);
1816 : }
1817 : }
1818 : }
1819 : }
1820 :
1821 : /************************************************************************/
1822 : /* ~HFARasterBand() */
1823 : /************************************************************************/
1824 :
1825 1340 : HFARasterBand::~HFARasterBand()
1826 :
1827 : {
1828 670 : FlushCache(true);
1829 :
1830 723 : for (int iOvIndex = 0; iOvIndex < nOverviews; iOvIndex++)
1831 : {
1832 53 : delete papoOverviewBands[iOvIndex];
1833 : }
1834 670 : CPLFree(papoOverviewBands);
1835 :
1836 670 : if (poCT != nullptr)
1837 9 : delete poCT;
1838 :
1839 670 : if (poDefaultRAT)
1840 616 : delete poDefaultRAT;
1841 1340 : }
1842 :
1843 : /************************************************************************/
1844 : /* ReadAuxMetadata() */
1845 : /************************************************************************/
1846 :
1847 616 : void HFARasterBand::ReadAuxMetadata()
1848 :
1849 : {
1850 : // Only load metadata for full resolution layer.
1851 616 : if (nThisOverview != -1)
1852 0 : return;
1853 :
1854 616 : HFABand *poBand = hHFA->papoBand[nBand - 1];
1855 :
1856 616 : const char *const *pszAuxMetaData = GetHFAAuxMetaDataList();
1857 9240 : for (int i = 0; pszAuxMetaData[i] != nullptr; i += 4)
1858 : {
1859 : HFAEntry *poEntry;
1860 8624 : if (strlen(pszAuxMetaData[i]) > 0)
1861 : {
1862 8008 : poEntry = poBand->poNode->GetNamedChild(pszAuxMetaData[i]);
1863 8008 : if (poEntry == nullptr)
1864 7078 : continue;
1865 : }
1866 : else
1867 : {
1868 616 : poEntry = poBand->poNode;
1869 616 : assert(poEntry);
1870 : }
1871 :
1872 1546 : const char *pszFieldName = pszAuxMetaData[i + 1] + 1;
1873 :
1874 1546 : switch (pszAuxMetaData[i + 1][0])
1875 : {
1876 692 : case 'd':
1877 : {
1878 1384 : CPLString osValueList;
1879 :
1880 692 : CPLErr eErr = CE_None;
1881 692 : int nCount = poEntry->GetFieldCount(pszFieldName, &eErr);
1882 692 : if (nCount > 65536)
1883 : {
1884 0 : nCount = 65536;
1885 0 : CPLDebug("HFA", "Limiting %s to %d entries",
1886 0 : pszAuxMetaData[i + 2], nCount);
1887 : }
1888 1326 : for (int iValue = 0; eErr == CE_None && iValue < nCount;
1889 : iValue++)
1890 : {
1891 660 : CPLString osSubFieldName;
1892 660 : osSubFieldName.Printf("%s[%d]", pszFieldName, iValue);
1893 : const double dfValue =
1894 660 : poEntry->GetDoubleField(osSubFieldName, &eErr);
1895 660 : if (eErr != CE_None)
1896 26 : break;
1897 :
1898 634 : char szValueAsString[100] = {};
1899 634 : CPLsnprintf(szValueAsString, sizeof(szValueAsString),
1900 : "%.14g", dfValue);
1901 :
1902 634 : if (iValue > 0)
1903 2 : osValueList += ",";
1904 634 : osValueList += szValueAsString;
1905 : }
1906 692 : if (eErr == CE_None)
1907 666 : SetMetadataItem(pszAuxMetaData[i + 2], osValueList);
1908 : }
1909 692 : break;
1910 229 : case 'i':
1911 : case 'l':
1912 : {
1913 458 : CPLString osValueList;
1914 :
1915 229 : CPLErr eErr = CE_None;
1916 229 : int nCount = poEntry->GetFieldCount(pszFieldName, &eErr);
1917 229 : if (nCount > 65536)
1918 : {
1919 0 : nCount = 65536;
1920 0 : CPLDebug("HFA", "Limiting %s to %d entries",
1921 0 : pszAuxMetaData[i + 2], nCount);
1922 : }
1923 445 : for (int iValue = 0; eErr == CE_None && iValue < nCount;
1924 : iValue++)
1925 : {
1926 229 : CPLString osSubFieldName;
1927 229 : osSubFieldName.Printf("%s[%d]", pszFieldName, iValue);
1928 229 : int nValue = poEntry->GetIntField(osSubFieldName, &eErr);
1929 229 : if (eErr != CE_None)
1930 13 : break;
1931 :
1932 216 : char szValueAsString[100] = {};
1933 216 : snprintf(szValueAsString, sizeof(szValueAsString), "%d",
1934 : nValue);
1935 :
1936 216 : if (iValue > 0)
1937 0 : osValueList += ",";
1938 216 : osValueList += szValueAsString;
1939 : }
1940 229 : if (eErr == CE_None)
1941 216 : SetMetadataItem(pszAuxMetaData[i + 2], osValueList);
1942 : }
1943 229 : break;
1944 625 : case 's':
1945 : case 'e':
1946 : {
1947 625 : CPLErr eErr = CE_None;
1948 : const char *pszValue =
1949 625 : poEntry->GetStringField(pszFieldName, &eErr);
1950 625 : if (eErr == CE_None)
1951 625 : SetMetadataItem(pszAuxMetaData[i + 2], pszValue);
1952 : }
1953 625 : break;
1954 0 : default:
1955 0 : CPLAssert(false);
1956 : }
1957 : }
1958 :
1959 : /* if we have a default RAT we can now set its thematic/athematic state
1960 : from the metadata we just read in */
1961 616 : if (GetDefaultRAT())
1962 : {
1963 616 : const char *psLayerType = GetMetadataItem("LAYER_TYPE", "");
1964 616 : if (psLayerType)
1965 : {
1966 1232 : GetDefaultRAT()->SetTableType(EQUALN(psLayerType, "athematic", 9)
1967 : ? GRTT_ATHEMATIC
1968 616 : : GRTT_THEMATIC);
1969 : }
1970 : }
1971 : }
1972 :
1973 : /************************************************************************/
1974 : /* ReadHistogramMetadata() */
1975 : /************************************************************************/
1976 :
1977 616 : void HFARasterBand::ReadHistogramMetadata()
1978 :
1979 : {
1980 : // Only load metadata for full resolution layer.
1981 616 : if (nThisOverview != -1)
1982 0 : return;
1983 :
1984 616 : HFABand *poBand = hHFA->papoBand[nBand - 1];
1985 :
1986 : HFAEntry *poEntry =
1987 616 : poBand->poNode->GetNamedChild("Descriptor_Table.Histogram");
1988 616 : if (poEntry == nullptr)
1989 539 : return;
1990 :
1991 77 : int nNumBins = poEntry->GetIntField("numRows");
1992 77 : if (nNumBins < 0)
1993 0 : return;
1994 : // TODO(schwehr): Can we do a better/tighter check?
1995 77 : if (nNumBins > 1000000)
1996 : {
1997 0 : CPLError(CE_Failure, CPLE_FileIO, "Unreasonably large histogram: %d",
1998 : nNumBins);
1999 0 : return;
2000 : }
2001 :
2002 : // Fetch the histogram values.
2003 77 : const int nOffset = poEntry->GetIntField("columnDataPtr");
2004 77 : const char *pszType = poEntry->GetStringField("dataType");
2005 77 : int nBinSize = 4;
2006 :
2007 77 : if (pszType != nullptr && STARTS_WITH_CI(pszType, "real"))
2008 75 : nBinSize = 8;
2009 :
2010 : GUIntBig *panHistValues = static_cast<GUIntBig *>(
2011 77 : VSI_MALLOC2_VERBOSE(sizeof(GUIntBig), nNumBins));
2012 : GByte *pabyWorkBuf =
2013 77 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nBinSize, nNumBins));
2014 :
2015 77 : if (panHistValues == nullptr || pabyWorkBuf == nullptr)
2016 : {
2017 0 : VSIFree(panHistValues);
2018 0 : VSIFree(pabyWorkBuf);
2019 0 : return;
2020 : }
2021 :
2022 154 : if (VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) != 0 ||
2023 : static_cast<int>(
2024 77 : VSIFReadL(pabyWorkBuf, nBinSize, nNumBins, hHFA->fp)) != nNumBins)
2025 : {
2026 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot read histogram values.");
2027 0 : CPLFree(panHistValues);
2028 0 : CPLFree(pabyWorkBuf);
2029 0 : return;
2030 : }
2031 :
2032 : // Swap into local order.
2033 16620 : for (int i = 0; i < nNumBins; i++)
2034 : HFAStandard(nBinSize, pabyWorkBuf + i * nBinSize);
2035 :
2036 77 : if (nBinSize == 8) // Source is doubles.
2037 : {
2038 75 : const double *padfWorkBuf = reinterpret_cast<double *>(pabyWorkBuf);
2039 16106 : for (int i = 0; i < nNumBins; i++)
2040 : {
2041 16031 : const double dfNumber = padfWorkBuf[i];
2042 16031 : if (dfNumber >=
2043 16031 : static_cast<double>(std::numeric_limits<GUIntBig>::max()) ||
2044 : dfNumber <
2045 32062 : static_cast<double>(std::numeric_limits<GUIntBig>::min()) ||
2046 16031 : std::isnan(dfNumber))
2047 : {
2048 0 : CPLError(CE_Failure, CPLE_FileIO, "Out of range hist vals.");
2049 0 : CPLFree(panHistValues);
2050 0 : CPLFree(pabyWorkBuf);
2051 0 : return;
2052 : }
2053 16031 : panHistValues[i] = static_cast<GUIntBig>(dfNumber);
2054 : }
2055 : }
2056 : else // Source is 32bit integers.
2057 : {
2058 2 : const int *panWorkBuf = reinterpret_cast<int *>(pabyWorkBuf);
2059 514 : for (int i = 0; i < nNumBins; i++)
2060 : {
2061 512 : const int nNumber = panWorkBuf[i];
2062 : // Positive int should always fit.
2063 512 : if (nNumber < 0)
2064 : {
2065 0 : CPLError(CE_Failure, CPLE_FileIO, "Out of range hist vals.");
2066 0 : CPLFree(panHistValues);
2067 0 : CPLFree(pabyWorkBuf);
2068 0 : return;
2069 : }
2070 512 : panHistValues[i] = static_cast<GUIntBig>(nNumber);
2071 : }
2072 : }
2073 :
2074 77 : CPLFree(pabyWorkBuf);
2075 77 : pabyWorkBuf = nullptr;
2076 :
2077 : // Do we have unique values for the bins?
2078 77 : double *padfBinValues = nullptr;
2079 : HFAEntry *poBinEntry =
2080 77 : poBand->poNode->GetNamedChild("Descriptor_Table.#Bin_Function840#");
2081 :
2082 90 : if (poBinEntry != nullptr &&
2083 13 : EQUAL(poBinEntry->GetType(), "Edsc_BinFunction840"))
2084 : {
2085 : const char *pszValue =
2086 13 : poBinEntry->GetStringField("binFunction.type.string");
2087 13 : if (pszValue && EQUAL(pszValue, "BFUnique"))
2088 13 : padfBinValues = HFAReadBFUniqueBins(poBinEntry, nNumBins);
2089 : }
2090 :
2091 77 : if (padfBinValues)
2092 : {
2093 13 : int nMaxValue = 0;
2094 13 : int nMinValue = 1000000;
2095 :
2096 1469 : for (int i = 0; i < nNumBins; i++)
2097 : {
2098 1456 : const double dfCurrent = padfBinValues[i];
2099 :
2100 1456 : if (dfCurrent != floor(dfCurrent) || /* not an integer value */
2101 1456 : dfCurrent < 0.0 || dfCurrent > 1000.0)
2102 : {
2103 0 : CPLFree(padfBinValues);
2104 0 : CPLFree(panHistValues);
2105 0 : CPLDebug("HFA",
2106 : "Unable to offer histogram because unique values "
2107 : "list is not convenient to reform as HISTOBINVALUES.");
2108 0 : return;
2109 : }
2110 :
2111 1456 : nMaxValue = std::max(nMaxValue, static_cast<int>(dfCurrent));
2112 1456 : nMinValue = std::min(nMinValue, static_cast<int>(dfCurrent));
2113 : }
2114 :
2115 13 : const int nNewBins = nMaxValue + 1;
2116 : GUIntBig *panNewHistValues =
2117 13 : static_cast<GUIntBig *>(CPLCalloc(sizeof(GUIntBig), nNewBins));
2118 :
2119 1469 : for (int i = 0; i < nNumBins; i++)
2120 1456 : panNewHistValues[static_cast<int>(padfBinValues[i])] =
2121 1456 : panHistValues[i];
2122 :
2123 13 : CPLFree(panHistValues);
2124 13 : panHistValues = panNewHistValues;
2125 13 : nNumBins = nNewBins;
2126 :
2127 13 : SetMetadataItem("STATISTICS_HISTOMIN", "0");
2128 13 : SetMetadataItem("STATISTICS_HISTOMAX",
2129 26 : CPLString().Printf("%d", nMaxValue));
2130 13 : SetMetadataItem("STATISTICS_HISTONUMBINS",
2131 26 : CPLString().Printf("%d", nMaxValue + 1));
2132 :
2133 13 : CPLFree(padfBinValues);
2134 13 : padfBinValues = nullptr;
2135 : }
2136 :
2137 : // Format into HISTOBINVALUES text format.
2138 77 : unsigned int nBufSize = 1024;
2139 77 : char *pszBinValues = static_cast<char *>(CPLMalloc(nBufSize));
2140 77 : pszBinValues[0] = 0;
2141 77 : int nBinValuesLen = 0;
2142 :
2143 18103 : for (int nBin = 0; nBin < nNumBins; ++nBin)
2144 : {
2145 18026 : char szBuf[32] = {};
2146 18026 : snprintf(szBuf, 31, CPL_FRMT_GUIB, panHistValues[nBin]);
2147 18026 : if ((nBinValuesLen + strlen(szBuf) + 2) > nBufSize)
2148 : {
2149 6 : nBufSize *= 2;
2150 : char *pszNewBinValues = static_cast<char *>(
2151 6 : VSI_REALLOC_VERBOSE(pszBinValues, nBufSize));
2152 6 : if (pszNewBinValues == nullptr)
2153 : {
2154 0 : break;
2155 : }
2156 6 : pszBinValues = pszNewBinValues;
2157 : }
2158 18026 : strcat(pszBinValues + nBinValuesLen, szBuf);
2159 18026 : strcat(pszBinValues + nBinValuesLen, "|");
2160 18026 : nBinValuesLen += static_cast<int>(strlen(pszBinValues + nBinValuesLen));
2161 : }
2162 :
2163 77 : SetMetadataItem("STATISTICS_HISTOBINVALUES", pszBinValues);
2164 77 : CPLFree(panHistValues);
2165 77 : CPLFree(pszBinValues);
2166 : }
2167 :
2168 : /************************************************************************/
2169 : /* GetNoDataValue() */
2170 : /************************************************************************/
2171 :
2172 193 : double HFARasterBand::GetNoDataValue(int *pbSuccess)
2173 :
2174 : {
2175 193 : double dfNoData = 0.0;
2176 :
2177 193 : if (HFAGetBandNoData(hHFA, nBand, &dfNoData))
2178 : {
2179 21 : if (pbSuccess)
2180 17 : *pbSuccess = TRUE;
2181 21 : return dfNoData;
2182 : }
2183 :
2184 172 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
2185 : }
2186 :
2187 : /************************************************************************/
2188 : /* SetNoDataValue() */
2189 : /************************************************************************/
2190 :
2191 9 : CPLErr HFARasterBand::SetNoDataValue(double dfValue)
2192 : {
2193 9 : return HFASetBandNoData(hHFA, nBand, dfValue);
2194 : }
2195 :
2196 : /************************************************************************/
2197 : /* GetMinimum() */
2198 : /************************************************************************/
2199 :
2200 4 : double HFARasterBand::GetMinimum(int *pbSuccess)
2201 :
2202 : {
2203 4 : const char *pszValue = GetMetadataItem("STATISTICS_MINIMUM");
2204 :
2205 4 : if (pszValue != nullptr)
2206 : {
2207 3 : if (pbSuccess)
2208 3 : *pbSuccess = TRUE;
2209 3 : return CPLAtofM(pszValue);
2210 : }
2211 :
2212 1 : return GDALRasterBand::GetMinimum(pbSuccess);
2213 : }
2214 :
2215 : /************************************************************************/
2216 : /* GetMaximum() */
2217 : /************************************************************************/
2218 :
2219 4 : double HFARasterBand::GetMaximum(int *pbSuccess)
2220 :
2221 : {
2222 4 : const char *pszValue = GetMetadataItem("STATISTICS_MAXIMUM");
2223 :
2224 4 : if (pszValue != nullptr)
2225 : {
2226 3 : if (pbSuccess)
2227 3 : *pbSuccess = TRUE;
2228 3 : return CPLAtofM(pszValue);
2229 : }
2230 :
2231 1 : return GDALRasterBand::GetMaximum(pbSuccess);
2232 : }
2233 :
2234 : /************************************************************************/
2235 : /* EstablishOverviews() */
2236 : /* */
2237 : /* Delayed population of overview information. */
2238 : /************************************************************************/
2239 :
2240 277 : void HFARasterBand::EstablishOverviews()
2241 :
2242 : {
2243 277 : if (nOverviews != -1)
2244 155 : return;
2245 :
2246 122 : nOverviews = HFAGetOverviewCount(hHFA, nBand);
2247 122 : if (nOverviews > 0)
2248 : {
2249 29 : papoOverviewBands = static_cast<HFARasterBand **>(
2250 29 : CPLMalloc(sizeof(void *) * nOverviews));
2251 :
2252 72 : for (int iOvIndex = 0; iOvIndex < nOverviews; iOvIndex++)
2253 : {
2254 43 : papoOverviewBands[iOvIndex] = new HFARasterBand(
2255 43 : cpl::down_cast<HFADataset *>(poDS), nBand, iOvIndex);
2256 43 : if (papoOverviewBands[iOvIndex]->GetXSize() == 0)
2257 : {
2258 0 : delete papoOverviewBands[iOvIndex];
2259 0 : papoOverviewBands[iOvIndex] = nullptr;
2260 : }
2261 : }
2262 : }
2263 : }
2264 :
2265 : /************************************************************************/
2266 : /* GetOverviewCount() */
2267 : /************************************************************************/
2268 :
2269 185 : int HFARasterBand::GetOverviewCount()
2270 :
2271 : {
2272 185 : EstablishOverviews();
2273 :
2274 185 : if (nOverviews == 0)
2275 91 : return GDALRasterBand::GetOverviewCount();
2276 :
2277 94 : return nOverviews;
2278 : }
2279 :
2280 : /************************************************************************/
2281 : /* GetOverview() */
2282 : /************************************************************************/
2283 :
2284 81 : GDALRasterBand *HFARasterBand::GetOverview(int i)
2285 :
2286 : {
2287 81 : EstablishOverviews();
2288 :
2289 81 : if (nOverviews == 0)
2290 0 : return GDALRasterBand::GetOverview(i);
2291 81 : else if (i < 0 || i >= nOverviews)
2292 0 : return nullptr;
2293 : else
2294 81 : return papoOverviewBands[i];
2295 : }
2296 :
2297 : /************************************************************************/
2298 : /* IReadBlock() */
2299 : /************************************************************************/
2300 :
2301 1354 : CPLErr HFARasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
2302 :
2303 : {
2304 1354 : CPLErr eErr = CE_None;
2305 :
2306 1354 : if (nThisOverview == -1)
2307 1330 : eErr = HFAGetRasterBlockEx(hHFA, nBand, nBlockXOff, nBlockYOff, pImage,
2308 1330 : nBlockXSize * nBlockYSize *
2309 1330 : GDALGetDataTypeSizeBytes(eDataType));
2310 : else
2311 24 : eErr = HFAGetOverviewRasterBlockEx(
2312 : hHFA, nBand, nThisOverview, nBlockXOff, nBlockYOff, pImage,
2313 24 : nBlockXSize * nBlockYSize * GDALGetDataTypeSizeBytes(eDataType));
2314 :
2315 1354 : if (eErr == CE_None && eHFADataType == EPT_u4)
2316 : {
2317 2 : GByte *pabyData = static_cast<GByte *>(pImage);
2318 :
2319 4098 : for (int ii = nBlockXSize * nBlockYSize - 2; ii >= 0; ii -= 2)
2320 : {
2321 4096 : int k = ii >> 1;
2322 4096 : pabyData[ii + 1] = (pabyData[k] >> 4) & 0xf;
2323 4096 : pabyData[ii] = (pabyData[k]) & 0xf;
2324 : }
2325 : }
2326 1354 : if (eErr == CE_None && eHFADataType == EPT_u2)
2327 : {
2328 6 : GByte *pabyData = static_cast<GByte *>(pImage);
2329 :
2330 6150 : for (int ii = nBlockXSize * nBlockYSize - 4; ii >= 0; ii -= 4)
2331 : {
2332 6144 : int k = ii >> 2;
2333 6144 : pabyData[ii + 3] = (pabyData[k] >> 6) & 0x3;
2334 6144 : pabyData[ii + 2] = (pabyData[k] >> 4) & 0x3;
2335 6144 : pabyData[ii + 1] = (pabyData[k] >> 2) & 0x3;
2336 6144 : pabyData[ii] = (pabyData[k]) & 0x3;
2337 : }
2338 : }
2339 1354 : if (eErr == CE_None && eHFADataType == EPT_u1)
2340 : {
2341 52 : GByte *pabyData = static_cast<GByte *>(pImage);
2342 :
2343 213044 : for (int ii = nBlockXSize * nBlockYSize - 1; ii >= 0; ii--)
2344 : {
2345 212992 : if ((pabyData[ii >> 3] & (1 << (ii & 0x7))))
2346 160690 : pabyData[ii] = 1;
2347 : else
2348 52302 : pabyData[ii] = 0;
2349 : }
2350 : }
2351 :
2352 1354 : return eErr;
2353 : }
2354 :
2355 : /************************************************************************/
2356 : /* IWriteBlock() */
2357 : /************************************************************************/
2358 :
2359 133 : CPLErr HFARasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)
2360 :
2361 : {
2362 133 : GByte *pabyOutBuf = static_cast<GByte *>(pImage);
2363 :
2364 : // Do we need to pack 1/2/4 bit data?
2365 : // TODO(schwehr): Make symbolic constants with explanations.
2366 133 : if (eHFADataType == EPT_u1 || eHFADataType == EPT_u2 ||
2367 129 : eHFADataType == EPT_u4)
2368 : {
2369 6 : const int nPixCount = nBlockXSize * nBlockYSize;
2370 6 : pabyOutBuf = static_cast<GByte *>(VSIMalloc2(nBlockXSize, nBlockYSize));
2371 6 : if (pabyOutBuf == nullptr)
2372 0 : return CE_Failure;
2373 :
2374 6 : if (eHFADataType == EPT_u1)
2375 : {
2376 1026 : for (int ii = 0; ii < nPixCount - 7; ii += 8)
2377 : {
2378 1024 : const int k = ii >> 3;
2379 : // TODO(schwehr): Create a temp for (GByte *)pImage.
2380 1024 : pabyOutBuf[k] = (((GByte *)pImage)[ii] & 0x1) |
2381 1024 : ((((GByte *)pImage)[ii + 1] & 0x1) << 1) |
2382 1024 : ((((GByte *)pImage)[ii + 2] & 0x1) << 2) |
2383 1024 : ((((GByte *)pImage)[ii + 3] & 0x1) << 3) |
2384 1024 : ((((GByte *)pImage)[ii + 4] & 0x1) << 4) |
2385 1024 : ((((GByte *)pImage)[ii + 5] & 0x1) << 5) |
2386 1024 : ((((GByte *)pImage)[ii + 6] & 0x1) << 6) |
2387 1024 : ((((GByte *)pImage)[ii + 7] & 0x1) << 7);
2388 : }
2389 : }
2390 4 : else if (eHFADataType == EPT_u2)
2391 : {
2392 2050 : for (int ii = 0; ii < nPixCount - 3; ii += 4)
2393 : {
2394 2048 : const int k = ii >> 2;
2395 2048 : pabyOutBuf[k] = (((GByte *)pImage)[ii] & 0x3) |
2396 2048 : ((((GByte *)pImage)[ii + 1] & 0x3) << 2) |
2397 2048 : ((((GByte *)pImage)[ii + 2] & 0x3) << 4) |
2398 2048 : ((((GByte *)pImage)[ii + 3] & 0x3) << 6);
2399 : }
2400 : }
2401 2 : else if (eHFADataType == EPT_u4)
2402 : {
2403 4098 : for (int ii = 0; ii < nPixCount - 1; ii += 2)
2404 : {
2405 4096 : const int k = ii >> 1;
2406 4096 : pabyOutBuf[k] = (((GByte *)pImage)[ii] & 0xf) |
2407 4096 : ((((GByte *)pImage)[ii + 1] & 0xf) << 4);
2408 : }
2409 : }
2410 : }
2411 :
2412 : // Actually write out.
2413 : const CPLErr nRetCode =
2414 133 : nThisOverview == -1
2415 133 : ? HFASetRasterBlock(hHFA, nBand, nBlockXOff, nBlockYOff, pabyOutBuf)
2416 28 : : HFASetOverviewRasterBlock(hHFA, nBand, nThisOverview, nBlockXOff,
2417 133 : nBlockYOff, pabyOutBuf);
2418 :
2419 133 : if (pabyOutBuf != pImage)
2420 6 : CPLFree(pabyOutBuf);
2421 :
2422 133 : return nRetCode;
2423 : }
2424 :
2425 : /************************************************************************/
2426 : /* GetDescription() */
2427 : /************************************************************************/
2428 :
2429 213 : const char *HFARasterBand::GetDescription() const
2430 : {
2431 213 : const char *pszName = HFAGetBandName(hHFA, nBand);
2432 :
2433 213 : if (pszName == nullptr)
2434 0 : return GDALPamRasterBand::GetDescription();
2435 :
2436 213 : return pszName;
2437 : }
2438 :
2439 : /************************************************************************/
2440 : /* SetDescription() */
2441 : /************************************************************************/
2442 7 : void HFARasterBand::SetDescription(const char *pszName)
2443 : {
2444 7 : if (strlen(pszName) > 0)
2445 7 : HFASetBandName(hHFA, nBand, pszName);
2446 7 : }
2447 :
2448 : /************************************************************************/
2449 : /* GetColorInterpretation() */
2450 : /************************************************************************/
2451 :
2452 55 : GDALColorInterp HFARasterBand::GetColorInterpretation()
2453 :
2454 : {
2455 55 : if (poCT != nullptr)
2456 4 : return GCI_PaletteIndex;
2457 :
2458 51 : return GCI_Undefined;
2459 : }
2460 :
2461 : /************************************************************************/
2462 : /* GetColorTable() */
2463 : /************************************************************************/
2464 :
2465 40 : GDALColorTable *HFARasterBand::GetColorTable()
2466 : {
2467 40 : return poCT;
2468 : }
2469 :
2470 : /************************************************************************/
2471 : /* SetColorTable() */
2472 : /************************************************************************/
2473 :
2474 3 : CPLErr HFARasterBand::SetColorTable(GDALColorTable *poCTable)
2475 :
2476 : {
2477 3 : if (GetAccess() == GA_ReadOnly)
2478 : {
2479 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
2480 : "Unable to set color table on read-only file.");
2481 0 : return CE_Failure;
2482 : }
2483 :
2484 : // Special case if we are clearing the color table.
2485 3 : if (poCTable == nullptr)
2486 : {
2487 2 : delete poCT;
2488 2 : poCT = nullptr;
2489 :
2490 2 : HFASetPCT(hHFA, nBand, 0, nullptr, nullptr, nullptr, nullptr);
2491 :
2492 2 : return CE_None;
2493 : }
2494 :
2495 : // Write out the colortable, and update the configuration.
2496 1 : int nColors = poCTable->GetColorEntryCount();
2497 :
2498 : /* -------------------------------------------------------------------- */
2499 : /* If we already have a non-empty RAT set and it's smaller than */
2500 : /* the colour table, and all the trailing CT entries are the same, */
2501 : /* truncate the colour table. Helps when RATs travel via GTiff. */
2502 : /* -------------------------------------------------------------------- */
2503 1 : const GDALRasterAttributeTable *poRAT = GetDefaultRAT();
2504 1 : if (poRAT != nullptr && poRAT->GetRowCount() > 0 &&
2505 0 : poRAT->GetRowCount() < nColors)
2506 : {
2507 0 : bool match = true;
2508 : const GDALColorEntry *color1 =
2509 0 : poCTable->GetColorEntry(poRAT->GetRowCount());
2510 0 : for (int i = poRAT->GetRowCount() + 1; match && i < nColors; i++)
2511 : {
2512 0 : const GDALColorEntry *color2 = poCTable->GetColorEntry(i);
2513 0 : match = (color1->c1 == color2->c1 && color1->c2 == color2->c2 &&
2514 0 : color1->c3 == color2->c3 && color1->c4 == color2->c4);
2515 : }
2516 0 : if (match)
2517 : {
2518 0 : CPLDebug("HFA",
2519 : "SetColorTable: Truncating PCT size (%d) to RAT size (%d)",
2520 0 : nColors, poRAT->GetRowCount());
2521 0 : nColors = poRAT->GetRowCount();
2522 : }
2523 : }
2524 :
2525 : double *padfRed =
2526 1 : static_cast<double *>(CPLMalloc(sizeof(double) * nColors));
2527 : double *padfGreen =
2528 1 : static_cast<double *>(CPLMalloc(sizeof(double) * nColors));
2529 : double *padfBlue =
2530 1 : static_cast<double *>(CPLMalloc(sizeof(double) * nColors));
2531 : double *padfAlpha =
2532 1 : static_cast<double *>(CPLMalloc(sizeof(double) * nColors));
2533 :
2534 257 : for (int iColor = 0; iColor < nColors; iColor++)
2535 : {
2536 : GDALColorEntry sRGB;
2537 :
2538 256 : poCTable->GetColorEntryAsRGB(iColor, &sRGB);
2539 :
2540 256 : padfRed[iColor] = sRGB.c1 / 255.0;
2541 256 : padfGreen[iColor] = sRGB.c2 / 255.0;
2542 256 : padfBlue[iColor] = sRGB.c3 / 255.0;
2543 256 : padfAlpha[iColor] = sRGB.c4 / 255.0;
2544 : }
2545 :
2546 1 : HFASetPCT(hHFA, nBand, nColors, padfRed, padfGreen, padfBlue, padfAlpha);
2547 :
2548 1 : CPLFree(padfRed);
2549 1 : CPLFree(padfGreen);
2550 1 : CPLFree(padfBlue);
2551 1 : CPLFree(padfAlpha);
2552 :
2553 1 : if (poCT)
2554 0 : delete poCT;
2555 :
2556 1 : poCT = poCTable->Clone();
2557 :
2558 1 : return CE_None;
2559 : }
2560 :
2561 : /************************************************************************/
2562 : /* SetMetadata() */
2563 : /************************************************************************/
2564 :
2565 21 : CPLErr HFARasterBand::SetMetadata(char **papszMDIn, const char *pszDomain)
2566 :
2567 : {
2568 21 : bMetadataDirty = true;
2569 :
2570 21 : return GDALPamRasterBand::SetMetadata(papszMDIn, pszDomain);
2571 : }
2572 :
2573 : /************************************************************************/
2574 : /* SetMetadata() */
2575 : /************************************************************************/
2576 :
2577 1640 : CPLErr HFARasterBand::SetMetadataItem(const char *pszTag, const char *pszValue,
2578 : const char *pszDomain)
2579 :
2580 : {
2581 1640 : bMetadataDirty = true;
2582 :
2583 1640 : return GDALPamRasterBand::SetMetadataItem(pszTag, pszValue, pszDomain);
2584 : }
2585 :
2586 : /************************************************************************/
2587 : /* CleanOverviews() */
2588 : /************************************************************************/
2589 :
2590 1 : CPLErr HFARasterBand::CleanOverviews()
2591 :
2592 : {
2593 1 : if (nOverviews == 0)
2594 0 : return CE_None;
2595 :
2596 : // Clear our reference to overviews as bands.
2597 2 : for (int iOverview = 0; iOverview < nOverviews; iOverview++)
2598 1 : delete papoOverviewBands[iOverview];
2599 :
2600 1 : CPLFree(papoOverviewBands);
2601 1 : papoOverviewBands = nullptr;
2602 1 : nOverviews = 0;
2603 :
2604 : // Search for any RRDNamesList and destroy it.
2605 1 : HFABand *poBand = hHFA->papoBand[nBand - 1];
2606 1 : HFAEntry *poEntry = poBand->poNode->GetNamedChild("RRDNamesList");
2607 1 : if (poEntry != nullptr)
2608 : {
2609 1 : poEntry->RemoveAndDestroy();
2610 : }
2611 :
2612 : // Destroy and subsample layers under our band.
2613 5 : for (HFAEntry *poChild = poBand->poNode->GetChild(); poChild != nullptr;)
2614 : {
2615 4 : HFAEntry *poNext = poChild->GetNext();
2616 :
2617 4 : if (EQUAL(poChild->GetType(), "Eimg_Layer_SubSample"))
2618 0 : poChild->RemoveAndDestroy();
2619 :
2620 4 : poChild = poNext;
2621 : }
2622 :
2623 : // Clean up dependent file if we are the last band under the
2624 : // assumption there will be nothing else referencing it after
2625 : // this.
2626 1 : if (hHFA->psDependent != hHFA && hHFA->psDependent != nullptr)
2627 : {
2628 : const CPLString osFilename =
2629 0 : CPLFormFilenameSafe(hHFA->psDependent->pszPath,
2630 2 : hHFA->psDependent->pszFilename, nullptr);
2631 :
2632 1 : CPL_IGNORE_RET_VAL(HFAClose(hHFA->psDependent));
2633 1 : hHFA->psDependent = nullptr;
2634 :
2635 1 : CPLDebug("HFA", "Unlink(%s)", osFilename.c_str());
2636 1 : VSIUnlink(osFilename);
2637 : }
2638 :
2639 1 : return CE_None;
2640 : }
2641 :
2642 : /************************************************************************/
2643 : /* BuildOverviews() */
2644 : /************************************************************************/
2645 :
2646 11 : CPLErr HFARasterBand::BuildOverviews(const char *pszResampling,
2647 : int nReqOverviews,
2648 : const int *panOverviewList,
2649 : GDALProgressFunc pfnProgress,
2650 : void *pProgressData,
2651 : CSLConstList papszOptions)
2652 :
2653 : {
2654 11 : EstablishOverviews();
2655 :
2656 11 : if (nThisOverview != -1)
2657 : {
2658 0 : CPLError(CE_Failure, CPLE_AppDefined,
2659 : "Attempt to build overviews on an overview layer.");
2660 :
2661 0 : return CE_Failure;
2662 : }
2663 :
2664 11 : if (nReqOverviews == 0)
2665 1 : return CleanOverviews();
2666 :
2667 : GDALRasterBand **papoOvBands = static_cast<GDALRasterBand **>(
2668 10 : CPLCalloc(sizeof(void *), nReqOverviews));
2669 :
2670 : const bool bRegenerate =
2671 10 : CPLTestBool(CSLFetchNameValueDef(papszOptions, "REGENERATE", "YES"));
2672 :
2673 : // Loop over overview levels requested.
2674 22 : for (int iOverview = 0; iOverview < nReqOverviews; iOverview++)
2675 : {
2676 : // Find this overview level.
2677 12 : const int nReqOvLevel = GDALOvLevelAdjust2(panOverviewList[iOverview],
2678 : nRasterXSize, nRasterYSize);
2679 :
2680 20 : for (int i = 0; i < nOverviews && papoOvBands[iOverview] == nullptr;
2681 : i++)
2682 : {
2683 8 : if (papoOverviewBands[i] == nullptr)
2684 : {
2685 0 : CPLDebug("HFA", "Shouldn't happen happened at line %d",
2686 : __LINE__);
2687 0 : continue;
2688 : }
2689 :
2690 24 : const int nThisOvLevel = GDALComputeOvFactor(
2691 8 : papoOverviewBands[i]->GetXSize(), GetXSize(),
2692 8 : papoOverviewBands[i]->GetYSize(), GetYSize());
2693 :
2694 8 : if (nReqOvLevel == nThisOvLevel)
2695 1 : papoOvBands[iOverview] = papoOverviewBands[i];
2696 : }
2697 :
2698 : // If this overview level does not yet exist, create it now.
2699 12 : if (papoOvBands[iOverview] == nullptr)
2700 : {
2701 22 : const int iResult = HFACreateOverview(
2702 11 : hHFA, nBand, panOverviewList[iOverview], pszResampling);
2703 11 : if (iResult < 0)
2704 : {
2705 0 : CPLFree(papoOvBands);
2706 0 : return CE_Failure;
2707 : }
2708 :
2709 11 : if (papoOverviewBands == nullptr && nOverviews == 0 && iResult > 0)
2710 : {
2711 0 : CPLDebug("HFA", "Shouldn't happen happened at line %d",
2712 : __LINE__);
2713 0 : papoOverviewBands = static_cast<HFARasterBand **>(
2714 0 : CPLCalloc(sizeof(void *), iResult));
2715 : }
2716 :
2717 11 : nOverviews = iResult + 1;
2718 11 : papoOverviewBands = static_cast<HFARasterBand **>(
2719 11 : CPLRealloc(papoOverviewBands, sizeof(void *) * nOverviews));
2720 11 : papoOverviewBands[iResult] = new HFARasterBand(
2721 11 : cpl::down_cast<HFADataset *>(poDS), nBand, iResult);
2722 :
2723 11 : papoOvBands[iOverview] = papoOverviewBands[iResult];
2724 : }
2725 : }
2726 :
2727 10 : CPLErr eErr = CE_None;
2728 :
2729 10 : if (bRegenerate)
2730 5 : eErr = GDALRegenerateOverviewsEx((GDALRasterBandH)this, nReqOverviews,
2731 : (GDALRasterBandH *)papoOvBands,
2732 : pszResampling, pfnProgress,
2733 : pProgressData, papszOptions);
2734 :
2735 10 : CPLFree(papoOvBands);
2736 :
2737 10 : return eErr;
2738 : }
2739 :
2740 : /************************************************************************/
2741 : /* GetDefaultHistogram() */
2742 : /************************************************************************/
2743 :
2744 12 : CPLErr HFARasterBand::GetDefaultHistogram(double *pdfMin, double *pdfMax,
2745 : int *pnBuckets,
2746 : GUIntBig **ppanHistogram, int bForce,
2747 : GDALProgressFunc pfnProgress,
2748 : void *pProgressData)
2749 :
2750 : {
2751 12 : if (GetMetadataItem("STATISTICS_HISTOBINVALUES") != nullptr &&
2752 20 : GetMetadataItem("STATISTICS_HISTOMIN") != nullptr &&
2753 8 : GetMetadataItem("STATISTICS_HISTOMAX") != nullptr)
2754 : {
2755 8 : const char *pszBinValues = GetMetadataItem("STATISTICS_HISTOBINVALUES");
2756 :
2757 8 : *pdfMin = CPLAtof(GetMetadataItem("STATISTICS_HISTOMIN"));
2758 8 : *pdfMax = CPLAtof(GetMetadataItem("STATISTICS_HISTOMAX"));
2759 :
2760 8 : *pnBuckets = 0;
2761 5402 : for (int i = 0; pszBinValues[i] != '\0'; i++)
2762 : {
2763 5394 : if (pszBinValues[i] == '|')
2764 1976 : (*pnBuckets)++;
2765 : }
2766 :
2767 8 : *ppanHistogram =
2768 8 : static_cast<GUIntBig *>(CPLCalloc(sizeof(GUIntBig), *pnBuckets));
2769 :
2770 8 : const char *pszNextBin = pszBinValues;
2771 1984 : for (int i = 0; i < *pnBuckets; i++)
2772 : {
2773 1976 : (*ppanHistogram)[i] =
2774 1976 : static_cast<GUIntBig>(CPLAtoGIntBig(pszNextBin));
2775 :
2776 5394 : while (*pszNextBin != '|' && *pszNextBin != '\0')
2777 3418 : pszNextBin++;
2778 1976 : if (*pszNextBin == '|')
2779 1976 : pszNextBin++;
2780 : }
2781 :
2782 : // Adjust min/max to reflect outer edges of buckets.
2783 8 : double dfBucketWidth = (*pdfMax - *pdfMin) / (*pnBuckets - 1);
2784 8 : *pdfMax += 0.5 * dfBucketWidth;
2785 8 : *pdfMin -= 0.5 * dfBucketWidth;
2786 :
2787 8 : return CE_None;
2788 : }
2789 :
2790 4 : return GDALPamRasterBand::GetDefaultHistogram(pdfMin, pdfMax, pnBuckets,
2791 : ppanHistogram, bForce,
2792 4 : pfnProgress, pProgressData);
2793 : }
2794 :
2795 : /************************************************************************/
2796 : /* SetDefaultRAT() */
2797 : /************************************************************************/
2798 :
2799 6 : CPLErr HFARasterBand::SetDefaultRAT(const GDALRasterAttributeTable *poRAT)
2800 :
2801 : {
2802 6 : if (poRAT == nullptr)
2803 0 : return CE_Failure;
2804 :
2805 6 : delete poDefaultRAT;
2806 6 : poDefaultRAT = nullptr;
2807 :
2808 6 : CPLErr r = WriteNamedRAT("Descriptor_Table", poRAT);
2809 6 : if (!r)
2810 6 : GetDefaultRAT();
2811 :
2812 6 : return r;
2813 : }
2814 :
2815 : /************************************************************************/
2816 : /* GetDefaultRAT() */
2817 : /************************************************************************/
2818 :
2819 1309 : GDALRasterAttributeTable *HFARasterBand::GetDefaultRAT()
2820 :
2821 : {
2822 1309 : if (poDefaultRAT == nullptr)
2823 622 : poDefaultRAT = new HFARasterAttributeTable(this, "Descriptor_Table");
2824 :
2825 1309 : return poDefaultRAT;
2826 : }
2827 :
2828 : /************************************************************************/
2829 : /* WriteNamedRAT() */
2830 : /************************************************************************/
2831 :
2832 6 : CPLErr HFARasterBand::WriteNamedRAT(const char * /*pszName*/,
2833 : const GDALRasterAttributeTable *poRAT)
2834 : {
2835 : // Find the requested table.
2836 : HFAEntry *poDT =
2837 6 : hHFA->papoBand[nBand - 1]->poNode->GetNamedChild("Descriptor_Table");
2838 6 : if (poDT == nullptr || !EQUAL(poDT->GetType(), "Edsc_Table"))
2839 : poDT =
2840 6 : HFAEntry::New(hHFA->papoBand[nBand - 1]->psInfo, "Descriptor_Table",
2841 6 : "Edsc_Table", hHFA->papoBand[nBand - 1]->poNode);
2842 :
2843 6 : const int nRowCount = poRAT->GetRowCount();
2844 :
2845 6 : poDT->SetIntField("numrows", nRowCount);
2846 : // Check if binning is set on this RAT.
2847 6 : double dfBinSize = 0.0;
2848 6 : double dfRow0Min = 0.0;
2849 6 : if (poRAT->GetLinearBinning(&dfRow0Min, &dfBinSize))
2850 : {
2851 : // Then it should have an Edsc_BinFunction.
2852 4 : HFAEntry *poBinFunction = poDT->GetNamedChild("#Bin_Function#");
2853 4 : if (poBinFunction == nullptr ||
2854 0 : !EQUAL(poBinFunction->GetType(), "Edsc_BinFunction"))
2855 : {
2856 : poBinFunction =
2857 4 : HFAEntry::New(hHFA->papoBand[nBand - 1]->psInfo,
2858 : "#Bin_Function#", "Edsc_BinFunction", poDT);
2859 : }
2860 :
2861 : // direct for thematic layers, linear otherwise
2862 : const char *pszLayerType =
2863 4 : hHFA->papoBand[nBand - 1]->poNode->GetStringField("layerType");
2864 4 : if (pszLayerType == nullptr || STARTS_WITH_CI(pszLayerType, "thematic"))
2865 0 : poBinFunction->SetStringField("binFunctionType", "direct");
2866 : else
2867 4 : poBinFunction->SetStringField("binFunctionType", "linear");
2868 :
2869 4 : poBinFunction->SetDoubleField("minLimit", dfRow0Min);
2870 4 : poBinFunction->SetDoubleField("maxLimit",
2871 4 : (nRowCount - 1) * dfBinSize + dfRow0Min);
2872 4 : poBinFunction->SetIntField("numBins", nRowCount);
2873 : }
2874 :
2875 : // Loop through each column in the RAT.
2876 18 : for (int col = 0; col < poRAT->GetColumnCount(); col++)
2877 : {
2878 12 : const char *pszName = nullptr;
2879 :
2880 12 : if (poRAT->GetUsageOfCol(col) == GFU_Red)
2881 : {
2882 1 : pszName = "Red";
2883 : }
2884 11 : else if (poRAT->GetUsageOfCol(col) == GFU_Green)
2885 : {
2886 1 : pszName = "Green";
2887 : }
2888 10 : else if (poRAT->GetUsageOfCol(col) == GFU_Blue)
2889 : {
2890 1 : pszName = "Blue";
2891 : }
2892 9 : else if (poRAT->GetUsageOfCol(col) == GFU_Alpha)
2893 : {
2894 1 : pszName = "Opacity";
2895 : }
2896 8 : else if (poRAT->GetUsageOfCol(col) == GFU_PixelCount)
2897 : {
2898 6 : pszName = "Histogram";
2899 : }
2900 2 : else if (poRAT->GetUsageOfCol(col) == GFU_Name)
2901 : {
2902 0 : pszName = "Class_Names";
2903 : }
2904 : else
2905 : {
2906 2 : pszName = poRAT->GetNameOfCol(col);
2907 : }
2908 :
2909 : // Check to see if a column with pszName exists and create if
2910 : // if necessary.
2911 12 : HFAEntry *poColumn = poDT->GetNamedChild(pszName);
2912 :
2913 12 : if (poColumn == nullptr || !EQUAL(poColumn->GetType(), "Edsc_Column"))
2914 12 : poColumn = HFAEntry::New(hHFA->papoBand[nBand - 1]->psInfo, pszName,
2915 : "Edsc_Column", poDT);
2916 :
2917 12 : poColumn->SetIntField("numRows", nRowCount);
2918 : // Color cols which are integer in GDAL are written as floats in HFA.
2919 12 : bool bIsColorCol = false;
2920 12 : if (poRAT->GetUsageOfCol(col) == GFU_Red ||
2921 11 : poRAT->GetUsageOfCol(col) == GFU_Green ||
2922 32 : poRAT->GetUsageOfCol(col) == GFU_Blue ||
2923 9 : poRAT->GetUsageOfCol(col) == GFU_Alpha)
2924 : {
2925 4 : bIsColorCol = true;
2926 : }
2927 :
2928 : // Write float also if a color column or histogram.
2929 12 : if (poRAT->GetTypeOfCol(col) == GFT_Real || bIsColorCol ||
2930 0 : poRAT->GetUsageOfCol(col) == GFU_PixelCount)
2931 : {
2932 : const int nOffset =
2933 24 : HFAAllocateSpace(hHFA->papoBand[nBand - 1]->psInfo,
2934 12 : static_cast<GUInt32>(nRowCount) *
2935 12 : static_cast<GUInt32>(sizeof(double)));
2936 12 : poColumn->SetIntField("columnDataPtr", nOffset);
2937 12 : poColumn->SetStringField("dataType", "real");
2938 :
2939 : double *padfColData =
2940 12 : static_cast<double *>(CPLCalloc(nRowCount, sizeof(double)));
2941 1716 : for (int i = 0; i < nRowCount; i++)
2942 : {
2943 1704 : if (bIsColorCol)
2944 : // Stored 0..1
2945 300 : padfColData[i] = poRAT->GetValueAsInt(i, col) / 255.0;
2946 : else
2947 1404 : padfColData[i] = poRAT->GetValueAsDouble(i, col);
2948 : }
2949 : #ifdef CPL_MSB
2950 : GDALSwapWords(padfColData, 8, nRowCount, 8);
2951 : #endif
2952 24 : if (VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) != 0 ||
2953 12 : VSIFWriteL(padfColData, nRowCount, sizeof(double), hHFA->fp) !=
2954 : sizeof(double))
2955 : {
2956 0 : CPLError(CE_Failure, CPLE_FileIO, "WriteNamedRAT() failed");
2957 0 : CPLFree(padfColData);
2958 0 : return CE_Failure;
2959 : }
2960 12 : CPLFree(padfColData);
2961 : }
2962 0 : else if (poRAT->GetTypeOfCol(col) == GFT_String)
2963 : {
2964 0 : unsigned int nMaxNumChars = 0;
2965 : // Find the length of the longest string.
2966 0 : for (int i = 0; i < nRowCount; i++)
2967 : {
2968 : // Include terminating byte.
2969 : const unsigned int nNumChars = static_cast<unsigned int>(
2970 0 : strlen(poRAT->GetValueAsString(i, col)) + 1);
2971 0 : if (nMaxNumChars < nNumChars)
2972 : {
2973 0 : nMaxNumChars = nNumChars;
2974 : }
2975 : }
2976 :
2977 : const int nOffset =
2978 0 : HFAAllocateSpace(hHFA->papoBand[nBand - 1]->psInfo,
2979 0 : (nRowCount + 1) * nMaxNumChars);
2980 0 : poColumn->SetIntField("columnDataPtr", nOffset);
2981 0 : poColumn->SetStringField("dataType", "string");
2982 0 : poColumn->SetIntField("maxNumChars", nMaxNumChars);
2983 :
2984 : char *pachColData =
2985 0 : static_cast<char *>(CPLCalloc(nRowCount + 1, nMaxNumChars));
2986 0 : for (int i = 0; i < nRowCount; i++)
2987 : {
2988 0 : strcpy(&pachColData[nMaxNumChars * i],
2989 0 : poRAT->GetValueAsString(i, col));
2990 : }
2991 0 : if (VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) != 0 ||
2992 0 : VSIFWriteL(pachColData, nRowCount, nMaxNumChars, hHFA->fp) !=
2993 0 : nMaxNumChars)
2994 : {
2995 0 : CPLError(CE_Failure, CPLE_FileIO, "WriteNamedRAT() failed");
2996 0 : CPLFree(pachColData);
2997 0 : return CE_Failure;
2998 : }
2999 0 : CPLFree(pachColData);
3000 : }
3001 0 : else if (poRAT->GetTypeOfCol(col) == GFT_Integer)
3002 : {
3003 0 : const int nOffset = HFAAllocateSpace(
3004 0 : hHFA->papoBand[nBand - 1]->psInfo,
3005 0 : static_cast<GUInt32>(nRowCount) * (GUInt32)sizeof(GInt32));
3006 0 : poColumn->SetIntField("columnDataPtr", nOffset);
3007 0 : poColumn->SetStringField("dataType", "integer");
3008 :
3009 : GInt32 *panColData =
3010 0 : static_cast<GInt32 *>(CPLCalloc(nRowCount, sizeof(GInt32)));
3011 0 : for (int i = 0; i < nRowCount; i++)
3012 : {
3013 0 : panColData[i] = poRAT->GetValueAsInt(i, col);
3014 : }
3015 : #ifdef CPL_MSB
3016 : GDALSwapWords(panColData, 4, nRowCount, 4);
3017 : #endif
3018 0 : if (VSIFSeekL(hHFA->fp, nOffset, SEEK_SET) != 0 ||
3019 0 : VSIFWriteL(panColData, nRowCount, sizeof(GInt32), hHFA->fp) !=
3020 : sizeof(GInt32))
3021 : {
3022 0 : CPLError(CE_Failure, CPLE_FileIO, "WriteNamedRAT() failed");
3023 0 : CPLFree(panColData);
3024 0 : return CE_Failure;
3025 : }
3026 0 : CPLFree(panColData);
3027 : }
3028 : else
3029 : {
3030 : // Can't deal with any of the others yet.
3031 0 : CPLError(CE_Failure, CPLE_NotSupported,
3032 : "Writing this data type in a column is not supported "
3033 : "for this Raster Attribute Table.");
3034 : }
3035 : }
3036 :
3037 6 : return CE_None;
3038 : }
3039 :
3040 : /************************************************************************/
3041 : /* ==================================================================== */
3042 : /* HFADataset */
3043 : /* ==================================================================== */
3044 : /************************************************************************/
3045 :
3046 : /************************************************************************/
3047 : /* HFADataset() */
3048 : /************************************************************************/
3049 :
3050 540 : HFADataset::HFADataset()
3051 : {
3052 540 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3053 540 : }
3054 :
3055 : /************************************************************************/
3056 : /* ~HFADataset() */
3057 : /************************************************************************/
3058 :
3059 1080 : HFADataset::~HFADataset()
3060 :
3061 : {
3062 540 : HFADataset::FlushCache(true);
3063 :
3064 : // Destroy the raster bands if they exist. We forcibly clean
3065 : // them up now to avoid any effort to write to them after the
3066 : // file is closed.
3067 1156 : for (int i = 0; i < nBands && papoBands != nullptr; i++)
3068 : {
3069 616 : if (papoBands[i] != nullptr)
3070 616 : delete papoBands[i];
3071 : }
3072 :
3073 540 : CPLFree(papoBands);
3074 540 : papoBands = nullptr;
3075 :
3076 : // Close the file.
3077 540 : if (hHFA != nullptr)
3078 : {
3079 540 : if (HFAClose(hHFA) != 0)
3080 : {
3081 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
3082 : }
3083 540 : hHFA = nullptr;
3084 : }
3085 1080 : }
3086 :
3087 : /************************************************************************/
3088 : /* FlushCache() */
3089 : /************************************************************************/
3090 :
3091 572 : CPLErr HFADataset::FlushCache(bool bAtClosing)
3092 :
3093 : {
3094 572 : CPLErr eErr = GDALPamDataset::FlushCache(bAtClosing);
3095 :
3096 572 : if (eAccess != GA_Update)
3097 354 : return eErr;
3098 :
3099 218 : if (bGeoDirty)
3100 146 : WriteProjection();
3101 :
3102 218 : if (bMetadataDirty && GetMetadata() != nullptr)
3103 : {
3104 50 : HFASetMetadata(hHFA, 0, GetMetadata());
3105 50 : bMetadataDirty = false;
3106 : }
3107 :
3108 472 : for (int iBand = 0; iBand < nBands; iBand++)
3109 : {
3110 : HFARasterBand *poBand =
3111 254 : static_cast<HFARasterBand *>(GetRasterBand(iBand + 1));
3112 254 : if (poBand->bMetadataDirty && poBand->GetMetadata() != nullptr)
3113 : {
3114 14 : HFASetMetadata(hHFA, iBand + 1, poBand->GetMetadata());
3115 14 : poBand->bMetadataDirty = false;
3116 : }
3117 : }
3118 :
3119 218 : return eErr;
3120 : }
3121 :
3122 : /************************************************************************/
3123 : /* WriteProjection() */
3124 : /************************************************************************/
3125 :
3126 146 : CPLErr HFADataset::WriteProjection()
3127 :
3128 : {
3129 146 : bool bPEStringStored = false;
3130 :
3131 146 : bGeoDirty = false;
3132 :
3133 146 : const OGRSpatialReference &oSRS = m_oSRS;
3134 146 : const bool bHaveSRS = !oSRS.IsEmpty();
3135 :
3136 : // Initialize projection and datum.
3137 : Eprj_Datum sDatum;
3138 : Eprj_ProParameters sPro;
3139 : Eprj_MapInfo sMapInfo;
3140 146 : memset(&sPro, 0, sizeof(sPro));
3141 146 : memset(&sDatum, 0, sizeof(sDatum));
3142 146 : memset(&sMapInfo, 0, sizeof(sMapInfo));
3143 :
3144 : // Collect datum information.
3145 146 : OGRSpatialReference *poGeogSRS = bHaveSRS ? oSRS.CloneGeogCS() : nullptr;
3146 :
3147 146 : if (poGeogSRS)
3148 : {
3149 133 : sDatum.datumname =
3150 133 : const_cast<char *>(poGeogSRS->GetAttrValue("GEOGCS|DATUM"));
3151 133 : if (sDatum.datumname == nullptr)
3152 0 : sDatum.datumname = const_cast<char *>("");
3153 :
3154 : // WKT to Imagine translation.
3155 133 : const char *const *papszDatumMap = HFAGetDatumMap();
3156 563 : for (int i = 0; papszDatumMap[i] != nullptr; i += 2)
3157 : {
3158 518 : if (EQUAL(sDatum.datumname, papszDatumMap[i + 1]))
3159 : {
3160 88 : sDatum.datumname = (char *)papszDatumMap[i];
3161 88 : break;
3162 : }
3163 : }
3164 :
3165 : // Map some EPSG datum codes directly to Imagine names.
3166 133 : const int nGCS = poGeogSRS->GetEPSGGeogCS();
3167 :
3168 133 : if (nGCS == 4326)
3169 7 : sDatum.datumname = const_cast<char *>("WGS 84");
3170 126 : else if (nGCS == 4322)
3171 0 : sDatum.datumname = const_cast<char *>("WGS 1972");
3172 126 : else if (nGCS == 4267)
3173 30 : sDatum.datumname = const_cast<char *>("NAD27");
3174 96 : else if (nGCS == 4269)
3175 2 : sDatum.datumname = const_cast<char *>("NAD83");
3176 94 : else if (nGCS == 4283)
3177 0 : sDatum.datumname = const_cast<char *>("GDA94");
3178 94 : else if (nGCS == 4284)
3179 0 : sDatum.datumname = const_cast<char *>("Pulkovo 1942");
3180 94 : else if (nGCS == 4272)
3181 1 : sDatum.datumname = const_cast<char *>("Geodetic Datum 1949");
3182 :
3183 133 : if (poGeogSRS->GetTOWGS84(sDatum.params) == OGRERR_NONE)
3184 : {
3185 2 : sDatum.type = EPRJ_DATUM_PARAMETRIC;
3186 2 : sDatum.params[3] *= -ARCSEC2RAD;
3187 2 : sDatum.params[4] *= -ARCSEC2RAD;
3188 2 : sDatum.params[5] *= -ARCSEC2RAD;
3189 2 : sDatum.params[6] *= 1e-6;
3190 : }
3191 131 : else if (EQUAL(sDatum.datumname, "NAD27"))
3192 : {
3193 30 : sDatum.type = EPRJ_DATUM_GRID;
3194 30 : sDatum.gridname = const_cast<char *>("nadcon.dat");
3195 : }
3196 : else
3197 : {
3198 : // We will default to this (effectively WGS84) for now.
3199 101 : sDatum.type = EPRJ_DATUM_PARAMETRIC;
3200 : }
3201 :
3202 : // Verify if we need to write a ESRI PE string.
3203 133 : if (!bDisablePEString)
3204 132 : bPEStringStored = CPL_TO_BOOL(WritePeStringIfNeeded(&oSRS, hHFA));
3205 :
3206 133 : sPro.proSpheroid.sphereName =
3207 133 : (char *)poGeogSRS->GetAttrValue("GEOGCS|DATUM|SPHEROID");
3208 133 : sPro.proSpheroid.a = poGeogSRS->GetSemiMajor();
3209 133 : sPro.proSpheroid.b = poGeogSRS->GetSemiMinor();
3210 133 : sPro.proSpheroid.radius = sPro.proSpheroid.a;
3211 :
3212 133 : const double a2 = sPro.proSpheroid.a * sPro.proSpheroid.a;
3213 133 : const double b2 = sPro.proSpheroid.b * sPro.proSpheroid.b;
3214 :
3215 : // a2 == 0 is non sensical of course. Just to please fuzzers
3216 133 : sPro.proSpheroid.eSquared = (a2 == 0.0) ? 0.0 : (a2 - b2) / a2;
3217 : }
3218 :
3219 146 : if (sDatum.datumname == nullptr)
3220 13 : sDatum.datumname = const_cast<char *>("");
3221 146 : if (sPro.proSpheroid.sphereName == nullptr)
3222 13 : sPro.proSpheroid.sphereName = const_cast<char *>("");
3223 :
3224 : // Recognise various projections.
3225 146 : const char *pszProjName = nullptr;
3226 :
3227 146 : if (bHaveSRS)
3228 133 : pszProjName = oSRS.GetAttrValue("PROJCS|PROJECTION");
3229 :
3230 146 : if (bForceToPEString && !bPEStringStored)
3231 : {
3232 0 : char *pszPEString = nullptr;
3233 0 : const char *const apszOptions[] = {"FORMAT=WKT1_ESRI", nullptr};
3234 0 : oSRS.exportToWkt(&pszPEString, apszOptions);
3235 : // Need to transform this into ESRI format.
3236 0 : HFASetPEString(hHFA, pszPEString);
3237 0 : CPLFree(pszPEString);
3238 :
3239 0 : bPEStringStored = true;
3240 : }
3241 146 : else if (pszProjName == nullptr)
3242 : {
3243 54 : if (bHaveSRS && oSRS.IsGeographic())
3244 : {
3245 41 : sPro.proNumber = EPRJ_LATLONG;
3246 41 : sPro.proName = const_cast<char *>("Geographic (Lat/Lon)");
3247 : }
3248 : }
3249 : // TODO: Add State Plane.
3250 92 : else if (!bIgnoreUTM && oSRS.GetUTMZone(nullptr) != 0)
3251 : {
3252 32 : int bNorth = FALSE;
3253 32 : const int nZone = oSRS.GetUTMZone(&bNorth);
3254 32 : sPro.proNumber = EPRJ_UTM;
3255 32 : sPro.proName = const_cast<char *>("UTM");
3256 32 : sPro.proZone = nZone;
3257 32 : if (bNorth)
3258 32 : sPro.proParams[3] = 1.0;
3259 : else
3260 0 : sPro.proParams[3] = -1.0;
3261 : }
3262 60 : else if (EQUAL(pszProjName, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
3263 : {
3264 1 : sPro.proNumber = EPRJ_ALBERS_CONIC_EQUAL_AREA;
3265 1 : sPro.proName = const_cast<char *>("Albers Conical Equal Area");
3266 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3267 1 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_2) * D2R;
3268 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3269 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3270 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3271 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3272 : }
3273 59 : else if (EQUAL(pszProjName, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
3274 : {
3275 : // Not sure if Imagine has a mapping of LCC_1SP. In the mean time
3276 : // convert it to LCC_2SP
3277 : auto poTmpSRS = std::unique_ptr<OGRSpatialReference>(
3278 2 : oSRS.convertToOtherProjection(SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP));
3279 1 : if (poTmpSRS)
3280 : {
3281 1 : sPro.proNumber = EPRJ_LAMBERT_CONFORMAL_CONIC;
3282 1 : sPro.proName = const_cast<char *>("Lambert Conformal Conic");
3283 1 : sPro.proParams[2] =
3284 1 : poTmpSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3285 1 : sPro.proParams[3] =
3286 1 : poTmpSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_2) * D2R;
3287 1 : sPro.proParams[4] =
3288 1 : poTmpSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3289 1 : sPro.proParams[5] =
3290 1 : poTmpSRS->GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3291 1 : sPro.proParams[6] = poTmpSRS->GetProjParm(SRS_PP_FALSE_EASTING);
3292 1 : sPro.proParams[7] = poTmpSRS->GetProjParm(SRS_PP_FALSE_NORTHING);
3293 : }
3294 : }
3295 58 : else if (EQUAL(pszProjName, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
3296 : {
3297 3 : sPro.proNumber = EPRJ_LAMBERT_CONFORMAL_CONIC;
3298 3 : sPro.proName = const_cast<char *>("Lambert Conformal Conic");
3299 3 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3300 3 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_2) * D2R;
3301 3 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3302 3 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3303 3 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3304 3 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3305 : }
3306 57 : else if (EQUAL(pszProjName, SRS_PT_MERCATOR_1SP) &&
3307 2 : oSRS.GetProjParm(SRS_PP_SCALE_FACTOR) == 1.0)
3308 : {
3309 1 : sPro.proNumber = EPRJ_MERCATOR;
3310 1 : sPro.proName = const_cast<char *>("Mercator");
3311 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3312 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3313 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3314 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3315 : }
3316 54 : else if (EQUAL(pszProjName, SRS_PT_MERCATOR_1SP))
3317 : {
3318 1 : sPro.proNumber = EPRJ_MERCATOR_VARIANT_A;
3319 1 : sPro.proName = const_cast<char *>("Mercator (Variant A)");
3320 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3321 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3322 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR);
3323 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3324 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3325 : }
3326 53 : else if (EQUAL(pszProjName, SRS_PT_MERCATOR_2SP))
3327 : {
3328 : // Not sure if Imagine has a mapping of Mercator_2SP. In the mean time
3329 : // convert it to Mercator_1SP
3330 : auto poTmpSRS = std::unique_ptr<OGRSpatialReference>(
3331 2 : oSRS.convertToOtherProjection(SRS_PT_MERCATOR_1SP));
3332 1 : if (poTmpSRS)
3333 : {
3334 1 : sPro.proNumber = EPRJ_MERCATOR_VARIANT_A;
3335 1 : sPro.proName = const_cast<char *>("Mercator (Variant A)");
3336 1 : sPro.proParams[4] =
3337 1 : poTmpSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3338 1 : sPro.proParams[5] =
3339 1 : poTmpSRS->GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3340 1 : sPro.proParams[2] = poTmpSRS->GetProjParm(SRS_PP_SCALE_FACTOR);
3341 1 : sPro.proParams[6] = poTmpSRS->GetProjParm(SRS_PP_FALSE_EASTING);
3342 1 : sPro.proParams[7] = poTmpSRS->GetProjParm(SRS_PP_FALSE_NORTHING);
3343 : }
3344 : }
3345 52 : else if (EQUAL(pszProjName, SRS_PT_KROVAK))
3346 : {
3347 1 : sPro.proNumber = EPRJ_KROVAK;
3348 1 : sPro.proName = const_cast<char *>("Krovak");
3349 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR);
3350 1 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_AZIMUTH) * D2R;
3351 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3352 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3353 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3354 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3355 1 : sPro.proParams[9] = oSRS.GetProjParm(SRS_PP_PSEUDO_STD_PARALLEL_1);
3356 :
3357 1 : sPro.proParams[8] = 0.0; // XY plane rotation
3358 1 : sPro.proParams[10] = 1.0; // X scale
3359 1 : sPro.proParams[11] = 1.0; // Y scale
3360 : }
3361 51 : else if (EQUAL(pszProjName, SRS_PT_POLAR_STEREOGRAPHIC))
3362 : {
3363 2 : sPro.proNumber = EPRJ_POLAR_STEREOGRAPHIC;
3364 2 : sPro.proName = const_cast<char *>("Polar Stereographic");
3365 2 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3366 2 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3367 : // Hopefully the scale factor is 1.0!
3368 2 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3369 2 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3370 : }
3371 49 : else if (EQUAL(pszProjName, SRS_PT_POLYCONIC))
3372 : {
3373 2 : sPro.proNumber = EPRJ_POLYCONIC;
3374 2 : sPro.proName = const_cast<char *>("Polyconic");
3375 2 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3376 2 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3377 2 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3378 2 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3379 : }
3380 47 : else if (EQUAL(pszProjName, SRS_PT_EQUIDISTANT_CONIC))
3381 : {
3382 1 : sPro.proNumber = EPRJ_EQUIDISTANT_CONIC;
3383 1 : sPro.proName = const_cast<char *>("Equidistant Conic");
3384 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3385 1 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_2) * D2R;
3386 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3387 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3388 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3389 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3390 1 : sPro.proParams[8] = 1.0;
3391 : }
3392 46 : else if (EQUAL(pszProjName, SRS_PT_TRANSVERSE_MERCATOR))
3393 : {
3394 1 : sPro.proNumber = EPRJ_TRANSVERSE_MERCATOR;
3395 1 : sPro.proName = const_cast<char *>("Transverse Mercator");
3396 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3397 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3398 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3399 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3400 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3401 : }
3402 45 : else if (EQUAL(pszProjName, SRS_PT_STEREOGRAPHIC))
3403 : {
3404 0 : sPro.proNumber = EPRJ_STEREOGRAPHIC_EXTENDED;
3405 0 : sPro.proName = const_cast<char *>("Stereographic (Extended)");
3406 0 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3407 0 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3408 0 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3409 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3410 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3411 : }
3412 45 : else if (EQUAL(pszProjName, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
3413 : {
3414 1 : sPro.proNumber = EPRJ_LAMBERT_AZIMUTHAL_EQUAL_AREA;
3415 1 : sPro.proName = const_cast<char *>("Lambert Azimuthal Equal-area");
3416 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3417 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3418 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3419 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3420 : }
3421 44 : else if (EQUAL(pszProjName, SRS_PT_AZIMUTHAL_EQUIDISTANT))
3422 : {
3423 1 : sPro.proNumber = EPRJ_AZIMUTHAL_EQUIDISTANT;
3424 1 : sPro.proName = const_cast<char *>("Azimuthal Equidistant");
3425 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3426 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3427 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3428 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3429 : }
3430 43 : else if (EQUAL(pszProjName, SRS_PT_GNOMONIC))
3431 : {
3432 1 : sPro.proNumber = EPRJ_GNOMONIC;
3433 1 : sPro.proName = const_cast<char *>("Gnomonic");
3434 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3435 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3436 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3437 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3438 : }
3439 42 : else if (EQUAL(pszProjName, SRS_PT_ORTHOGRAPHIC))
3440 : {
3441 2 : sPro.proNumber = EPRJ_ORTHOGRAPHIC;
3442 2 : sPro.proName = const_cast<char *>("Orthographic");
3443 2 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3444 2 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3445 2 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3446 2 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3447 : }
3448 40 : else if (EQUAL(pszProjName, SRS_PT_SINUSOIDAL))
3449 : {
3450 1 : sPro.proNumber = EPRJ_SINUSOIDAL;
3451 1 : sPro.proName = const_cast<char *>("Sinusoidal");
3452 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3453 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3454 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3455 : }
3456 39 : else if (EQUAL(pszProjName, SRS_PT_EQUIRECTANGULAR))
3457 : {
3458 3 : sPro.proNumber = EPRJ_EQUIRECTANGULAR;
3459 3 : sPro.proName = const_cast<char *>("Equirectangular");
3460 3 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3461 3 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3462 3 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3463 3 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3464 : }
3465 36 : else if (EQUAL(pszProjName, SRS_PT_MILLER_CYLINDRICAL))
3466 : {
3467 1 : sPro.proNumber = EPRJ_MILLER_CYLINDRICAL;
3468 1 : sPro.proName = const_cast<char *>("Miller Cylindrical");
3469 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3470 : // Hopefully the latitude is zero!
3471 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3472 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3473 : }
3474 35 : else if (EQUAL(pszProjName, SRS_PT_VANDERGRINTEN))
3475 : {
3476 1 : sPro.proNumber = EPRJ_VANDERGRINTEN;
3477 1 : sPro.proName = const_cast<char *>("Van der Grinten");
3478 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3479 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3480 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3481 : }
3482 34 : else if (EQUAL(pszProjName, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
3483 : {
3484 2 : if (oSRS.GetProjParm(SRS_PP_RECTIFIED_GRID_ANGLE) == 0.0)
3485 : {
3486 0 : sPro.proNumber = EPRJ_HOTINE_OBLIQUE_MERCATOR;
3487 0 : sPro.proName = const_cast<char *>("Oblique Mercator (Hotine)");
3488 0 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3489 0 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_AZIMUTH) * D2R;
3490 0 : sPro.proParams[4] =
3491 0 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3492 0 : sPro.proParams[5] =
3493 0 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3494 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3495 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3496 0 : sPro.proParams[12] = 1.0;
3497 : }
3498 : else
3499 : {
3500 2 : sPro.proNumber = EPRJ_HOTINE_OBLIQUE_MERCATOR_VARIANT_A;
3501 2 : sPro.proName =
3502 : const_cast<char *>("Hotine Oblique Mercator (Variant A)");
3503 2 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3504 2 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_AZIMUTH) * D2R;
3505 2 : sPro.proParams[4] =
3506 2 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3507 2 : sPro.proParams[5] =
3508 2 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3509 2 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3510 2 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3511 2 : sPro.proParams[8] =
3512 2 : oSRS.GetProjParm(SRS_PP_RECTIFIED_GRID_ANGLE) * D2R;
3513 : }
3514 : }
3515 32 : else if (EQUAL(pszProjName, SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
3516 : {
3517 2 : sPro.proNumber = EPRJ_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER;
3518 2 : sPro.proName =
3519 : const_cast<char *>("Hotine Oblique Mercator Azimuth Center");
3520 2 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3521 2 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_AZIMUTH) * D2R;
3522 2 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3523 2 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3524 2 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3525 2 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3526 2 : sPro.proParams[12] = 1.0;
3527 : }
3528 30 : else if (EQUAL(pszProjName, SRS_PT_ROBINSON))
3529 : {
3530 1 : sPro.proNumber = EPRJ_ROBINSON;
3531 1 : sPro.proName = const_cast<char *>("Robinson");
3532 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3533 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3534 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3535 : }
3536 29 : else if (EQUAL(pszProjName, SRS_PT_MOLLWEIDE))
3537 : {
3538 1 : sPro.proNumber = EPRJ_MOLLWEIDE;
3539 1 : sPro.proName = const_cast<char *>("Mollweide");
3540 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3541 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3542 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3543 : }
3544 28 : else if (EQUAL(pszProjName, SRS_PT_ECKERT_I))
3545 : {
3546 1 : sPro.proNumber = EPRJ_ECKERT_I;
3547 1 : sPro.proName = const_cast<char *>("Eckert I");
3548 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3549 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3550 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3551 : }
3552 27 : else if (EQUAL(pszProjName, SRS_PT_ECKERT_II))
3553 : {
3554 1 : sPro.proNumber = EPRJ_ECKERT_II;
3555 1 : sPro.proName = const_cast<char *>("Eckert II");
3556 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3557 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3558 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3559 : }
3560 26 : else if (EQUAL(pszProjName, SRS_PT_ECKERT_III))
3561 : {
3562 1 : sPro.proNumber = EPRJ_ECKERT_III;
3563 1 : sPro.proName = const_cast<char *>("Eckert III");
3564 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3565 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3566 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3567 : }
3568 25 : else if (EQUAL(pszProjName, SRS_PT_ECKERT_IV))
3569 : {
3570 1 : sPro.proNumber = EPRJ_ECKERT_IV;
3571 1 : sPro.proName = const_cast<char *>("Eckert IV");
3572 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3573 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3574 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3575 : }
3576 24 : else if (EQUAL(pszProjName, SRS_PT_ECKERT_V))
3577 : {
3578 1 : sPro.proNumber = EPRJ_ECKERT_V;
3579 1 : sPro.proName = const_cast<char *>("Eckert V");
3580 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3581 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3582 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3583 : }
3584 23 : else if (EQUAL(pszProjName, SRS_PT_ECKERT_VI))
3585 : {
3586 1 : sPro.proNumber = EPRJ_ECKERT_VI;
3587 1 : sPro.proName = const_cast<char *>("Eckert VI");
3588 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3589 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3590 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3591 : }
3592 22 : else if (EQUAL(pszProjName, SRS_PT_GALL_STEREOGRAPHIC))
3593 : {
3594 1 : sPro.proNumber = EPRJ_GALL_STEREOGRAPHIC;
3595 1 : sPro.proName = const_cast<char *>("Gall Stereographic");
3596 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3597 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3598 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3599 : }
3600 21 : else if (EQUAL(pszProjName, SRS_PT_CASSINI_SOLDNER))
3601 : {
3602 2 : sPro.proNumber = EPRJ_CASSINI;
3603 2 : sPro.proName = const_cast<char *>("Cassini");
3604 2 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3605 2 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3606 2 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3607 2 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3608 : }
3609 19 : else if (EQUAL(pszProjName, SRS_PT_TWO_POINT_EQUIDISTANT))
3610 : {
3611 1 : sPro.proNumber = EPRJ_TWO_POINT_EQUIDISTANT;
3612 1 : sPro.proName = const_cast<char *>("Two_Point_Equidistant");
3613 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3614 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3615 1 : sPro.proParams[8] =
3616 1 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_POINT_1, 0.0) * D2R;
3617 1 : sPro.proParams[9] =
3618 1 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_POINT_1, 0.0) * D2R;
3619 1 : sPro.proParams[10] =
3620 1 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_POINT_2, 60.0) * D2R;
3621 1 : sPro.proParams[11] =
3622 1 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_POINT_2, 60.0) * D2R;
3623 : }
3624 18 : else if (EQUAL(pszProjName, SRS_PT_BONNE))
3625 : {
3626 1 : sPro.proNumber = EPRJ_BONNE;
3627 1 : sPro.proName = const_cast<char *>("Bonne");
3628 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3629 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3630 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3631 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3632 : }
3633 17 : else if (EQUAL(pszProjName, "Loximuthal"))
3634 : {
3635 1 : sPro.proNumber = EPRJ_LOXIMUTHAL;
3636 1 : sPro.proName = const_cast<char *>("Loximuthal");
3637 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3638 1 : sPro.proParams[5] = oSRS.GetProjParm("latitude_of_origin") * D2R;
3639 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3640 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3641 : }
3642 16 : else if (EQUAL(pszProjName, "Quartic_Authalic"))
3643 : {
3644 1 : sPro.proNumber = EPRJ_QUARTIC_AUTHALIC;
3645 1 : sPro.proName = const_cast<char *>("Quartic Authalic");
3646 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3647 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3648 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3649 : }
3650 15 : else if (EQUAL(pszProjName, "Winkel_I"))
3651 : {
3652 1 : sPro.proNumber = EPRJ_WINKEL_I;
3653 1 : sPro.proName = const_cast<char *>("Winkel I");
3654 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3655 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3656 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3657 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3658 : }
3659 14 : else if (EQUAL(pszProjName, "Winkel_II"))
3660 : {
3661 1 : sPro.proNumber = EPRJ_WINKEL_II;
3662 1 : sPro.proName = const_cast<char *>("Winkel II");
3663 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3664 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3665 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3666 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3667 : }
3668 13 : else if (EQUAL(pszProjName, "Behrmann"))
3669 : {
3670 : // Mapped to Lambert Cylindrical Equal Area in recent PROJ versions
3671 0 : sPro.proNumber = EPRJ_BEHRMANN;
3672 0 : sPro.proName = const_cast<char *>("Behrmann");
3673 0 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3674 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3675 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3676 : }
3677 13 : else if (EQUAL(pszProjName, "Equidistant_Cylindrical"))
3678 : {
3679 : // Dead code path. Mapped to Equirectangular
3680 0 : sPro.proNumber = EPRJ_EQUIDISTANT_CYLINDRICAL;
3681 0 : sPro.proName = const_cast<char *>("Equidistant_Cylindrical");
3682 0 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3683 0 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3684 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3685 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3686 : }
3687 13 : else if (EQUAL(pszProjName, SRS_PT_KROVAK))
3688 : {
3689 0 : sPro.proNumber = EPRJ_KROVAK;
3690 0 : sPro.proName = const_cast<char *>("Krovak");
3691 0 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3692 0 : sPro.proParams[3] = oSRS.GetProjParm(SRS_PP_AZIMUTH) * D2R;
3693 0 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER) * D2R;
3694 0 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER) * D2R;
3695 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3696 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3697 0 : sPro.proParams[8] = oSRS.GetProjParm("XY_Plane_Rotation", 0.0) * D2R;
3698 0 : sPro.proParams[9] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3699 0 : sPro.proParams[10] = oSRS.GetProjParm("X_Scale", 1.0);
3700 0 : sPro.proParams[11] = oSRS.GetProjParm("Y_Scale", 1.0);
3701 : }
3702 13 : else if (EQUAL(pszProjName, "Double_Stereographic"))
3703 : {
3704 0 : sPro.proNumber = EPRJ_DOUBLE_STEREOGRAPHIC;
3705 0 : sPro.proName = const_cast<char *>("Double_Stereographic");
3706 0 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3707 0 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3708 0 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3709 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3710 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3711 : }
3712 13 : else if (EQUAL(pszProjName, "Aitoff"))
3713 : {
3714 1 : sPro.proNumber = EPRJ_AITOFF;
3715 1 : sPro.proName = const_cast<char *>("Aitoff");
3716 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3717 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3718 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3719 : }
3720 12 : else if (EQUAL(pszProjName, "Craster_Parabolic"))
3721 : {
3722 1 : sPro.proNumber = EPRJ_CRASTER_PARABOLIC;
3723 1 : sPro.proName = const_cast<char *>("Craster_Parabolic");
3724 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3725 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3726 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3727 : }
3728 11 : else if (EQUAL(pszProjName, SRS_PT_CYLINDRICAL_EQUAL_AREA))
3729 : {
3730 1 : sPro.proNumber = EPRJ_CYLINDRICAL_EQUAL_AREA;
3731 1 : sPro.proName = const_cast<char *>("Cylindrical_Equal_Area");
3732 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3733 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3734 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3735 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3736 : }
3737 10 : else if (EQUAL(pszProjName, "Flat_Polar_Quartic"))
3738 : {
3739 1 : sPro.proNumber = EPRJ_FLAT_POLAR_QUARTIC;
3740 1 : sPro.proName = const_cast<char *>("Flat_Polar_Quartic");
3741 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3742 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3743 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3744 : }
3745 9 : else if (EQUAL(pszProjName, "Times"))
3746 : {
3747 1 : sPro.proNumber = EPRJ_TIMES;
3748 1 : sPro.proName = const_cast<char *>("Times");
3749 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3750 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3751 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3752 : }
3753 8 : else if (EQUAL(pszProjName, "Winkel_Tripel"))
3754 : {
3755 1 : sPro.proNumber = EPRJ_WINKEL_TRIPEL;
3756 1 : sPro.proName = const_cast<char *>("Winkel_Tripel");
3757 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_STANDARD_PARALLEL_1) * D2R;
3758 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3759 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3760 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3761 : }
3762 7 : else if (EQUAL(pszProjName, "Hammer_Aitoff"))
3763 : {
3764 1 : sPro.proNumber = EPRJ_HAMMER_AITOFF;
3765 1 : sPro.proName = const_cast<char *>("Hammer_Aitoff");
3766 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3767 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3768 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3769 : }
3770 6 : else if (
3771 6 : EQUAL(pszProjName,
3772 : "Vertical_Near_Side_Perspective")) // ESRI WKT, before PROJ 6.3.0
3773 : {
3774 1 : sPro.proNumber = EPRJ_VERTICAL_NEAR_SIDE_PERSPECTIVE;
3775 1 : sPro.proName = const_cast<char *>("Vertical_Near_Side_Perspective");
3776 1 : sPro.proParams[2] = oSRS.GetProjParm("Height");
3777 1 : sPro.proParams[4] =
3778 1 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_CENTER, 75.0) * D2R;
3779 1 : sPro.proParams[5] =
3780 1 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER, 40.0) * D2R;
3781 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3782 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3783 : }
3784 5 : else if (EQUAL(pszProjName,
3785 : "Vertical Perspective")) // WKT2, starting with PROJ 6.3.0
3786 : {
3787 0 : sPro.proNumber = EPRJ_VERTICAL_NEAR_SIDE_PERSPECTIVE;
3788 0 : sPro.proName = const_cast<char *>("Vertical_Near_Side_Perspective");
3789 0 : sPro.proParams[2] = oSRS.GetProjParm("Viewpoint height");
3790 0 : sPro.proParams[4] =
3791 0 : oSRS.GetProjParm("Longitude of topocentric origin", 75.0) * D2R;
3792 0 : sPro.proParams[5] =
3793 0 : oSRS.GetProjParm("Latitude of topocentric origin", 40.0) * D2R;
3794 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3795 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3796 : }
3797 5 : else if (EQUAL(pszProjName, "Hotine_Oblique_Mercator_Two_Point_Center"))
3798 : {
3799 0 : sPro.proNumber = EPRJ_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_CENTER;
3800 0 : sPro.proName =
3801 : const_cast<char *>("Hotine_Oblique_Mercator_Two_Point_Center");
3802 0 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3803 0 : sPro.proParams[5] =
3804 0 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER, 40.0) * D2R;
3805 0 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3806 0 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3807 0 : sPro.proParams[8] =
3808 0 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_POINT_1, 0.0) * D2R;
3809 0 : sPro.proParams[9] =
3810 0 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_POINT_1, 0.0) * D2R;
3811 0 : sPro.proParams[10] =
3812 0 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_POINT_2, 60.0) * D2R;
3813 0 : sPro.proParams[11] =
3814 0 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_POINT_2, 60.0) * D2R;
3815 : }
3816 5 : else if (EQUAL(pszProjName,
3817 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN))
3818 : {
3819 1 : sPro.proNumber = EPRJ_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN;
3820 1 : sPro.proName = const_cast<char *>(
3821 : "Hotine_Oblique_Mercator_Two_Point_Natural_Origin");
3822 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3823 1 : sPro.proParams[5] =
3824 1 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_CENTER, 40.0) * D2R;
3825 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3826 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3827 1 : sPro.proParams[8] =
3828 1 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_POINT_1, 0.0) * D2R;
3829 1 : sPro.proParams[9] =
3830 1 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_POINT_1, 0.0) * D2R;
3831 1 : sPro.proParams[10] =
3832 1 : oSRS.GetProjParm(SRS_PP_LONGITUDE_OF_POINT_2, 60.0) * D2R;
3833 1 : sPro.proParams[11] =
3834 1 : oSRS.GetProjParm(SRS_PP_LATITUDE_OF_POINT_2, 60.0) * D2R;
3835 : }
3836 4 : else if (EQUAL(pszProjName, "New_Zealand_Map_Grid"))
3837 : {
3838 1 : sPro.proType = EPRJ_EXTERNAL;
3839 1 : sPro.proNumber = 0;
3840 1 : sPro.proExeName = const_cast<char *>(EPRJ_EXTERNAL_NZMG);
3841 1 : sPro.proName = const_cast<char *>("New Zealand Map Grid");
3842 1 : sPro.proZone = 0;
3843 1 : sPro.proParams[0] = 0; // False easting etc not stored in .img it seems
3844 1 : sPro.proParams[1] = 0; // always fixed by definition.
3845 1 : sPro.proParams[2] = 0;
3846 1 : sPro.proParams[3] = 0;
3847 1 : sPro.proParams[4] = 0;
3848 1 : sPro.proParams[5] = 0;
3849 1 : sPro.proParams[6] = 0;
3850 1 : sPro.proParams[7] = 0;
3851 : }
3852 3 : else if (EQUAL(pszProjName, SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED))
3853 : {
3854 1 : sPro.proNumber = EPRJ_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED;
3855 1 : sPro.proName =
3856 : const_cast<char *>("Transverse Mercator (South Orientated)");
3857 1 : sPro.proParams[4] = oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN) * D2R;
3858 1 : sPro.proParams[5] = oSRS.GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN) * D2R;
3859 1 : sPro.proParams[2] = oSRS.GetProjParm(SRS_PP_SCALE_FACTOR, 1.0);
3860 1 : sPro.proParams[6] = oSRS.GetProjParm(SRS_PP_FALSE_EASTING);
3861 1 : sPro.proParams[7] = oSRS.GetProjParm(SRS_PP_FALSE_NORTHING);
3862 : }
3863 :
3864 : // Anything we can't map, we store as an ESRI PE_STRING.
3865 2 : else if (oSRS.IsProjected() || oSRS.IsGeographic())
3866 : {
3867 2 : if (!bPEStringStored)
3868 : {
3869 0 : char *pszPEString = nullptr;
3870 0 : const char *const apszOptions[] = {"FORMAT=WKT1_ESRI", nullptr};
3871 0 : oSRS.exportToWkt(&pszPEString, apszOptions);
3872 : // Need to transform this into ESRI format.
3873 0 : HFASetPEString(hHFA, pszPEString);
3874 0 : CPLFree(pszPEString);
3875 0 : bPEStringStored = true;
3876 : }
3877 : }
3878 : else
3879 : {
3880 0 : CPLError(CE_Warning, CPLE_NotSupported,
3881 : "Projection %s not supported for translation to Imagine.",
3882 : pszProjName);
3883 : }
3884 :
3885 : // MapInfo
3886 146 : const char *pszPROJCS = oSRS.GetAttrValue("PROJCS");
3887 :
3888 146 : if (pszPROJCS)
3889 92 : sMapInfo.proName = (char *)pszPROJCS;
3890 54 : else if (bHaveSRS && sPro.proName != nullptr)
3891 41 : sMapInfo.proName = sPro.proName;
3892 : else
3893 13 : sMapInfo.proName = const_cast<char *>("Unknown");
3894 :
3895 146 : sMapInfo.upperLeftCenter.x = m_gt[0] + m_gt[1] * 0.5;
3896 146 : sMapInfo.upperLeftCenter.y = m_gt[3] + m_gt[5] * 0.5;
3897 :
3898 146 : sMapInfo.lowerRightCenter.x = m_gt[0] + m_gt[1] * (GetRasterXSize() - 0.5);
3899 146 : sMapInfo.lowerRightCenter.y = m_gt[3] + m_gt[5] * (GetRasterYSize() - 0.5);
3900 :
3901 146 : sMapInfo.pixelSize.width = std::abs(m_gt[1]);
3902 146 : sMapInfo.pixelSize.height = std::abs(m_gt[5]);
3903 :
3904 : // Handle units. Try to match up with a known name.
3905 146 : sMapInfo.units = const_cast<char *>("meters");
3906 :
3907 146 : if (bHaveSRS && oSRS.IsGeographic())
3908 41 : sMapInfo.units = const_cast<char *>("dd");
3909 105 : else if (bHaveSRS && oSRS.GetLinearUnits() != 1.0)
3910 : {
3911 5 : double dfClosestDiff = 100.0;
3912 5 : int iClosest = -1;
3913 5 : const char *pszUnitName = nullptr;
3914 5 : const double dfActualSize = oSRS.GetLinearUnits(&pszUnitName);
3915 :
3916 5 : const char *const *papszUnitMap = HFAGetUnitMap();
3917 180 : for (int iUnit = 0; papszUnitMap[iUnit] != nullptr; iUnit += 2)
3918 : {
3919 175 : if (fabs(CPLAtof(papszUnitMap[iUnit + 1]) - dfActualSize) <
3920 : dfClosestDiff)
3921 : {
3922 17 : iClosest = iUnit;
3923 17 : dfClosestDiff =
3924 17 : fabs(CPLAtof(papszUnitMap[iUnit + 1]) - dfActualSize);
3925 : }
3926 : }
3927 :
3928 5 : if (iClosest == -1 || fabs(dfClosestDiff / dfActualSize) > 0.0001)
3929 : {
3930 1 : CPLError(CE_Warning, CPLE_NotSupported,
3931 : "Unable to identify Erdas units matching %s/%gm, "
3932 : "output units will be wrong.",
3933 : pszUnitName, dfActualSize);
3934 : }
3935 : else
3936 : {
3937 4 : sMapInfo.units = (char *)papszUnitMap[iClosest];
3938 : }
3939 :
3940 : // We need to convert false easting and northing to meters.
3941 5 : sPro.proParams[6] *= dfActualSize;
3942 5 : sPro.proParams[7] *= dfActualSize;
3943 : }
3944 :
3945 : // Write out definitions.
3946 146 : if (m_gt[2] == 0.0 && m_gt[4] == 0.0)
3947 : {
3948 145 : HFASetMapInfo(hHFA, &sMapInfo);
3949 : }
3950 : else
3951 : {
3952 1 : HFASetGeoTransform(hHFA, sMapInfo.proName, sMapInfo.units, m_gt.data());
3953 : }
3954 :
3955 146 : if (bHaveSRS && sPro.proName != nullptr)
3956 : {
3957 131 : HFASetProParameters(hHFA, &sPro);
3958 131 : HFASetDatum(hHFA, &sDatum);
3959 :
3960 131 : if (!bPEStringStored)
3961 1 : HFASetPEString(hHFA, "");
3962 : }
3963 15 : else if (!bPEStringStored)
3964 : {
3965 13 : ClearSR(hHFA);
3966 : }
3967 :
3968 146 : if (poGeogSRS != nullptr)
3969 133 : delete poGeogSRS;
3970 :
3971 146 : return CE_None;
3972 : }
3973 :
3974 : /************************************************************************/
3975 : /* WritePeStringIfNeeded() */
3976 : /************************************************************************/
3977 132 : int WritePeStringIfNeeded(const OGRSpatialReference *poSRS, HFAHandle hHFA)
3978 : {
3979 132 : if (!poSRS || !hHFA)
3980 0 : return FALSE;
3981 :
3982 132 : const char *pszGEOGCS = poSRS->GetAttrValue("GEOGCS");
3983 132 : if (pszGEOGCS == nullptr)
3984 0 : pszGEOGCS = "";
3985 :
3986 132 : const char *pszDatum = poSRS->GetAttrValue("DATUM");
3987 132 : if (pszDatum == nullptr)
3988 0 : pszDatum = "";
3989 :
3990 : // The strlen() checks are just there to make Coverity happy because it
3991 : // doesn't seem to realize that STARTS_WITH() success implies them.
3992 132 : const size_t gcsNameOffset =
3993 132 : (strlen(pszGEOGCS) > strlen("GCS_") && STARTS_WITH(pszGEOGCS, "GCS_"))
3994 264 : ? strlen("GCS_")
3995 : : 0;
3996 :
3997 132 : const size_t datumNameOffset =
3998 132 : (strlen(pszDatum) > strlen("D_") && STARTS_WITH(pszDatum, "D_"))
3999 264 : ? strlen("D_")
4000 : : 0;
4001 :
4002 132 : bool ret = false;
4003 132 : if (CPLString(pszGEOGCS + gcsNameOffset).replaceAll(' ', '_').tolower() !=
4004 264 : CPLString(pszDatum + datumNameOffset).replaceAll(' ', '_').tolower())
4005 : {
4006 123 : ret = true;
4007 : }
4008 : else
4009 : {
4010 9 : const char *name = poSRS->GetAttrValue("PRIMEM");
4011 9 : if (name && !EQUAL(name, "Greenwich"))
4012 0 : ret = true;
4013 :
4014 9 : if (!ret)
4015 : {
4016 9 : const OGR_SRSNode *poAUnits = poSRS->GetAttrNode("GEOGCS|UNIT");
4017 : const OGR_SRSNode *poChild =
4018 9 : poAUnits == nullptr ? nullptr : poAUnits->GetChild(0);
4019 9 : name = poChild == nullptr ? nullptr : poChild->GetValue();
4020 9 : if (name && !EQUAL(name, "Degree"))
4021 0 : ret = true;
4022 : }
4023 9 : if (!ret)
4024 : {
4025 9 : name = poSRS->GetAttrValue("UNIT");
4026 9 : if (name)
4027 : {
4028 9 : ret = true;
4029 9 : const char *const *papszUnitMap = HFAGetUnitMap();
4030 324 : for (int i = 0; papszUnitMap[i] != nullptr; i += 2)
4031 315 : if (EQUAL(name, papszUnitMap[i]))
4032 0 : ret = false;
4033 : }
4034 : }
4035 9 : if (!ret)
4036 : {
4037 0 : const int nGCS = poSRS->GetEPSGGeogCS();
4038 0 : switch (nGCS)
4039 : {
4040 0 : case 4326:
4041 0 : if (!EQUAL(pszDatum + datumNameOffset, "WGS_84"))
4042 0 : ret = true;
4043 0 : break;
4044 0 : case 4322:
4045 0 : if (!EQUAL(pszDatum + datumNameOffset, "WGS_72"))
4046 0 : ret = true;
4047 0 : break;
4048 0 : case 4267:
4049 0 : if (!EQUAL(pszDatum + datumNameOffset,
4050 : "North_America_1927"))
4051 0 : ret = true;
4052 0 : break;
4053 0 : case 4269:
4054 0 : if (!EQUAL(pszDatum + datumNameOffset,
4055 : "North_America_1983"))
4056 0 : ret = true;
4057 0 : break;
4058 : }
4059 : }
4060 : }
4061 132 : if (ret)
4062 : {
4063 132 : char *pszPEString = nullptr;
4064 264 : OGRSpatialReference oSRSForESRI(*poSRS);
4065 132 : oSRSForESRI.morphToESRI();
4066 132 : oSRSForESRI.exportToWkt(&pszPEString);
4067 132 : HFASetPEString(hHFA, pszPEString);
4068 132 : CPLFree(pszPEString);
4069 : }
4070 :
4071 132 : return ret;
4072 : }
4073 :
4074 : /************************************************************************/
4075 : /* ClearSR() */
4076 : /************************************************************************/
4077 13 : void ClearSR(HFAHandle hHFA)
4078 : {
4079 26 : for (int iBand = 0; iBand < hHFA->nBands; iBand++)
4080 : {
4081 13 : HFAEntry *poMIEntry = nullptr;
4082 26 : if (hHFA->papoBand[iBand]->poNode &&
4083 13 : (poMIEntry = hHFA->papoBand[iBand]->poNode->GetNamedChild(
4084 : "Projection")) != nullptr)
4085 : {
4086 0 : poMIEntry->MarkDirty();
4087 0 : poMIEntry->SetIntField("proType", 0);
4088 0 : poMIEntry->SetIntField("proNumber", 0);
4089 0 : poMIEntry->SetStringField("proExeName", "");
4090 0 : poMIEntry->SetStringField("proName", "");
4091 0 : poMIEntry->SetIntField("proZone", 0);
4092 0 : poMIEntry->SetDoubleField("proParams[0]", 0.0);
4093 0 : poMIEntry->SetDoubleField("proParams[1]", 0.0);
4094 0 : poMIEntry->SetDoubleField("proParams[2]", 0.0);
4095 0 : poMIEntry->SetDoubleField("proParams[3]", 0.0);
4096 0 : poMIEntry->SetDoubleField("proParams[4]", 0.0);
4097 0 : poMIEntry->SetDoubleField("proParams[5]", 0.0);
4098 0 : poMIEntry->SetDoubleField("proParams[6]", 0.0);
4099 0 : poMIEntry->SetDoubleField("proParams[7]", 0.0);
4100 0 : poMIEntry->SetDoubleField("proParams[8]", 0.0);
4101 0 : poMIEntry->SetDoubleField("proParams[9]", 0.0);
4102 0 : poMIEntry->SetDoubleField("proParams[10]", 0.0);
4103 0 : poMIEntry->SetDoubleField("proParams[11]", 0.0);
4104 0 : poMIEntry->SetDoubleField("proParams[12]", 0.0);
4105 0 : poMIEntry->SetDoubleField("proParams[13]", 0.0);
4106 0 : poMIEntry->SetDoubleField("proParams[14]", 0.0);
4107 0 : poMIEntry->SetStringField("proSpheroid.sphereName", "");
4108 0 : poMIEntry->SetDoubleField("proSpheroid.a", 0.0);
4109 0 : poMIEntry->SetDoubleField("proSpheroid.b", 0.0);
4110 0 : poMIEntry->SetDoubleField("proSpheroid.eSquared", 0.0);
4111 0 : poMIEntry->SetDoubleField("proSpheroid.radius", 0.0);
4112 0 : HFAEntry *poDatumEntry = poMIEntry->GetNamedChild("Datum");
4113 0 : if (poDatumEntry != nullptr)
4114 : {
4115 0 : poDatumEntry->MarkDirty();
4116 0 : poDatumEntry->SetStringField("datumname", "");
4117 0 : poDatumEntry->SetIntField("type", 0);
4118 0 : poDatumEntry->SetDoubleField("params[0]", 0.0);
4119 0 : poDatumEntry->SetDoubleField("params[1]", 0.0);
4120 0 : poDatumEntry->SetDoubleField("params[2]", 0.0);
4121 0 : poDatumEntry->SetDoubleField("params[3]", 0.0);
4122 0 : poDatumEntry->SetDoubleField("params[4]", 0.0);
4123 0 : poDatumEntry->SetDoubleField("params[5]", 0.0);
4124 0 : poDatumEntry->SetDoubleField("params[6]", 0.0);
4125 0 : poDatumEntry->SetStringField("gridname", "");
4126 : }
4127 0 : poMIEntry->FlushToDisk();
4128 0 : char *peStr = HFAGetPEString(hHFA);
4129 0 : if (peStr != nullptr && strlen(peStr) > 0)
4130 0 : HFASetPEString(hHFA, "");
4131 : }
4132 : }
4133 13 : }
4134 :
4135 : /************************************************************************/
4136 : /* ReadProjection() */
4137 : /************************************************************************/
4138 :
4139 536 : CPLErr HFADataset::ReadProjection()
4140 :
4141 : {
4142 : // General case for Erdas style projections.
4143 : //
4144 : // We make a particular effort to adapt the mapinfo->proname as
4145 : // the PROJCS[] name per #2422.
4146 536 : const Eprj_Datum *psDatum = HFAGetDatum(hHFA);
4147 536 : const Eprj_ProParameters *psPro = HFAGetProParameters(hHFA);
4148 536 : const Eprj_MapInfo *psMapInfo = HFAGetMapInfo(hHFA);
4149 :
4150 536 : HFAEntry *poMapInformation = nullptr;
4151 536 : if (psMapInfo == nullptr)
4152 : {
4153 290 : HFABand *poBand = hHFA->papoBand[0];
4154 290 : poMapInformation = poBand->poNode->GetNamedChild("MapInformation");
4155 : }
4156 :
4157 536 : m_oSRS.Clear();
4158 :
4159 536 : if (psMapInfo == nullptr && poMapInformation == nullptr)
4160 : {
4161 284 : return CE_None;
4162 : }
4163 252 : else if (((!psDatum || strlen(psDatum->datumname) == 0 ||
4164 252 : EQUAL(psDatum->datumname, "Unknown")) &&
4165 0 : (!psPro || strlen(psPro->proName) == 0 ||
4166 49 : EQUAL(psPro->proName, "Unknown")) &&
4167 49 : (psMapInfo && (strlen(psMapInfo->proName) == 0 ||
4168 49 : EQUAL(psMapInfo->proName, "Unknown"))) &&
4169 0 : (!psPro || psPro->proZone == 0)))
4170 : {
4171 : // It is not clear if Erdas Imagine would recognize a ESRI_PE string
4172 : // alone, but versions of GDAL between 3.0 and 3.6.3 have written CRS
4173 : // using for example the Vertical Projection with a ESRI_PE string only.
4174 43 : char *pszPE_COORDSYS = HFAGetPEString(hHFA);
4175 43 : OGRSpatialReference oSRSFromPE;
4176 43 : oSRSFromPE.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4177 1 : if (pszPE_COORDSYS != nullptr && strlen(pszPE_COORDSYS) > 0 &&
4178 : // Config option for testing purposes only/mostly
4179 45 : CPLTestBool(CPLGetConfigOption("HFA_USE_ESRI_PE_STRING", "YES")) &&
4180 1 : oSRSFromPE.importFromWkt(pszPE_COORDSYS) == OGRERR_NONE)
4181 : {
4182 : const char *pszProjName =
4183 1 : oSRSFromPE.GetAttrValue("PROJCS|PROJECTION");
4184 2 : if (pszProjName &&
4185 1 : (EQUAL(pszProjName, "Vertical Perspective") ||
4186 3 : EQUAL(pszProjName, "Vertical_Near_Side_Perspective")) &&
4187 1 : CPLTestBool(CPLGetConfigOption(
4188 : "HFA_SHOW_ESRI_PE_STRING_ONLY_WARNING", "YES")))
4189 : {
4190 1 : CPLError(CE_Warning, CPLE_AppDefined,
4191 : "A ESRI_PE string encoding a CRS has been found for "
4192 : "projection method %s, but no corresponding "
4193 : "Eprj_ProParameters are present. This file has likely "
4194 : "been generated by GDAL >= 3.0 and <= 3.6.2. It is "
4195 : "recommended to recreate it, e.g with gdal_translate, "
4196 : "with GDAL >= 3.6.3. This warning can be suppressed "
4197 : "by setting the HFA_SHOW_ESRI_PE_STRING_ONLY_WARNING "
4198 : "configuration option to NO.",
4199 : pszProjName);
4200 : }
4201 1 : m_oSRS = std::move(oSRSFromPE);
4202 : }
4203 43 : CPLFree(pszPE_COORDSYS);
4204 43 : return m_oSRS.IsEmpty() ? CE_Failure : CE_None;
4205 : }
4206 :
4207 418 : auto poSRS = HFAPCSStructToOSR(psDatum, psPro, psMapInfo, poMapInformation);
4208 209 : if (poSRS)
4209 209 : m_oSRS = *poSRS;
4210 :
4211 : // If we got a valid projection and managed to identify a EPSG code,
4212 : // then do not use the ESRI PE String.
4213 : const bool bTryReadingPEString =
4214 209 : poSRS == nullptr || poSRS->GetAuthorityCode(nullptr) == nullptr;
4215 :
4216 : // Special logic for PE string in ProjectionX node.
4217 209 : char *pszPE_COORDSYS = nullptr;
4218 209 : if (bTryReadingPEString)
4219 38 : pszPE_COORDSYS = HFAGetPEString(hHFA);
4220 :
4221 209 : OGRSpatialReference oSRSFromPE;
4222 209 : oSRSFromPE.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4223 18 : if (pszPE_COORDSYS != nullptr && strlen(pszPE_COORDSYS) > 0 &&
4224 : // Config option for testing purposes only/mostly
4225 240 : CPLTestBool(CPLGetConfigOption("HFA_USE_ESRI_PE_STRING", "YES")) &&
4226 13 : oSRSFromPE.importFromWkt(pszPE_COORDSYS) == OGRERR_NONE)
4227 : {
4228 13 : m_oSRS = std::move(oSRSFromPE);
4229 :
4230 : // Copy TOWGS84 clause from HFA SRS to PE SRS.
4231 13 : if (poSRS != nullptr)
4232 : {
4233 : double adfCoeffs[7];
4234 : double adfCoeffsUnused[7];
4235 13 : if (poSRS->GetTOWGS84(adfCoeffs, 7) == OGRERR_NONE &&
4236 0 : m_oSRS.GetTOWGS84(adfCoeffsUnused, 7) == OGRERR_FAILURE)
4237 : {
4238 0 : m_oSRS.SetTOWGS84(adfCoeffs[0], adfCoeffs[1], adfCoeffs[2],
4239 : adfCoeffs[3], adfCoeffs[4], adfCoeffs[5],
4240 : adfCoeffs[6]);
4241 : }
4242 : }
4243 : }
4244 :
4245 209 : CPLFree(pszPE_COORDSYS);
4246 :
4247 209 : return m_oSRS.IsEmpty() ? CE_Failure : CE_None;
4248 : }
4249 :
4250 : /************************************************************************/
4251 : /* IBuildOverviews() */
4252 : /************************************************************************/
4253 :
4254 11 : CPLErr HFADataset::IBuildOverviews(const char *pszResampling, int nOverviews,
4255 : const int *panOverviewList, int nListBands,
4256 : const int *panBandList,
4257 : GDALProgressFunc pfnProgress,
4258 : void *pProgressData,
4259 : CSLConstList papszOptions)
4260 :
4261 : {
4262 11 : if (GetAccess() == GA_ReadOnly)
4263 : {
4264 0 : for (int i = 0; i < nListBands; i++)
4265 : {
4266 0 : if (HFAGetOverviewCount(hHFA, panBandList[i]) > 0)
4267 : {
4268 0 : CPLError(CE_Failure, CPLE_NotSupported,
4269 : "Cannot add external overviews when there are already "
4270 : "internal overviews");
4271 0 : return CE_Failure;
4272 : }
4273 : }
4274 :
4275 0 : return GDALDataset::IBuildOverviews(
4276 : pszResampling, nOverviews, panOverviewList, nListBands, panBandList,
4277 0 : pfnProgress, pProgressData, papszOptions);
4278 : }
4279 :
4280 22 : for (int i = 0; i < nListBands; i++)
4281 : {
4282 22 : void *pScaledProgressData = GDALCreateScaledProgress(
4283 11 : i * 1.0 / nListBands, (i + 1) * 1.0 / nListBands, pfnProgress,
4284 : pProgressData);
4285 :
4286 11 : GDALRasterBand *poBand = GetRasterBand(panBandList[i]);
4287 :
4288 : // GetRasterBand can return NULL.
4289 11 : if (poBand == nullptr)
4290 : {
4291 0 : CPLError(CE_Failure, CPLE_ObjectNull, "GetRasterBand failed");
4292 0 : GDALDestroyScaledProgress(pScaledProgressData);
4293 0 : return CE_Failure;
4294 : }
4295 :
4296 22 : const CPLErr eErr = poBand->BuildOverviews(
4297 : pszResampling, nOverviews, panOverviewList, GDALScaledProgress,
4298 11 : pScaledProgressData, papszOptions);
4299 :
4300 11 : GDALDestroyScaledProgress(pScaledProgressData);
4301 :
4302 11 : if (eErr != CE_None)
4303 0 : return eErr;
4304 : }
4305 :
4306 11 : return CE_None;
4307 : }
4308 :
4309 : /************************************************************************/
4310 : /* Identify() */
4311 : /************************************************************************/
4312 :
4313 64079 : int HFADataset::Identify(GDALOpenInfo *poOpenInfo)
4314 :
4315 : {
4316 : // Verify that this is a HFA file.
4317 64079 : if (poOpenInfo->nHeaderBytes < 15 ||
4318 8832 : !STARTS_WITH_CI((char *)poOpenInfo->pabyHeader, "EHFA_HEADER_TAG"))
4319 62990 : return FALSE;
4320 :
4321 1089 : return TRUE;
4322 : }
4323 :
4324 : /************************************************************************/
4325 : /* Open() */
4326 : /************************************************************************/
4327 :
4328 540 : GDALDataset *HFADataset::Open(GDALOpenInfo *poOpenInfo)
4329 :
4330 : {
4331 : // Verify that this is a HFA file.
4332 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
4333 : // During fuzzing, do not use Identify to reject crazy content.
4334 540 : if (!Identify(poOpenInfo))
4335 0 : return nullptr;
4336 : #endif
4337 :
4338 : // Open the file.
4339 540 : HFAHandle hHFA = HFAOpen(poOpenInfo->pszFilename,
4340 540 : (poOpenInfo->eAccess == GA_Update ? "r+" : "r"));
4341 :
4342 540 : if (hHFA == nullptr)
4343 0 : return nullptr;
4344 :
4345 : // Create a corresponding GDALDataset.
4346 540 : HFADataset *poDS = new HFADataset();
4347 :
4348 540 : poDS->hHFA = hHFA;
4349 540 : poDS->eAccess = poOpenInfo->eAccess;
4350 :
4351 : // Establish raster info.
4352 540 : HFAGetRasterInfo(hHFA, &poDS->nRasterXSize, &poDS->nRasterYSize,
4353 : &poDS->nBands);
4354 :
4355 540 : if (poDS->nBands == 0)
4356 : {
4357 4 : delete poDS;
4358 4 : CPLError(CE_Failure, CPLE_AppDefined,
4359 : "Unable to open %s, it has zero usable bands.",
4360 : poOpenInfo->pszFilename);
4361 4 : return nullptr;
4362 : }
4363 :
4364 536 : if (poDS->nRasterXSize == 0 || poDS->nRasterYSize == 0)
4365 : {
4366 0 : delete poDS;
4367 0 : CPLError(CE_Failure, CPLE_AppDefined,
4368 : "Unable to open %s, it has no pixels.",
4369 : poOpenInfo->pszFilename);
4370 0 : return nullptr;
4371 : }
4372 :
4373 : // Get geotransform, or if that fails, try to find XForms to
4374 : // build gcps, and metadata.
4375 536 : if (!HFAGetGeoTransform(hHFA, poDS->m_gt.data()))
4376 : {
4377 286 : Efga_Polynomial *pasPolyListForward = nullptr;
4378 286 : Efga_Polynomial *pasPolyListReverse = nullptr;
4379 : const int nStepCount =
4380 286 : HFAReadXFormStack(hHFA, &pasPolyListForward, &pasPolyListReverse);
4381 :
4382 286 : if (nStepCount > 0)
4383 : {
4384 1 : poDS->UseXFormStack(nStepCount, pasPolyListForward,
4385 : pasPolyListReverse);
4386 1 : CPLFree(pasPolyListForward);
4387 1 : CPLFree(pasPolyListReverse);
4388 : }
4389 : }
4390 :
4391 536 : poDS->ReadProjection();
4392 :
4393 536 : char **papszCM = HFAReadCameraModel(hHFA);
4394 :
4395 536 : if (papszCM != nullptr)
4396 : {
4397 1 : poDS->SetMetadata(papszCM, "CAMERA_MODEL");
4398 1 : CSLDestroy(papszCM);
4399 : }
4400 :
4401 1152 : for (int i = 0; i < poDS->nBands; i++)
4402 : {
4403 616 : poDS->SetBand(i + 1, new HFARasterBand(poDS, i + 1, -1));
4404 : }
4405 :
4406 : // Collect GDAL custom Metadata, and "auxiliary" metadata from
4407 : // well known HFA structures for the bands. We defer this till
4408 : // now to ensure that the bands are properly setup before
4409 : // interacting with PAM.
4410 1152 : for (int i = 0; i < poDS->nBands; i++)
4411 : {
4412 : HFARasterBand *poBand =
4413 616 : static_cast<HFARasterBand *>(poDS->GetRasterBand(i + 1));
4414 :
4415 616 : char **papszMD = HFAGetMetadata(hHFA, i + 1);
4416 616 : if (papszMD != nullptr)
4417 : {
4418 10 : poBand->SetMetadata(papszMD);
4419 10 : CSLDestroy(papszMD);
4420 : }
4421 :
4422 616 : poBand->ReadAuxMetadata();
4423 616 : poBand->ReadHistogramMetadata();
4424 : }
4425 :
4426 : // Check for GDAL style metadata.
4427 536 : char **papszMD = HFAGetMetadata(hHFA, 0);
4428 536 : if (papszMD != nullptr)
4429 : {
4430 85 : poDS->SetMetadata(papszMD);
4431 85 : CSLDestroy(papszMD);
4432 : }
4433 :
4434 : // Read the elevation metadata, if present.
4435 1152 : for (int iBand = 0; iBand < poDS->nBands; iBand++)
4436 : {
4437 : HFARasterBand *poBand =
4438 616 : static_cast<HFARasterBand *>(poDS->GetRasterBand(iBand + 1));
4439 616 : const char *pszEU = HFAReadElevationUnit(hHFA, iBand);
4440 :
4441 616 : if (pszEU != nullptr)
4442 : {
4443 3 : poBand->SetUnitType(pszEU);
4444 3 : if (poDS->nBands == 1)
4445 : {
4446 3 : poDS->SetMetadataItem("ELEVATION_UNITS", pszEU);
4447 : }
4448 : }
4449 : }
4450 :
4451 : // Check for dependent dataset value.
4452 536 : HFAInfo_t *psInfo = hHFA;
4453 536 : HFAEntry *poEntry = psInfo->poRoot->GetNamedChild("DependentFile");
4454 536 : if (poEntry != nullptr)
4455 : {
4456 30 : poDS->SetMetadataItem("HFA_DEPENDENT_FILE",
4457 : poEntry->GetStringField("dependent.string"),
4458 : "HFA");
4459 : }
4460 :
4461 : // Initialize any PAM information.
4462 536 : poDS->SetDescription(poOpenInfo->pszFilename);
4463 536 : poDS->TryLoadXML();
4464 :
4465 : // Check for external overviews.
4466 536 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
4467 :
4468 : // Clear dirty metadata flags.
4469 1152 : for (int i = 0; i < poDS->nBands; i++)
4470 : {
4471 : HFARasterBand *poBand =
4472 616 : static_cast<HFARasterBand *>(poDS->GetRasterBand(i + 1));
4473 616 : poBand->bMetadataDirty = false;
4474 : }
4475 536 : poDS->bMetadataDirty = false;
4476 :
4477 536 : return poDS;
4478 : }
4479 :
4480 : /************************************************************************/
4481 : /* GetSpatialRef() */
4482 : /************************************************************************/
4483 :
4484 154 : const OGRSpatialReference *HFADataset::GetSpatialRef() const
4485 : {
4486 154 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
4487 : }
4488 :
4489 : /************************************************************************/
4490 : /* SetSpatialRef() */
4491 : /************************************************************************/
4492 :
4493 133 : CPLErr HFADataset::SetSpatialRef(const OGRSpatialReference *poSRS)
4494 :
4495 : {
4496 133 : m_oSRS.Clear();
4497 133 : if (poSRS)
4498 133 : m_oSRS = *poSRS;
4499 133 : bGeoDirty = true;
4500 :
4501 133 : return CE_None;
4502 : }
4503 :
4504 : /************************************************************************/
4505 : /* SetMetadata() */
4506 : /************************************************************************/
4507 :
4508 136 : CPLErr HFADataset::SetMetadata(char **papszMDIn, const char *pszDomain)
4509 :
4510 : {
4511 136 : bMetadataDirty = true;
4512 :
4513 136 : return GDALPamDataset::SetMetadata(papszMDIn, pszDomain);
4514 : }
4515 :
4516 : /************************************************************************/
4517 : /* SetMetadata() */
4518 : /************************************************************************/
4519 :
4520 33 : CPLErr HFADataset::SetMetadataItem(const char *pszTag, const char *pszValue,
4521 : const char *pszDomain)
4522 :
4523 : {
4524 33 : bMetadataDirty = true;
4525 :
4526 33 : return GDALPamDataset::SetMetadataItem(pszTag, pszValue, pszDomain);
4527 : }
4528 :
4529 : /************************************************************************/
4530 : /* GetGeoTransform() */
4531 : /************************************************************************/
4532 :
4533 142 : CPLErr HFADataset::GetGeoTransform(GDALGeoTransform >) const
4534 :
4535 : {
4536 154 : if (m_gt[0] != 0.0 || m_gt[1] != 1.0 || m_gt[2] != 0.0 || m_gt[3] != 0.0 ||
4537 154 : m_gt[4] != 0.0 || m_gt[5] != 1.0)
4538 : {
4539 136 : gt = m_gt;
4540 136 : return CE_None;
4541 : }
4542 :
4543 6 : return GDALPamDataset::GetGeoTransform(gt);
4544 : }
4545 :
4546 : /************************************************************************/
4547 : /* SetGeoTransform() */
4548 : /************************************************************************/
4549 :
4550 86 : CPLErr HFADataset::SetGeoTransform(const GDALGeoTransform >)
4551 :
4552 : {
4553 86 : m_gt = gt;
4554 86 : bGeoDirty = true;
4555 :
4556 86 : return CE_None;
4557 : }
4558 :
4559 : /************************************************************************/
4560 : /* IRasterIO() */
4561 : /* */
4562 : /* Multi-band raster io handler. Here we ensure that the block */
4563 : /* based loading is used for spill file rasters. That is */
4564 : /* because they are effectively pixel interleaved, so */
4565 : /* processing all bands for a given block together avoid extra */
4566 : /* seeks. */
4567 : /************************************************************************/
4568 :
4569 80 : CPLErr HFADataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
4570 : int nXSize, int nYSize, void *pData, int nBufXSize,
4571 : int nBufYSize, GDALDataType eBufType,
4572 : int nBandCount, BANDMAP_TYPE panBandMap,
4573 : GSpacing nPixelSpace, GSpacing nLineSpace,
4574 : GSpacing nBandSpace,
4575 : GDALRasterIOExtraArg *psExtraArg)
4576 :
4577 : {
4578 80 : if (hHFA->papoBand[panBandMap[0] - 1]->fpExternal != nullptr &&
4579 : nBandCount > 1)
4580 0 : return GDALDataset::BlockBasedRasterIO(
4581 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
4582 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
4583 0 : nBandSpace, psExtraArg);
4584 :
4585 80 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
4586 : nBufXSize, nBufYSize, eBufType, nBandCount,
4587 : panBandMap, nPixelSpace, nLineSpace,
4588 80 : nBandSpace, psExtraArg);
4589 : }
4590 :
4591 : /************************************************************************/
4592 : /* UseXFormStack() */
4593 : /************************************************************************/
4594 :
4595 1 : void HFADataset::UseXFormStack(int nStepCount, Efga_Polynomial *pasPLForward,
4596 : Efga_Polynomial *pasPLReverse)
4597 :
4598 : {
4599 : // Generate GCPs using the transform.
4600 7 : for (double dfYRatio = 0.0; dfYRatio < 1.001; dfYRatio += 0.2)
4601 : {
4602 42 : for (double dfXRatio = 0.0; dfXRatio < 1.001; dfXRatio += 0.2)
4603 : {
4604 36 : const double dfLine = 0.5 + (GetRasterYSize() - 1) * dfYRatio;
4605 36 : const double dfPixel = 0.5 + (GetRasterXSize() - 1) * dfXRatio;
4606 :
4607 : gdal::GCP gcp("", "", dfPixel, dfLine,
4608 : /* X = */ dfPixel,
4609 72 : /* Y = */ dfLine);
4610 36 : if (HFAEvaluateXFormStack(nStepCount, FALSE, pasPLReverse,
4611 72 : &(gcp.X()), &(gcp.Y())))
4612 : {
4613 36 : m_aoGCPs.emplace_back(std::move(gcp));
4614 : }
4615 : }
4616 : }
4617 :
4618 : // Store the transform as metadata.
4619 1 : GDALMajorObject::SetMetadataItem(
4620 2 : "XFORM_STEPS", CPLString().Printf("%d", nStepCount), "XFORMS");
4621 :
4622 3 : for (int iStep = 0; iStep < nStepCount; iStep++)
4623 : {
4624 2 : GDALMajorObject::SetMetadataItem(
4625 4 : CPLString().Printf("XFORM%d_ORDER", iStep),
4626 4 : CPLString().Printf("%d", pasPLForward[iStep].order), "XFORMS");
4627 :
4628 2 : if (pasPLForward[iStep].order == 1)
4629 : {
4630 5 : for (int i = 0; i < 4; i++)
4631 4 : GDALMajorObject::SetMetadataItem(
4632 8 : CPLString().Printf("XFORM%d_POLYCOEFMTX[%d]", iStep, i),
4633 8 : CPLString().Printf("%.15g",
4634 4 : pasPLForward[iStep].polycoefmtx[i]),
4635 : "XFORMS");
4636 :
4637 3 : for (int i = 0; i < 2; i++)
4638 2 : GDALMajorObject::SetMetadataItem(
4639 4 : CPLString().Printf("XFORM%d_POLYCOEFVECTOR[%d]", iStep, i),
4640 4 : CPLString().Printf("%.15g",
4641 2 : pasPLForward[iStep].polycoefvector[i]),
4642 : "XFORMS");
4643 :
4644 1 : continue;
4645 : }
4646 :
4647 1 : int nCoefCount = 10;
4648 :
4649 1 : if (pasPLForward[iStep].order != 2)
4650 : {
4651 1 : CPLAssert(pasPLForward[iStep].order == 3);
4652 1 : nCoefCount = 18;
4653 : }
4654 :
4655 19 : for (int i = 0; i < nCoefCount; i++)
4656 18 : GDALMajorObject::SetMetadataItem(
4657 36 : CPLString().Printf("XFORM%d_FWD_POLYCOEFMTX[%d]", iStep, i),
4658 36 : CPLString().Printf("%.15g", pasPLForward[iStep].polycoefmtx[i]),
4659 : "XFORMS");
4660 :
4661 3 : for (int i = 0; i < 2; i++)
4662 2 : GDALMajorObject::SetMetadataItem(
4663 4 : CPLString().Printf("XFORM%d_FWD_POLYCOEFVECTOR[%d]", iStep, i),
4664 4 : CPLString().Printf("%.15g",
4665 2 : pasPLForward[iStep].polycoefvector[i]),
4666 : "XFORMS");
4667 :
4668 19 : for (int i = 0; i < nCoefCount; i++)
4669 18 : GDALMajorObject::SetMetadataItem(
4670 36 : CPLString().Printf("XFORM%d_REV_POLYCOEFMTX[%d]", iStep, i),
4671 36 : CPLString().Printf("%.15g", pasPLReverse[iStep].polycoefmtx[i]),
4672 : "XFORMS");
4673 :
4674 3 : for (int i = 0; i < 2; i++)
4675 2 : GDALMajorObject::SetMetadataItem(
4676 4 : CPLString().Printf("XFORM%d_REV_POLYCOEFVECTOR[%d]", iStep, i),
4677 4 : CPLString().Printf("%.15g",
4678 2 : pasPLReverse[iStep].polycoefvector[i]),
4679 : "XFORMS");
4680 : }
4681 1 : }
4682 :
4683 : /************************************************************************/
4684 : /* GetGCPCount() */
4685 : /************************************************************************/
4686 :
4687 25 : int HFADataset::GetGCPCount()
4688 : {
4689 25 : const int nPAMCount = GDALPamDataset::GetGCPCount();
4690 25 : return nPAMCount > 0 ? nPAMCount : static_cast<int>(m_aoGCPs.size());
4691 : }
4692 :
4693 : /************************************************************************/
4694 : /* GetGCPSpatialRef() */
4695 : /************************************************************************/
4696 :
4697 1 : const OGRSpatialReference *HFADataset::GetGCPSpatialRef() const
4698 :
4699 : {
4700 1 : const OGRSpatialReference *poSRS = GDALPamDataset::GetGCPSpatialRef();
4701 1 : if (poSRS)
4702 1 : return poSRS;
4703 0 : return !m_aoGCPs.empty() && !m_oSRS.IsEmpty() ? &m_oSRS : nullptr;
4704 : }
4705 :
4706 : /************************************************************************/
4707 : /* GetGCPs() */
4708 : /************************************************************************/
4709 :
4710 2 : const GDAL_GCP *HFADataset::GetGCPs()
4711 : {
4712 2 : const GDAL_GCP *psPAMGCPs = GDALPamDataset::GetGCPs();
4713 2 : if (psPAMGCPs)
4714 1 : return psPAMGCPs;
4715 1 : return gdal::GCP::c_ptr(m_aoGCPs);
4716 : }
4717 :
4718 : /************************************************************************/
4719 : /* GetFileList() */
4720 : /************************************************************************/
4721 :
4722 92 : char **HFADataset::GetFileList()
4723 :
4724 : {
4725 184 : CPLStringList oFileList(GDALPamDataset::GetFileList());
4726 :
4727 184 : const std::string osIGEFilename = HFAGetIGEFilename(hHFA);
4728 92 : if (!osIGEFilename.empty())
4729 : {
4730 14 : oFileList.push_back(osIGEFilename);
4731 : }
4732 :
4733 : // Request an overview to force opening of dependent overview files.
4734 92 : if (nBands > 0 && GetRasterBand(1)->GetOverviewCount() > 0)
4735 12 : GetRasterBand(1)->GetOverview(0);
4736 :
4737 92 : if (hHFA->psDependent != nullptr)
4738 : {
4739 6 : HFAInfo_t *psDep = hHFA->psDependent;
4740 :
4741 6 : oFileList.push_back(
4742 12 : CPLFormFilenameSafe(psDep->pszPath, psDep->pszFilename, nullptr));
4743 :
4744 12 : const std::string osIGEFilenameDep = HFAGetIGEFilename(psDep);
4745 6 : if (!osIGEFilenameDep.empty())
4746 5 : oFileList.push_back(osIGEFilenameDep);
4747 : }
4748 :
4749 184 : return oFileList.StealList();
4750 : }
4751 :
4752 : /************************************************************************/
4753 : /* Create() */
4754 : /************************************************************************/
4755 :
4756 212 : GDALDataset *HFADataset::Create(const char *pszFilenameIn, int nXSize,
4757 : int nYSize, int nBandsIn, GDALDataType eType,
4758 : char **papszParamList)
4759 :
4760 : {
4761 212 : const int nBits = CSLFetchNameValue(papszParamList, "NBITS") != nullptr
4762 212 : ? atoi(CSLFetchNameValue(papszParamList, "NBITS"))
4763 212 : : 0;
4764 :
4765 212 : const char *pszPixelType = CSLFetchNameValue(papszParamList, "PIXELTYPE");
4766 212 : if (pszPixelType == nullptr)
4767 212 : pszPixelType = "";
4768 :
4769 : // Translate the data type.
4770 : EPTType eHfaDataType;
4771 212 : switch (eType)
4772 : {
4773 131 : case GDT_Byte:
4774 131 : if (nBits == 1)
4775 2 : eHfaDataType = EPT_u1;
4776 129 : else if (nBits == 2)
4777 2 : eHfaDataType = EPT_u2;
4778 127 : else if (nBits == 4)
4779 2 : eHfaDataType = EPT_u4;
4780 125 : else if (EQUAL(pszPixelType, "SIGNEDBYTE"))
4781 0 : eHfaDataType = EPT_s8;
4782 : else
4783 125 : eHfaDataType = EPT_u8;
4784 131 : break;
4785 :
4786 2 : case GDT_Int8:
4787 2 : eHfaDataType = EPT_s8;
4788 2 : break;
4789 :
4790 12 : case GDT_UInt16:
4791 12 : eHfaDataType = EPT_u16;
4792 12 : break;
4793 :
4794 7 : case GDT_Int16:
4795 7 : eHfaDataType = EPT_s16;
4796 7 : break;
4797 :
4798 8 : case GDT_Int32:
4799 8 : eHfaDataType = EPT_s32;
4800 8 : break;
4801 :
4802 8 : case GDT_UInt32:
4803 8 : eHfaDataType = EPT_u32;
4804 8 : break;
4805 :
4806 9 : case GDT_Float32:
4807 9 : eHfaDataType = EPT_f32;
4808 9 : break;
4809 :
4810 10 : case GDT_Float64:
4811 10 : eHfaDataType = EPT_f64;
4812 10 : break;
4813 :
4814 7 : case GDT_CFloat32:
4815 7 : eHfaDataType = EPT_c64;
4816 7 : break;
4817 :
4818 7 : case GDT_CFloat64:
4819 7 : eHfaDataType = EPT_c128;
4820 7 : break;
4821 :
4822 11 : default:
4823 11 : CPLError(
4824 : CE_Failure, CPLE_NotSupported,
4825 : "Data type %s not supported by Erdas Imagine (HFA) format.",
4826 : GDALGetDataTypeName(eType));
4827 11 : return nullptr;
4828 : }
4829 :
4830 : const bool bForceToPEString =
4831 201 : CPLFetchBool(papszParamList, "FORCETOPESTRING", false);
4832 : const bool bDisablePEString =
4833 201 : CPLFetchBool(papszParamList, "DISABLEPESTRING", false);
4834 201 : if (bForceToPEString && bDisablePEString)
4835 : {
4836 0 : CPLError(CE_Failure, CPLE_AppDefined,
4837 : "FORCETOPESTRING and DISABLEPESTRING are mutually exclusive");
4838 0 : return nullptr;
4839 : }
4840 :
4841 : // Create the new file.
4842 201 : HFAHandle hHFA = HFACreate(pszFilenameIn, nXSize, nYSize, nBandsIn,
4843 : eHfaDataType, papszParamList);
4844 201 : if (hHFA == nullptr)
4845 12 : return nullptr;
4846 :
4847 189 : if (HFAClose(hHFA) != 0)
4848 : {
4849 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
4850 0 : return nullptr;
4851 : }
4852 :
4853 : // Open the dataset normally.
4854 189 : HFADataset *poDS = (HFADataset *)GDALOpen(pszFilenameIn, GA_Update);
4855 :
4856 : // Special creation option to disable checking for UTM
4857 : // parameters when writing the projection. This is a special
4858 : // hack for sam.gillingham@nrm.qld.gov.au.
4859 189 : if (poDS != nullptr)
4860 : {
4861 187 : poDS->bIgnoreUTM = CPLFetchBool(papszParamList, "IGNOREUTM", false);
4862 : }
4863 :
4864 : // Sometimes we can improve ArcGIS compatibility by forcing
4865 : // generation of a PEString instead of traditional Imagine
4866 : // coordinate system descriptions.
4867 189 : if (poDS != nullptr)
4868 : {
4869 187 : poDS->bForceToPEString = bForceToPEString;
4870 187 : poDS->bDisablePEString = bDisablePEString;
4871 : }
4872 :
4873 189 : return poDS;
4874 : }
4875 :
4876 : /************************************************************************/
4877 : /* Rename() */
4878 : /* */
4879 : /* Custom Rename() implementation that knows how to update */
4880 : /* filename references in .img and .aux files. */
4881 : /************************************************************************/
4882 :
4883 1 : CPLErr HFADataset::Rename(const char *pszNewName, const char *pszOldName)
4884 :
4885 : {
4886 : // Rename all the files at the filesystem level.
4887 1 : CPLErr eErr = GDALDriver::DefaultRename(pszNewName, pszOldName);
4888 1 : if (eErr != CE_None)
4889 0 : return eErr;
4890 :
4891 : // Now try to go into the .img file and update RRDNames[] lists.
4892 2 : CPLString osOldBasename = CPLGetBasenameSafe(pszOldName);
4893 1 : CPLString osNewBasename = CPLGetBasenameSafe(pszNewName);
4894 :
4895 1 : if (osOldBasename != osNewBasename)
4896 : {
4897 1 : HFAHandle hHFA = HFAOpen(pszNewName, "r+");
4898 :
4899 1 : if (hHFA != nullptr)
4900 : {
4901 1 : eErr = HFARenameReferences(hHFA, osNewBasename, osOldBasename);
4902 :
4903 1 : HFAGetOverviewCount(hHFA, 1);
4904 :
4905 1 : if (hHFA->psDependent != nullptr)
4906 1 : HFARenameReferences(hHFA->psDependent, osNewBasename,
4907 : osOldBasename);
4908 :
4909 1 : if (HFAClose(hHFA) != 0)
4910 0 : eErr = CE_Failure;
4911 : }
4912 : }
4913 :
4914 1 : return eErr;
4915 : }
4916 :
4917 : /************************************************************************/
4918 : /* CopyFiles() */
4919 : /* */
4920 : /* Custom CopyFiles() implementation that knows how to update */
4921 : /* filename references in .img and .aux files. */
4922 : /************************************************************************/
4923 :
4924 1 : CPLErr HFADataset::CopyFiles(const char *pszNewName, const char *pszOldName)
4925 :
4926 : {
4927 : // Rename all the files at the filesystem level.
4928 1 : CPLErr eErr = GDALDriver::DefaultCopyFiles(pszNewName, pszOldName);
4929 :
4930 1 : if (eErr != CE_None)
4931 0 : return eErr;
4932 :
4933 : // Now try to go into the .img file and update RRDNames[] lists.
4934 2 : CPLString osOldBasename = CPLGetBasenameSafe(pszOldName);
4935 1 : CPLString osNewBasename = CPLGetBasenameSafe(pszNewName);
4936 :
4937 1 : if (osOldBasename != osNewBasename)
4938 : {
4939 1 : HFAHandle hHFA = HFAOpen(pszNewName, "r+");
4940 :
4941 1 : if (hHFA != nullptr)
4942 : {
4943 1 : eErr = HFARenameReferences(hHFA, osNewBasename, osOldBasename);
4944 :
4945 1 : HFAGetOverviewCount(hHFA, 1);
4946 :
4947 1 : if (hHFA->psDependent != nullptr)
4948 1 : HFARenameReferences(hHFA->psDependent, osNewBasename,
4949 : osOldBasename);
4950 :
4951 1 : if (HFAClose(hHFA) != 0)
4952 0 : eErr = CE_Failure;
4953 : }
4954 : }
4955 :
4956 1 : return eErr;
4957 : }
4958 :
4959 : /************************************************************************/
4960 : /* CreateCopy() */
4961 : /************************************************************************/
4962 :
4963 66 : GDALDataset *HFADataset::CreateCopy(const char *pszFilename,
4964 : GDALDataset *poSrcDS, int /* bStrict */,
4965 : char **papszOptions,
4966 : GDALProgressFunc pfnProgress,
4967 : void *pProgressData)
4968 : {
4969 : // Do we really just want to create an .aux file?
4970 66 : const bool bCreateAux = CPLFetchBool(papszOptions, "AUX", false);
4971 :
4972 : // Establish a representative data type to use.
4973 66 : char **papszModOptions = CSLDuplicate(papszOptions);
4974 66 : if (!pfnProgress(0.0, nullptr, pProgressData))
4975 : {
4976 0 : CSLDestroy(papszModOptions);
4977 0 : return nullptr;
4978 : }
4979 :
4980 66 : const int nBandCount = poSrcDS->GetRasterCount();
4981 66 : GDALDataType eType = GDT_Unknown;
4982 :
4983 141 : for (int iBand = 0; iBand < nBandCount; iBand++)
4984 : {
4985 75 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(iBand + 1);
4986 75 : if (iBand == 0)
4987 65 : eType = poBand->GetRasterDataType();
4988 : else
4989 10 : eType = GDALDataTypeUnion(eType, poBand->GetRasterDataType());
4990 : }
4991 :
4992 : // If we have PIXELTYPE metadata in the source, pass it
4993 : // through as a creation option.
4994 132 : if (CSLFetchNameValue(papszOptions, "PIXELTYPE") == nullptr &&
4995 132 : nBandCount > 0 && eType == GDT_Byte)
4996 : {
4997 42 : auto poSrcBand = poSrcDS->GetRasterBand(1);
4998 42 : poSrcBand->EnablePixelTypeSignedByteWarning(false);
4999 : const char *pszPixelType =
5000 42 : poSrcBand->GetMetadataItem("PIXELTYPE", "IMAGE_STRUCTURE");
5001 42 : poSrcBand->EnablePixelTypeSignedByteWarning(true);
5002 42 : if (pszPixelType)
5003 : {
5004 : papszModOptions =
5005 0 : CSLSetNameValue(papszModOptions, "PIXELTYPE", pszPixelType);
5006 : }
5007 : }
5008 :
5009 66 : HFADataset *poDS = cpl::down_cast<HFADataset *>(
5010 : Create(pszFilename, poSrcDS->GetRasterXSize(),
5011 : poSrcDS->GetRasterYSize(), nBandCount, eType, papszModOptions));
5012 :
5013 66 : CSLDestroy(papszModOptions);
5014 :
5015 66 : if (poDS == nullptr)
5016 16 : return nullptr;
5017 :
5018 : // Does the source have a PCT or RAT for any of the bands? If so, copy it
5019 : // over.
5020 110 : for (int iBand = 0; iBand < nBandCount; iBand++)
5021 : {
5022 60 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(iBand + 1);
5023 :
5024 60 : GDALColorTable *poCT = poBand->GetColorTable();
5025 60 : if (poCT != nullptr)
5026 : {
5027 1 : poDS->GetRasterBand(iBand + 1)->SetColorTable(poCT);
5028 : }
5029 :
5030 60 : if (poBand->GetDefaultRAT() != nullptr)
5031 5 : poDS->GetRasterBand(iBand + 1)->SetDefaultRAT(
5032 5 : poBand->GetDefaultRAT());
5033 : }
5034 :
5035 : // Do we have metadata for any of the bands or the dataset as a whole?
5036 50 : if (poSrcDS->GetMetadata() != nullptr)
5037 39 : poDS->SetMetadata(poSrcDS->GetMetadata());
5038 :
5039 110 : for (int iBand = 0; iBand < nBandCount; iBand++)
5040 : {
5041 60 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand + 1);
5042 60 : GDALRasterBand *poDstBand = poDS->GetRasterBand(iBand + 1);
5043 :
5044 60 : if (poSrcBand->GetMetadata() != nullptr)
5045 7 : poDstBand->SetMetadata(poSrcBand->GetMetadata());
5046 :
5047 60 : if (strlen(poSrcBand->GetDescription()) > 0)
5048 6 : poDstBand->SetDescription(poSrcBand->GetDescription());
5049 :
5050 60 : int bSuccess = FALSE;
5051 60 : const double dfNoDataValue = poSrcBand->GetNoDataValue(&bSuccess);
5052 60 : if (bSuccess)
5053 2 : poDstBand->SetNoDataValue(dfNoDataValue);
5054 : }
5055 :
5056 : // Copy projection information.
5057 50 : GDALGeoTransform gt;
5058 50 : if (poSrcDS->GetGeoTransform(gt) == CE_None)
5059 48 : poDS->SetGeoTransform(gt);
5060 :
5061 50 : const OGRSpatialReference *poSrcSRS = poSrcDS->GetSpatialRef();
5062 50 : if (poSrcSRS)
5063 46 : poDS->SetSpatialRef(poSrcSRS);
5064 :
5065 : // Copy the imagery.
5066 50 : if (!bCreateAux)
5067 : {
5068 50 : const CPLErr eErr = GDALDatasetCopyWholeRaster(
5069 : (GDALDatasetH)poSrcDS, (GDALDatasetH)poDS, nullptr, pfnProgress,
5070 : pProgressData);
5071 :
5072 50 : if (eErr != CE_None)
5073 : {
5074 0 : delete poDS;
5075 0 : return nullptr;
5076 : }
5077 : }
5078 :
5079 : // Do we want to generate statistics and a histogram?
5080 50 : if (CPLFetchBool(papszOptions, "STATISTICS", false))
5081 : {
5082 2 : for (int iBand = 0; iBand < nBandCount; iBand++)
5083 : {
5084 1 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand + 1);
5085 1 : double dfMin = 0.0;
5086 1 : double dfMax = 0.0;
5087 1 : double dfMean = 0.0;
5088 1 : double dfStdDev = 0.0;
5089 1 : char **papszStatsMD = nullptr;
5090 :
5091 : // Statistics
5092 3 : if (poSrcBand->GetStatistics(TRUE, FALSE, &dfMin, &dfMax, &dfMean,
5093 2 : &dfStdDev) == CE_None ||
5094 1 : poSrcBand->ComputeStatistics(TRUE, &dfMin, &dfMax, &dfMean,
5095 : &dfStdDev, pfnProgress,
5096 1 : pProgressData) == CE_None)
5097 : {
5098 1 : CPLString osValue;
5099 :
5100 : papszStatsMD =
5101 1 : CSLSetNameValue(papszStatsMD, "STATISTICS_MINIMUM",
5102 1 : osValue.Printf("%.15g", dfMin));
5103 : papszStatsMD =
5104 1 : CSLSetNameValue(papszStatsMD, "STATISTICS_MAXIMUM",
5105 1 : osValue.Printf("%.15g", dfMax));
5106 1 : papszStatsMD = CSLSetNameValue(papszStatsMD, "STATISTICS_MEAN",
5107 1 : osValue.Printf("%.15g", dfMean));
5108 : papszStatsMD =
5109 1 : CSLSetNameValue(papszStatsMD, "STATISTICS_STDDEV",
5110 1 : osValue.Printf("%.15g", dfStdDev));
5111 : }
5112 :
5113 : // Histogram
5114 1 : int nBuckets = 0;
5115 1 : GUIntBig *panHistogram = nullptr;
5116 :
5117 2 : if (poSrcBand->GetDefaultHistogram(&dfMin, &dfMax, &nBuckets,
5118 : &panHistogram, TRUE, pfnProgress,
5119 1 : pProgressData) == CE_None)
5120 : {
5121 2 : CPLString osValue;
5122 1 : const double dfBinWidth = (dfMax - dfMin) / nBuckets;
5123 :
5124 1 : papszStatsMD = CSLSetNameValue(
5125 : papszStatsMD, "STATISTICS_HISTOMIN",
5126 1 : osValue.Printf("%.15g", dfMin + dfBinWidth * 0.5));
5127 1 : papszStatsMD = CSLSetNameValue(
5128 : papszStatsMD, "STATISTICS_HISTOMAX",
5129 1 : osValue.Printf("%.15g", dfMax - dfBinWidth * 0.5));
5130 : papszStatsMD =
5131 1 : CSLSetNameValue(papszStatsMD, "STATISTICS_HISTONUMBINS",
5132 1 : osValue.Printf("%d", nBuckets));
5133 :
5134 1 : int nBinValuesLen = 0;
5135 : char *pszBinValues =
5136 1 : static_cast<char *>(CPLCalloc(20, nBuckets + 1));
5137 257 : for (int iBin = 0; iBin < nBuckets; iBin++)
5138 : {
5139 :
5140 512 : strcat(pszBinValues + nBinValuesLen,
5141 256 : osValue.Printf(CPL_FRMT_GUIB, panHistogram[iBin]));
5142 256 : strcat(pszBinValues + nBinValuesLen, "|");
5143 256 : nBinValuesLen +=
5144 256 : static_cast<int>(strlen(pszBinValues + nBinValuesLen));
5145 : }
5146 1 : papszStatsMD = CSLSetNameValue(
5147 : papszStatsMD, "STATISTICS_HISTOBINVALUES", pszBinValues);
5148 1 : CPLFree(pszBinValues);
5149 : }
5150 :
5151 1 : CPLFree(panHistogram);
5152 :
5153 1 : if (CSLCount(papszStatsMD) > 0)
5154 1 : HFASetMetadata(poDS->hHFA, iBand + 1, papszStatsMD);
5155 :
5156 1 : CSLDestroy(papszStatsMD);
5157 : }
5158 : }
5159 :
5160 : // All report completion.
5161 50 : if (!pfnProgress(1.0, nullptr, pProgressData))
5162 : {
5163 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
5164 0 : delete poDS;
5165 :
5166 0 : GDALDriver *poHFADriver = (GDALDriver *)GDALGetDriverByName("HFA");
5167 0 : poHFADriver->Delete(pszFilename);
5168 0 : return nullptr;
5169 : }
5170 :
5171 50 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
5172 :
5173 50 : return poDS;
5174 : }
5175 :
5176 : /************************************************************************/
5177 : /* GDALRegister_HFA() */
5178 : /************************************************************************/
5179 :
5180 1935 : void GDALRegister_HFA()
5181 :
5182 : {
5183 1935 : if (GDALGetDriverByName("HFA") != nullptr)
5184 282 : return;
5185 :
5186 1653 : GDALDriver *poDriver = new GDALDriver();
5187 :
5188 1653 : poDriver->SetDescription("HFA");
5189 1653 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
5190 1653 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Erdas Imagine Images (.img)");
5191 1653 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/hfa.html");
5192 1653 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "img");
5193 1653 : poDriver->SetMetadataItem(
5194 : GDAL_DMD_CREATIONDATATYPES,
5195 : "Byte Int8 Int16 UInt16 Int32 UInt32 Float32 Float64 "
5196 1653 : "CFloat32 CFloat64");
5197 :
5198 1653 : poDriver->SetMetadataItem(
5199 : GDAL_DMD_CREATIONOPTIONLIST,
5200 : "<CreationOptionList>"
5201 : " <Option name='BLOCKSIZE' type='integer' description='tile "
5202 : "width/height (32-2048)' default='64'/>"
5203 : " <Option name='USE_SPILL' type='boolean' description='Force use of "
5204 : "spill file'/>"
5205 : " <Option name='COMPRESSED' alias='COMPRESS' type='boolean' "
5206 : "description='compress blocks'/>"
5207 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
5208 : "use Int8) By setting this to SIGNEDBYTE, a new Byte file can be "
5209 : "forced to be written as signed byte'/>"
5210 : " <Option name='AUX' type='boolean' description='Create an .aux "
5211 : "file'/>"
5212 : " <Option name='IGNOREUTM' type='boolean' description='Ignore UTM "
5213 : "when selecting coordinate system - will use Transverse Mercator. Only "
5214 : "used for Create() method'/>"
5215 : " <Option name='NBITS' type='integer' description='Create file with "
5216 : "special sub-byte data type (1/2/4)'/>"
5217 : " <Option name='STATISTICS' type='boolean' description='Generate "
5218 : "statistics and a histogram'/>"
5219 : " <Option name='DEPENDENT_FILE' type='string' description='Name of "
5220 : "dependent file (must not have absolute path)'/>"
5221 : " <Option name='FORCETOPESTRING' type='boolean' description='Force "
5222 : "use of ArcGIS PE String in file instead of Imagine coordinate system "
5223 : "format' default='NO'/>"
5224 : " <Option name='DISABLEPESTRING' type='boolean' description='Disable "
5225 : "use of ArcGIS PE String' default='NO'/>"
5226 1653 : "</CreationOptionList>");
5227 :
5228 1653 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
5229 :
5230 1653 : poDriver->SetMetadataItem(GDAL_DCAP_UPDATE, "YES");
5231 1653 : poDriver->SetMetadataItem(GDAL_DMD_UPDATE_ITEMS,
5232 : "GeoTransform SRS NoData "
5233 : "RasterValues "
5234 1653 : "DatasetMetadata BandMetadata");
5235 :
5236 1653 : poDriver->pfnOpen = HFADataset::Open;
5237 1653 : poDriver->pfnCreate = HFADataset::Create;
5238 1653 : poDriver->pfnCreateCopy = HFADataset::CreateCopy;
5239 1653 : poDriver->pfnIdentify = HFADataset::Identify;
5240 1653 : poDriver->pfnRename = HFADataset::Rename;
5241 1653 : poDriver->pfnCopyFiles = HFADataset::CopyFiles;
5242 :
5243 1653 : GetGDALDriverManager()->RegisterDriver(poDriver);
5244 : }
|
