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