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