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