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