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