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