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