1 /******************************************************************************
2 * $Id: hfaopen.cpp 28831 2015-04-01 16:46:05Z rouault $
3 *
4 * Project: Erdas Imagine (.img) Translator
5 * Purpose: Supporting functions for HFA (.img) ... main (C callable) API
6 * that is not dependent on GDAL (just CPL).
7 * Author: Frank Warmerdam, warmerdam@pobox.com
8 *
9 ******************************************************************************
10 * Copyright (c) 1999, Intergraph Corporation
11 * Copyright (c) 2007-2011, 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 * hfaopen.cpp
33 *
34 * Supporting routines for reading Erdas Imagine (.imf) Heirarchical
35 * File Architecture files. This is intended to be a library independent
36 * of the GDAL core, but dependent on the Common Portability Library.
37 *
38 */
39
40 #include "hfa_p.h"
41 #include "cpl_conv.h"
42 #include <limits.h>
43 #include <vector>
44
45 CPL_CVSID("$Id: hfaopen.cpp 28831 2015-04-01 16:46:05Z rouault $");
46
47
48 static const char *apszAuxMetadataItems[] = {
49
50 // node/entry field_name metadata_key type
51
52 "Statistics", "dminimum", "STATISTICS_MINIMUM", "Esta_Statistics",
53 "Statistics", "dmaximum", "STATISTICS_MAXIMUM", "Esta_Statistics",
54 "Statistics", "dmean", "STATISTICS_MEAN", "Esta_Statistics",
55 "Statistics", "dmedian", "STATISTICS_MEDIAN", "Esta_Statistics",
56 "Statistics", "dmode", "STATISTICS_MODE", "Esta_Statistics",
57 "Statistics", "dstddev", "STATISTICS_STDDEV", "Esta_Statistics",
58 "HistogramParameters", "lBinFunction.numBins", "STATISTICS_HISTONUMBINS","Eimg_StatisticsParameters830",
59 "HistogramParameters", "dBinFunction.minLimit", "STATISTICS_HISTOMIN", "Eimg_StatisticsParameters830",
60 "HistogramParameters", "dBinFunction.maxLimit", "STATISTICS_HISTOMAX", "Eimg_StatisticsParameters830",
61 "StatisticsParameters", "lSkipFactorX", "STATISTICS_SKIPFACTORX", "",
62 "StatisticsParameters", "lSkipFactorY", "STATISTICS_SKIPFACTORY", "",
63 "StatisticsParameters", "dExcludedValues", "STATISTICS_EXCLUDEDVALUES","",
64 "", "elayerType", "LAYER_TYPE", "",
65 NULL
66 };
67
68
GetHFAAuxMetaDataList()69 const char ** GetHFAAuxMetaDataList()
70 {
71 return apszAuxMetadataItems;
72 }
73
74
75 /************************************************************************/
76 /* HFAGetDictionary() */
77 /************************************************************************/
78
HFAGetDictionary(HFAHandle hHFA)79 static char * HFAGetDictionary( HFAHandle hHFA )
80
81 {
82 int nDictMax = 100;
83 char *pszDictionary = (char *) CPLMalloc(nDictMax);
84 int nDictSize = 0;
85
86 VSIFSeekL( hHFA->fp, hHFA->nDictionaryPos, SEEK_SET );
87
88 while( TRUE )
89 {
90 if( nDictSize >= nDictMax-1 )
91 {
92 nDictMax = nDictSize * 2 + 100;
93 pszDictionary = (char *) CPLRealloc(pszDictionary, nDictMax );
94 }
95
96 if( VSIFReadL( pszDictionary + nDictSize, 1, 1, hHFA->fp ) < 1
97 || pszDictionary[nDictSize] == '\0'
98 || (nDictSize > 2 && pszDictionary[nDictSize-2] == ','
99 && pszDictionary[nDictSize-1] == '.') )
100 break;
101
102 nDictSize++;
103 }
104
105 pszDictionary[nDictSize] = '\0';
106
107
108 return( pszDictionary );
109 }
110
111 /************************************************************************/
112 /* HFAOpen() */
113 /************************************************************************/
114
HFAOpen(const char * pszFilename,const char * pszAccess)115 HFAHandle HFAOpen( const char * pszFilename, const char * pszAccess )
116
117 {
118 VSILFILE *fp;
119 char szHeader[16];
120 HFAInfo_t *psInfo;
121 GUInt32 nHeaderPos;
122
123 /* -------------------------------------------------------------------- */
124 /* Open the file. */
125 /* -------------------------------------------------------------------- */
126 if( EQUAL(pszAccess,"r") || EQUAL(pszAccess,"rb" ) )
127 fp = VSIFOpenL( pszFilename, "rb" );
128 else
129 fp = VSIFOpenL( pszFilename, "r+b" );
130
131 /* should this be changed to use some sort of CPLFOpen() which will
132 set the error? */
133 if( fp == NULL )
134 {
135 CPLError( CE_Failure, CPLE_OpenFailed,
136 "File open of %s failed.",
137 pszFilename );
138
139 return NULL;
140 }
141
142 /* -------------------------------------------------------------------- */
143 /* Read and verify the header. */
144 /* -------------------------------------------------------------------- */
145 if( VSIFReadL( szHeader, 16, 1, fp ) < 1 )
146 {
147 CPLError( CE_Failure, CPLE_AppDefined,
148 "Attempt to read 16 byte header failed for\n%s.",
149 pszFilename );
150
151 return NULL;
152 }
153
154 if( !EQUALN(szHeader,"EHFA_HEADER_TAG",15) )
155 {
156 CPLError( CE_Failure, CPLE_AppDefined,
157 "File %s is not an Imagine HFA file ... header wrong.",
158 pszFilename );
159
160 return NULL;
161 }
162
163 /* -------------------------------------------------------------------- */
164 /* Create the HFAInfo_t */
165 /* -------------------------------------------------------------------- */
166 psInfo = (HFAInfo_t *) CPLCalloc(sizeof(HFAInfo_t),1);
167
168 psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
169 psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
170 psInfo->fp = fp;
171 if( EQUAL(pszAccess,"r") || EQUAL(pszAccess,"rb" ) )
172 psInfo->eAccess = HFA_ReadOnly;
173 else
174 psInfo->eAccess = HFA_Update;
175 psInfo->bTreeDirty = FALSE;
176
177 /* -------------------------------------------------------------------- */
178 /* Where is the header? */
179 /* -------------------------------------------------------------------- */
180 VSIFReadL( &nHeaderPos, sizeof(GInt32), 1, fp );
181 HFAStandard( 4, &nHeaderPos );
182
183 /* -------------------------------------------------------------------- */
184 /* Read the header. */
185 /* -------------------------------------------------------------------- */
186 VSIFSeekL( fp, nHeaderPos, SEEK_SET );
187
188 VSIFReadL( &(psInfo->nVersion), sizeof(GInt32), 1, fp );
189 HFAStandard( 4, &(psInfo->nVersion) );
190
191 VSIFReadL( szHeader, 4, 1, fp ); /* skip freeList */
192
193 VSIFReadL( &(psInfo->nRootPos), sizeof(GInt32), 1, fp );
194 HFAStandard( 4, &(psInfo->nRootPos) );
195
196 VSIFReadL( &(psInfo->nEntryHeaderLength), sizeof(GInt16), 1, fp );
197 HFAStandard( 2, &(psInfo->nEntryHeaderLength) );
198
199 VSIFReadL( &(psInfo->nDictionaryPos), sizeof(GInt32), 1, fp );
200 HFAStandard( 4, &(psInfo->nDictionaryPos) );
201
202 /* -------------------------------------------------------------------- */
203 /* Collect file size. */
204 /* -------------------------------------------------------------------- */
205 VSIFSeekL( fp, 0, SEEK_END );
206 psInfo->nEndOfFile = (GUInt32) VSIFTellL( fp );
207
208 /* -------------------------------------------------------------------- */
209 /* Instantiate the root entry. */
210 /* -------------------------------------------------------------------- */
211 psInfo->poRoot = HFAEntry::New( psInfo, psInfo->nRootPos, NULL, NULL );
212 if( psInfo->poRoot == NULL )
213 {
214 CPLFree(psInfo);
215 return NULL;
216 }
217
218 /* -------------------------------------------------------------------- */
219 /* Read the dictionary */
220 /* -------------------------------------------------------------------- */
221 psInfo->pszDictionary = HFAGetDictionary( psInfo );
222 psInfo->poDictionary = new HFADictionary( psInfo->pszDictionary );
223
224 /* -------------------------------------------------------------------- */
225 /* Collect band definitions. */
226 /* -------------------------------------------------------------------- */
227 HFAParseBandInfo( psInfo );
228
229 return psInfo;
230 }
231
232 /************************************************************************/
233 /* HFACreateDependent() */
234 /* */
235 /* Create a .rrd file for the named file if it does not exist, */
236 /* or return the existing dependent if it already exists. */
237 /************************************************************************/
238
HFACreateDependent(HFAInfo_t * psBase)239 HFAInfo_t *HFACreateDependent( HFAInfo_t *psBase )
240
241 {
242 if( psBase->psDependent != NULL )
243 return psBase->psDependent;
244
245 /* -------------------------------------------------------------------- */
246 /* Create desired RRD filename. */
247 /* -------------------------------------------------------------------- */
248 CPLString oBasename = CPLGetBasename( psBase->pszFilename );
249 CPLString oRRDFilename =
250 CPLFormFilename( psBase->pszPath, oBasename, "rrd" );
251
252 /* -------------------------------------------------------------------- */
253 /* Does this file already exist? If so, re-use it. */
254 /* -------------------------------------------------------------------- */
255 VSILFILE *fp = VSIFOpenL( oRRDFilename, "rb" );
256 if( fp != NULL )
257 {
258 VSIFCloseL( fp );
259 psBase->psDependent = HFAOpen( oRRDFilename, "rb" );
260 }
261
262 /* -------------------------------------------------------------------- */
263 /* Otherwise create it now. */
264 /* -------------------------------------------------------------------- */
265 HFAInfo_t *psDep;
266 psDep = psBase->psDependent = HFACreateLL( oRRDFilename );
267
268 /* -------------------------------------------------------------------- */
269 /* Add the DependentFile node with the pointer back to the */
270 /* parent. When working from an .aux file we really want the */
271 /* .rrd to point back to the original file, not the .aux file. */
272 /* -------------------------------------------------------------------- */
273 HFAEntry *poEntry = psBase->poRoot->GetNamedChild("DependentFile");
274 const char *pszDependentFile = NULL;
275 if( poEntry != NULL )
276 pszDependentFile = poEntry->GetStringField( "dependent.string" );
277 if( pszDependentFile == NULL )
278 pszDependentFile = psBase->pszFilename;
279
280 HFAEntry *poDF = new HFAEntry( psDep, "DependentFile",
281 "Eimg_DependentFile", psDep->poRoot );
282
283 poDF->MakeData( strlen(pszDependentFile) + 50 );
284 poDF->SetPosition();
285 poDF->SetStringField( "dependent.string", pszDependentFile );
286
287 return psDep;
288 }
289
290 /************************************************************************/
291 /* HFAGetDependent() */
292 /************************************************************************/
293
HFAGetDependent(HFAInfo_t * psBase,const char * pszFilename)294 HFAInfo_t *HFAGetDependent( HFAInfo_t *psBase, const char *pszFilename )
295
296 {
297 if( EQUAL(pszFilename,psBase->pszFilename) )
298 return psBase;
299
300 if( psBase->psDependent != NULL )
301 {
302 if( EQUAL(pszFilename,psBase->psDependent->pszFilename) )
303 return psBase->psDependent;
304 else
305 return NULL;
306 }
307
308 /* -------------------------------------------------------------------- */
309 /* Try to open the dependent file. */
310 /* -------------------------------------------------------------------- */
311 char *pszDependent;
312 VSILFILE *fp;
313 const char* pszMode = psBase->eAccess == HFA_Update ? "r+b" : "rb";
314
315 pszDependent = CPLStrdup(
316 CPLFormFilename( psBase->pszPath, pszFilename, NULL ) );
317
318 fp = VSIFOpenL( pszDependent, pszMode );
319 if( fp != NULL )
320 {
321 VSIFCloseL( fp );
322 psBase->psDependent = HFAOpen( pszDependent, pszMode );
323 }
324
325 CPLFree( pszDependent );
326
327 return psBase->psDependent;
328 }
329
330
331 /************************************************************************/
332 /* HFAParseBandInfo() */
333 /* */
334 /* This is used by HFAOpen() and HFACreate() to initialize the */
335 /* band structures. */
336 /************************************************************************/
337
HFAParseBandInfo(HFAInfo_t * psInfo)338 CPLErr HFAParseBandInfo( HFAInfo_t *psInfo )
339
340 {
341 HFAEntry *poNode;
342
343 /* -------------------------------------------------------------------- */
344 /* Find the first band node. */
345 /* -------------------------------------------------------------------- */
346 psInfo->nBands = 0;
347 poNode = psInfo->poRoot->GetChild();
348 while( poNode != NULL )
349 {
350 if( EQUAL(poNode->GetType(),"Eimg_Layer")
351 && poNode->GetIntField("width") > 0
352 && poNode->GetIntField("height") > 0 )
353 {
354 if( psInfo->nBands == 0 )
355 {
356 psInfo->nXSize = poNode->GetIntField("width");
357 psInfo->nYSize = poNode->GetIntField("height");
358 }
359 else if( poNode->GetIntField("width") != psInfo->nXSize
360 || poNode->GetIntField("height") != psInfo->nYSize )
361 {
362 return CE_Failure;
363 }
364
365 psInfo->papoBand = (HFABand **)
366 CPLRealloc(psInfo->papoBand,
367 sizeof(HFABand *) * (psInfo->nBands+1));
368 psInfo->papoBand[psInfo->nBands] = new HFABand( psInfo, poNode );
369 if (psInfo->papoBand[psInfo->nBands]->nWidth == 0)
370 {
371 delete psInfo->papoBand[psInfo->nBands];
372 return CE_Failure;
373 }
374 psInfo->nBands++;
375 }
376
377 poNode = poNode->GetNext();
378 }
379
380 return CE_None;
381 }
382
383 /************************************************************************/
384 /* HFAClose() */
385 /************************************************************************/
386
HFAClose(HFAHandle hHFA)387 void HFAClose( HFAHandle hHFA )
388
389 {
390 int i;
391
392 if( hHFA->eAccess == HFA_Update && (hHFA->bTreeDirty || hHFA->poDictionary->bDictionaryTextDirty) )
393 HFAFlush( hHFA );
394
395 if( hHFA->psDependent != NULL )
396 HFAClose( hHFA->psDependent );
397
398 delete hHFA->poRoot;
399
400 VSIFCloseL( hHFA->fp );
401
402 if( hHFA->poDictionary != NULL )
403 delete hHFA->poDictionary;
404
405 CPLFree( hHFA->pszDictionary );
406 CPLFree( hHFA->pszFilename );
407 CPLFree( hHFA->pszIGEFilename );
408 CPLFree( hHFA->pszPath );
409
410 for( i = 0; i < hHFA->nBands; i++ )
411 {
412 delete hHFA->papoBand[i];
413 }
414
415 CPLFree( hHFA->papoBand );
416
417 if( hHFA->pProParameters != NULL )
418 {
419 Eprj_ProParameters *psProParms = (Eprj_ProParameters *)
420 hHFA->pProParameters;
421
422 CPLFree( psProParms->proExeName );
423 CPLFree( psProParms->proName );
424 CPLFree( psProParms->proSpheroid.sphereName );
425
426 CPLFree( psProParms );
427 }
428
429 if( hHFA->pDatum != NULL )
430 {
431 CPLFree( ((Eprj_Datum *) hHFA->pDatum)->datumname );
432 CPLFree( ((Eprj_Datum *) hHFA->pDatum)->gridname );
433 CPLFree( hHFA->pDatum );
434 }
435
436 if( hHFA->pMapInfo != NULL )
437 {
438 CPLFree( ((Eprj_MapInfo *) hHFA->pMapInfo)->proName );
439 CPLFree( ((Eprj_MapInfo *) hHFA->pMapInfo)->units );
440 CPLFree( hHFA->pMapInfo );
441 }
442
443 CPLFree( hHFA );
444 }
445
446 /************************************************************************/
447 /* HFARemove() */
448 /* Used from HFADelete() function. */
449 /************************************************************************/
450
HFARemove(const char * pszFilename)451 CPLErr HFARemove( const char *pszFilename )
452
453 {
454 VSIStatBufL sStat;
455
456 if( VSIStatL( pszFilename, &sStat ) == 0 && VSI_ISREG( sStat.st_mode ) )
457 {
458 if( VSIUnlink( pszFilename ) == 0 )
459 return CE_None;
460 else
461 {
462 CPLError( CE_Failure, CPLE_AppDefined,
463 "Attempt to unlink %s failed.\n", pszFilename );
464 return CE_Failure;
465 }
466 }
467 else
468 {
469 CPLError( CE_Failure, CPLE_AppDefined,
470 "Unable to delete %s, not a file.\n", pszFilename );
471 return CE_Failure;
472 }
473 }
474
475 /************************************************************************/
476 /* HFADelete() */
477 /************************************************************************/
478
HFADelete(const char * pszFilename)479 CPLErr HFADelete( const char *pszFilename )
480
481 {
482 HFAInfo_t *psInfo = HFAOpen( pszFilename, "rb" );
483 HFAEntry *poDMS = NULL;
484 HFAEntry *poLayer = NULL;
485 HFAEntry *poNode = NULL;
486
487 if( psInfo != NULL )
488 {
489 poNode = psInfo->poRoot->GetChild();
490 while( ( poNode != NULL ) && ( poLayer == NULL ) )
491 {
492 if( EQUAL(poNode->GetType(),"Eimg_Layer") )
493 {
494 poLayer = poNode;
495 }
496 poNode = poNode->GetNext();
497 }
498
499 if( poLayer != NULL )
500 poDMS = poLayer->GetNamedChild( "ExternalRasterDMS" );
501
502 if ( poDMS )
503 {
504 const char *pszRawFilename =
505 poDMS->GetStringField( "fileName.string" );
506
507 if( pszRawFilename != NULL )
508 HFARemove( CPLFormFilename( psInfo->pszPath,
509 pszRawFilename, NULL ) );
510 }
511
512 HFAClose( psInfo );
513 }
514 return HFARemove( pszFilename );
515 }
516
517 /************************************************************************/
518 /* HFAGetRasterInfo() */
519 /************************************************************************/
520
HFAGetRasterInfo(HFAHandle hHFA,int * pnXSize,int * pnYSize,int * pnBands)521 CPLErr HFAGetRasterInfo( HFAHandle hHFA, int * pnXSize, int * pnYSize,
522 int * pnBands )
523
524 {
525 if( pnXSize != NULL )
526 *pnXSize = hHFA->nXSize;
527 if( pnYSize != NULL )
528 *pnYSize = hHFA->nYSize;
529 if( pnBands != NULL )
530 *pnBands = hHFA->nBands;
531 return CE_None;
532 }
533
534 /************************************************************************/
535 /* HFAGetBandInfo() */
536 /************************************************************************/
537
HFAGetBandInfo(HFAHandle hHFA,int nBand,int * pnDataType,int * pnBlockXSize,int * pnBlockYSize,int * pnCompressionType)538 CPLErr HFAGetBandInfo( HFAHandle hHFA, int nBand, int * pnDataType,
539 int * pnBlockXSize, int * pnBlockYSize,
540 int *pnCompressionType )
541
542 {
543 if( nBand < 0 || nBand > hHFA->nBands )
544 {
545 CPLAssert( FALSE );
546 return CE_Failure;
547 }
548
549 HFABand *poBand = hHFA->papoBand[nBand-1];
550
551 if( pnDataType != NULL )
552 *pnDataType = poBand->nDataType;
553
554 if( pnBlockXSize != NULL )
555 *pnBlockXSize = poBand->nBlockXSize;
556
557 if( pnBlockYSize != NULL )
558 *pnBlockYSize = poBand->nBlockYSize;
559
560 /* -------------------------------------------------------------------- */
561 /* Get compression code from RasterDMS. */
562 /* -------------------------------------------------------------------- */
563 if( pnCompressionType != NULL )
564 {
565 HFAEntry *poDMS;
566
567 *pnCompressionType = 0;
568
569 poDMS = poBand->poNode->GetNamedChild( "RasterDMS" );
570
571 if( poDMS != NULL )
572 *pnCompressionType = poDMS->GetIntField( "compressionType" );
573 }
574
575 return( CE_None );
576 }
577
578 /************************************************************************/
579 /* HFAGetBandNoData() */
580 /* */
581 /* returns TRUE if value is set, otherwise FALSE. */
582 /************************************************************************/
583
HFAGetBandNoData(HFAHandle hHFA,int nBand,double * pdfNoData)584 int HFAGetBandNoData( HFAHandle hHFA, int nBand, double *pdfNoData )
585
586 {
587 if( nBand < 0 || nBand > hHFA->nBands )
588 {
589 CPLAssert( FALSE );
590 return CE_Failure;
591 }
592
593 HFABand *poBand = hHFA->papoBand[nBand-1];
594
595 if( !poBand->bNoDataSet && poBand->nOverviews > 0 )
596 {
597 poBand = poBand->papoOverviews[0];
598 if( poBand == NULL )
599 return FALSE;
600 }
601
602 *pdfNoData = poBand->dfNoData;
603 return poBand->bNoDataSet;
604 }
605
606 /************************************************************************/
607 /* HFASetBandNoData() */
608 /* */
609 /* attempts to set a no-data value on the given band */
610 /************************************************************************/
611
HFASetBandNoData(HFAHandle hHFA,int nBand,double dfValue)612 CPLErr HFASetBandNoData( HFAHandle hHFA, int nBand, double dfValue )
613
614 {
615 if ( nBand < 0 || nBand > hHFA->nBands )
616 {
617 CPLAssert( FALSE );
618 return CE_Failure;
619 }
620
621 HFABand *poBand = hHFA->papoBand[nBand - 1];
622
623 return poBand->SetNoDataValue( dfValue );
624 }
625
626 /************************************************************************/
627 /* HFAGetOverviewCount() */
628 /************************************************************************/
629
HFAGetOverviewCount(HFAHandle hHFA,int nBand)630 int HFAGetOverviewCount( HFAHandle hHFA, int nBand )
631
632 {
633 HFABand *poBand;
634
635 if( nBand < 0 || nBand > hHFA->nBands )
636 {
637 CPLAssert( FALSE );
638 return CE_Failure;
639 }
640
641 poBand = hHFA->papoBand[nBand-1];
642 poBand->LoadOverviews();
643
644 return poBand->nOverviews;
645 }
646
647 /************************************************************************/
648 /* HFAGetOverviewInfo() */
649 /************************************************************************/
650
HFAGetOverviewInfo(HFAHandle hHFA,int nBand,int iOverview,int * pnXSize,int * pnYSize,int * pnBlockXSize,int * pnBlockYSize,int * pnHFADataType)651 CPLErr HFAGetOverviewInfo( HFAHandle hHFA, int nBand, int iOverview,
652 int * pnXSize, int * pnYSize,
653 int * pnBlockXSize, int * pnBlockYSize,
654 int * pnHFADataType )
655
656 {
657 HFABand *poBand;
658
659 if( nBand < 0 || nBand > hHFA->nBands )
660 {
661 CPLAssert( FALSE );
662 return CE_Failure;
663 }
664
665 poBand = hHFA->papoBand[nBand-1];
666 poBand->LoadOverviews();
667
668 if( iOverview < 0 || iOverview >= poBand->nOverviews )
669 {
670 CPLAssert( FALSE );
671 return CE_Failure;
672 }
673 poBand = poBand->papoOverviews[iOverview];
674 if( poBand == NULL )
675 {
676 return CE_Failure;
677 }
678
679 if( pnXSize != NULL )
680 *pnXSize = poBand->nWidth;
681
682 if( pnYSize != NULL )
683 *pnYSize = poBand->nHeight;
684
685 if( pnBlockXSize != NULL )
686 *pnBlockXSize = poBand->nBlockXSize;
687
688 if( pnBlockYSize != NULL )
689 *pnBlockYSize = poBand->nBlockYSize;
690
691 if( pnHFADataType != NULL )
692 *pnHFADataType = poBand->nDataType;
693
694 return( CE_None );
695 }
696
697 /************************************************************************/
698 /* HFAGetRasterBlock() */
699 /************************************************************************/
700
HFAGetRasterBlock(HFAHandle hHFA,int nBand,int nXBlock,int nYBlock,void * pData)701 CPLErr HFAGetRasterBlock( HFAHandle hHFA, int nBand,
702 int nXBlock, int nYBlock, void * pData )
703
704 {
705 return HFAGetRasterBlockEx(hHFA, nBand, nXBlock, nYBlock, pData, -1);
706 }
707
708 /************************************************************************/
709 /* HFAGetRasterBlockEx() */
710 /************************************************************************/
711
HFAGetRasterBlockEx(HFAHandle hHFA,int nBand,int nXBlock,int nYBlock,void * pData,int nDataSize)712 CPLErr HFAGetRasterBlockEx( HFAHandle hHFA, int nBand,
713 int nXBlock, int nYBlock, void * pData, int nDataSize )
714
715 {
716 if( nBand < 1 || nBand > hHFA->nBands )
717 return CE_Failure;
718
719 return( hHFA->papoBand[nBand-1]->GetRasterBlock(nXBlock,nYBlock,pData,nDataSize) );
720 }
721
722 /************************************************************************/
723 /* HFAGetOverviewRasterBlock() */
724 /************************************************************************/
725
HFAGetOverviewRasterBlock(HFAHandle hHFA,int nBand,int iOverview,int nXBlock,int nYBlock,void * pData)726 CPLErr HFAGetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
727 int nXBlock, int nYBlock, void * pData )
728
729 {
730 return HFAGetOverviewRasterBlockEx(hHFA, nBand, iOverview, nXBlock, nYBlock, pData, -1);
731 }
732
733 /************************************************************************/
734 /* HFAGetOverviewRasterBlockEx() */
735 /************************************************************************/
736
HFAGetOverviewRasterBlockEx(HFAHandle hHFA,int nBand,int iOverview,int nXBlock,int nYBlock,void * pData,int nDataSize)737 CPLErr HFAGetOverviewRasterBlockEx( HFAHandle hHFA, int nBand, int iOverview,
738 int nXBlock, int nYBlock, void * pData, int nDataSize )
739
740 {
741 if( nBand < 1 || nBand > hHFA->nBands )
742 return CE_Failure;
743
744 if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand-1]->nOverviews )
745 return CE_Failure;
746
747 return( hHFA->papoBand[nBand-1]->papoOverviews[iOverview]->
748 GetRasterBlock(nXBlock,nYBlock,pData, nDataSize) );
749 }
750
751 /************************************************************************/
752 /* HFASetRasterBlock() */
753 /************************************************************************/
754
HFASetRasterBlock(HFAHandle hHFA,int nBand,int nXBlock,int nYBlock,void * pData)755 CPLErr HFASetRasterBlock( HFAHandle hHFA, int nBand,
756 int nXBlock, int nYBlock, void * pData )
757
758 {
759 if( nBand < 1 || nBand > hHFA->nBands )
760 return CE_Failure;
761
762 return( hHFA->papoBand[nBand-1]->SetRasterBlock(nXBlock,nYBlock,pData) );
763 }
764
765 /************************************************************************/
766 /* HFASetRasterBlock() */
767 /************************************************************************/
768
HFASetOverviewRasterBlock(HFAHandle hHFA,int nBand,int iOverview,int nXBlock,int nYBlock,void * pData)769 CPLErr HFASetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
770 int nXBlock, int nYBlock, void * pData )
771
772 {
773 if( nBand < 1 || nBand > hHFA->nBands )
774 return CE_Failure;
775
776 if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand-1]->nOverviews )
777 return CE_Failure;
778
779 return( hHFA->papoBand[nBand-1]->papoOverviews[iOverview]->
780 SetRasterBlock(nXBlock,nYBlock,pData) );
781 }
782
783 /************************************************************************/
784 /* HFAGetBandName() */
785 /************************************************************************/
786
HFAGetBandName(HFAHandle hHFA,int nBand)787 const char * HFAGetBandName( HFAHandle hHFA, int nBand )
788 {
789 if( nBand < 1 || nBand > hHFA->nBands )
790 return "";
791
792 return( hHFA->papoBand[nBand-1]->GetBandName() );
793 }
794
795 /************************************************************************/
796 /* HFASetBandName() */
797 /************************************************************************/
798
HFASetBandName(HFAHandle hHFA,int nBand,const char * pszName)799 void HFASetBandName( HFAHandle hHFA, int nBand, const char *pszName )
800 {
801 if( nBand < 1 || nBand > hHFA->nBands )
802 return;
803
804 hHFA->papoBand[nBand-1]->SetBandName( pszName );
805 }
806
807 /************************************************************************/
808 /* HFAGetDataTypeBits() */
809 /************************************************************************/
810
HFAGetDataTypeBits(int nDataType)811 int HFAGetDataTypeBits( int nDataType )
812
813 {
814 switch( nDataType )
815 {
816 case EPT_u1:
817 return 1;
818
819 case EPT_u2:
820 return 2;
821
822 case EPT_u4:
823 return 4;
824
825 case EPT_u8:
826 case EPT_s8:
827 return 8;
828
829 case EPT_u16:
830 case EPT_s16:
831 return 16;
832
833 case EPT_u32:
834 case EPT_s32:
835 case EPT_f32:
836 return 32;
837
838 case EPT_f64:
839 case EPT_c64:
840 return 64;
841
842 case EPT_c128:
843 return 128;
844 }
845
846 return 0;
847 }
848
849 /************************************************************************/
850 /* HFAGetDataTypeName() */
851 /************************************************************************/
852
HFAGetDataTypeName(int nDataType)853 const char *HFAGetDataTypeName( int nDataType )
854
855 {
856 switch( nDataType )
857 {
858 case EPT_u1:
859 return "u1";
860
861 case EPT_u2:
862 return "u2";
863
864 case EPT_u4:
865 return "u4";
866
867 case EPT_u8:
868 return "u8";
869
870 case EPT_s8:
871 return "s8";
872
873 case EPT_u16:
874 return "u16";
875
876 case EPT_s16:
877 return "s16";
878
879 case EPT_u32:
880 return "u32";
881
882 case EPT_s32:
883 return "s32";
884
885 case EPT_f32:
886 return "f32";
887
888 case EPT_f64:
889 return "f64";
890
891 case EPT_c64:
892 return "c64";
893
894 case EPT_c128:
895 return "c128";
896
897 default:
898 return "unknown";
899 }
900 }
901
902 /************************************************************************/
903 /* HFAGetMapInfo() */
904 /************************************************************************/
905
HFAGetMapInfo(HFAHandle hHFA)906 const Eprj_MapInfo *HFAGetMapInfo( HFAHandle hHFA )
907
908 {
909 HFAEntry *poMIEntry;
910 Eprj_MapInfo *psMapInfo;
911
912 if( hHFA->nBands < 1 )
913 return NULL;
914
915 /* -------------------------------------------------------------------- */
916 /* Do we already have it? */
917 /* -------------------------------------------------------------------- */
918 if( hHFA->pMapInfo != NULL )
919 return( (Eprj_MapInfo *) hHFA->pMapInfo );
920
921 /* -------------------------------------------------------------------- */
922 /* Get the HFA node. If we don't find it under the usual name */
923 /* we search for any node of the right type (#3338). */
924 /* -------------------------------------------------------------------- */
925 poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild( "Map_Info" );
926 if( poMIEntry == NULL )
927 {
928 HFAEntry *poChild;
929 for( poChild = hHFA->papoBand[0]->poNode->GetChild();
930 poChild != NULL && poMIEntry == NULL;
931 poChild = poChild->GetNext() )
932 {
933 if( EQUAL(poChild->GetType(),"Eprj_MapInfo") )
934 poMIEntry = poChild;
935 }
936 }
937
938 if( poMIEntry == NULL )
939 {
940 return NULL;
941 }
942
943 /* -------------------------------------------------------------------- */
944 /* Allocate the structure. */
945 /* -------------------------------------------------------------------- */
946 psMapInfo = (Eprj_MapInfo *) CPLCalloc(sizeof(Eprj_MapInfo),1);
947
948 /* -------------------------------------------------------------------- */
949 /* Fetch the fields. */
950 /* -------------------------------------------------------------------- */
951 CPLErr eErr;
952
953 psMapInfo->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
954
955 psMapInfo->upperLeftCenter.x =
956 poMIEntry->GetDoubleField("upperLeftCenter.x");
957 psMapInfo->upperLeftCenter.y =
958 poMIEntry->GetDoubleField("upperLeftCenter.y");
959
960 psMapInfo->lowerRightCenter.x =
961 poMIEntry->GetDoubleField("lowerRightCenter.x");
962 psMapInfo->lowerRightCenter.y =
963 poMIEntry->GetDoubleField("lowerRightCenter.y");
964
965 psMapInfo->pixelSize.width =
966 poMIEntry->GetDoubleField("pixelSize.width",&eErr);
967 psMapInfo->pixelSize.height =
968 poMIEntry->GetDoubleField("pixelSize.height",&eErr);
969
970 // The following is basically a hack to get files with
971 // non-standard MapInfo's that misname the pixelSize fields. (#3338)
972 if( eErr != CE_None )
973 {
974 psMapInfo->pixelSize.width =
975 poMIEntry->GetDoubleField("pixelSize.x");
976 psMapInfo->pixelSize.height =
977 poMIEntry->GetDoubleField("pixelSize.y");
978 }
979
980 psMapInfo->units = CPLStrdup(poMIEntry->GetStringField("units"));
981
982 hHFA->pMapInfo = (void *) psMapInfo;
983
984 return psMapInfo;
985 }
986
987 /************************************************************************/
988 /* HFAInvGeoTransform() */
989 /************************************************************************/
990
HFAInvGeoTransform(double * gt_in,double * gt_out)991 static int HFAInvGeoTransform( double *gt_in, double *gt_out )
992
993 {
994 double det, inv_det;
995
996 /* we assume a 3rd row that is [1 0 0] */
997
998 /* Compute determinate */
999
1000 det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
1001
1002 if( fabs(det) < 0.000000000000001 )
1003 return 0;
1004
1005 inv_det = 1.0 / det;
1006
1007 /* compute adjoint, and devide by determinate */
1008
1009 gt_out[1] = gt_in[5] * inv_det;
1010 gt_out[4] = -gt_in[4] * inv_det;
1011
1012 gt_out[2] = -gt_in[2] * inv_det;
1013 gt_out[5] = gt_in[1] * inv_det;
1014
1015 gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
1016 gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
1017
1018 return 1;
1019 }
1020
1021 /************************************************************************/
1022 /* HFAGetGeoTransform() */
1023 /************************************************************************/
1024
HFAGetGeoTransform(HFAHandle hHFA,double * padfGeoTransform)1025 int HFAGetGeoTransform( HFAHandle hHFA, double *padfGeoTransform )
1026
1027 {
1028 const Eprj_MapInfo *psMapInfo = HFAGetMapInfo( hHFA );
1029
1030 padfGeoTransform[0] = 0.0;
1031 padfGeoTransform[1] = 1.0;
1032 padfGeoTransform[2] = 0.0;
1033 padfGeoTransform[3] = 0.0;
1034 padfGeoTransform[4] = 0.0;
1035 padfGeoTransform[5] = 1.0;
1036
1037 /* -------------------------------------------------------------------- */
1038 /* Simple (north up) MapInfo approach. */
1039 /* -------------------------------------------------------------------- */
1040 if( psMapInfo != NULL )
1041 {
1042 padfGeoTransform[0] = psMapInfo->upperLeftCenter.x
1043 - psMapInfo->pixelSize.width*0.5;
1044 padfGeoTransform[1] = psMapInfo->pixelSize.width;
1045 if(padfGeoTransform[1] == 0.0)
1046 padfGeoTransform[1] = 1.0;
1047 padfGeoTransform[2] = 0.0;
1048 if( psMapInfo->upperLeftCenter.y >= psMapInfo->lowerRightCenter.y )
1049 padfGeoTransform[5] = - psMapInfo->pixelSize.height;
1050 else
1051 padfGeoTransform[5] = psMapInfo->pixelSize.height;
1052 if(padfGeoTransform[5] == 0.0)
1053 padfGeoTransform[5] = 1.0;
1054
1055 padfGeoTransform[3] = psMapInfo->upperLeftCenter.y
1056 - padfGeoTransform[5]*0.5;
1057 padfGeoTransform[4] = 0.0;
1058
1059 // special logic to fixup odd angular units.
1060 if( EQUAL(psMapInfo->units,"ds") )
1061 {
1062 padfGeoTransform[0] /= 3600.0;
1063 padfGeoTransform[1] /= 3600.0;
1064 padfGeoTransform[2] /= 3600.0;
1065 padfGeoTransform[3] /= 3600.0;
1066 padfGeoTransform[4] /= 3600.0;
1067 padfGeoTransform[5] /= 3600.0;
1068 }
1069
1070 return TRUE;
1071 }
1072
1073 /* -------------------------------------------------------------------- */
1074 /* Try for a MapToPixelXForm affine polynomial supporting */
1075 /* rotated and sheared affine transformations. */
1076 /* -------------------------------------------------------------------- */
1077 if( hHFA->nBands == 0 )
1078 return FALSE;
1079
1080 HFAEntry *poXForm0 =
1081 hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm.XForm0" );
1082
1083 if( poXForm0 == NULL )
1084 return FALSE;
1085
1086 if( poXForm0->GetIntField( "order" ) != 1
1087 || poXForm0->GetIntField( "numdimtransform" ) != 2
1088 || poXForm0->GetIntField( "numdimpolynomial" ) != 2
1089 || poXForm0->GetIntField( "termcount" ) != 3 )
1090 return FALSE;
1091
1092 // Verify that there aren't any further xform steps.
1093 if( hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm.XForm1" )
1094 != NULL )
1095 return FALSE;
1096
1097 // we should check that the exponent list is 0 0 1 0 0 1 but
1098 // we don't because we are lazy
1099
1100 // fetch geotransform values.
1101 double adfXForm[6];
1102
1103 adfXForm[0] = poXForm0->GetDoubleField( "polycoefvector[0]" );
1104 adfXForm[1] = poXForm0->GetDoubleField( "polycoefmtx[0]" );
1105 adfXForm[4] = poXForm0->GetDoubleField( "polycoefmtx[1]" );
1106 adfXForm[3] = poXForm0->GetDoubleField( "polycoefvector[1]" );
1107 adfXForm[2] = poXForm0->GetDoubleField( "polycoefmtx[2]" );
1108 adfXForm[5] = poXForm0->GetDoubleField( "polycoefmtx[3]" );
1109
1110 // invert
1111
1112 HFAInvGeoTransform( adfXForm, padfGeoTransform );
1113
1114 // Adjust origin from center of top left pixel to top left corner
1115 // of top left pixel.
1116
1117 padfGeoTransform[0] -= padfGeoTransform[1] * 0.5;
1118 padfGeoTransform[0] -= padfGeoTransform[2] * 0.5;
1119 padfGeoTransform[3] -= padfGeoTransform[4] * 0.5;
1120 padfGeoTransform[3] -= padfGeoTransform[5] * 0.5;
1121
1122 return TRUE;
1123 }
1124
1125 /************************************************************************/
1126 /* HFASetMapInfo() */
1127 /************************************************************************/
1128
HFASetMapInfo(HFAHandle hHFA,const Eprj_MapInfo * poMapInfo)1129 CPLErr HFASetMapInfo( HFAHandle hHFA, const Eprj_MapInfo *poMapInfo )
1130
1131 {
1132 /* -------------------------------------------------------------------- */
1133 /* Loop over bands, setting information on each one. */
1134 /* -------------------------------------------------------------------- */
1135 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1136 {
1137 HFAEntry *poMIEntry;
1138
1139 /* -------------------------------------------------------------------- */
1140 /* Create a new Map_Info if there isn't one present already. */
1141 /* -------------------------------------------------------------------- */
1142 poMIEntry = hHFA->papoBand[iBand]->poNode->GetNamedChild( "Map_Info" );
1143 if( poMIEntry == NULL )
1144 {
1145 poMIEntry = new HFAEntry( hHFA, "Map_Info", "Eprj_MapInfo",
1146 hHFA->papoBand[iBand]->poNode );
1147 }
1148
1149 poMIEntry->MarkDirty();
1150
1151 /* -------------------------------------------------------------------- */
1152 /* Ensure we have enough space for all the data. */
1153 /* -------------------------------------------------------------------- */
1154 int nSize;
1155 GByte *pabyData;
1156
1157 nSize = 48 + 40
1158 + strlen(poMapInfo->proName) + 1
1159 + strlen(poMapInfo->units) + 1;
1160
1161 pabyData = poMIEntry->MakeData( nSize );
1162 memset( pabyData, 0, nSize );
1163
1164 poMIEntry->SetPosition();
1165
1166 /* -------------------------------------------------------------------- */
1167 /* Write the various fields. */
1168 /* -------------------------------------------------------------------- */
1169 poMIEntry->SetStringField( "proName", poMapInfo->proName );
1170
1171 poMIEntry->SetDoubleField( "upperLeftCenter.x",
1172 poMapInfo->upperLeftCenter.x );
1173 poMIEntry->SetDoubleField( "upperLeftCenter.y",
1174 poMapInfo->upperLeftCenter.y );
1175
1176 poMIEntry->SetDoubleField( "lowerRightCenter.x",
1177 poMapInfo->lowerRightCenter.x );
1178 poMIEntry->SetDoubleField( "lowerRightCenter.y",
1179 poMapInfo->lowerRightCenter.y );
1180
1181 poMIEntry->SetDoubleField( "pixelSize.width",
1182 poMapInfo->pixelSize.width );
1183 poMIEntry->SetDoubleField( "pixelSize.height",
1184 poMapInfo->pixelSize.height );
1185
1186 poMIEntry->SetStringField( "units", poMapInfo->units );
1187 }
1188
1189 return CE_None;
1190 }
1191
1192 /************************************************************************/
1193 /* HFAGetPEString() */
1194 /* */
1195 /* Some files have a ProjectionX node contining the ESRI style */
1196 /* PE_STRING. This function allows fetching from it. */
1197 /************************************************************************/
1198
HFAGetPEString(HFAHandle hHFA)1199 char *HFAGetPEString( HFAHandle hHFA )
1200
1201 {
1202 if( hHFA->nBands == 0 )
1203 return NULL;
1204
1205 /* -------------------------------------------------------------------- */
1206 /* Get the HFA node. */
1207 /* -------------------------------------------------------------------- */
1208 HFAEntry *poProX;
1209
1210 poProX = hHFA->papoBand[0]->poNode->GetNamedChild( "ProjectionX" );
1211 if( poProX == NULL )
1212 return NULL;
1213
1214 const char *pszType = poProX->GetStringField( "projection.type.string" );
1215 if( pszType == NULL || !EQUAL(pszType,"PE_COORDSYS") )
1216 return NULL;
1217
1218 /* -------------------------------------------------------------------- */
1219 /* Use a gross hack to scan ahead to the actual projection */
1220 /* string. We do it this way because we don't have general */
1221 /* handling for MIFObjects. */
1222 /* -------------------------------------------------------------------- */
1223 GByte *pabyData = poProX->GetData();
1224 int nDataSize = poProX->GetDataSize();
1225
1226 while( nDataSize > 10
1227 && !EQUALN((const char *) pabyData,"PE_COORDSYS,.",13) ) {
1228 pabyData++;
1229 nDataSize--;
1230 }
1231
1232 if( nDataSize < 31 )
1233 return NULL;
1234
1235 /* -------------------------------------------------------------------- */
1236 /* Skip ahead to the actual string. */
1237 /* -------------------------------------------------------------------- */
1238 pabyData += 30;
1239 nDataSize -= 30;
1240
1241 return CPLStrdup( (const char *) pabyData );
1242 }
1243
1244 /************************************************************************/
1245 /* HFASetPEString() */
1246 /************************************************************************/
1247
HFASetPEString(HFAHandle hHFA,const char * pszPEString)1248 CPLErr HFASetPEString( HFAHandle hHFA, const char *pszPEString )
1249
1250 {
1251 /* -------------------------------------------------------------------- */
1252 /* Loop over bands, setting information on each one. */
1253 /* -------------------------------------------------------------------- */
1254 int iBand;
1255
1256 for( iBand = 0; iBand < hHFA->nBands; iBand++ )
1257 {
1258 HFAEntry *poProX;
1259
1260 /* -------------------------------------------------------------------- */
1261 /* Verify we don't already have the node, since update-in-place */
1262 /* is likely to be more complicated. */
1263 /* -------------------------------------------------------------------- */
1264 poProX = hHFA->papoBand[iBand]->poNode->GetNamedChild( "ProjectionX" );
1265
1266 /* -------------------------------------------------------------------- */
1267 /* If we are setting an empty string then a missing entry is */
1268 /* equivelent. */
1269 /* -------------------------------------------------------------------- */
1270 if( strlen(pszPEString) == 0 && poProX == NULL )
1271 continue;
1272
1273 /* -------------------------------------------------------------------- */
1274 /* Create the node. */
1275 /* -------------------------------------------------------------------- */
1276 if( poProX == NULL )
1277 {
1278 poProX = new HFAEntry( hHFA, "ProjectionX","Eprj_MapProjection842",
1279 hHFA->papoBand[iBand]->poNode );
1280 if( poProX == NULL || poProX->GetTypeObject() == NULL )
1281 return CE_Failure;
1282 }
1283
1284 /* -------------------------------------------------------------------- */
1285 /* Prepare the data area with some extra space just in case. */
1286 /* -------------------------------------------------------------------- */
1287 GByte *pabyData = poProX->MakeData( 700 + strlen(pszPEString) );
1288 if( !pabyData )
1289 return CE_Failure;
1290
1291 memset( pabyData, 0, 250+strlen(pszPEString) );
1292
1293 poProX->SetPosition();
1294
1295 poProX->SetStringField( "projection.type.string", "PE_COORDSYS" );
1296 poProX->SetStringField( "projection.MIFDictionary.string",
1297 "{0:pcstring,}Emif_String,{1:x{0:pcstring,}Emif_String,coordSys,}PE_COORDSYS,." );
1298
1299 /* -------------------------------------------------------------------- */
1300 /* Use a gross hack to scan ahead to the actual projection */
1301 /* string. We do it this way because we don't have general */
1302 /* handling for MIFObjects. */
1303 /* -------------------------------------------------------------------- */
1304 pabyData = poProX->GetData();
1305 int nDataSize = poProX->GetDataSize();
1306 GUInt32 iOffset = poProX->GetDataPos();
1307 GUInt32 nSize;
1308
1309 while( nDataSize > 10
1310 && !EQUALN((const char *) pabyData,"PE_COORDSYS,.",13) ) {
1311 pabyData++;
1312 nDataSize--;
1313 iOffset++;
1314 }
1315
1316 CPLAssert( nDataSize > (int) strlen(pszPEString) + 10 );
1317
1318 pabyData += 14;
1319 iOffset += 14;
1320
1321 /* -------------------------------------------------------------------- */
1322 /* Set the size and offset of the mifobject. */
1323 /* -------------------------------------------------------------------- */
1324 iOffset += 8;
1325
1326 nSize = strlen(pszPEString) + 9;
1327
1328 HFAStandard( 4, &nSize );
1329 memcpy( pabyData, &nSize, 4 );
1330 pabyData += 4;
1331
1332 HFAStandard( 4, &iOffset );
1333 memcpy( pabyData, &iOffset, 4 );
1334 pabyData += 4;
1335
1336 /* -------------------------------------------------------------------- */
1337 /* Set the size and offset of the string value. */
1338 /* -------------------------------------------------------------------- */
1339 nSize = strlen(pszPEString) + 1;
1340
1341 HFAStandard( 4, &nSize );
1342 memcpy( pabyData, &nSize, 4 );
1343 pabyData += 4;
1344
1345 iOffset = 8;
1346 HFAStandard( 4, &iOffset );
1347 memcpy( pabyData, &iOffset, 4 );
1348 pabyData += 4;
1349
1350 /* -------------------------------------------------------------------- */
1351 /* Place the string itself. */
1352 /* -------------------------------------------------------------------- */
1353 memcpy( pabyData, pszPEString, strlen(pszPEString)+1 );
1354
1355 poProX->SetStringField( "title.string", "PE" );
1356 }
1357
1358 return CE_None;
1359 }
1360
1361 /************************************************************************/
1362 /* HFAGetProParameters() */
1363 /************************************************************************/
1364
HFAGetProParameters(HFAHandle hHFA)1365 const Eprj_ProParameters *HFAGetProParameters( HFAHandle hHFA )
1366
1367 {
1368 HFAEntry *poMIEntry;
1369 Eprj_ProParameters *psProParms;
1370 int i;
1371
1372 if( hHFA->nBands < 1 )
1373 return NULL;
1374
1375 /* -------------------------------------------------------------------- */
1376 /* Do we already have it? */
1377 /* -------------------------------------------------------------------- */
1378 if( hHFA->pProParameters != NULL )
1379 return( (Eprj_ProParameters *) hHFA->pProParameters );
1380
1381 /* -------------------------------------------------------------------- */
1382 /* Get the HFA node. */
1383 /* -------------------------------------------------------------------- */
1384 poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild( "Projection" );
1385 if( poMIEntry == NULL )
1386 return NULL;
1387
1388 /* -------------------------------------------------------------------- */
1389 /* Allocate the structure. */
1390 /* -------------------------------------------------------------------- */
1391 psProParms = (Eprj_ProParameters *)CPLCalloc(sizeof(Eprj_ProParameters),1);
1392
1393 /* -------------------------------------------------------------------- */
1394 /* Fetch the fields. */
1395 /* -------------------------------------------------------------------- */
1396 psProParms->proType = (Eprj_ProType) poMIEntry->GetIntField("proType");
1397 psProParms->proNumber = poMIEntry->GetIntField("proNumber");
1398 psProParms->proExeName =CPLStrdup(poMIEntry->GetStringField("proExeName"));
1399 psProParms->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
1400 psProParms->proZone = poMIEntry->GetIntField("proZone");
1401
1402 for( i = 0; i < 15; i++ )
1403 {
1404 char szFieldName[40];
1405
1406 sprintf( szFieldName, "proParams[%d]", i );
1407 psProParms->proParams[i] = poMIEntry->GetDoubleField(szFieldName);
1408 }
1409
1410 psProParms->proSpheroid.sphereName =
1411 CPLStrdup(poMIEntry->GetStringField("proSpheroid.sphereName"));
1412 psProParms->proSpheroid.a = poMIEntry->GetDoubleField("proSpheroid.a");
1413 psProParms->proSpheroid.b = poMIEntry->GetDoubleField("proSpheroid.b");
1414 psProParms->proSpheroid.eSquared =
1415 poMIEntry->GetDoubleField("proSpheroid.eSquared");
1416 psProParms->proSpheroid.radius =
1417 poMIEntry->GetDoubleField("proSpheroid.radius");
1418
1419 hHFA->pProParameters = (void *) psProParms;
1420
1421 return psProParms;
1422 }
1423
1424 /************************************************************************/
1425 /* HFASetProParameters() */
1426 /************************************************************************/
1427
HFASetProParameters(HFAHandle hHFA,const Eprj_ProParameters * poPro)1428 CPLErr HFASetProParameters( HFAHandle hHFA, const Eprj_ProParameters *poPro )
1429
1430 {
1431 /* -------------------------------------------------------------------- */
1432 /* Loop over bands, setting information on each one. */
1433 /* -------------------------------------------------------------------- */
1434 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1435 {
1436 HFAEntry *poMIEntry;
1437
1438 /* -------------------------------------------------------------------- */
1439 /* Create a new Projection if there isn't one present already. */
1440 /* -------------------------------------------------------------------- */
1441 poMIEntry = hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1442 if( poMIEntry == NULL )
1443 {
1444 poMIEntry = new HFAEntry( hHFA, "Projection","Eprj_ProParameters",
1445 hHFA->papoBand[iBand]->poNode );
1446 }
1447
1448 poMIEntry->MarkDirty();
1449
1450 /* -------------------------------------------------------------------- */
1451 /* Ensure we have enough space for all the data. */
1452 /* -------------------------------------------------------------------- */
1453 int nSize;
1454 GByte *pabyData;
1455
1456 nSize = 34 + 15 * 8
1457 + 8 + strlen(poPro->proName) + 1
1458 + 32 + 8 + strlen(poPro->proSpheroid.sphereName) + 1;
1459
1460 if( poPro->proExeName != NULL )
1461 nSize += strlen(poPro->proExeName) + 1;
1462
1463 pabyData = poMIEntry->MakeData( nSize );
1464 if(!pabyData)
1465 return CE_Failure;
1466
1467 poMIEntry->SetPosition();
1468
1469 // Initialize the whole thing to zeros for a clean start.
1470 memset( poMIEntry->GetData(), 0, poMIEntry->GetDataSize() );
1471
1472 /* -------------------------------------------------------------------- */
1473 /* Write the various fields. */
1474 /* -------------------------------------------------------------------- */
1475 poMIEntry->SetIntField( "proType", poPro->proType );
1476
1477 poMIEntry->SetIntField( "proNumber", poPro->proNumber );
1478
1479 poMIEntry->SetStringField( "proExeName", poPro->proExeName );
1480 poMIEntry->SetStringField( "proName", poPro->proName );
1481 poMIEntry->SetIntField( "proZone", poPro->proZone );
1482 poMIEntry->SetDoubleField( "proParams[0]", poPro->proParams[0] );
1483 poMIEntry->SetDoubleField( "proParams[1]", poPro->proParams[1] );
1484 poMIEntry->SetDoubleField( "proParams[2]", poPro->proParams[2] );
1485 poMIEntry->SetDoubleField( "proParams[3]", poPro->proParams[3] );
1486 poMIEntry->SetDoubleField( "proParams[4]", poPro->proParams[4] );
1487 poMIEntry->SetDoubleField( "proParams[5]", poPro->proParams[5] );
1488 poMIEntry->SetDoubleField( "proParams[6]", poPro->proParams[6] );
1489 poMIEntry->SetDoubleField( "proParams[7]", poPro->proParams[7] );
1490 poMIEntry->SetDoubleField( "proParams[8]", poPro->proParams[8] );
1491 poMIEntry->SetDoubleField( "proParams[9]", poPro->proParams[9] );
1492 poMIEntry->SetDoubleField( "proParams[10]", poPro->proParams[10] );
1493 poMIEntry->SetDoubleField( "proParams[11]", poPro->proParams[11] );
1494 poMIEntry->SetDoubleField( "proParams[12]", poPro->proParams[12] );
1495 poMIEntry->SetDoubleField( "proParams[13]", poPro->proParams[13] );
1496 poMIEntry->SetDoubleField( "proParams[14]", poPro->proParams[14] );
1497 poMIEntry->SetStringField( "proSpheroid.sphereName",
1498 poPro->proSpheroid.sphereName );
1499 poMIEntry->SetDoubleField( "proSpheroid.a",
1500 poPro->proSpheroid.a );
1501 poMIEntry->SetDoubleField( "proSpheroid.b",
1502 poPro->proSpheroid.b );
1503 poMIEntry->SetDoubleField( "proSpheroid.eSquared",
1504 poPro->proSpheroid.eSquared );
1505 poMIEntry->SetDoubleField( "proSpheroid.radius",
1506 poPro->proSpheroid.radius );
1507 }
1508
1509 return CE_None;
1510 }
1511
1512 /************************************************************************/
1513 /* HFAGetDatum() */
1514 /************************************************************************/
1515
HFAGetDatum(HFAHandle hHFA)1516 const Eprj_Datum *HFAGetDatum( HFAHandle hHFA )
1517
1518 {
1519 HFAEntry *poMIEntry;
1520 Eprj_Datum *psDatum;
1521 int i;
1522
1523 if( hHFA->nBands < 1 )
1524 return NULL;
1525
1526 /* -------------------------------------------------------------------- */
1527 /* Do we already have it? */
1528 /* -------------------------------------------------------------------- */
1529 if( hHFA->pDatum != NULL )
1530 return( (Eprj_Datum *) hHFA->pDatum );
1531
1532 /* -------------------------------------------------------------------- */
1533 /* Get the HFA node. */
1534 /* -------------------------------------------------------------------- */
1535 poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild( "Projection.Datum" );
1536 if( poMIEntry == NULL )
1537 return NULL;
1538
1539 /* -------------------------------------------------------------------- */
1540 /* Allocate the structure. */
1541 /* -------------------------------------------------------------------- */
1542 psDatum = (Eprj_Datum *) CPLCalloc(sizeof(Eprj_Datum),1);
1543
1544 /* -------------------------------------------------------------------- */
1545 /* Fetch the fields. */
1546 /* -------------------------------------------------------------------- */
1547 psDatum->datumname = CPLStrdup(poMIEntry->GetStringField("datumname"));
1548 psDatum->type = (Eprj_DatumType) poMIEntry->GetIntField("type");
1549
1550 for( i = 0; i < 7; i++ )
1551 {
1552 char szFieldName[30];
1553
1554 sprintf( szFieldName, "params[%d]", i );
1555 psDatum->params[i] = poMIEntry->GetDoubleField(szFieldName);
1556 }
1557
1558 psDatum->gridname = CPLStrdup(poMIEntry->GetStringField("gridname"));
1559
1560 hHFA->pDatum = (void *) psDatum;
1561
1562 return psDatum;
1563 }
1564
1565 /************************************************************************/
1566 /* HFASetDatum() */
1567 /************************************************************************/
1568
HFASetDatum(HFAHandle hHFA,const Eprj_Datum * poDatum)1569 CPLErr HFASetDatum( HFAHandle hHFA, const Eprj_Datum *poDatum )
1570
1571 {
1572 /* -------------------------------------------------------------------- */
1573 /* Loop over bands, setting information on each one. */
1574 /* -------------------------------------------------------------------- */
1575 for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1576 {
1577 HFAEntry *poDatumEntry=NULL, *poProParms;
1578
1579 /* -------------------------------------------------------------------- */
1580 /* Create a new Projection if there isn't one present already. */
1581 /* -------------------------------------------------------------------- */
1582 poProParms =
1583 hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1584 if( poProParms == NULL )
1585 {
1586 CPLError( CE_Failure, CPLE_AppDefined,
1587 "Can't add Eprj_Datum with no Eprj_ProjParameters." );
1588 return CE_Failure;
1589 }
1590
1591 poDatumEntry = poProParms->GetNamedChild("Datum");
1592 if( poDatumEntry == NULL )
1593 {
1594 poDatumEntry = new HFAEntry( hHFA, "Datum","Eprj_Datum",
1595 poProParms );
1596 }
1597
1598 poDatumEntry->MarkDirty();
1599
1600 /* -------------------------------------------------------------------- */
1601 /* Ensure we have enough space for all the data. */
1602 /* -------------------------------------------------------------------- */
1603 int nSize;
1604 GByte *pabyData;
1605
1606 nSize = 26 + strlen(poDatum->datumname) + 1 + 7*8;
1607
1608 if( poDatum->gridname != NULL )
1609 nSize += strlen(poDatum->gridname) + 1;
1610
1611 pabyData = poDatumEntry->MakeData( nSize );
1612 if(!pabyData)
1613 return CE_Failure;
1614
1615 poDatumEntry->SetPosition();
1616
1617 // Initialize the whole thing to zeros for a clean start.
1618 memset( poDatumEntry->GetData(), 0, poDatumEntry->GetDataSize() );
1619
1620 /* -------------------------------------------------------------------- */
1621 /* Write the various fields. */
1622 /* -------------------------------------------------------------------- */
1623 poDatumEntry->SetStringField( "datumname", poDatum->datumname );
1624 poDatumEntry->SetIntField( "type", poDatum->type );
1625
1626 poDatumEntry->SetDoubleField( "params[0]", poDatum->params[0] );
1627 poDatumEntry->SetDoubleField( "params[1]", poDatum->params[1] );
1628 poDatumEntry->SetDoubleField( "params[2]", poDatum->params[2] );
1629 poDatumEntry->SetDoubleField( "params[3]", poDatum->params[3] );
1630 poDatumEntry->SetDoubleField( "params[4]", poDatum->params[4] );
1631 poDatumEntry->SetDoubleField( "params[5]", poDatum->params[5] );
1632 poDatumEntry->SetDoubleField( "params[6]", poDatum->params[6] );
1633
1634 poDatumEntry->SetStringField( "gridname", poDatum->gridname );
1635 }
1636
1637 return CE_None;
1638 }
1639
1640 /************************************************************************/
1641 /* HFAGetPCT() */
1642 /* */
1643 /* Read the PCT from a band, if it has one. */
1644 /************************************************************************/
1645
HFAGetPCT(HFAHandle hHFA,int nBand,int * pnColors,double ** ppadfRed,double ** ppadfGreen,double ** ppadfBlue,double ** ppadfAlpha,double ** ppadfBins)1646 CPLErr HFAGetPCT( HFAHandle hHFA, int nBand, int *pnColors,
1647 double **ppadfRed, double **ppadfGreen,
1648 double **ppadfBlue , double **ppadfAlpha,
1649 double **ppadfBins )
1650
1651 {
1652 if( nBand < 1 || nBand > hHFA->nBands )
1653 return CE_Failure;
1654
1655 return( hHFA->papoBand[nBand-1]->GetPCT( pnColors, ppadfRed,
1656 ppadfGreen, ppadfBlue,
1657 ppadfAlpha, ppadfBins ) );
1658 }
1659
1660 /************************************************************************/
1661 /* HFASetPCT() */
1662 /* */
1663 /* Set the PCT on a band. */
1664 /************************************************************************/
1665
HFASetPCT(HFAHandle hHFA,int nBand,int nColors,double * padfRed,double * padfGreen,double * padfBlue,double * padfAlpha)1666 CPLErr HFASetPCT( HFAHandle hHFA, int nBand, int nColors,
1667 double *padfRed, double *padfGreen, double *padfBlue,
1668 double *padfAlpha )
1669
1670 {
1671 if( nBand < 1 || nBand > hHFA->nBands )
1672 return CE_Failure;
1673
1674 return( hHFA->papoBand[nBand-1]->SetPCT( nColors, padfRed,
1675 padfGreen, padfBlue, padfAlpha ) );
1676 }
1677
1678 /************************************************************************/
1679 /* HFAGetDataRange() */
1680 /************************************************************************/
1681
HFAGetDataRange(HFAHandle hHFA,int nBand,double * pdfMin,double * pdfMax)1682 CPLErr HFAGetDataRange( HFAHandle hHFA, int nBand,
1683 double * pdfMin, double *pdfMax )
1684
1685 {
1686 HFAEntry *poBinInfo;
1687
1688 if( nBand < 1 || nBand > hHFA->nBands )
1689 return CE_Failure;
1690
1691 poBinInfo = hHFA->papoBand[nBand-1]->poNode->GetNamedChild("Statistics" );
1692
1693 if( poBinInfo == NULL )
1694 return( CE_Failure );
1695
1696 *pdfMin = poBinInfo->GetDoubleField( "minimum" );
1697 *pdfMax = poBinInfo->GetDoubleField( "maximum" );
1698
1699 if( *pdfMax > *pdfMin )
1700 return CE_None;
1701 else
1702 return CE_Failure;
1703 }
1704
1705 /************************************************************************/
1706 /* HFADumpNode() */
1707 /************************************************************************/
1708
HFADumpNode(HFAEntry * poEntry,int nIndent,int bVerbose,FILE * fp)1709 static void HFADumpNode( HFAEntry *poEntry, int nIndent, int bVerbose,
1710 FILE * fp )
1711
1712 {
1713 static char szSpaces[256];
1714 int i;
1715
1716 for( i = 0; i < nIndent*2; i++ )
1717 szSpaces[i] = ' ';
1718 szSpaces[nIndent*2] = '\0';
1719
1720 fprintf( fp, "%s%s(%s) @ %d + %d @ %d\n", szSpaces,
1721 poEntry->GetName(), poEntry->GetType(),
1722 poEntry->GetFilePos(),
1723 poEntry->GetDataSize(), poEntry->GetDataPos() );
1724
1725 if( bVerbose )
1726 {
1727 strcat( szSpaces, "+ " );
1728 poEntry->DumpFieldValues( fp, szSpaces );
1729 fprintf( fp, "\n" );
1730 }
1731
1732 if( poEntry->GetChild() != NULL )
1733 HFADumpNode( poEntry->GetChild(), nIndent+1, bVerbose, fp );
1734
1735 if( poEntry->GetNext() != NULL )
1736 HFADumpNode( poEntry->GetNext(), nIndent, bVerbose, fp );
1737 }
1738
1739 /************************************************************************/
1740 /* HFADumpTree() */
1741 /* */
1742 /* Dump the tree of information in a HFA file. */
1743 /************************************************************************/
1744
HFADumpTree(HFAHandle hHFA,FILE * fpOut)1745 void HFADumpTree( HFAHandle hHFA, FILE * fpOut )
1746
1747 {
1748 HFADumpNode( hHFA->poRoot, 0, TRUE, fpOut );
1749 }
1750
1751 /************************************************************************/
1752 /* HFADumpDictionary() */
1753 /* */
1754 /* Dump the dictionary (in raw, and parsed form) to the named */
1755 /* device. */
1756 /************************************************************************/
1757
HFADumpDictionary(HFAHandle hHFA,FILE * fpOut)1758 void HFADumpDictionary( HFAHandle hHFA, FILE * fpOut )
1759
1760 {
1761 fprintf( fpOut, "%s\n", hHFA->pszDictionary );
1762
1763 hHFA->poDictionary->Dump( fpOut );
1764 }
1765
1766 /************************************************************************/
1767 /* HFAStandard() */
1768 /* */
1769 /* Swap byte order on MSB systems. */
1770 /************************************************************************/
1771
1772 #ifdef CPL_MSB
HFAStandard(int nBytes,void * pData)1773 void HFAStandard( int nBytes, void * pData )
1774
1775 {
1776 int i;
1777 GByte *pabyData = (GByte *) pData;
1778
1779 for( i = nBytes/2-1; i >= 0; i-- )
1780 {
1781 GByte byTemp;
1782
1783 byTemp = pabyData[i];
1784 pabyData[i] = pabyData[nBytes-i-1];
1785 pabyData[nBytes-i-1] = byTemp;
1786 }
1787 }
1788 #endif
1789
1790 /* ==================================================================== */
1791 /* Default data dictionary. Emitted verbatim into the imagine */
1792 /* file. */
1793 /* ==================================================================== */
1794
1795 static const char *aszDefaultDD[] = {
1796 "{1:lversion,1:LfreeList,1:LrootEntryPtr,1:sentryHeaderLength,1:LdictionaryPtr,}Ehfa_File,{1:Lnext,1:Lprev,1:Lparent,1:Lchild,1:Ldata,1:ldataSize,64:cname,32:ctype,1:tmodTime,}Ehfa_Entry,{16:clabel,1:LheaderPtr,}Ehfa_HeaderTag,{1:LfreeList,1:lfreeSize,}Ehfa_FreeListNode,{1:lsize,1:Lptr,}Ehfa_Data,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,",
1797 "1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer_SubSample,{1:e2:raster,vector,type,1:LdictionaryPtr,}Ehfa_Layer,{1:LspaceUsedForRasterData,}ImgFormatInfo831,{1:sfileCode,1:Loffset,1:lsize,1:e2:false,true,logvalid,",
1798 "1:e2:no compression,ESRI GRID compression,compressionType,}Edms_VirtualBlockInfo,{1:lmin,1:lmax,}Edms_FreeIDList,{1:lnumvirtualblocks,1:lnumobjectsperblock,1:lnextobjectnum,1:e2:no compression,RLC compression,compressionType,0:poEdms_VirtualBlockInfo,blockinfo,0:poEdms_FreeIDList,freelist,1:tmodTime,}Edms_State,{0:pcstring,}Emif_String,{1:oEmif_String,fileName,2:LlayerStackValidFlagsOffset,2:LlayerStackDataOffset,1:LlayerStackCount,1:LlayerStackIndex,}ImgExternalRaster,{1:oEmif_String,algorithm,0:poEmif_String,nameList,}Eimg_RRDNamesList,{1:oEmif_String,projection,1:oEmif_String,units,}Eimg_MapInformation,",
1799 "{1:oEmif_String,dependent,}Eimg_DependentFile,{1:oEmif_String,ImageLayerName,}Eimg_DependentLayerName,{1:lnumrows,1:lnumcolumns,1:e13:EGDA_TYPE_U1,EGDA_TYPE_U2,EGDA_TYPE_U4,EGDA_TYPE_U8,EGDA_TYPE_S8,EGDA_TYPE_U16,EGDA_TYPE_S16,EGDA_TYPE_U32,EGDA_TYPE_S32,EGDA_TYPE_F32,EGDA_TYPE_F64,EGDA_TYPE_C64,EGDA_TYPE_C128,datatype,1:e4:EGDA_SCALAR_OBJECT,EGDA_TABLE_OBJECT,EGDA_MATRIX_OBJECT,EGDA_RASTER_OBJECT,objecttype,}Egda_BaseData,{1:*bvalueBD,}Eimg_NonInitializedValue,{1:dx,1:dy,}Eprj_Coordinate,{1:dwidth,1:dheight,}Eprj_Size,{0:pcproName,1:*oEprj_Coordinate,upperLeftCenter,",
1800 "1:*oEprj_Coordinate,lowerRightCenter,1:*oEprj_Size,pixelSize,0:pcunits,}Eprj_MapInfo,{0:pcdatumname,1:e3:EPRJ_DATUM_PARAMETRIC,EPRJ_DATUM_GRID,EPRJ_DATUM_REGRESSION,type,0:pdparams,0:pcgridname,}Eprj_Datum,{0:pcsphereName,1:da,1:db,1:deSquared,1:dradius,}Eprj_Spheroid,{1:e2:EPRJ_INTERNAL,EPRJ_EXTERNAL,proType,1:lproNumber,0:pcproExeName,0:pcproName,1:lproZone,0:pdproParams,1:*oEprj_Spheroid,proSpheroid,}Eprj_ProParameters,{1:dminimum,1:dmaximum,1:dmean,1:dmedian,1:dmode,1:dstddev,}Esta_Statistics,{1:lnumBins,1:e4:direct,linear,logarithmic,explicit,binFunctionType,1:dminLimit,1:dmaxLimit,1:*bbinLimits,}Edsc_BinFunction,{0:poEmif_String,LayerNames,1:*bExcludedValues,1:oEmif_String,AOIname,",
1801 "1:lSkipFactorX,1:lSkipFactorY,1:*oEdsc_BinFunction,BinFunction,}Eimg_StatisticsParameters830,{1:lnumrows,}Edsc_Table,{1:lnumRows,1:LcolumnDataPtr,1:e4:integer,real,complex,string,dataType,1:lmaxNumChars,}Edsc_Column,{1:lposition,0:pcname,1:e2:EMSC_FALSE,EMSC_TRUE,editable,1:e3:LEFT,CENTER,RIGHT,alignment,0:pcformat,1:e3:DEFAULT,APPLY,AUTO-APPLY,formulamode,0:pcformula,1:dcolumnwidth,0:pcunits,1:e5:NO_COLOR,RED,GREEN,BLUE,COLOR,colorflag,0:pcgreenname,0:pcbluename,}Eded_ColumnAttributes_1,{1:lversion,1:lnumobjects,1:e2:EAOI_UNION,EAOI_INTERSECTION,operation,}Eaoi_AreaOfInterest,",
1802 "{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,",
1803 "{1:x{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,projection,1:x{0:pcstring,}Emif_String,title,}Eprj_MapProjection842,",
1804 "{0:poEmif_String,titleList,}Exfr_GenericXFormHeader,{1:lorder,1:lnumdimtransform,1:lnumdimpolynomial,1:ltermcount,0:plexponentlist,1:*bpolycoefmtx,1:*bpolycoefvector,}Efga_Polynomial,",
1805 ".",
1806 NULL
1807 };
1808
1809
1810
1811 /************************************************************************/
1812 /* HFACreateLL() */
1813 /* */
1814 /* Low level creation of an Imagine file. Writes out the */
1815 /* Ehfa_HeaderTag, dictionary and Ehfa_File. */
1816 /************************************************************************/
1817
HFACreateLL(const char * pszFilename)1818 HFAHandle HFACreateLL( const char * pszFilename )
1819
1820 {
1821 VSILFILE *fp;
1822 HFAInfo_t *psInfo;
1823
1824 /* -------------------------------------------------------------------- */
1825 /* Create the file in the file system. */
1826 /* -------------------------------------------------------------------- */
1827 fp = VSIFOpenL( pszFilename, "w+b" );
1828 if( fp == NULL )
1829 {
1830 CPLError( CE_Failure, CPLE_OpenFailed,
1831 "Creation of file %s failed.",
1832 pszFilename );
1833 return NULL;
1834 }
1835
1836 /* -------------------------------------------------------------------- */
1837 /* Create the HFAInfo_t */
1838 /* -------------------------------------------------------------------- */
1839 psInfo = (HFAInfo_t *) CPLCalloc(sizeof(HFAInfo_t),1);
1840
1841 psInfo->fp = fp;
1842 psInfo->eAccess = HFA_Update;
1843 psInfo->nXSize = 0;
1844 psInfo->nYSize = 0;
1845 psInfo->nBands = 0;
1846 psInfo->papoBand = NULL;
1847 psInfo->pMapInfo = NULL;
1848 psInfo->pDatum = NULL;
1849 psInfo->pProParameters = NULL;
1850 psInfo->bTreeDirty = FALSE;
1851 psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
1852 psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
1853
1854 /* -------------------------------------------------------------------- */
1855 /* Write out the Ehfa_HeaderTag */
1856 /* -------------------------------------------------------------------- */
1857 GInt32 nHeaderPos;
1858
1859 VSIFWriteL( (void *) "EHFA_HEADER_TAG", 1, 16, fp );
1860
1861 nHeaderPos = 20;
1862 HFAStandard( 4, &nHeaderPos );
1863 VSIFWriteL( &nHeaderPos, 4, 1, fp );
1864
1865 /* -------------------------------------------------------------------- */
1866 /* Write the Ehfa_File node, locked in at offset 20. */
1867 /* -------------------------------------------------------------------- */
1868 GInt32 nVersion = 1, nFreeList = 0, nRootEntry = 0;
1869 GInt16 nEntryHeaderLength = 128;
1870 GInt32 nDictionaryPtr = 38;
1871
1872 psInfo->nEntryHeaderLength = nEntryHeaderLength;
1873 psInfo->nRootPos = 0;
1874 psInfo->nDictionaryPos = nDictionaryPtr;
1875 psInfo->nVersion = nVersion;
1876
1877 HFAStandard( 4, &nVersion );
1878 HFAStandard( 4, &nFreeList );
1879 HFAStandard( 4, &nRootEntry );
1880 HFAStandard( 2, &nEntryHeaderLength );
1881 HFAStandard( 4, &nDictionaryPtr );
1882
1883 VSIFWriteL( &nVersion, 4, 1, fp );
1884 VSIFWriteL( &nFreeList, 4, 1, fp );
1885 VSIFWriteL( &nRootEntry, 4, 1, fp );
1886 VSIFWriteL( &nEntryHeaderLength, 2, 1, fp );
1887 VSIFWriteL( &nDictionaryPtr, 4, 1, fp );
1888
1889 /* -------------------------------------------------------------------- */
1890 /* Write the dictionary, locked in at location 38. Note that */
1891 /* we jump through a bunch of hoops to operate on the */
1892 /* dictionary in chunks because some compiles (such as VC++) */
1893 /* don't allow particularly large static strings. */
1894 /* -------------------------------------------------------------------- */
1895 int nDictLen = 0, iChunk;
1896
1897 for( iChunk = 0; aszDefaultDD[iChunk] != NULL; iChunk++ )
1898 nDictLen += strlen(aszDefaultDD[iChunk]);
1899
1900 psInfo->pszDictionary = (char *) CPLMalloc(nDictLen+1);
1901 psInfo->pszDictionary[0] = '\0';
1902
1903 for( iChunk = 0; aszDefaultDD[iChunk] != NULL; iChunk++ )
1904 strcat( psInfo->pszDictionary, aszDefaultDD[iChunk] );
1905
1906 VSIFWriteL( (void *) psInfo->pszDictionary, 1,
1907 strlen(psInfo->pszDictionary)+1, fp );
1908
1909 psInfo->poDictionary = new HFADictionary( psInfo->pszDictionary );
1910
1911 psInfo->nEndOfFile = (GUInt32) VSIFTellL( fp );
1912
1913 /* -------------------------------------------------------------------- */
1914 /* Create a root entry. */
1915 /* -------------------------------------------------------------------- */
1916 psInfo->poRoot = new HFAEntry( psInfo, "root", "root", NULL );
1917
1918 /* -------------------------------------------------------------------- */
1919 /* If an .ige or .rrd file exists with the same base name, */
1920 /* delete them. (#1784) */
1921 /* -------------------------------------------------------------------- */
1922 CPLString osExtension = CPLGetExtension(pszFilename);
1923 if( !EQUAL(osExtension,"rrd") && !EQUAL(osExtension,"aux") )
1924 {
1925 CPLString osPath = CPLGetPath( pszFilename );
1926 CPLString osBasename = CPLGetBasename( pszFilename );
1927 VSIStatBufL sStatBuf;
1928 CPLString osSupFile = CPLFormCIFilename( osPath, osBasename, "rrd" );
1929
1930 if( VSIStatL( osSupFile, &sStatBuf ) == 0 )
1931 VSIUnlink( osSupFile );
1932
1933 osSupFile = CPLFormCIFilename( osPath, osBasename, "ige" );
1934
1935 if( VSIStatL( osSupFile, &sStatBuf ) == 0 )
1936 VSIUnlink( osSupFile );
1937 }
1938
1939 return psInfo;
1940 }
1941
1942 /************************************************************************/
1943 /* HFAAllocateSpace() */
1944 /* */
1945 /* Return an area in the file to the caller to write the */
1946 /* requested number of bytes. Currently this is always at the */
1947 /* end of the file, but eventually we might actually keep track */
1948 /* of free space. The HFAInfo_t's concept of file size is */
1949 /* updated, even if nothing ever gets written to this region. */
1950 /* */
1951 /* Returns the offset to the requested space, or zero one */
1952 /* failure. */
1953 /************************************************************************/
1954
HFAAllocateSpace(HFAInfo_t * psInfo,GUInt32 nBytes)1955 GUInt32 HFAAllocateSpace( HFAInfo_t *psInfo, GUInt32 nBytes )
1956
1957 {
1958 /* should check if this will wrap over 2GB limit */
1959
1960 psInfo->nEndOfFile += nBytes;
1961 return psInfo->nEndOfFile - nBytes;
1962 }
1963
1964 /************************************************************************/
1965 /* HFAFlush() */
1966 /* */
1967 /* Write out any dirty tree information to disk, putting the */
1968 /* disk file in a consistent state. */
1969 /************************************************************************/
1970
HFAFlush(HFAHandle hHFA)1971 CPLErr HFAFlush( HFAHandle hHFA )
1972
1973 {
1974 CPLErr eErr;
1975
1976 if( !hHFA->bTreeDirty && !hHFA->poDictionary->bDictionaryTextDirty )
1977 return CE_None;
1978
1979 CPLAssert( hHFA->poRoot != NULL );
1980
1981 /* -------------------------------------------------------------------- */
1982 /* Flush HFAEntry tree to disk. */
1983 /* -------------------------------------------------------------------- */
1984 if( hHFA->bTreeDirty )
1985 {
1986 eErr = hHFA->poRoot->FlushToDisk();
1987 if( eErr != CE_None )
1988 return eErr;
1989
1990 hHFA->bTreeDirty = FALSE;
1991 }
1992
1993 /* -------------------------------------------------------------------- */
1994 /* Flush Dictionary to disk. */
1995 /* -------------------------------------------------------------------- */
1996 GUInt32 nNewDictionaryPos = hHFA->nDictionaryPos;
1997
1998 if( hHFA->poDictionary->bDictionaryTextDirty )
1999 {
2000 VSIFSeekL( hHFA->fp, 0, SEEK_END );
2001 nNewDictionaryPos = (GUInt32) VSIFTellL( hHFA->fp );
2002 VSIFWriteL( hHFA->poDictionary->osDictionaryText.c_str(),
2003 strlen(hHFA->poDictionary->osDictionaryText.c_str()) + 1,
2004 1, hHFA->fp );
2005 hHFA->poDictionary->bDictionaryTextDirty = FALSE;
2006 }
2007
2008 /* -------------------------------------------------------------------- */
2009 /* do we need to update the Ehfa_File pointer to the root node? */
2010 /* -------------------------------------------------------------------- */
2011 if( hHFA->nRootPos != hHFA->poRoot->GetFilePos()
2012 || nNewDictionaryPos != hHFA->nDictionaryPos )
2013 {
2014 GUInt32 nOffset;
2015 GUInt32 nHeaderPos;
2016
2017 VSIFSeekL( hHFA->fp, 16, SEEK_SET );
2018 VSIFReadL( &nHeaderPos, sizeof(GInt32), 1, hHFA->fp );
2019 HFAStandard( 4, &nHeaderPos );
2020
2021 nOffset = hHFA->nRootPos = hHFA->poRoot->GetFilePos();
2022 HFAStandard( 4, &nOffset );
2023 VSIFSeekL( hHFA->fp, nHeaderPos+8, SEEK_SET );
2024 VSIFWriteL( &nOffset, 4, 1, hHFA->fp );
2025
2026 nOffset = hHFA->nDictionaryPos = nNewDictionaryPos;
2027 HFAStandard( 4, &nOffset );
2028 VSIFSeekL( hHFA->fp, nHeaderPos+14, SEEK_SET );
2029 VSIFWriteL( &nOffset, 4, 1, hHFA->fp );
2030 }
2031
2032 return CE_None;
2033 }
2034
2035 /************************************************************************/
2036 /* HFACreateLayer() */
2037 /* */
2038 /* Create a layer object, and corresponding RasterDMS. */
2039 /* Suitable for use with primary layers, and overviews. */
2040 /************************************************************************/
2041
2042 int
HFACreateLayer(HFAHandle psInfo,HFAEntry * poParent,const char * pszLayerName,int bOverview,int nBlockSize,int bCreateCompressed,int bCreateLargeRaster,int bDependentLayer,int nXSize,int nYSize,int nDataType,CPL_UNUSED char ** papszOptions,GIntBig nStackValidFlagsOffset,GIntBig nStackDataOffset,int nStackCount,int nStackIndex)2043 HFACreateLayer( HFAHandle psInfo, HFAEntry *poParent,
2044 const char *pszLayerName,
2045 int bOverview, int nBlockSize,
2046 int bCreateCompressed, int bCreateLargeRaster,
2047 int bDependentLayer,
2048 int nXSize, int nYSize, int nDataType,
2049 CPL_UNUSED char **papszOptions,
2050
2051 // these are only related to external (large) files
2052 GIntBig nStackValidFlagsOffset,
2053 GIntBig nStackDataOffset,
2054 int nStackCount, int nStackIndex )
2055
2056 {
2057
2058 HFAEntry *poEimg_Layer;
2059 const char *pszLayerType;
2060
2061 if( bOverview )
2062 pszLayerType = "Eimg_Layer_SubSample";
2063 else
2064 pszLayerType = "Eimg_Layer";
2065
2066 if (nBlockSize <= 0)
2067 {
2068 CPLError(CE_Failure, CPLE_IllegalArg, "HFACreateLayer : nBlockXSize < 0");
2069 return FALSE;
2070 }
2071
2072 /* -------------------------------------------------------------------- */
2073 /* Work out some details about the tiling scheme. */
2074 /* -------------------------------------------------------------------- */
2075 int nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock;
2076
2077 nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2078 nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2079 nBlocks = nBlocksPerRow * nBlocksPerColumn;
2080 nBytesPerBlock = (nBlockSize * nBlockSize
2081 * HFAGetDataTypeBits(nDataType) + 7) / 8;
2082
2083 /* -------------------------------------------------------------------- */
2084 /* Create the Eimg_Layer for the band. */
2085 /* -------------------------------------------------------------------- */
2086 poEimg_Layer =
2087 new HFAEntry( psInfo, pszLayerName, pszLayerType, poParent );
2088
2089 poEimg_Layer->SetIntField( "width", nXSize );
2090 poEimg_Layer->SetIntField( "height", nYSize );
2091 poEimg_Layer->SetStringField( "layerType", "athematic" );
2092 poEimg_Layer->SetIntField( "pixelType", nDataType );
2093 poEimg_Layer->SetIntField( "blockWidth", nBlockSize );
2094 poEimg_Layer->SetIntField( "blockHeight", nBlockSize );
2095
2096 /* -------------------------------------------------------------------- */
2097 /* Create the RasterDMS (block list). This is a complex type */
2098 /* with pointers, and variable size. We set the superstructure */
2099 /* ourselves rather than trying to have the HFA type management */
2100 /* system do it for us (since this would be hard to implement). */
2101 /* -------------------------------------------------------------------- */
2102 if ( !bCreateLargeRaster && !bDependentLayer )
2103 {
2104 int nDmsSize;
2105 HFAEntry *poEdms_State;
2106 GByte *pabyData;
2107
2108 poEdms_State =
2109 new HFAEntry( psInfo, "RasterDMS", "Edms_State", poEimg_Layer );
2110
2111 nDmsSize = 14 * nBlocks + 38;
2112 pabyData = poEdms_State->MakeData( nDmsSize );
2113
2114 /* set some simple values */
2115 poEdms_State->SetIntField( "numvirtualblocks", nBlocks );
2116 poEdms_State->SetIntField( "numobjectsperblock",
2117 nBlockSize*nBlockSize );
2118 poEdms_State->SetIntField( "nextobjectnum",
2119 nBlockSize*nBlockSize*nBlocks );
2120
2121 /* Is file compressed or not? */
2122 if( bCreateCompressed )
2123 {
2124 poEdms_State->SetStringField( "compressionType", "RLC compression" );
2125 }
2126 else
2127 {
2128 poEdms_State->SetStringField( "compressionType", "no compression" );
2129 }
2130
2131 /* we need to hardcode file offset into the data, so locate it now */
2132 poEdms_State->SetPosition();
2133
2134 /* Set block info headers */
2135 GUInt32 nValue;
2136
2137 /* blockinfo count */
2138 nValue = nBlocks;
2139 HFAStandard( 4, &nValue );
2140 memcpy( pabyData + 14, &nValue, 4 );
2141
2142 /* blockinfo position */
2143 nValue = poEdms_State->GetDataPos() + 22;
2144 HFAStandard( 4, &nValue );
2145 memcpy( pabyData + 18, &nValue, 4 );
2146
2147 /* Set each blockinfo */
2148 for( int iBlock = 0; iBlock < nBlocks; iBlock++ )
2149 {
2150 GInt16 nValue16;
2151 int nOffset = 22 + 14 * iBlock;
2152
2153 /* fileCode */
2154 nValue16 = 0;
2155 HFAStandard( 2, &nValue16 );
2156 memcpy( pabyData + nOffset, &nValue16, 2 );
2157
2158 /* offset */
2159 if( bCreateCompressed )
2160 {
2161 /* flag it with zero offset - will allocate space when we compress it */
2162 nValue = 0;
2163 }
2164 else
2165 {
2166 nValue = HFAAllocateSpace( psInfo, nBytesPerBlock );
2167 }
2168 HFAStandard( 4, &nValue );
2169 memcpy( pabyData + nOffset + 2, &nValue, 4 );
2170
2171 /* size */
2172 if( bCreateCompressed )
2173 {
2174 /* flag it with zero size - don't know until we compress it */
2175 nValue = 0;
2176 }
2177 else
2178 {
2179 nValue = nBytesPerBlock;
2180 }
2181 HFAStandard( 4, &nValue );
2182 memcpy( pabyData + nOffset + 6, &nValue, 4 );
2183
2184 /* logValid (false) */
2185 nValue16 = 0;
2186 HFAStandard( 2, &nValue16 );
2187 memcpy( pabyData + nOffset + 10, &nValue16, 2 );
2188
2189 /* compressionType */
2190 if( bCreateCompressed )
2191 nValue16 = 1;
2192 else
2193 nValue16 = 0;
2194
2195 HFAStandard( 2, &nValue16 );
2196 memcpy( pabyData + nOffset + 12, &nValue16, 2 );
2197 }
2198
2199 }
2200 /* -------------------------------------------------------------------- */
2201 /* Create ExternalRasterDMS object. */
2202 /* -------------------------------------------------------------------- */
2203 else if( bCreateLargeRaster )
2204 {
2205 HFAEntry *poEdms_State;
2206
2207 poEdms_State =
2208 new HFAEntry( psInfo, "ExternalRasterDMS",
2209 "ImgExternalRaster", poEimg_Layer );
2210 poEdms_State->MakeData( 8 + strlen(psInfo->pszIGEFilename) + 1 + 6 * 4 );
2211
2212 poEdms_State->SetStringField( "fileName.string",
2213 psInfo->pszIGEFilename );
2214
2215 poEdms_State->SetIntField( "layerStackValidFlagsOffset[0]",
2216 (int) (nStackValidFlagsOffset & 0xFFFFFFFF));
2217 poEdms_State->SetIntField( "layerStackValidFlagsOffset[1]",
2218 (int) (nStackValidFlagsOffset >> 32) );
2219
2220 poEdms_State->SetIntField( "layerStackDataOffset[0]",
2221 (int) (nStackDataOffset & 0xFFFFFFFF) );
2222 poEdms_State->SetIntField( "layerStackDataOffset[1]",
2223 (int) (nStackDataOffset >> 32 ) );
2224 poEdms_State->SetIntField( "layerStackCount", nStackCount );
2225 poEdms_State->SetIntField( "layerStackIndex", nStackIndex );
2226 }
2227
2228 /* -------------------------------------------------------------------- */
2229 /* Dependent... */
2230 /* -------------------------------------------------------------------- */
2231 else if( bDependentLayer )
2232 {
2233 HFAEntry *poDepLayerName;
2234
2235 poDepLayerName =
2236 new HFAEntry( psInfo, "DependentLayerName",
2237 "Eimg_DependentLayerName", poEimg_Layer );
2238 poDepLayerName->MakeData( 8 + strlen(pszLayerName) + 2 );
2239
2240 poDepLayerName->SetStringField( "ImageLayerName.string",
2241 pszLayerName );
2242 }
2243
2244 /* -------------------------------------------------------------------- */
2245 /* Create the Ehfa_Layer. */
2246 /* -------------------------------------------------------------------- */
2247 HFAEntry *poEhfa_Layer;
2248 GUInt32 nLDict;
2249 char szLDict[128], chBandType;
2250
2251 if( nDataType == EPT_u1 )
2252 chBandType = '1';
2253 else if( nDataType == EPT_u2 )
2254 chBandType = '2';
2255 else if( nDataType == EPT_u4 )
2256 chBandType = '4';
2257 else if( nDataType == EPT_u8 )
2258 chBandType = 'c';
2259 else if( nDataType == EPT_s8 )
2260 chBandType = 'C';
2261 else if( nDataType == EPT_u16 )
2262 chBandType = 's';
2263 else if( nDataType == EPT_s16 )
2264 chBandType = 'S';
2265 else if( nDataType == EPT_u32 )
2266 // for some reason erdas imagine expects an L for unsinged 32 bit ints
2267 // otherwise it gives strange "out of memory errors"
2268 chBandType = 'L';
2269 else if( nDataType == EPT_s32 )
2270 chBandType = 'L';
2271 else if( nDataType == EPT_f32 )
2272 chBandType = 'f';
2273 else if( nDataType == EPT_f64 )
2274 chBandType = 'd';
2275 else if( nDataType == EPT_c64 )
2276 chBandType = 'm';
2277 else if( nDataType == EPT_c128 )
2278 chBandType = 'M';
2279 else
2280 {
2281 CPLAssert( FALSE );
2282 chBandType = 'c';
2283 }
2284
2285 // the first value in the entry below gives the number of pixels within a block
2286 sprintf( szLDict, "{%d:%cdata,}RasterDMS,.", nBlockSize*nBlockSize, chBandType );
2287
2288 poEhfa_Layer = new HFAEntry( psInfo, "Ehfa_Layer", "Ehfa_Layer",
2289 poEimg_Layer );
2290 poEhfa_Layer->MakeData();
2291 poEhfa_Layer->SetPosition();
2292 nLDict = HFAAllocateSpace( psInfo, strlen(szLDict) + 1 );
2293
2294 poEhfa_Layer->SetStringField( "type", "raster" );
2295 poEhfa_Layer->SetIntField( "dictionaryPtr", nLDict );
2296
2297 VSIFSeekL( psInfo->fp, nLDict, SEEK_SET );
2298 VSIFWriteL( (void *) szLDict, strlen(szLDict) + 1, 1, psInfo->fp );
2299
2300 return TRUE;
2301 }
2302
2303
2304 /************************************************************************/
2305 /* HFACreate() */
2306 /************************************************************************/
2307
HFACreate(const char * pszFilename,int nXSize,int nYSize,int nBands,int nDataType,char ** papszOptions)2308 HFAHandle HFACreate( const char * pszFilename,
2309 int nXSize, int nYSize, int nBands,
2310 int nDataType, char ** papszOptions )
2311
2312 {
2313 HFAHandle psInfo;
2314 int nBlockSize = 64;
2315 const char * pszValue = CSLFetchNameValue( papszOptions, "BLOCKSIZE" );
2316
2317 if ( pszValue != NULL )
2318 {
2319 nBlockSize = atoi( pszValue );
2320 // check for sane values
2321 if ( (( nBlockSize < 32 ) || (nBlockSize > 2048))
2322 && !CSLTestBoolean(CPLGetConfigOption("FORCE_BLOCKSIZE", "NO")) )
2323 {
2324 nBlockSize = 64;
2325 }
2326 }
2327 int bCreateLargeRaster = CSLFetchBoolean(papszOptions,"USE_SPILL",
2328 FALSE);
2329 int bCreateCompressed =
2330 CSLFetchBoolean(papszOptions,"COMPRESS", FALSE)
2331 || CSLFetchBoolean(papszOptions,"COMPRESSED", FALSE);
2332 int bCreateAux = CSLFetchBoolean(papszOptions,"AUX", FALSE);
2333
2334 char *pszFullFilename = NULL, *pszRawFilename = NULL;
2335
2336 /* -------------------------------------------------------------------- */
2337 /* Create the low level structure. */
2338 /* -------------------------------------------------------------------- */
2339 psInfo = HFACreateLL( pszFilename );
2340 if( psInfo == NULL )
2341 return NULL;
2342
2343 /* -------------------------------------------------------------------- */
2344 /* Create the DependentFile node if requested. */
2345 /* -------------------------------------------------------------------- */
2346 const char *pszDependentFile =
2347 CSLFetchNameValue( papszOptions, "DEPENDENT_FILE" );
2348
2349 if( pszDependentFile != NULL )
2350 {
2351 HFAEntry *poDF = new HFAEntry( psInfo, "DependentFile",
2352 "Eimg_DependentFile", psInfo->poRoot );
2353
2354 poDF->MakeData( strlen(pszDependentFile) + 50 );
2355 poDF->SetPosition();
2356 poDF->SetStringField( "dependent.string", pszDependentFile );
2357 }
2358
2359 /* -------------------------------------------------------------------- */
2360 /* Work out some details about the tiling scheme. */
2361 /* -------------------------------------------------------------------- */
2362 int nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock;
2363
2364 nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2365 nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2366 nBlocks = nBlocksPerRow * nBlocksPerColumn;
2367 nBytesPerBlock = (nBlockSize * nBlockSize
2368 * HFAGetDataTypeBits(nDataType) + 7) / 8;
2369
2370 CPLDebug( "HFACreate", "Blocks per row %d, blocks per column %d, "
2371 "total number of blocks %d, bytes per block %d.",
2372 nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock );
2373
2374 /* -------------------------------------------------------------------- */
2375 /* Check whether we should create external large file with */
2376 /* image. We create a spill file if the amount of imagery is */
2377 /* close to 2GB. We don't check the amount of auxiliary */
2378 /* information, so in theory if there were an awful lot of */
2379 /* non-imagery data our approximate size could be smaller than */
2380 /* the file will actually we be. We leave room for 10MB of */
2381 /* auxiliary data. */
2382 /* We can also force spill file creation using option */
2383 /* SPILL_FILE=YES. */
2384 /* -------------------------------------------------------------------- */
2385 double dfApproxSize = (double)nBytesPerBlock * (double)nBlocks *
2386 (double)nBands + 10000000.0;
2387
2388 if( dfApproxSize > 2147483648.0 && !bCreateAux )
2389 bCreateLargeRaster = TRUE;
2390
2391 // erdas imagine creates this entry even if an external spill file is used
2392 if( !bCreateAux )
2393 {
2394 HFAEntry *poImgFormat;
2395 poImgFormat = new HFAEntry( psInfo, "IMGFormatInfo",
2396 "ImgFormatInfo831", psInfo->poRoot );
2397 poImgFormat->MakeData();
2398 if ( bCreateLargeRaster )
2399 {
2400 poImgFormat->SetIntField( "spaceUsedForRasterData", 0 );
2401 bCreateCompressed = FALSE; // Can't be compressed if we are creating a spillfile
2402 }
2403 else
2404 {
2405 poImgFormat->SetIntField( "spaceUsedForRasterData",
2406 nBytesPerBlock*nBlocks*nBands );
2407 }
2408 }
2409
2410 /* -------------------------------------------------------------------- */
2411 /* Create external file and write its header. */
2412 /* -------------------------------------------------------------------- */
2413 GIntBig nValidFlagsOffset = 0, nDataOffset = 0;
2414
2415 if( bCreateLargeRaster )
2416 {
2417 if( !HFACreateSpillStack( psInfo, nXSize, nYSize, nBands,
2418 nBlockSize, nDataType,
2419 &nValidFlagsOffset, &nDataOffset ) )
2420 {
2421 CPLFree( pszRawFilename );
2422 CPLFree( pszFullFilename );
2423 return NULL;
2424 }
2425 }
2426
2427 /* ==================================================================== */
2428 /* Create each band (layer) */
2429 /* ==================================================================== */
2430 int iBand;
2431
2432 for( iBand = 0; iBand < nBands; iBand++ )
2433 {
2434 char szName[128];
2435
2436 sprintf( szName, "Layer_%d", iBand + 1 );
2437
2438 if( !HFACreateLayer( psInfo, psInfo->poRoot, szName, FALSE, nBlockSize,
2439 bCreateCompressed, bCreateLargeRaster, bCreateAux,
2440 nXSize, nYSize, nDataType, papszOptions,
2441 nValidFlagsOffset, nDataOffset,
2442 nBands, iBand ) )
2443 {
2444 HFAClose( psInfo );
2445 return NULL;
2446 }
2447 }
2448
2449 /* -------------------------------------------------------------------- */
2450 /* Initialize the band information. */
2451 /* -------------------------------------------------------------------- */
2452 HFAParseBandInfo( psInfo );
2453
2454 return psInfo;
2455 }
2456
2457 /************************************************************************/
2458 /* HFACreateOverview() */
2459 /* */
2460 /* Create an overview layer object for a band. */
2461 /************************************************************************/
2462
HFACreateOverview(HFAHandle hHFA,int nBand,int nOverviewLevel,const char * pszResampling)2463 int HFACreateOverview( HFAHandle hHFA, int nBand, int nOverviewLevel,
2464 const char *pszResampling )
2465
2466 {
2467 if( nBand < 1 || nBand > hHFA->nBands )
2468 return -1;
2469 else
2470 {
2471 HFABand *poBand = hHFA->papoBand[nBand-1];
2472 return poBand->CreateOverview( nOverviewLevel, pszResampling );
2473 }
2474 }
2475
2476 /************************************************************************/
2477 /* HFAGetMetadata() */
2478 /* */
2479 /* Read metadata structured in a table called GDAL_MetaData. */
2480 /************************************************************************/
2481
HFAGetMetadata(HFAHandle hHFA,int nBand)2482 char ** HFAGetMetadata( HFAHandle hHFA, int nBand )
2483
2484 {
2485 HFAEntry *poTable;
2486
2487 if( nBand > 0 && nBand <= hHFA->nBands )
2488 poTable = hHFA->papoBand[nBand - 1]->poNode->GetChild();
2489 else if( nBand == 0 )
2490 poTable = hHFA->poRoot->GetChild();
2491 else
2492 return NULL;
2493
2494 for( ; poTable != NULL && !EQUAL(poTable->GetName(),"GDAL_MetaData");
2495 poTable = poTable->GetNext() ) {}
2496
2497 if( poTable == NULL || !EQUAL(poTable->GetType(),"Edsc_Table") )
2498 return NULL;
2499
2500 if( poTable->GetIntField( "numRows" ) != 1 )
2501 {
2502 CPLDebug( "HFADataset", "GDAL_MetaData.numRows = %d, expected 1!",
2503 poTable->GetIntField( "numRows" ) );
2504 return NULL;
2505 }
2506
2507 /* -------------------------------------------------------------------- */
2508 /* Loop over each column. Each column will be one metadata */
2509 /* entry, with the title being the key, and the row value being */
2510 /* the value. There is only ever one row in GDAL_MetaData */
2511 /* tables. */
2512 /* -------------------------------------------------------------------- */
2513 HFAEntry *poColumn;
2514 char **papszMD = NULL;
2515
2516 for( poColumn = poTable->GetChild();
2517 poColumn != NULL;
2518 poColumn = poColumn->GetNext() )
2519 {
2520 const char *pszValue;
2521 int columnDataPtr;
2522
2523 // Skip the #Bin_Function# entry.
2524 if( EQUALN(poColumn->GetName(),"#",1) )
2525 continue;
2526
2527 pszValue = poColumn->GetStringField( "dataType" );
2528 if( pszValue == NULL || !EQUAL(pszValue,"string") )
2529 continue;
2530
2531 columnDataPtr = poColumn->GetIntField( "columnDataPtr" );
2532 if( columnDataPtr == 0 )
2533 continue;
2534
2535 /* -------------------------------------------------------------------- */
2536 /* read up to nMaxNumChars bytes from the indicated location. */
2537 /* allocate required space temporarily */
2538 /* nMaxNumChars should have been set by GDAL orginally so we should*/
2539 /* trust it, but who knows... */
2540 /* -------------------------------------------------------------------- */
2541 int nMaxNumChars = poColumn->GetIntField( "maxNumChars" );
2542
2543 if( nMaxNumChars == 0 )
2544 {
2545 papszMD = CSLSetNameValue( papszMD, poColumn->GetName(), "" );
2546 }
2547 else
2548 {
2549 char *pszMDValue = (char*) VSIMalloc(nMaxNumChars);
2550 if (pszMDValue == NULL)
2551 {
2552 CPLError(CE_Failure, CPLE_OutOfMemory,
2553 "HFAGetMetadata : Out of memory while allocating %d bytes", nMaxNumChars);
2554 continue;
2555 }
2556
2557 if( VSIFSeekL( hHFA->fp, columnDataPtr, SEEK_SET ) != 0 )
2558 continue;
2559
2560 int nMDBytes = VSIFReadL( pszMDValue, 1, nMaxNumChars, hHFA->fp );
2561 if( nMDBytes == 0 )
2562 {
2563 CPLFree( pszMDValue );
2564 continue;
2565 }
2566
2567 pszMDValue[nMaxNumChars-1] = '\0';
2568
2569 papszMD = CSLSetNameValue( papszMD, poColumn->GetName(),
2570 pszMDValue );
2571 CPLFree( pszMDValue );
2572 }
2573 }
2574
2575 return papszMD;
2576 }
2577
2578 /************************************************************************/
2579 /* HFASetGDALMetadata() */
2580 /* */
2581 /* This function is used to set metadata in a table called */
2582 /* GDAL_MetaData. It is called by HFASetMetadata() for all */
2583 /* metadata items that aren't some specific supported */
2584 /* information (like histogram or stats info). */
2585 /************************************************************************/
2586
2587 static CPLErr
HFASetGDALMetadata(HFAHandle hHFA,int nBand,char ** papszMD)2588 HFASetGDALMetadata( HFAHandle hHFA, int nBand, char **papszMD )
2589
2590 {
2591 if( papszMD == NULL )
2592 return CE_None;
2593
2594 HFAEntry *poNode;
2595
2596 if( nBand > 0 && nBand <= hHFA->nBands )
2597 poNode = hHFA->papoBand[nBand - 1]->poNode;
2598 else if( nBand == 0 )
2599 poNode = hHFA->poRoot;
2600 else
2601 return CE_Failure;
2602
2603 /* -------------------------------------------------------------------- */
2604 /* Create the Descriptor table. */
2605 /* Check we have no table with this name already */
2606 /* -------------------------------------------------------------------- */
2607 HFAEntry *poEdsc_Table = poNode->GetNamedChild( "GDAL_MetaData" );
2608
2609 if( poEdsc_Table == NULL || !EQUAL(poEdsc_Table->GetType(),"Edsc_Table") )
2610 poEdsc_Table = new HFAEntry( hHFA, "GDAL_MetaData", "Edsc_Table",
2611 poNode );
2612
2613 poEdsc_Table->SetIntField( "numrows", 1 );
2614
2615 /* -------------------------------------------------------------------- */
2616 /* Create the Binning function node. I am not sure that we */
2617 /* really need this though. */
2618 /* Check it doesn't exist already */
2619 /* -------------------------------------------------------------------- */
2620 HFAEntry *poEdsc_BinFunction =
2621 poEdsc_Table->GetNamedChild( "#Bin_Function#" );
2622
2623 if( poEdsc_BinFunction == NULL
2624 || !EQUAL(poEdsc_BinFunction->GetType(),"Edsc_BinFunction") )
2625 poEdsc_BinFunction = new HFAEntry( hHFA, "#Bin_Function#",
2626 "Edsc_BinFunction", poEdsc_Table );
2627
2628 // Because of the BaseData we have to hardcode the size.
2629 poEdsc_BinFunction->MakeData( 30 );
2630
2631 poEdsc_BinFunction->SetIntField( "numBins", 1 );
2632 poEdsc_BinFunction->SetStringField( "binFunction", "direct" );
2633 poEdsc_BinFunction->SetDoubleField( "minLimit", 0.0 );
2634 poEdsc_BinFunction->SetDoubleField( "maxLimit", 0.0 );
2635
2636 /* -------------------------------------------------------------------- */
2637 /* Process each metadata item as a separate column. */
2638 /* -------------------------------------------------------------------- */
2639 for( int iColumn = 0; papszMD[iColumn] != NULL; iColumn++ )
2640 {
2641 HFAEntry *poEdsc_Column;
2642 char *pszKey = NULL;
2643 const char *pszValue;
2644
2645 pszValue = CPLParseNameValue( papszMD[iColumn], &pszKey );
2646 if( pszValue == NULL )
2647 continue;
2648
2649 /* -------------------------------------------------------------------- */
2650 /* Create the Edsc_Column. */
2651 /* Check it doesn't exist already */
2652 /* -------------------------------------------------------------------- */
2653 poEdsc_Column = poEdsc_Table->GetNamedChild(pszKey);
2654
2655 if( poEdsc_Column == NULL
2656 || !EQUAL(poEdsc_Column->GetType(),"Edsc_Column") )
2657 poEdsc_Column = new HFAEntry( hHFA, pszKey, "Edsc_Column",
2658 poEdsc_Table );
2659
2660 poEdsc_Column->SetIntField( "numRows", 1 );
2661 poEdsc_Column->SetStringField( "dataType", "string" );
2662 poEdsc_Column->SetIntField( "maxNumChars", strlen(pszValue)+1 );
2663
2664 /* -------------------------------------------------------------------- */
2665 /* Write the data out. */
2666 /* -------------------------------------------------------------------- */
2667 int nOffset = HFAAllocateSpace( hHFA, strlen(pszValue)+1);
2668
2669 poEdsc_Column->SetIntField( "columnDataPtr", nOffset );
2670
2671 VSIFSeekL( hHFA->fp, nOffset, SEEK_SET );
2672 VSIFWriteL( (void *) pszValue, 1, strlen(pszValue)+1, hHFA->fp );
2673
2674 CPLFree( pszKey );
2675 }
2676
2677 return CE_Failure;
2678 }
2679
2680 /************************************************************************/
2681 /* HFASetMetadata() */
2682 /************************************************************************/
2683
HFASetMetadata(HFAHandle hHFA,int nBand,char ** papszMD)2684 CPLErr HFASetMetadata( HFAHandle hHFA, int nBand, char **papszMD )
2685
2686 {
2687 char **papszGDALMD = NULL;
2688
2689 if( CSLCount(papszMD) == 0 )
2690 return CE_None;
2691
2692 HFAEntry *poNode;
2693
2694 if( nBand > 0 && nBand <= hHFA->nBands )
2695 poNode = hHFA->papoBand[nBand - 1]->poNode;
2696 else if( nBand == 0 )
2697 poNode = hHFA->poRoot;
2698 else
2699 return CE_Failure;
2700
2701 /* -------------------------------------------------------------------- */
2702 /* Check if the Metadata is an "known" entity which should be */
2703 /* stored in a better place. */
2704 /* -------------------------------------------------------------------- */
2705 char * pszBinValues = NULL;
2706 int bCreatedHistogramParameters = FALSE;
2707 int bCreatedStatistics = FALSE;
2708 const char ** pszAuxMetaData = GetHFAAuxMetaDataList();
2709 // check each metadata item
2710 for( int iColumn = 0; papszMD[iColumn] != NULL; iColumn++ )
2711 {
2712 char *pszKey = NULL;
2713 const char *pszValue;
2714
2715 pszValue = CPLParseNameValue( papszMD[iColumn], &pszKey );
2716 if( pszValue == NULL )
2717 continue;
2718
2719 // know look if its known
2720 int i;
2721 for( i = 0; pszAuxMetaData[i] != NULL; i += 4 )
2722 {
2723 if ( EQUALN( pszAuxMetaData[i + 2], pszKey, strlen(pszKey) ) )
2724 break;
2725 }
2726 if ( pszAuxMetaData[i] != NULL )
2727 {
2728 // found one, get the right entry
2729 HFAEntry *poEntry;
2730
2731 if( strlen(pszAuxMetaData[i]) > 0 )
2732 poEntry = poNode->GetNamedChild( pszAuxMetaData[i] );
2733 else
2734 poEntry = poNode;
2735
2736 if( poEntry == NULL && strlen(pszAuxMetaData[i+3]) > 0 )
2737 {
2738 // child does not yet exist --> create it
2739 poEntry = new HFAEntry( hHFA, pszAuxMetaData[i], pszAuxMetaData[i+3],
2740 poNode );
2741
2742 if ( EQUALN( "Statistics", pszAuxMetaData[i], 10 ) )
2743 bCreatedStatistics = TRUE;
2744
2745 if( EQUALN( "HistogramParameters", pszAuxMetaData[i], 19 ) )
2746 {
2747 // this is a bit nasty I need to set the string field for the object
2748 // first because the SetStringField sets the count for the object
2749 // BinFunction to the length of the string
2750 poEntry->MakeData( 70 );
2751 poEntry->SetStringField( "BinFunction.binFunctionType", "direct" );
2752
2753 bCreatedHistogramParameters = TRUE;
2754 }
2755 }
2756 if ( poEntry == NULL )
2757 {
2758 CPLFree( pszKey );
2759 continue;
2760 }
2761
2762 const char *pszFieldName = pszAuxMetaData[i+1] + 1;
2763 switch( pszAuxMetaData[i+1][0] )
2764 {
2765 case 'd':
2766 {
2767 double dfValue = CPLAtof( pszValue );
2768 poEntry->SetDoubleField( pszFieldName, dfValue );
2769 }
2770 break;
2771 case 'i':
2772 case 'l':
2773 {
2774 int nValue = atoi( pszValue );
2775 poEntry->SetIntField( pszFieldName, nValue );
2776 }
2777 break;
2778 case 's':
2779 case 'e':
2780 {
2781 poEntry->SetStringField( pszFieldName, pszValue );
2782 }
2783 break;
2784 default:
2785 CPLAssert( FALSE );
2786 }
2787 }
2788 else if ( EQUALN( "STATISTICS_HISTOBINVALUES", pszKey, strlen(pszKey) ) )
2789 {
2790 pszBinValues = CPLStrdup( pszValue );
2791 }
2792 else
2793 papszGDALMD = CSLAddString( papszGDALMD, papszMD[iColumn] );
2794
2795 CPLFree( pszKey );
2796 }
2797
2798 /* -------------------------------------------------------------------- */
2799 /* Special case to write out the histogram. */
2800 /* -------------------------------------------------------------------- */
2801 if ( pszBinValues != NULL )
2802 {
2803 HFAEntry * poEntry = poNode->GetNamedChild( "HistogramParameters" );
2804 if ( poEntry != NULL && bCreatedHistogramParameters )
2805 {
2806 // if this node exists we have added Histogram data -- complete with some defaults
2807 poEntry->SetIntField( "SkipFactorX", 1 );
2808 poEntry->SetIntField( "SkipFactorY", 1 );
2809
2810 int nNumBins = poEntry->GetIntField( "BinFunction.numBins" );
2811 double dMinLimit = poEntry->GetDoubleField( "BinFunction.minLimit" );
2812 double dMaxLimit = poEntry->GetDoubleField( "BinFunction.maxLimit" );
2813
2814 // fill the descriptor table - check it isn't there already
2815 poEntry = poNode->GetNamedChild( "Descriptor_Table" );
2816 if( poEntry == NULL || !EQUAL(poEntry->GetType(),"Edsc_Table") )
2817 poEntry = new HFAEntry( hHFA, "Descriptor_Table", "Edsc_Table", poNode );
2818
2819 poEntry->SetIntField( "numRows", nNumBins );
2820
2821 // bin function
2822 HFAEntry * poBinFunc = poEntry->GetNamedChild( "#Bin_Function#" );
2823 if( poBinFunc == NULL || !EQUAL(poBinFunc->GetType(),"Edsc_BinFunction") )
2824 poBinFunc = new HFAEntry( hHFA, "#Bin_Function#", "Edsc_BinFunction", poEntry );
2825
2826 poBinFunc->MakeData( 30 );
2827 poBinFunc->SetIntField( "numBins", nNumBins );
2828 poBinFunc->SetDoubleField( "minLimit", dMinLimit );
2829 poBinFunc->SetDoubleField( "maxLimit", dMaxLimit );
2830 // direct for thematic layers, linear otherwise
2831 if ( EQUALN ( poNode->GetStringField("layerType"), "thematic", 8) )
2832 poBinFunc->SetStringField( "binFunctionType", "direct" );
2833 else
2834 poBinFunc->SetStringField( "binFunctionType", "linear" );
2835
2836 // we need a child named histogram
2837 HFAEntry * poHisto = poEntry->GetNamedChild( "Histogram" );
2838 if( poHisto == NULL || !EQUAL(poHisto->GetType(),"Edsc_Column") )
2839 poHisto = new HFAEntry( hHFA, "Histogram", "Edsc_Column", poEntry );
2840
2841 poHisto->SetIntField( "numRows", nNumBins );
2842 // allocate space for the bin values
2843 GUInt32 nOffset = HFAAllocateSpace( hHFA, nNumBins*8 );
2844 poHisto->SetIntField( "columnDataPtr", nOffset );
2845 poHisto->SetStringField( "dataType", "real" );
2846 poHisto->SetIntField( "maxNumChars", 0 );
2847 // write out histogram data
2848 char * pszWork = pszBinValues;
2849 for ( int nBin = 0; nBin < nNumBins; ++nBin )
2850 {
2851 char * pszEnd = strchr( pszWork, '|' );
2852 if ( pszEnd != NULL )
2853 {
2854 *pszEnd = 0;
2855 VSIFSeekL( hHFA->fp, nOffset + 8*nBin, SEEK_SET );
2856 double nValue = CPLAtof( pszWork );
2857 HFAStandard( 8, &nValue );
2858
2859 VSIFWriteL( (void *)&nValue, 1, 8, hHFA->fp );
2860 pszWork = pszEnd + 1;
2861 }
2862 }
2863 }
2864 else if ( poEntry != NULL )
2865 {
2866 // In this case, there are HistogramParameters present, but we did not
2867 // create them. However, we might be modifying them, in the case where
2868 // the data has changed and the histogram counts need to be updated. It could
2869 // be worse than that, but that is all we are going to cope with for now.
2870 // We are assuming that we did not change any of the other stuff, like
2871 // skip factors and so forth. The main need for this case is for programs
2872 // (such as Imagine itself) which will happily modify the pixel values
2873 // without re-calculating the histogram counts.
2874 int nNumBins = poEntry->GetIntField( "BinFunction.numBins" );
2875 HFAEntry *poEntryDescrTbl = poNode->GetNamedChild( "Descriptor_Table" );
2876 HFAEntry *poHisto = NULL;
2877 if ( poEntryDescrTbl != NULL) {
2878 poHisto = poEntryDescrTbl->GetNamedChild( "Histogram" );
2879 }
2880 if ( poHisto != NULL ) {
2881 int nOffset = poHisto->GetIntField( "columnDataPtr" );
2882 // write out histogram data
2883 char * pszWork = pszBinValues;
2884
2885 // Check whether histogram counts were written as int or double
2886 bool bCountIsInt = TRUE;
2887 const char *pszDataType = poHisto->GetStringField("dataType");
2888 if ( EQUALN(pszDataType, "real", strlen(pszDataType)) )
2889 {
2890 bCountIsInt = FALSE;
2891 }
2892 for ( int nBin = 0; nBin < nNumBins; ++nBin )
2893 {
2894 char * pszEnd = strchr( pszWork, '|' );
2895 if ( pszEnd != NULL )
2896 {
2897 *pszEnd = 0;
2898 if ( bCountIsInt ) {
2899 // Histogram counts were written as ints, so re-write them the same way
2900 VSIFSeekL( hHFA->fp, nOffset + 4*nBin, SEEK_SET );
2901 int nValue = atoi( pszWork );
2902 HFAStandard( 4, &nValue );
2903 VSIFWriteL( (void *)&nValue, 1, 4, hHFA->fp );
2904 } else {
2905 // Histogram were written as doubles, as is now the default behaviour
2906 VSIFSeekL( hHFA->fp, nOffset + 8*nBin, SEEK_SET );
2907 double nValue = CPLAtof( pszWork );
2908 HFAStandard( 8, &nValue );
2909 VSIFWriteL( (void *)&nValue, 1, 8, hHFA->fp );
2910 }
2911 pszWork = pszEnd + 1;
2912 }
2913 }
2914 }
2915 }
2916 CPLFree( pszBinValues );
2917 }
2918
2919 /* -------------------------------------------------------------------- */
2920 /* If we created a statistics node then try to create a */
2921 /* StatisticsParameters node too. */
2922 /* -------------------------------------------------------------------- */
2923 if( bCreatedStatistics )
2924 {
2925 HFAEntry *poEntry =
2926 new HFAEntry( hHFA, "StatisticsParameters",
2927 "Eimg_StatisticsParameters830", poNode );
2928
2929 poEntry->MakeData( 70 );
2930 //poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
2931
2932 poEntry->SetIntField( "SkipFactorX", 1 );
2933 poEntry->SetIntField( "SkipFactorY", 1 );
2934 }
2935
2936 /* -------------------------------------------------------------------- */
2937 /* Write out metadata items without a special place. */
2938 /* -------------------------------------------------------------------- */
2939 if( CSLCount( papszGDALMD) != 0 )
2940 {
2941 CPLErr eErr = HFASetGDALMetadata( hHFA, nBand, papszGDALMD );
2942
2943 CSLDestroy( papszGDALMD );
2944 return eErr;
2945 }
2946 else
2947 return CE_Failure;
2948 }
2949
2950 /************************************************************************/
2951 /* HFAGetIGEFilename() */
2952 /* */
2953 /* Returns the .ige filename if one is associated with this */
2954 /* object. For files not newly created we need to scan the */
2955 /* bands for spill files. Presumably there will only be one. */
2956 /* */
2957 /* NOTE: Returns full path, not just the filename portion. */
2958 /************************************************************************/
2959
HFAGetIGEFilename(HFAHandle hHFA)2960 const char *HFAGetIGEFilename( HFAHandle hHFA )
2961
2962 {
2963 if( hHFA->pszIGEFilename == NULL )
2964 {
2965 HFAEntry *poDMS = NULL;
2966 std::vector<HFAEntry*> apoDMSList =
2967 hHFA->poRoot->FindChildren( NULL, "ImgExternalRaster" );
2968
2969 if( apoDMSList.size() > 0 )
2970 poDMS = apoDMSList[0];
2971
2972 /* -------------------------------------------------------------------- */
2973 /* Get the IGE filename from if we have an ExternalRasterDMS */
2974 /* -------------------------------------------------------------------- */
2975 if( poDMS )
2976 {
2977 const char *pszRawFilename =
2978 poDMS->GetStringField( "fileName.string" );
2979
2980 if( pszRawFilename != NULL )
2981 {
2982 VSIStatBufL sStatBuf;
2983 CPLString osFullFilename =
2984 CPLFormFilename( hHFA->pszPath, pszRawFilename, NULL );
2985
2986 if( VSIStatL( osFullFilename, &sStatBuf ) != 0 )
2987 {
2988 CPLString osExtension = CPLGetExtension(pszRawFilename);
2989 CPLString osBasename = CPLGetBasename(hHFA->pszFilename);
2990 CPLString osFullFilename =
2991 CPLFormFilename( hHFA->pszPath, osBasename,
2992 osExtension );
2993
2994 if( VSIStatL( osFullFilename, &sStatBuf ) == 0 )
2995 hHFA->pszIGEFilename =
2996 CPLStrdup(
2997 CPLFormFilename( NULL, osBasename,
2998 osExtension ) );
2999 else
3000 hHFA->pszIGEFilename = CPLStrdup( pszRawFilename );
3001 }
3002 else
3003 hHFA->pszIGEFilename = CPLStrdup( pszRawFilename );
3004 }
3005 }
3006 }
3007
3008 /* -------------------------------------------------------------------- */
3009 /* Return the full filename. */
3010 /* -------------------------------------------------------------------- */
3011 if( hHFA->pszIGEFilename )
3012 return CPLFormFilename( hHFA->pszPath, hHFA->pszIGEFilename, NULL );
3013 else
3014 return NULL;
3015 }
3016
3017 /************************************************************************/
3018 /* HFACreateSpillStack() */
3019 /* */
3020 /* Create a new stack of raster layers in the spill (.ige) */
3021 /* file. Create the spill file if it didn't exist before. */
3022 /************************************************************************/
3023
HFACreateSpillStack(HFAInfo_t * psInfo,int nXSize,int nYSize,int nLayers,int nBlockSize,int nDataType,GIntBig * pnValidFlagsOffset,GIntBig * pnDataOffset)3024 int HFACreateSpillStack( HFAInfo_t *psInfo, int nXSize, int nYSize,
3025 int nLayers, int nBlockSize, int nDataType,
3026 GIntBig *pnValidFlagsOffset,
3027 GIntBig *pnDataOffset )
3028
3029 {
3030 /* -------------------------------------------------------------------- */
3031 /* Form .ige filename. */
3032 /* -------------------------------------------------------------------- */
3033 char *pszFullFilename;
3034
3035 if (nBlockSize <= 0)
3036 {
3037 CPLError(CE_Failure, CPLE_IllegalArg, "HFACreateSpillStack : nBlockXSize < 0");
3038 return FALSE;
3039 }
3040
3041 if( psInfo->pszIGEFilename == NULL )
3042 {
3043 if( EQUAL(CPLGetExtension(psInfo->pszFilename),"rrd") )
3044 psInfo->pszIGEFilename =
3045 CPLStrdup( CPLResetExtension( psInfo->pszFilename, "rde" ) );
3046 else if( EQUAL(CPLGetExtension(psInfo->pszFilename),"aux") )
3047 psInfo->pszIGEFilename =
3048 CPLStrdup( CPLResetExtension( psInfo->pszFilename, "axe" ) );
3049 else
3050 psInfo->pszIGEFilename =
3051 CPLStrdup( CPLResetExtension( psInfo->pszFilename, "ige" ) );
3052 }
3053
3054 pszFullFilename =
3055 CPLStrdup( CPLFormFilename( psInfo->pszPath, psInfo->pszIGEFilename, NULL ) );
3056
3057 /* -------------------------------------------------------------------- */
3058 /* Try and open it. If we fail, create it and write the magic */
3059 /* header. */
3060 /* -------------------------------------------------------------------- */
3061 static const char *pszMagick = "ERDAS_IMG_EXTERNAL_RASTER";
3062 VSILFILE *fpVSIL;
3063
3064 fpVSIL = VSIFOpenL( pszFullFilename, "r+b" );
3065 if( fpVSIL == NULL )
3066 {
3067 fpVSIL = VSIFOpenL( pszFullFilename, "w+" );
3068 if( fpVSIL == NULL )
3069 {
3070 CPLError( CE_Failure, CPLE_OpenFailed,
3071 "Failed to create spill file %s.\n%s",
3072 psInfo->pszIGEFilename, VSIStrerror( errno ) );
3073 return FALSE;
3074 }
3075
3076 VSIFWriteL( (void *) pszMagick, 1, strlen(pszMagick)+1, fpVSIL );
3077 }
3078
3079 CPLFree( pszFullFilename );
3080
3081 /* -------------------------------------------------------------------- */
3082 /* Work out some details about the tiling scheme. */
3083 /* -------------------------------------------------------------------- */
3084 int nBlocksPerRow, nBlocksPerColumn, /* nBlocks, */ nBytesPerBlock;
3085 int nBytesPerRow, nBlockMapSize /* , iFlagsSize */;
3086
3087 nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
3088 nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
3089 /* nBlocks = nBlocksPerRow * nBlocksPerColumn; */
3090 nBytesPerBlock = (nBlockSize * nBlockSize
3091 * HFAGetDataTypeBits(nDataType) + 7) / 8;
3092
3093 nBytesPerRow = ( nBlocksPerRow + 7 ) / 8;
3094 nBlockMapSize = nBytesPerRow * nBlocksPerColumn;
3095 /* iFlagsSize = nBlockMapSize + 20; */
3096
3097 /* -------------------------------------------------------------------- */
3098 /* Write stack prefix information. */
3099 /* -------------------------------------------------------------------- */
3100 GByte bUnknown;
3101 GInt32 nValue32;
3102
3103 VSIFSeekL( fpVSIL, 0, SEEK_END );
3104
3105 bUnknown = 1;
3106 VSIFWriteL( &bUnknown, 1, 1, fpVSIL );
3107 nValue32 = nLayers;
3108 HFAStandard( 4, &nValue32 );
3109 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3110 nValue32 = nXSize;
3111 HFAStandard( 4, &nValue32 );
3112 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3113 nValue32 = nYSize;
3114 HFAStandard( 4, &nValue32 );
3115 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3116 nValue32 = nBlockSize;
3117 HFAStandard( 4, &nValue32 );
3118 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3119 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3120 bUnknown = 3;
3121 VSIFWriteL( &bUnknown, 1, 1, fpVSIL );
3122 bUnknown = 0;
3123 VSIFWriteL( &bUnknown, 1, 1, fpVSIL );
3124
3125 /* -------------------------------------------------------------------- */
3126 /* Write out ValidFlags section(s). */
3127 /* -------------------------------------------------------------------- */
3128 unsigned char *pabyBlockMap;
3129 int iBand;
3130
3131 *pnValidFlagsOffset = VSIFTellL( fpVSIL );
3132
3133 pabyBlockMap = (unsigned char *) VSIMalloc( nBlockMapSize );
3134 if (pabyBlockMap == NULL)
3135 {
3136 CPLError(CE_Failure, CPLE_OutOfMemory, "HFACreateSpillStack : Out of memory");
3137 VSIFCloseL( fpVSIL );
3138 return FALSE;
3139 }
3140
3141 memset( pabyBlockMap, 0xff, nBlockMapSize );
3142 for ( iBand = 0; iBand < nLayers; iBand++ )
3143 {
3144 int i, iRemainder;
3145
3146 nValue32 = 1; // Unknown
3147 HFAStandard( 4, &nValue32 );
3148 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3149 nValue32 = 0; // Unknown
3150 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3151 nValue32 = nBlocksPerColumn;
3152 HFAStandard( 4, &nValue32 );
3153 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3154 nValue32 = nBlocksPerRow;
3155 HFAStandard( 4, &nValue32 );
3156 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3157 nValue32 = 0x30000; // Unknown
3158 HFAStandard( 4, &nValue32 );
3159 VSIFWriteL( &nValue32, 4, 1, fpVSIL );
3160
3161 iRemainder = nBlocksPerRow % 8;
3162 CPLDebug( "HFACreate",
3163 "Block map size %d, bytes per row %d, remainder %d.",
3164 nBlockMapSize, nBytesPerRow, iRemainder );
3165 if ( iRemainder )
3166 {
3167 for ( i = nBytesPerRow - 1; i < nBlockMapSize; i+=nBytesPerRow )
3168 pabyBlockMap[i] = (GByte) ((1<<iRemainder) - 1);
3169 }
3170
3171 VSIFWriteL( pabyBlockMap, 1, nBlockMapSize, fpVSIL );
3172 }
3173 CPLFree(pabyBlockMap);
3174 pabyBlockMap = NULL;
3175
3176 /* -------------------------------------------------------------------- */
3177 /* Extend the file to account for all the imagery space. */
3178 /* -------------------------------------------------------------------- */
3179 GIntBig nTileDataSize = ((GIntBig) nBytesPerBlock)
3180 * nBlocksPerRow * nBlocksPerColumn * nLayers;
3181
3182 *pnDataOffset = VSIFTellL( fpVSIL );
3183
3184 if( VSIFSeekL( fpVSIL, nTileDataSize - 1 + *pnDataOffset, SEEK_SET ) != 0
3185 || VSIFWriteL( (void *) "", 1, 1, fpVSIL ) != 1 )
3186 {
3187 CPLError( CE_Failure, CPLE_FileIO,
3188 "Failed to extend %s to full size (%g bytes),\n"
3189 "likely out of disk space.\n%s",
3190 psInfo->pszIGEFilename,
3191 (double) nTileDataSize - 1 + *pnDataOffset,
3192 VSIStrerror( errno ) );
3193
3194 VSIFCloseL( fpVSIL );
3195 return FALSE;
3196 }
3197
3198 VSIFCloseL( fpVSIL );
3199
3200 return TRUE;
3201 }
3202
3203 /************************************************************************/
3204 /* HFAReadAndValidatePoly() */
3205 /************************************************************************/
3206
HFAReadAndValidatePoly(HFAEntry * poTarget,const char * pszName,Efga_Polynomial * psRetPoly)3207 static int HFAReadAndValidatePoly( HFAEntry *poTarget,
3208 const char *pszName,
3209 Efga_Polynomial *psRetPoly )
3210
3211 {
3212 CPLString osFldName;
3213
3214 memset( psRetPoly, 0, sizeof(Efga_Polynomial) );
3215
3216 osFldName.Printf( "%sorder", pszName );
3217 psRetPoly->order = poTarget->GetIntField(osFldName);
3218
3219 if( psRetPoly->order < 1 || psRetPoly->order > 3 )
3220 return FALSE;
3221
3222 /* -------------------------------------------------------------------- */
3223 /* Validate that things are in a "well known" form. */
3224 /* -------------------------------------------------------------------- */
3225 int numdimtransform, numdimpolynomial, termcount;
3226
3227 osFldName.Printf( "%snumdimtransform", pszName );
3228 numdimtransform = poTarget->GetIntField(osFldName);
3229
3230 osFldName.Printf( "%snumdimpolynomial", pszName );
3231 numdimpolynomial = poTarget->GetIntField(osFldName);
3232
3233 osFldName.Printf( "%stermcount", pszName );
3234 termcount = poTarget->GetIntField(osFldName);
3235
3236 if( numdimtransform != 2 || numdimpolynomial != 2 )
3237 return FALSE;
3238
3239 if( (psRetPoly->order == 1 && termcount != 3)
3240 || (psRetPoly->order == 2 && termcount != 6)
3241 || (psRetPoly->order == 3 && termcount != 10) )
3242 return FALSE;
3243
3244 // we don't check the exponent organization for now. Hopefully
3245 // it is always standard.
3246
3247 /* -------------------------------------------------------------------- */
3248 /* Get coefficients. */
3249 /* -------------------------------------------------------------------- */
3250 int i;
3251
3252 for( i = 0; i < termcount*2 - 2; i++ )
3253 {
3254 osFldName.Printf( "%spolycoefmtx[%d]", pszName, i );
3255 psRetPoly->polycoefmtx[i] = poTarget->GetDoubleField(osFldName);
3256 }
3257
3258 for( i = 0; i < 2; i++ )
3259 {
3260 osFldName.Printf( "%spolycoefvector[%d]", pszName, i );
3261 psRetPoly->polycoefvector[i] = poTarget->GetDoubleField(osFldName);
3262 }
3263
3264 return TRUE;
3265 }
3266
3267 /************************************************************************/
3268 /* HFAReadXFormStack() */
3269 /************************************************************************/
3270
3271
HFAReadXFormStack(HFAHandle hHFA,Efga_Polynomial ** ppasPolyListForward,Efga_Polynomial ** ppasPolyListReverse)3272 int HFAReadXFormStack( HFAHandle hHFA,
3273 Efga_Polynomial **ppasPolyListForward,
3274 Efga_Polynomial **ppasPolyListReverse )
3275
3276 {
3277 if( hHFA->nBands == 0 )
3278 return 0;
3279
3280 /* -------------------------------------------------------------------- */
3281 /* Get the HFA node. */
3282 /* -------------------------------------------------------------------- */
3283 HFAEntry *poXFormHeader;
3284
3285 poXFormHeader = hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm" );
3286 if( poXFormHeader == NULL )
3287 return 0;
3288
3289 /* -------------------------------------------------------------------- */
3290 /* Loop over children, collecting XForms. */
3291 /* -------------------------------------------------------------------- */
3292 HFAEntry *poXForm;
3293 int nStepCount = 0;
3294 *ppasPolyListForward = NULL;
3295 *ppasPolyListReverse = NULL;
3296
3297 for( poXForm = poXFormHeader->GetChild();
3298 poXForm != NULL;
3299 poXForm = poXForm->GetNext() )
3300 {
3301 int bSuccess = FALSE;
3302 Efga_Polynomial sForward, sReverse;
3303 memset( &sForward, 0, sizeof(sForward) );
3304 memset( &sReverse, 0, sizeof(sReverse) );
3305
3306 if( EQUAL(poXForm->GetType(),"Efga_Polynomial") )
3307 {
3308 bSuccess =
3309 HFAReadAndValidatePoly( poXForm, "", &sForward );
3310
3311 if( bSuccess )
3312 {
3313 double adfGT[6], adfInvGT[6];
3314
3315 adfGT[0] = sForward.polycoefvector[0];
3316 adfGT[1] = sForward.polycoefmtx[0];
3317 adfGT[2] = sForward.polycoefmtx[2];
3318 adfGT[3] = sForward.polycoefvector[1];
3319 adfGT[4] = sForward.polycoefmtx[1];
3320 adfGT[5] = sForward.polycoefmtx[3];
3321
3322 bSuccess = HFAInvGeoTransform( adfGT, adfInvGT );
3323
3324 sReverse.order = sForward.order;
3325 sReverse.polycoefvector[0] = adfInvGT[0];
3326 sReverse.polycoefmtx[0] = adfInvGT[1];
3327 sReverse.polycoefmtx[2] = adfInvGT[2];
3328 sReverse.polycoefvector[1] = adfInvGT[3];
3329 sReverse.polycoefmtx[1] = adfInvGT[4];
3330 sReverse.polycoefmtx[3] = adfInvGT[5];
3331 }
3332 }
3333 else if( EQUAL(poXForm->GetType(),"GM_PolyPair") )
3334 {
3335 bSuccess =
3336 HFAReadAndValidatePoly( poXForm, "forward.", &sForward );
3337 bSuccess = bSuccess &&
3338 HFAReadAndValidatePoly( poXForm, "reverse.", &sReverse );
3339 }
3340
3341 if( bSuccess )
3342 {
3343 nStepCount++;
3344 *ppasPolyListForward = (Efga_Polynomial *)
3345 CPLRealloc( *ppasPolyListForward,
3346 sizeof(Efga_Polynomial) * nStepCount);
3347 memcpy( *ppasPolyListForward + nStepCount - 1,
3348 &sForward, sizeof(sForward) );
3349
3350 *ppasPolyListReverse = (Efga_Polynomial *)
3351 CPLRealloc( *ppasPolyListReverse,
3352 sizeof(Efga_Polynomial) * nStepCount);
3353 memcpy( *ppasPolyListReverse + nStepCount - 1,
3354 &sReverse, sizeof(sReverse) );
3355 }
3356 }
3357
3358 return nStepCount;
3359 }
3360
3361 /************************************************************************/
3362 /* HFAEvaluateXFormStack() */
3363 /************************************************************************/
3364
HFAEvaluateXFormStack(int nStepCount,int bForward,Efga_Polynomial * pasPolyList,double * pdfX,double * pdfY)3365 int HFAEvaluateXFormStack( int nStepCount, int bForward,
3366 Efga_Polynomial *pasPolyList,
3367 double *pdfX, double *pdfY )
3368
3369 {
3370 int iStep;
3371
3372 for( iStep = 0; iStep < nStepCount; iStep++ )
3373 {
3374 double dfXOut, dfYOut;
3375 Efga_Polynomial *psStep;
3376
3377 if( bForward )
3378 psStep = pasPolyList + iStep;
3379 else
3380 psStep = pasPolyList + nStepCount - iStep - 1;
3381
3382 if( psStep->order == 1 )
3383 {
3384 dfXOut = psStep->polycoefvector[0]
3385 + psStep->polycoefmtx[0] * *pdfX
3386 + psStep->polycoefmtx[2] * *pdfY;
3387
3388 dfYOut = psStep->polycoefvector[1]
3389 + psStep->polycoefmtx[1] * *pdfX
3390 + psStep->polycoefmtx[3] * *pdfY;
3391
3392 *pdfX = dfXOut;
3393 *pdfY = dfYOut;
3394 }
3395 else if( psStep->order == 2 )
3396 {
3397 dfXOut = psStep->polycoefvector[0]
3398 + psStep->polycoefmtx[0] * *pdfX
3399 + psStep->polycoefmtx[2] * *pdfY
3400 + psStep->polycoefmtx[4] * *pdfX * *pdfX
3401 + psStep->polycoefmtx[6] * *pdfX * *pdfY
3402 + psStep->polycoefmtx[8] * *pdfY * *pdfY;
3403 dfYOut = psStep->polycoefvector[1]
3404 + psStep->polycoefmtx[1] * *pdfX
3405 + psStep->polycoefmtx[3] * *pdfY
3406 + psStep->polycoefmtx[5] * *pdfX * *pdfX
3407 + psStep->polycoefmtx[7] * *pdfX * *pdfY
3408 + psStep->polycoefmtx[9] * *pdfY * *pdfY;
3409
3410 *pdfX = dfXOut;
3411 *pdfY = dfYOut;
3412 }
3413 else if( psStep->order == 3 )
3414 {
3415 dfXOut = psStep->polycoefvector[0]
3416 + psStep->polycoefmtx[ 0] * *pdfX
3417 + psStep->polycoefmtx[ 2] * *pdfY
3418 + psStep->polycoefmtx[ 4] * *pdfX * *pdfX
3419 + psStep->polycoefmtx[ 6] * *pdfX * *pdfY
3420 + psStep->polycoefmtx[ 8] * *pdfY * *pdfY
3421 + psStep->polycoefmtx[10] * *pdfX * *pdfX * *pdfX
3422 + psStep->polycoefmtx[12] * *pdfX * *pdfX * *pdfY
3423 + psStep->polycoefmtx[14] * *pdfX * *pdfY * *pdfY
3424 + psStep->polycoefmtx[16] * *pdfY * *pdfY * *pdfY;
3425 dfYOut = psStep->polycoefvector[1]
3426 + psStep->polycoefmtx[ 1] * *pdfX
3427 + psStep->polycoefmtx[ 3] * *pdfY
3428 + psStep->polycoefmtx[ 5] * *pdfX * *pdfX
3429 + psStep->polycoefmtx[ 7] * *pdfX * *pdfY
3430 + psStep->polycoefmtx[ 9] * *pdfY * *pdfY
3431 + psStep->polycoefmtx[11] * *pdfX * *pdfX * *pdfX
3432 + psStep->polycoefmtx[13] * *pdfX * *pdfX * *pdfY
3433 + psStep->polycoefmtx[15] * *pdfX * *pdfY * *pdfY
3434 + psStep->polycoefmtx[17] * *pdfY * *pdfY * *pdfY;
3435
3436 *pdfX = dfXOut;
3437 *pdfY = dfYOut;
3438 }
3439 else
3440 return FALSE;
3441 }
3442
3443 return TRUE;
3444 }
3445
3446 /************************************************************************/
3447 /* HFAWriteXFormStack() */
3448 /************************************************************************/
3449
HFAWriteXFormStack(HFAHandle hHFA,int nBand,int nXFormCount,Efga_Polynomial ** ppasPolyListForward,Efga_Polynomial ** ppasPolyListReverse)3450 CPLErr HFAWriteXFormStack( HFAHandle hHFA, int nBand, int nXFormCount,
3451 Efga_Polynomial **ppasPolyListForward,
3452 Efga_Polynomial **ppasPolyListReverse )
3453
3454 {
3455 if( nXFormCount == 0 )
3456 return CE_None;
3457
3458 if( ppasPolyListForward[0]->order != 1 )
3459 {
3460 CPLError( CE_Failure, CPLE_AppDefined,
3461 "For now HFAWriteXFormStack() only supports order 1 polynomials" );
3462 return CE_Failure;
3463 }
3464
3465 if( nBand < 0 || nBand > hHFA->nBands )
3466 return CE_Failure;
3467
3468 /* -------------------------------------------------------------------- */
3469 /* If no band number is provided, operate on all bands. */
3470 /* -------------------------------------------------------------------- */
3471 if( nBand == 0 )
3472 {
3473 CPLErr eErr = CE_None;
3474
3475 for( nBand = 1; nBand <= hHFA->nBands; nBand++ )
3476 {
3477 eErr = HFAWriteXFormStack( hHFA, nBand, nXFormCount,
3478 ppasPolyListForward,
3479 ppasPolyListReverse );
3480 if( eErr != CE_None )
3481 return eErr;
3482 }
3483
3484 return eErr;
3485 }
3486
3487 /* -------------------------------------------------------------------- */
3488 /* Fetch our band node. */
3489 /* -------------------------------------------------------------------- */
3490 HFAEntry *poBandNode = hHFA->papoBand[nBand-1]->poNode;
3491 HFAEntry *poXFormHeader;
3492
3493 poXFormHeader = poBandNode->GetNamedChild( "MapToPixelXForm" );
3494 if( poXFormHeader == NULL )
3495 {
3496 poXFormHeader = new HFAEntry( hHFA, "MapToPixelXForm",
3497 "Exfr_GenericXFormHeader", poBandNode );
3498 poXFormHeader->MakeData( 23 );
3499 poXFormHeader->SetPosition();
3500 poXFormHeader->SetStringField( "titleList.string", "Affine" );
3501 }
3502
3503 /* -------------------------------------------------------------------- */
3504 /* Loop over XForms. */
3505 /* -------------------------------------------------------------------- */
3506 for( int iXForm = 0; iXForm < nXFormCount; iXForm++ )
3507 {
3508 Efga_Polynomial *psForward = *ppasPolyListForward + iXForm;
3509 CPLString osXFormName;
3510 osXFormName.Printf( "XForm%d", iXForm );
3511
3512 HFAEntry *poXForm = poXFormHeader->GetNamedChild( osXFormName );
3513
3514 if( poXForm == NULL )
3515 {
3516 poXForm = new HFAEntry( hHFA, osXFormName, "Efga_Polynomial",
3517 poXFormHeader );
3518 poXForm->MakeData( 136 );
3519 poXForm->SetPosition();
3520 }
3521
3522 poXForm->SetIntField( "order", 1 );
3523 poXForm->SetIntField( "numdimtransform", 2 );
3524 poXForm->SetIntField( "numdimpolynomial", 2 );
3525 poXForm->SetIntField( "termcount", 3 );
3526 poXForm->SetIntField( "exponentlist[0]", 0 );
3527 poXForm->SetIntField( "exponentlist[1]", 0 );
3528 poXForm->SetIntField( "exponentlist[2]", 1 );
3529 poXForm->SetIntField( "exponentlist[3]", 0 );
3530 poXForm->SetIntField( "exponentlist[4]", 0 );
3531 poXForm->SetIntField( "exponentlist[5]", 1 );
3532
3533 poXForm->SetIntField( "polycoefmtx[-3]", EPT_f64 );
3534 poXForm->SetIntField( "polycoefmtx[-2]", 2 );
3535 poXForm->SetIntField( "polycoefmtx[-1]", 2 );
3536 poXForm->SetDoubleField( "polycoefmtx[0]",
3537 psForward->polycoefmtx[0] );
3538 poXForm->SetDoubleField( "polycoefmtx[1]",
3539 psForward->polycoefmtx[1] );
3540 poXForm->SetDoubleField( "polycoefmtx[2]",
3541 psForward->polycoefmtx[2] );
3542 poXForm->SetDoubleField( "polycoefmtx[3]",
3543 psForward->polycoefmtx[3] );
3544
3545 poXForm->SetIntField( "polycoefvector[-3]", EPT_f64 );
3546 poXForm->SetIntField( "polycoefvector[-2]", 1 );
3547 poXForm->SetIntField( "polycoefvector[-1]", 2 );
3548 poXForm->SetDoubleField( "polycoefvector[0]",
3549 psForward->polycoefvector[0] );
3550 poXForm->SetDoubleField( "polycoefvector[1]",
3551 psForward->polycoefvector[1] );
3552 }
3553
3554 return CE_None;
3555 }
3556
3557 /************************************************************************/
3558 /* HFAReadCameraModel() */
3559 /************************************************************************/
3560
HFAReadCameraModel(HFAHandle hHFA)3561 char **HFAReadCameraModel( HFAHandle hHFA )
3562
3563 {
3564 if( hHFA->nBands == 0 )
3565 return NULL;
3566
3567 /* -------------------------------------------------------------------- */
3568 /* Get the camera model node, and confirm it's type. */
3569 /* -------------------------------------------------------------------- */
3570 HFAEntry *poXForm;
3571
3572 poXForm =
3573 hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm.XForm0" );
3574 if( poXForm == NULL )
3575 return NULL;
3576
3577 if( !EQUAL(poXForm->GetType(),"Camera_ModelX") )
3578 return NULL;
3579
3580 /* -------------------------------------------------------------------- */
3581 /* Convert the values to metadata. */
3582 /* -------------------------------------------------------------------- */
3583 const char *pszValue;
3584 int i;
3585 char **papszMD = NULL;
3586 static const char *apszFields[] = {
3587 "direction", "refType", "demsource", "PhotoDirection", "RotationSystem",
3588 "demfilename", "demzunits",
3589 "forSrcAffine[0]", "forSrcAffine[1]", "forSrcAffine[2]",
3590 "forSrcAffine[3]", "forSrcAffine[4]", "forSrcAffine[5]",
3591 "forDstAffine[0]", "forDstAffine[1]", "forDstAffine[2]",
3592 "forDstAffine[3]", "forDstAffine[4]", "forDstAffine[5]",
3593 "invSrcAffine[0]", "invSrcAffine[1]", "invSrcAffine[2]",
3594 "invSrcAffine[3]", "invSrcAffine[4]", "invSrcAffine[5]",
3595 "invDstAffine[0]", "invDstAffine[1]", "invDstAffine[2]",
3596 "invDstAffine[3]", "invDstAffine[4]", "invDstAffine[5]",
3597 "z_mean", "lat0", "lon0",
3598 "coeffs[0]", "coeffs[1]", "coeffs[2]",
3599 "coeffs[3]", "coeffs[4]", "coeffs[5]",
3600 "coeffs[6]", "coeffs[7]", "coeffs[8]",
3601 "LensDistortion[0]", "LensDistortion[1]", "LensDistortion[2]",
3602 NULL };
3603
3604 for( i = 0; apszFields[i] != NULL; i++ )
3605 {
3606 pszValue = poXForm->GetStringField( apszFields[i] );
3607 if( pszValue == NULL )
3608 pszValue = "";
3609
3610 papszMD = CSLSetNameValue( papszMD, apszFields[i], pszValue );
3611 }
3612
3613 /* -------------------------------------------------------------------- */
3614 /* Create a pseudo-entry for the MIFObject with the */
3615 /* outputProjection. */
3616 /* -------------------------------------------------------------------- */
3617 HFAEntry *poProjInfo = HFAEntry::BuildEntryFromMIFObject( poXForm, "outputProjection" );
3618 if (poProjInfo)
3619 {
3620 /* -------------------------------------------------------------------- */
3621 /* Fetch the datum. */
3622 /* -------------------------------------------------------------------- */
3623 Eprj_Datum sDatum;
3624
3625 memset( &sDatum, 0, sizeof(sDatum));
3626
3627 sDatum.datumname =
3628 (char *) poProjInfo->GetStringField("earthModel.datum.datumname");
3629 sDatum.type = (Eprj_DatumType) poProjInfo->GetIntField(
3630 "earthModel.datum.type");
3631
3632 for( i = 0; i < 7; i++ )
3633 {
3634 char szFieldName[60];
3635
3636 sprintf( szFieldName, "earthModel.datum.params[%d]", i );
3637 sDatum.params[i] = poProjInfo->GetDoubleField(szFieldName);
3638 }
3639
3640 sDatum.gridname = (char *)
3641 poProjInfo->GetStringField("earthModel.datum.gridname");
3642
3643 /* -------------------------------------------------------------------- */
3644 /* Fetch the projection parameters. */
3645 /* -------------------------------------------------------------------- */
3646 Eprj_ProParameters sPro;
3647
3648 memset( &sPro, 0, sizeof(sPro) );
3649
3650 sPro.proType = (Eprj_ProType) poProjInfo->GetIntField("projectionObject.proType");
3651 sPro.proNumber = poProjInfo->GetIntField("projectionObject.proNumber");
3652 sPro.proExeName = (char *) poProjInfo->GetStringField("projectionObject.proExeName");
3653 sPro.proName = (char *) poProjInfo->GetStringField("projectionObject.proName");
3654 sPro.proZone = poProjInfo->GetIntField("projectionObject.proZone");
3655
3656 for( i = 0; i < 15; i++ )
3657 {
3658 char szFieldName[40];
3659
3660 sprintf( szFieldName, "projectionObject.proParams[%d]", i );
3661 sPro.proParams[i] = poProjInfo->GetDoubleField(szFieldName);
3662 }
3663
3664 /* -------------------------------------------------------------------- */
3665 /* Fetch the spheroid. */
3666 /* -------------------------------------------------------------------- */
3667 sPro.proSpheroid.sphereName = (char *)
3668 poProjInfo->GetStringField("earthModel.proSpheroid.sphereName");
3669 sPro.proSpheroid.a = poProjInfo->GetDoubleField("earthModel.proSpheroid.a");
3670 sPro.proSpheroid.b = poProjInfo->GetDoubleField("earthModel.proSpheroid.b");
3671 sPro.proSpheroid.eSquared =
3672 poProjInfo->GetDoubleField("earthModel.proSpheroid.eSquared");
3673 sPro.proSpheroid.radius =
3674 poProjInfo->GetDoubleField("earthModel.proSpheroid.radius");
3675
3676 /* -------------------------------------------------------------------- */
3677 /* Fetch the projection info. */
3678 /* -------------------------------------------------------------------- */
3679 char *pszProjection;
3680
3681 // poProjInfo->DumpFieldValues( stdout, "" );
3682
3683 pszProjection = HFAPCSStructToWKT( &sDatum, &sPro, NULL, NULL );
3684
3685 if( pszProjection )
3686 {
3687 papszMD =
3688 CSLSetNameValue( papszMD, "outputProjection", pszProjection );
3689 CPLFree( pszProjection );
3690 }
3691
3692 delete poProjInfo;
3693 }
3694
3695 /* -------------------------------------------------------------------- */
3696 /* Fetch the horizontal units. */
3697 /* -------------------------------------------------------------------- */
3698 pszValue = poXForm->GetStringField( "outputHorizontalUnits.string" );
3699 if( pszValue == NULL )
3700 pszValue = "";
3701
3702 papszMD = CSLSetNameValue( papszMD, "outputHorizontalUnits", pszValue );
3703
3704 /* -------------------------------------------------------------------- */
3705 /* Fetch the elevationinfo. */
3706 /* -------------------------------------------------------------------- */
3707 HFAEntry *poElevInfo = HFAEntry::BuildEntryFromMIFObject( poXForm, "outputElevationInfo" );
3708 if ( poElevInfo )
3709 {
3710 //poElevInfo->DumpFieldValues( stdout, "" );
3711
3712 if( poElevInfo->GetDataSize() != 0 )
3713 {
3714 static const char *apszEFields[] = {
3715 "verticalDatum.datumname",
3716 "verticalDatum.type",
3717 "elevationUnit",
3718 "elevationType",
3719 NULL };
3720
3721 for( i = 0; apszEFields[i] != NULL; i++ )
3722 {
3723 pszValue = poElevInfo->GetStringField( apszEFields[i] );
3724 if( pszValue == NULL )
3725 pszValue = "";
3726
3727 papszMD = CSLSetNameValue( papszMD, apszEFields[i], pszValue );
3728 }
3729 }
3730
3731 delete poElevInfo;
3732 }
3733
3734 return papszMD;
3735 }
3736
3737 /************************************************************************/
3738 /* HFASetGeoTransform() */
3739 /* */
3740 /* Set a MapInformation and XForm block. Allows for rotated */
3741 /* and shared geotransforms. */
3742 /************************************************************************/
3743
HFASetGeoTransform(HFAHandle hHFA,const char * pszProName,const char * pszUnits,double * padfGeoTransform)3744 CPLErr HFASetGeoTransform( HFAHandle hHFA,
3745 const char *pszProName,
3746 const char *pszUnits,
3747 double *padfGeoTransform )
3748
3749 {
3750 /* -------------------------------------------------------------------- */
3751 /* Write MapInformation. */
3752 /* -------------------------------------------------------------------- */
3753 int nBand;
3754
3755 for( nBand = 1; nBand <= hHFA->nBands; nBand++ )
3756 {
3757 HFAEntry *poBandNode = hHFA->papoBand[nBand-1]->poNode;
3758
3759 HFAEntry *poMI = poBandNode->GetNamedChild( "MapInformation" );
3760 if( poMI == NULL )
3761 {
3762 poMI = new HFAEntry( hHFA, "MapInformation",
3763 "Eimg_MapInformation", poBandNode );
3764 poMI->MakeData( 18 + strlen(pszProName) + strlen(pszUnits) );
3765 poMI->SetPosition();
3766 }
3767
3768 poMI->SetStringField( "projection.string", pszProName );
3769 poMI->SetStringField( "units.string", pszUnits );
3770 }
3771
3772 /* -------------------------------------------------------------------- */
3773 /* Write XForm. */
3774 /* -------------------------------------------------------------------- */
3775 Efga_Polynomial sForward, sReverse;
3776 double adfAdjTransform[6], adfRevTransform[6];
3777
3778 // Offset by half pixel.
3779
3780 memcpy( adfAdjTransform, padfGeoTransform, sizeof(double) * 6 );
3781 adfAdjTransform[0] += adfAdjTransform[1] * 0.5;
3782 adfAdjTransform[0] += adfAdjTransform[2] * 0.5;
3783 adfAdjTransform[3] += adfAdjTransform[4] * 0.5;
3784 adfAdjTransform[3] += adfAdjTransform[5] * 0.5;
3785
3786 // Invert
3787 HFAInvGeoTransform( adfAdjTransform, adfRevTransform );
3788
3789 // Assign to polynomial object.
3790
3791 sForward.order = 1;
3792 sForward.polycoefvector[0] = adfRevTransform[0];
3793 sForward.polycoefmtx[0] = adfRevTransform[1];
3794 sForward.polycoefmtx[1] = adfRevTransform[4];
3795 sForward.polycoefvector[1] = adfRevTransform[3];
3796 sForward.polycoefmtx[2] = adfRevTransform[2];
3797 sForward.polycoefmtx[3] = adfRevTransform[5];
3798
3799 sReverse = sForward;
3800 Efga_Polynomial *psForward=&sForward, *psReverse=&sReverse;
3801
3802 return HFAWriteXFormStack( hHFA, 0, 1, &psForward, &psReverse );
3803 }
3804
3805 /************************************************************************/
3806 /* HFARenameReferences() */
3807 /* */
3808 /* Rename references in this .img file from the old basename to */
3809 /* a new basename. This should be passed on to .aux and .rrd */
3810 /* files and should include references to .aux, .rrd and .ige. */
3811 /************************************************************************/
3812
HFARenameReferences(HFAHandle hHFA,const char * pszNewBase,const char * pszOldBase)3813 CPLErr HFARenameReferences( HFAHandle hHFA,
3814 const char *pszNewBase,
3815 const char *pszOldBase )
3816
3817 {
3818 /* -------------------------------------------------------------------- */
3819 /* Handle RRDNamesList updates. */
3820 /* -------------------------------------------------------------------- */
3821 size_t iNode;
3822 std::vector<HFAEntry*> apoNodeList =
3823 hHFA->poRoot->FindChildren( "RRDNamesList", NULL );
3824
3825 for( iNode = 0; iNode < apoNodeList.size(); iNode++ )
3826 {
3827 HFAEntry *poRRDNL = apoNodeList[iNode];
3828 std::vector<CPLString> aosNL;
3829
3830 // Collect all the existing names.
3831 int i, nNameCount = poRRDNL->GetFieldCount( "nameList" );
3832
3833 CPLString osAlgorithm = poRRDNL->GetStringField("algorithm.string");
3834 for( i = 0; i < nNameCount; i++ )
3835 {
3836 CPLString osFN;
3837 osFN.Printf( "nameList[%d].string", i );
3838 aosNL.push_back( poRRDNL->GetStringField(osFN) );
3839 }
3840
3841 // Adjust the names to the new form.
3842 for( i = 0; i < nNameCount; i++ )
3843 {
3844 if( strncmp(aosNL[i],pszOldBase,strlen(pszOldBase)) == 0 )
3845 {
3846 CPLString osNew = pszNewBase;
3847 osNew += aosNL[i].c_str() + strlen(pszOldBase);
3848 aosNL[i] = osNew;
3849 }
3850 }
3851
3852 // try to make sure the RRDNamesList is big enough to hold the
3853 // adjusted name list.
3854 if( strlen(pszNewBase) > strlen(pszOldBase) )
3855 {
3856 CPLDebug( "HFA", "Growing RRDNamesList to hold new names" );
3857 poRRDNL->MakeData( poRRDNL->GetDataSize()
3858 + nNameCount * (strlen(pszNewBase) - strlen(pszOldBase)) );
3859 }
3860
3861 // Initialize the whole thing to zeros for a clean start.
3862 memset( poRRDNL->GetData(), 0, poRRDNL->GetDataSize() );
3863
3864 // Write the updates back to the file.
3865 poRRDNL->SetStringField( "algorithm.string", osAlgorithm );
3866 for( i = 0; i < nNameCount; i++ )
3867 {
3868 CPLString osFN;
3869 osFN.Printf( "nameList[%d].string", i );
3870 poRRDNL->SetStringField( osFN, aosNL[i] );
3871 }
3872 }
3873
3874 /* -------------------------------------------------------------------- */
3875 /* spill file references. */
3876 /* -------------------------------------------------------------------- */
3877 apoNodeList =
3878 hHFA->poRoot->FindChildren( "ExternalRasterDMS", "ImgExternalRaster" );
3879
3880 for( iNode = 0; iNode < apoNodeList.size(); iNode++ )
3881 {
3882 HFAEntry *poERDMS = apoNodeList[iNode];
3883
3884 if( poERDMS == NULL )
3885 continue;
3886
3887 // Fetch all existing values.
3888 CPLString osFileName = poERDMS->GetStringField("fileName.string");
3889 GInt32 anValidFlagsOffset[2], anStackDataOffset[2];
3890 GInt32 nStackCount, nStackIndex;
3891
3892 anValidFlagsOffset[0] =
3893 poERDMS->GetIntField( "layerStackValidFlagsOffset[0]" );
3894 anValidFlagsOffset[1] =
3895 poERDMS->GetIntField( "layerStackValidFlagsOffset[1]" );
3896
3897 anStackDataOffset[0] =
3898 poERDMS->GetIntField( "layerStackDataOffset[0]" );
3899 anStackDataOffset[1] =
3900 poERDMS->GetIntField( "layerStackDataOffset[1]" );
3901
3902 nStackCount = poERDMS->GetIntField( "layerStackCount" );
3903 nStackIndex = poERDMS->GetIntField( "layerStackIndex" );
3904
3905 // Update the filename.
3906 if( strncmp(osFileName,pszOldBase,strlen(pszOldBase)) == 0 )
3907 {
3908 CPLString osNew = pszNewBase;
3909 osNew += osFileName.c_str() + strlen(pszOldBase);
3910 osFileName = osNew;
3911 }
3912
3913 // Grow the node if needed.
3914 if( strlen(pszNewBase) > strlen(pszOldBase) )
3915 {
3916 CPLDebug( "HFA", "Growing ExternalRasterDMS to hold new names" );
3917 poERDMS->MakeData( poERDMS->GetDataSize()
3918 + (strlen(pszNewBase) - strlen(pszOldBase)) );
3919 }
3920
3921 // Initialize the whole thing to zeros for a clean start.
3922 memset( poERDMS->GetData(), 0, poERDMS->GetDataSize() );
3923
3924 // Write it all out again, this may change the size of the node.
3925 poERDMS->SetStringField( "fileName.string", osFileName );
3926 poERDMS->SetIntField( "layerStackValidFlagsOffset[0]",
3927 anValidFlagsOffset[0] );
3928 poERDMS->SetIntField( "layerStackValidFlagsOffset[1]",
3929 anValidFlagsOffset[1] );
3930
3931 poERDMS->SetIntField( "layerStackDataOffset[0]",
3932 anStackDataOffset[0] );
3933 poERDMS->SetIntField( "layerStackDataOffset[1]",
3934 anStackDataOffset[1] );
3935
3936 poERDMS->SetIntField( "layerStackCount", nStackCount );
3937 poERDMS->SetIntField( "layerStackIndex", nStackIndex );
3938 }
3939
3940 /* -------------------------------------------------------------------- */
3941 /* DependentFile */
3942 /* -------------------------------------------------------------------- */
3943 apoNodeList =
3944 hHFA->poRoot->FindChildren( "DependentFile", "Eimg_DependentFile" );
3945
3946 for( iNode = 0; iNode < apoNodeList.size(); iNode++ )
3947 {
3948 CPLString osFileName = apoNodeList[iNode]->
3949 GetStringField("dependent.string");
3950
3951 // Grow the node if needed.
3952 if( strlen(pszNewBase) > strlen(pszOldBase) )
3953 {
3954 CPLDebug( "HFA", "Growing DependentFile to hold new names" );
3955 apoNodeList[iNode]->MakeData( apoNodeList[iNode]->GetDataSize()
3956 + (strlen(pszNewBase)
3957 - strlen(pszOldBase)) );
3958 }
3959
3960 // Update the filename.
3961 if( strncmp(osFileName,pszOldBase,strlen(pszOldBase)) == 0 )
3962 {
3963 CPLString osNew = pszNewBase;
3964 osNew += osFileName.c_str() + strlen(pszOldBase);
3965 osFileName = osNew;
3966 }
3967
3968 apoNodeList[iNode]->SetStringField( "dependent.string", osFileName );
3969 }
3970
3971 return CE_None;
3972 }
3973