1 : /******************************************************************************
2 : * $Id: ehdrdataset.cpp 21783 2011-02-21 22:32:11Z rouault $
3 : *
4 : * Project: ESRI .hdr Driver
5 : * Purpose: Implementation of EHdrDataset
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam <warmerdam@pobox.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "ogr_spatialref.h"
32 : #include "cpl_string.h"
33 :
34 : CPL_CVSID("$Id: ehdrdataset.cpp 21783 2011-02-21 22:32:11Z rouault $");
35 :
36 : CPL_C_START
37 : void GDALRegister_EHdr(void);
38 : CPL_C_END
39 :
40 : #define HAS_MIN_FLAG 0x1
41 : #define HAS_MAX_FLAG 0x2
42 : #define HAS_MEAN_FLAG 0x4
43 : #define HAS_STDDEV_FLAG 0x8
44 : #define HAS_ALL_FLAGS (HAS_MIN_FLAG | HAS_MAX_FLAG | HAS_MEAN_FLAG | HAS_STDDEV_FLAG)
45 :
46 : /************************************************************************/
47 : /* ==================================================================== */
48 : /* EHdrDataset */
49 : /* ==================================================================== */
50 : /************************************************************************/
51 :
52 : class EHdrRasterBand;
53 :
54 : class EHdrDataset : public RawDataset
55 : {
56 : friend class EHdrRasterBand;
57 :
58 : VSILFILE *fpImage; // image data file.
59 :
60 : CPLString osHeaderExt;
61 :
62 : int bGotTransform;
63 : double adfGeoTransform[6];
64 : char *pszProjection;
65 :
66 : int bHDRDirty;
67 : char **papszHDR;
68 :
69 : int bCLRDirty;
70 :
71 : CPLErr ReadSTX();
72 : CPLErr RewriteSTX();
73 : CPLErr RewriteHDR();
74 : void ResetKeyValue( const char *pszKey, const char *pszValue );
75 : const char *GetKeyValue( const char *pszKey, const char *pszDefault = "" );
76 : void RewriteColorTable( GDALColorTable * );
77 :
78 : public:
79 : EHdrDataset();
80 : ~EHdrDataset();
81 :
82 : virtual CPLErr GetGeoTransform( double * padfTransform );
83 : virtual CPLErr SetGeoTransform( double *padfTransform );
84 : virtual const char *GetProjectionRef(void);
85 : virtual CPLErr SetProjection( const char * );
86 :
87 : virtual char **GetFileList();
88 :
89 : static GDALDataset *Open( GDALOpenInfo * );
90 : static GDALDataset *Create( const char * pszFilename,
91 : int nXSize, int nYSize, int nBands,
92 : GDALDataType eType, char ** papszParmList );
93 : static GDALDataset *CreateCopy( const char * pszFilename,
94 : GDALDataset * poSrcDS,
95 : int bStrict, char ** papszOptions,
96 : GDALProgressFunc pfnProgress,
97 : void * pProgressData );
98 : static CPLString GetImageRepFilename(const char* pszFilename);
99 : };
100 :
101 : /************************************************************************/
102 : /* ==================================================================== */
103 : /* EHdrRasterBand */
104 : /* ==================================================================== */
105 : /************************************************************************/
106 :
107 : class EHdrRasterBand : public RawRasterBand
108 252 : {
109 : friend class EHdrDataset;
110 :
111 : int nBits;
112 : long nStartBit;
113 : int nPixelOffsetBits;
114 : int nLineOffsetBits;
115 :
116 : double dfMin;
117 : double dfMax;
118 : double dfMean;
119 : double dfStdDev;
120 :
121 : int minmaxmeanstddev;
122 :
123 : virtual CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
124 : void *, int, int, GDALDataType,
125 : int, int );
126 :
127 : public:
128 : EHdrRasterBand( GDALDataset *poDS, int nBand, VSILFILE * fpRaw,
129 : vsi_l_offset nImgOffset, int nPixelOffset,
130 : int nLineOffset,
131 : GDALDataType eDataType, int bNativeOrder,
132 : int nBits);
133 :
134 : virtual CPLErr IReadBlock( int, int, void * );
135 : virtual CPLErr IWriteBlock( int, int, void * );
136 :
137 : virtual double GetMinimum( int *pbSuccess = NULL );
138 : virtual double GetMaximum(int *pbSuccess = NULL );
139 : virtual CPLErr GetStatistics( int bApproxOK, int bForce,
140 : double *pdfMin, double *pdfMax,
141 : double *pdfMean, double *pdfStdDev );
142 : virtual CPLErr SetStatistics( double dfMin, double dfMax,
143 : double dfMean, double dfStdDev );
144 : virtual CPLErr SetColorTable( GDALColorTable *poNewCT );
145 :
146 : };
147 :
148 : /************************************************************************/
149 : /* EHdrRasterBand() */
150 : /************************************************************************/
151 :
152 252 : EHdrRasterBand::EHdrRasterBand( GDALDataset *poDS,
153 : int nBand, VSILFILE * fpRaw,
154 : vsi_l_offset nImgOffset, int nPixelOffset,
155 : int nLineOffset,
156 : GDALDataType eDataType, int bNativeOrder,
157 : int nBits)
158 : : RawRasterBand( poDS, nBand, fpRaw, nImgOffset, nPixelOffset, nLineOffset,
159 : eDataType, bNativeOrder, TRUE ),
160 : nBits(nBits),
161 : dfMin(0),
162 : dfMax(0),
163 252 : minmaxmeanstddev(0)
164 : {
165 252 : EHdrDataset* poEDS = (EHdrDataset*)poDS;
166 :
167 252 : if (nBits < 8)
168 : {
169 6 : nStartBit = atoi(poEDS->GetKeyValue("SKIPBYTES")) * 8;
170 6 : if (nBand >= 2)
171 : {
172 0 : long nRowBytes = atoi(poEDS->GetKeyValue("BANDROWBYTES"));
173 0 : if (nRowBytes == 0)
174 0 : nRowBytes = (nBits * poDS->GetRasterXSize() + 7) / 8;
175 :
176 0 : nStartBit += nRowBytes * (nBand-1) * 8;
177 : }
178 :
179 6 : nPixelOffsetBits = nBits;
180 6 : nLineOffsetBits = atoi(poEDS->GetKeyValue("TOTALROWBYTES")) * 8;
181 6 : if( nLineOffsetBits == 0 )
182 0 : nLineOffsetBits = nPixelOffsetBits * poDS->GetRasterXSize();
183 :
184 6 : nBlockXSize = poDS->GetRasterXSize();
185 6 : nBlockYSize = 1;
186 :
187 : SetMetadataItem( "NBITS",
188 : CPLString().Printf( "%d", nBits ),
189 6 : "IMAGE_STRUCTURE" );
190 : }
191 :
192 252 : if( eDataType == GDT_Byte
193 : && EQUAL(poEDS->GetKeyValue("PIXELTYPE",""),"SIGNEDINT") )
194 : SetMetadataItem( "PIXELTYPE", "SIGNEDBYTE",
195 4 : "IMAGE_STRUCTURE" );
196 252 : }
197 :
198 : /************************************************************************/
199 : /* IReadBlock() */
200 : /************************************************************************/
201 :
202 512 : CPLErr EHdrRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
203 : void * pImage )
204 :
205 : {
206 512 : if (nBits >= 8)
207 448 : return RawRasterBand::IReadBlock(nBlockXOff, nBlockYOff, pImage);
208 :
209 : vsi_l_offset nLineStart;
210 : unsigned int nLineBytes;
211 : int iBitOffset;
212 : GByte *pabyBuffer;
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* Establish desired position. */
216 : /* -------------------------------------------------------------------- */
217 64 : nLineBytes = (nPixelOffsetBits*nBlockXSize + 7)/8;
218 64 : nLineStart = (nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) / 8;
219 : iBitOffset = (int)
220 64 : ((nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) % 8);
221 :
222 : /* -------------------------------------------------------------------- */
223 : /* Read data into buffer. */
224 : /* -------------------------------------------------------------------- */
225 64 : pabyBuffer = (GByte *) CPLCalloc(nLineBytes,1);
226 :
227 64 : if( VSIFSeekL( GetFPL(), nLineStart, SEEK_SET ) != 0
228 : || VSIFReadL( pabyBuffer, 1, nLineBytes, GetFPL() ) != nLineBytes )
229 : {
230 : CPLError( CE_Failure, CPLE_FileIO,
231 : "Failed to read %u bytes at offset %lu.\n%s",
232 : nLineBytes, (unsigned long)nLineStart,
233 0 : VSIStrerror( errno ) );
234 0 : return CE_Failure;
235 : }
236 :
237 : /* -------------------------------------------------------------------- */
238 : /* Copy data, promoting to 8bit. */
239 : /* -------------------------------------------------------------------- */
240 64 : int iPixel = 0, iX;
241 :
242 2112 : for( iX = 0; iX < nBlockXSize; iX++ )
243 : {
244 2048 : int nOutWord = 0, iBit;
245 :
246 4096 : for( iBit = 0; iBit < nBits; iBit++ )
247 : {
248 2048 : if( pabyBuffer[iBitOffset>>3] & (0x80 >>(iBitOffset & 7)) )
249 400 : nOutWord |= (1 << (nBits - 1 - iBit));
250 2048 : iBitOffset++;
251 : }
252 :
253 2048 : iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
254 :
255 2048 : ((GByte *) pImage)[iPixel++] = (GByte) nOutWord;
256 : }
257 :
258 64 : CPLFree( pabyBuffer );
259 :
260 64 : return CE_None;
261 : }
262 :
263 : /************************************************************************/
264 : /* IWriteBlock() */
265 : /************************************************************************/
266 :
267 944 : CPLErr EHdrRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
268 : void * pImage )
269 :
270 : {
271 944 : if (nBits >= 8)
272 880 : return RawRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
273 :
274 : vsi_l_offset nLineStart;
275 : unsigned int nLineBytes;
276 : int iBitOffset;
277 : GByte *pabyBuffer;
278 :
279 : /* -------------------------------------------------------------------- */
280 : /* Establish desired position. */
281 : /* -------------------------------------------------------------------- */
282 64 : nLineBytes = (nPixelOffsetBits*nBlockXSize + 7)/8;
283 64 : nLineStart = (nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) / 8;
284 : iBitOffset = (int)
285 64 : ((nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) % 8);
286 :
287 : /* -------------------------------------------------------------------- */
288 : /* Read data into buffer. */
289 : /* -------------------------------------------------------------------- */
290 64 : pabyBuffer = (GByte *) CPLCalloc(nLineBytes,1);
291 :
292 64 : if( VSIFSeekL( GetFPL(), nLineStart, SEEK_SET ) != 0 )
293 : {
294 : CPLError( CE_Failure, CPLE_FileIO,
295 : "Failed to read %u bytes at offset %lu.\n%s",
296 : nLineBytes, (unsigned long)nLineStart,
297 0 : VSIStrerror( errno ) );
298 0 : return CE_Failure;
299 : }
300 :
301 64 : VSIFReadL( pabyBuffer, 1, nLineBytes, GetFPL() );
302 :
303 : /* -------------------------------------------------------------------- */
304 : /* Copy data, promoting to 8bit. */
305 : /* -------------------------------------------------------------------- */
306 64 : int iPixel = 0, iX;
307 :
308 2112 : for( iX = 0; iX < nBlockXSize; iX++ )
309 : {
310 : int iBit;
311 2048 : int nOutWord = ((GByte *) pImage)[iPixel++];
312 :
313 4096 : for( iBit = 0; iBit < nBits; iBit++ )
314 : {
315 2048 : if( nOutWord & (1 << (nBits - 1 - iBit)) )
316 400 : pabyBuffer[iBitOffset>>3] |= (0x80 >>(iBitOffset & 7));
317 : else
318 1648 : pabyBuffer[iBitOffset>>3] &= ~((0x80 >>(iBitOffset & 7)));
319 :
320 2048 : iBitOffset++;
321 : }
322 :
323 2048 : iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
324 : }
325 :
326 : /* -------------------------------------------------------------------- */
327 : /* Write the data back out. */
328 : /* -------------------------------------------------------------------- */
329 64 : if( VSIFSeekL( GetFPL(), nLineStart, SEEK_SET ) != 0
330 : || VSIFWriteL( pabyBuffer, 1, nLineBytes, GetFPL() ) != nLineBytes )
331 : {
332 : CPLError( CE_Failure, CPLE_FileIO,
333 : "Failed to write %u bytes at offset %lu.\n%s",
334 : nLineBytes, (unsigned long)nLineStart,
335 0 : VSIStrerror( errno ) );
336 0 : return CE_Failure;
337 : }
338 :
339 64 : CPLFree( pabyBuffer );
340 :
341 64 : return CE_None;
342 : }
343 :
344 : /************************************************************************/
345 : /* IRasterIO() */
346 : /************************************************************************/
347 :
348 816 : CPLErr EHdrRasterBand::IRasterIO( GDALRWFlag eRWFlag,
349 : int nXOff, int nYOff, int nXSize, int nYSize,
350 : void * pData, int nBufXSize, int nBufYSize,
351 : GDALDataType eBufType,
352 : int nPixelSpace, int nLineSpace )
353 :
354 : {
355 : // Defer to RawRasterBand
356 816 : if (nBits >= 8)
357 : return RawRasterBand::IRasterIO( eRWFlag,
358 : nXOff, nYOff, nXSize, nYSize,
359 : pData, nBufXSize, nBufYSize,
360 750 : eBufType, nPixelSpace, nLineSpace );
361 :
362 : // Force use of IReadBlock() and IWriteBlock()
363 : else
364 : return GDALRasterBand::IRasterIO( eRWFlag,
365 : nXOff, nYOff, nXSize, nYSize,
366 : pData, nBufXSize, nBufYSize,
367 66 : eBufType, nPixelSpace, nLineSpace );
368 : }
369 :
370 : /************************************************************************/
371 : /* OSR_GDS() */
372 : /************************************************************************/
373 :
374 32 : static const char*OSR_GDS( char* pszResult, int nResultLen,
375 : char **papszNV, const char * pszField,
376 : const char *pszDefaultValue )
377 :
378 : {
379 : int iLine;
380 :
381 32 : if( papszNV == NULL || papszNV[0] == NULL )
382 0 : return pszDefaultValue;
383 :
384 250 : for( iLine = 0;
385 124 : papszNV[iLine] != NULL &&
386 94 : !EQUALN(papszNV[iLine],pszField,strlen(pszField));
387 : iLine++ ) {}
388 :
389 32 : if( papszNV[iLine] == NULL )
390 30 : return pszDefaultValue;
391 : else
392 : {
393 : char **papszTokens;
394 :
395 2 : papszTokens = CSLTokenizeString(papszNV[iLine]);
396 :
397 2 : if( CSLCount(papszTokens) > 1 )
398 2 : strncpy( pszResult, papszTokens[1], nResultLen);
399 : else
400 0 : strncpy( pszResult, pszDefaultValue, nResultLen);
401 2 : pszResult[nResultLen-1] = '\0';
402 :
403 2 : CSLDestroy( papszTokens );
404 2 : return pszResult;
405 : }
406 : }
407 :
408 :
409 : /************************************************************************/
410 : /* ==================================================================== */
411 : /* EHdrDataset */
412 : /* ==================================================================== */
413 : /************************************************************************/
414 :
415 : /************************************************************************/
416 : /* EHdrDataset() */
417 : /************************************************************************/
418 :
419 152 : EHdrDataset::EHdrDataset()
420 : {
421 152 : fpImage = NULL;
422 152 : pszProjection = CPLStrdup("");
423 152 : bGotTransform = FALSE;
424 152 : adfGeoTransform[0] = 0.0;
425 152 : adfGeoTransform[1] = 1.0;
426 152 : adfGeoTransform[2] = 0.0;
427 152 : adfGeoTransform[3] = 0.0;
428 152 : adfGeoTransform[4] = 0.0;
429 152 : adfGeoTransform[5] = 1.0;
430 152 : papszHDR = NULL;
431 152 : bHDRDirty = FALSE;
432 152 : bCLRDirty = FALSE;
433 152 : osHeaderExt = "hdr";
434 152 : }
435 :
436 : /************************************************************************/
437 : /* ~EHdrDataset() */
438 : /************************************************************************/
439 :
440 152 : EHdrDataset::~EHdrDataset()
441 :
442 : {
443 152 : FlushCache();
444 :
445 152 : if( nBands > 0 && GetAccess() == GA_Update )
446 : {
447 : int bNoDataSet;
448 : double dfNoData;
449 68 : RawRasterBand *poBand = (RawRasterBand *) GetRasterBand( 1 );
450 :
451 68 : dfNoData = poBand->GetNoDataValue(&bNoDataSet);
452 68 : if( bNoDataSet )
453 : {
454 : ResetKeyValue( "NODATA",
455 2 : CPLString().Printf( "%.8g", dfNoData ) );
456 : }
457 :
458 68 : if( bCLRDirty )
459 4 : RewriteColorTable( poBand->GetColorTable() );
460 :
461 68 : if( bHDRDirty )
462 66 : RewriteHDR();
463 : }
464 :
465 152 : if( fpImage != NULL )
466 152 : VSIFCloseL( fpImage );
467 :
468 152 : CPLFree( pszProjection );
469 152 : CSLDestroy( papszHDR );
470 152 : }
471 :
472 : /************************************************************************/
473 : /* GetKeyValue() */
474 : /************************************************************************/
475 :
476 132 : const char *EHdrDataset::GetKeyValue( const char *pszKey,
477 : const char *pszDefault )
478 :
479 : {
480 : int i;
481 :
482 1204 : for( i = 0; papszHDR[i] != NULL; i++ )
483 : {
484 1324 : if( EQUALN(pszKey,papszHDR[i],strlen(pszKey))
485 126 : && isspace((unsigned char)papszHDR[i][strlen(pszKey)]) )
486 : {
487 126 : const char *pszValue = papszHDR[i] + strlen(pszKey);
488 984 : while( isspace((unsigned char)*pszValue) )
489 732 : pszValue++;
490 :
491 126 : return pszValue;
492 : }
493 : }
494 :
495 6 : return pszDefault;
496 : }
497 :
498 : /************************************************************************/
499 : /* ResetKeyValue() */
500 : /* */
501 : /* Replace or add the keyword with the indicated value in the */
502 : /* papszHDR list. */
503 : /************************************************************************/
504 :
505 258 : void EHdrDataset::ResetKeyValue( const char *pszKey, const char *pszValue )
506 :
507 : {
508 : int i;
509 : char szNewLine[82];
510 :
511 258 : if( strlen(pszValue) > 65 )
512 : {
513 0 : CPLAssert( strlen(pszValue) <= 65 );
514 0 : return;
515 : }
516 :
517 258 : sprintf( szNewLine, "%-15s%s", pszKey, pszValue );
518 :
519 2964 : for( i = CSLCount(papszHDR)-1; i >= 0; i-- )
520 : {
521 2706 : if( EQUALN(papszHDR[i],szNewLine,strlen(pszKey)+1 ) )
522 : {
523 0 : if( strcmp(papszHDR[i],szNewLine) != 0 )
524 : {
525 0 : CPLFree( papszHDR[i] );
526 0 : papszHDR[i] = CPLStrdup( szNewLine );
527 0 : bHDRDirty = TRUE;
528 : }
529 0 : return;
530 : }
531 : }
532 :
533 258 : bHDRDirty = TRUE;
534 258 : papszHDR = CSLAddString( papszHDR, szNewLine );
535 : }
536 :
537 : /************************************************************************/
538 : /* RewriteColorTable() */
539 : /************************************************************************/
540 :
541 4 : void EHdrDataset::RewriteColorTable( GDALColorTable *poTable )
542 :
543 : {
544 4 : CPLString osCLRFilename = CPLResetExtension( GetDescription(), "clr" );
545 4 : if( poTable )
546 : {
547 4 : VSILFILE *fp = VSIFOpenL( osCLRFilename, "wt" );
548 4 : if( fp != NULL )
549 : {
550 16 : for( int iColor = 0; iColor < poTable->GetColorEntryCount(); iColor++ )
551 : {
552 12 : CPLString oLine;
553 : GDALColorEntry sEntry;
554 :
555 12 : poTable->GetColorEntryAsRGB( iColor, &sEntry );
556 :
557 : // I wish we had a way to mark transparency.
558 : oLine.Printf( "%3d %3d %3d %3d\n",
559 12 : iColor, sEntry.c1, sEntry.c2, sEntry.c3 );
560 12 : VSIFWriteL( (void *) oLine.c_str(), 1, strlen(oLine), fp );
561 : }
562 4 : VSIFCloseL( fp );
563 : }
564 : else
565 : {
566 : CPLError( CE_Failure, CPLE_OpenFailed,
567 : "Unable to create color file %s.",
568 0 : osCLRFilename.c_str() );
569 : }
570 : }
571 : else
572 0 : VSIUnlink( osCLRFilename );
573 4 : }
574 :
575 : /************************************************************************/
576 : /* GetProjectionRef() */
577 : /************************************************************************/
578 :
579 32 : const char *EHdrDataset::GetProjectionRef()
580 :
581 : {
582 32 : if (pszProjection && strlen(pszProjection) > 0)
583 32 : return pszProjection;
584 :
585 0 : return GDALPamDataset::GetProjectionRef();
586 : }
587 :
588 : /************************************************************************/
589 : /* SetProjection() */
590 : /************************************************************************/
591 :
592 60 : CPLErr EHdrDataset::SetProjection( const char *pszSRS )
593 :
594 : {
595 : /* -------------------------------------------------------------------- */
596 : /* Reset coordinate system on the dataset. */
597 : /* -------------------------------------------------------------------- */
598 60 : CPLFree( pszProjection );
599 60 : pszProjection = CPLStrdup( pszSRS );
600 :
601 60 : if( strlen(pszSRS) == 0 )
602 0 : return CE_None;
603 :
604 : /* -------------------------------------------------------------------- */
605 : /* Convert to ESRI WKT. */
606 : /* -------------------------------------------------------------------- */
607 60 : OGRSpatialReference oSRS( pszSRS );
608 60 : char *pszESRI_SRS = NULL;
609 :
610 60 : oSRS.morphToESRI();
611 60 : oSRS.exportToWkt( &pszESRI_SRS );
612 :
613 : /* -------------------------------------------------------------------- */
614 : /* Write to .prj file. */
615 : /* -------------------------------------------------------------------- */
616 60 : CPLString osPrjFilename = CPLResetExtension( GetDescription(), "prj" );
617 : VSILFILE *fp;
618 :
619 60 : fp = VSIFOpenL( osPrjFilename.c_str(), "wt" );
620 60 : if( fp != NULL )
621 : {
622 60 : VSIFWriteL( pszESRI_SRS, 1, strlen(pszESRI_SRS), fp );
623 60 : VSIFWriteL( (void *) "\n", 1, 1, fp );
624 60 : VSIFCloseL( fp );
625 : }
626 :
627 60 : CPLFree( pszESRI_SRS );
628 :
629 60 : return CE_None;
630 : }
631 :
632 : /************************************************************************/
633 : /* GetGeoTransform() */
634 : /************************************************************************/
635 :
636 48 : CPLErr EHdrDataset::GetGeoTransform( double * padfTransform )
637 :
638 : {
639 48 : if( bGotTransform )
640 : {
641 48 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
642 48 : return CE_None;
643 : }
644 : else
645 : {
646 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
647 : }
648 : }
649 :
650 : /************************************************************************/
651 : /* SetGeoTransform() */
652 : /************************************************************************/
653 :
654 64 : CPLErr EHdrDataset::SetGeoTransform( double *padfGeoTransform )
655 :
656 : {
657 : /* -------------------------------------------------------------------- */
658 : /* We only support non-rotated images with info in the .HDR file. */
659 : /* -------------------------------------------------------------------- */
660 128 : if( padfGeoTransform[2] != 0.0
661 64 : || padfGeoTransform[4] != 0.0 )
662 : {
663 0 : return GDALPamDataset::SetGeoTransform( padfGeoTransform );
664 : }
665 :
666 : /* -------------------------------------------------------------------- */
667 : /* Record new geotransform. */
668 : /* -------------------------------------------------------------------- */
669 64 : bGotTransform = TRUE;
670 64 : memcpy( adfGeoTransform, padfGeoTransform, sizeof(double) * 6 );
671 :
672 : /* -------------------------------------------------------------------- */
673 : /* Strip out all old geotransform keywords from HDR records. */
674 : /* -------------------------------------------------------------------- */
675 : int i;
676 640 : for( i = CSLCount(papszHDR)-1; i >= 0; i-- )
677 : {
678 2304 : if( EQUALN(papszHDR[i],"ul",2)
679 576 : || EQUALN(papszHDR[i]+1,"ll",2)
680 576 : || EQUALN(papszHDR[i],"cell",4)
681 576 : || EQUALN(papszHDR[i]+1,"dim",3) )
682 : {
683 0 : papszHDR = CSLRemoveStrings( papszHDR, i, 1, NULL );
684 : }
685 : }
686 :
687 : /* -------------------------------------------------------------------- */
688 : /* Set the transformation information. */
689 : /* -------------------------------------------------------------------- */
690 64 : CPLString oValue;
691 :
692 64 : oValue.Printf( "%.15g", adfGeoTransform[0] + adfGeoTransform[1] * 0.5 );
693 64 : ResetKeyValue( "ULXMAP", oValue );
694 :
695 64 : oValue.Printf( "%.15g", adfGeoTransform[3] + adfGeoTransform[5] * 0.5 );
696 64 : ResetKeyValue( "ULYMAP", oValue );
697 :
698 64 : oValue.Printf( "%.15g", adfGeoTransform[1] );
699 64 : ResetKeyValue( "XDIM", oValue );
700 :
701 64 : oValue.Printf( "%.15g", fabs(adfGeoTransform[5]) );
702 64 : ResetKeyValue( "YDIM", oValue );
703 :
704 64 : return CE_None;
705 : }
706 :
707 : /************************************************************************/
708 : /* RewriteHDR() */
709 : /************************************************************************/
710 :
711 66 : CPLErr EHdrDataset::RewriteHDR()
712 :
713 : {
714 66 : CPLString osPath = CPLGetPath( GetDescription() );
715 66 : CPLString osName = CPLGetBasename( GetDescription() );
716 66 : CPLString osHDRFilename = CPLFormCIFilename( osPath, osName, osHeaderExt );
717 :
718 : /* -------------------------------------------------------------------- */
719 : /* Write .hdr file. */
720 : /* -------------------------------------------------------------------- */
721 : VSILFILE *fp;
722 : int i;
723 :
724 66 : fp = VSIFOpenL( osHDRFilename, "wt" );
725 :
726 66 : if( fp == NULL )
727 : {
728 : CPLError( CE_Failure, CPLE_OpenFailed,
729 : "Failed to rewrite .hdr file %s.",
730 0 : osHDRFilename.c_str() );
731 0 : return CE_Failure;
732 : }
733 :
734 918 : for( i = 0; papszHDR[i] != NULL; i++ )
735 : {
736 852 : VSIFWriteL( papszHDR[i], 1, strlen(papszHDR[i]), fp );
737 852 : VSIFWriteL( (void *) "\n", 1, 1, fp );
738 : }
739 :
740 66 : VSIFCloseL( fp );
741 :
742 66 : bHDRDirty = FALSE;
743 :
744 66 : return CE_None;
745 : }
746 :
747 : /************************************************************************/
748 : /* RewriteSTX() */
749 : /************************************************************************/
750 :
751 4 : CPLErr EHdrDataset::RewriteSTX()
752 : {
753 4 : CPLString osPath = CPLGetPath( GetDescription() );
754 4 : CPLString osName = CPLGetBasename( GetDescription() );
755 4 : CPLString osSTXFilename = CPLFormCIFilename( osPath, osName, "stx" );
756 :
757 : /* -------------------------------------------------------------------- */
758 : /* Write .stx file. */
759 : /* -------------------------------------------------------------------- */
760 : VSILFILE *fp;
761 4 : fp = VSIFOpenL( osSTXFilename, "wt" );
762 4 : if( fp == NULL )
763 : {
764 : CPLDebug( "EHDR", "Failed to rewrite .stx file %s.",
765 0 : osSTXFilename.c_str() );
766 0 : return CE_Failure;
767 : }
768 :
769 8 : for (int i = 0; i < nBands; ++i)
770 : {
771 4 : EHdrRasterBand* poBand = (EHdrRasterBand*)papoBands[i];
772 4 : VSIFPrintfL( fp, "%d %.10f %.10f ", i+1, poBand->dfMin, poBand->dfMax );
773 4 : if ( poBand->minmaxmeanstddev & HAS_MEAN_FLAG )
774 4 : VSIFPrintfL( fp, "%.10f ", poBand->dfMean);
775 : else
776 0 : VSIFPrintfL( fp, "# ");
777 :
778 4 : if ( poBand->minmaxmeanstddev & HAS_STDDEV_FLAG )
779 4 : VSIFPrintfL( fp, "%.10f\n", poBand->dfStdDev);
780 : else
781 0 : VSIFPrintfL( fp, "#\n");
782 : }
783 :
784 4 : VSIFCloseL( fp );
785 :
786 4 : return CE_None;
787 : }
788 :
789 : /************************************************************************/
790 : /* ReadSTX() */
791 : /************************************************************************/
792 :
793 152 : CPLErr EHdrDataset::ReadSTX()
794 : {
795 152 : CPLString osPath = CPLGetPath( GetDescription() );
796 152 : CPLString osName = CPLGetBasename( GetDescription() );
797 152 : CPLString osSTXFilename = CPLFormCIFilename( osPath, osName, "stx" );
798 :
799 : /* -------------------------------------------------------------------- */
800 : /* Read .stx file. */
801 : /* -------------------------------------------------------------------- */
802 : VSILFILE *fp;
803 152 : if ((fp = VSIFOpenL( osSTXFilename, "rt" )) != NULL)
804 : {
805 : const char * pszLine;
806 12 : while( (pszLine = CPLReadLineL( fp )) != NULL )
807 : {
808 : char **papszTokens;
809 4 : papszTokens = CSLTokenizeStringComplex( pszLine, " \t", TRUE, FALSE );
810 4 : int nTokens = CSLCount( papszTokens );
811 4 : if( nTokens >= 5 )
812 : {
813 4 : int i = atoi(papszTokens[0]);
814 4 : if (i > 0 && i <= nBands)
815 : {
816 4 : EHdrRasterBand* poBand = (EHdrRasterBand*)papoBands[i-1];
817 4 : poBand->dfMin = atof(papszTokens[1]);
818 4 : poBand->dfMax = atof(papszTokens[2]);
819 :
820 4 : int bNoDataSet = FALSE;
821 4 : double dfNoData = poBand->GetNoDataValue(&bNoDataSet);
822 4 : if (bNoDataSet && dfNoData == poBand->dfMin)
823 : {
824 : /* Triggered by /vsicurl/http://eros.usgs.gov/archive/nslrsda/GeoTowns/HongKong/srtm/n22e113.zip/n22e113.bil */
825 : CPLDebug("EHDr", "Ignoring .stx file where min == nodata. "
826 : "The nodata value shouldn't be taken into account "
827 0 : "in minimum value computation.");
828 0 : CSLDestroy( papszTokens );
829 0 : papszTokens = NULL;
830 0 : break;
831 : }
832 :
833 4 : poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
834 : // reads optional mean and stddev
835 4 : if ( !EQUAL(papszTokens[3], "#") )
836 : {
837 4 : poBand->dfMean = atof(papszTokens[3]);
838 4 : poBand->minmaxmeanstddev |= HAS_MEAN_FLAG;
839 : }
840 4 : if ( !EQUAL(papszTokens[4], "#") )
841 : {
842 4 : poBand->dfStdDev = atof(papszTokens[4]);
843 4 : poBand->minmaxmeanstddev |= HAS_STDDEV_FLAG;
844 : }
845 :
846 4 : if( nTokens >= 6 && !EQUAL(papszTokens[5], "#") )
847 0 : poBand->SetMetadataItem( "STRETCHMIN", papszTokens[5], "RENDERING_HINTS" );
848 :
849 4 : if( nTokens >= 7 && !EQUAL(papszTokens[6], "#") )
850 0 : poBand->SetMetadataItem( "STRETCHMAX", papszTokens[6], "RENDERING_HINTS" );
851 : }
852 : }
853 :
854 4 : CSLDestroy( papszTokens );
855 : }
856 :
857 4 : VSIFCloseL( fp );
858 : }
859 :
860 152 : return CE_None;
861 : }
862 :
863 :
864 : /************************************************************************/
865 : /* GetImageRepFilename() */
866 : /************************************************************************/
867 :
868 : /* -------------------------------------------------------------------- */
869 : /* Check for IMAGE.REP (Spatiocarte Defense 1.0) or name_of_image.rep */
870 : /* if it's a GIS-GeoSPOT image */
871 : /* For the specification of SPDF (in French), */
872 : /* see http://eden.ign.fr/download/pub/doc/emabgi/spdf10.pdf/download */
873 : /* -------------------------------------------------------------------- */
874 :
875 106 : CPLString EHdrDataset::GetImageRepFilename(const char* pszFilename)
876 : {
877 : VSIStatBufL sStatBuf;
878 :
879 106 : CPLString osPath = CPLGetPath( pszFilename );
880 106 : CPLString osName = CPLGetBasename( pszFilename );
881 : CPLString osREPFilename =
882 106 : CPLFormCIFilename( osPath, osName, "rep" );
883 106 : if( VSIStatExL( osREPFilename.c_str(), &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
884 0 : return osREPFilename;
885 :
886 106 : if (EQUAL(CPLGetFilename(pszFilename), "imspatio.bil") ||
887 : EQUAL(CPLGetFilename(pszFilename), "haspatio.bil"))
888 : {
889 0 : CPLString osImageRepFilename(CPLFormCIFilename( osPath, "image", "rep" ));
890 0 : if( VSIStatExL( osImageRepFilename.c_str(), &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
891 0 : return osImageRepFilename;
892 :
893 : /* Try in the upper directories if not found in the BIL image directory */
894 0 : CPLString dirName(CPLGetDirname(osPath));
895 0 : if (CPLIsFilenameRelative(osPath.c_str()))
896 : {
897 0 : char* cwd = CPLGetCurrentDir();
898 0 : if (cwd)
899 : {
900 0 : dirName = CPLFormFilename(cwd, dirName.c_str(), NULL);
901 0 : CPLFree(cwd);
902 : }
903 : }
904 0 : while (dirName[0] != 0 && EQUAL(dirName, ".") == FALSE && EQUAL(dirName, "/") == FALSE)
905 : {
906 0 : osImageRepFilename = CPLFormCIFilename( dirName.c_str(), "image", "rep" );
907 0 : if( VSIStatExL( osImageRepFilename.c_str(), &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
908 0 : return osImageRepFilename;
909 :
910 : /* Don't try to recurse above the 'image' subdirectory */
911 0 : if (EQUAL(dirName, "image"))
912 : {
913 0 : break;
914 : }
915 0 : dirName = CPLString(CPLGetDirname(dirName));
916 0 : }
917 : }
918 106 : return CPLString();
919 : }
920 :
921 : /************************************************************************/
922 : /* GetFileList() */
923 : /************************************************************************/
924 :
925 18 : char **EHdrDataset::GetFileList()
926 :
927 : {
928 : VSIStatBufL sStatBuf;
929 18 : CPLString osPath = CPLGetPath( GetDescription() );
930 18 : CPLString osName = CPLGetBasename( GetDescription() );
931 18 : char **papszFileList = NULL;
932 :
933 : // Main data file, etc.
934 18 : papszFileList = GDALPamDataset::GetFileList();
935 :
936 : // Header file.
937 18 : CPLString osFilename = CPLFormCIFilename( osPath, osName, osHeaderExt );
938 18 : papszFileList = CSLAddString( papszFileList, osFilename );
939 :
940 : // Statistics file
941 18 : osFilename = CPLFormCIFilename( osPath, osName, "stx" );
942 18 : if( VSIStatExL( osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
943 2 : papszFileList = CSLAddString( papszFileList, osFilename );
944 :
945 : // color table file.
946 18 : osFilename = CPLFormCIFilename( osPath, osName, "clr" );
947 18 : if( VSIStatExL( osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
948 4 : papszFileList = CSLAddString( papszFileList, osFilename );
949 :
950 : // projections file.
951 18 : osFilename = CPLFormCIFilename( osPath, osName, "prj" );
952 18 : if( VSIStatExL( osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
953 10 : papszFileList = CSLAddString( papszFileList, osFilename );
954 :
955 18 : CPLString imageRepFilename = GetImageRepFilename( GetDescription() );
956 18 : if (!imageRepFilename.empty())
957 0 : papszFileList = CSLAddString( papszFileList, imageRepFilename.c_str() );
958 :
959 18 : return papszFileList;
960 : }
961 :
962 : /************************************************************************/
963 : /* Open() */
964 : /************************************************************************/
965 :
966 23030 : GDALDataset *EHdrDataset::Open( GDALOpenInfo * poOpenInfo )
967 :
968 : {
969 : int i, bSelectedHDR;
970 :
971 : /* -------------------------------------------------------------------- */
972 : /* We assume the user is pointing to the binary (ie. .bil) file. */
973 : /* -------------------------------------------------------------------- */
974 23030 : if( poOpenInfo->nHeaderBytes < 2 )
975 21302 : return NULL;
976 :
977 : /* -------------------------------------------------------------------- */
978 : /* Now we need to tear apart the filename to form a .HDR */
979 : /* filename. */
980 : /* -------------------------------------------------------------------- */
981 1728 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
982 1728 : CPLString osName = CPLGetBasename( poOpenInfo->pszFilename );
983 1728 : CPLString osHDRFilename;
984 :
985 1728 : const char* pszHeaderExt = "hdr";
986 1728 : if( EQUAL( CPLGetExtension( poOpenInfo->pszFilename ), "SRC" ) &&
987 : osName.size() == 7 &&
988 : (osName[0] == 'e' || osName[0] == 'E' || osName[0] == 'w' || osName[0] == 'W') &&
989 : (osName[4] == 'n' || osName[4] == 'N' || osName[4] == 's' || osName[4] == 'S') )
990 : {
991 : /* It is a GTOPO30 or SRTM30 source file, whose header extension is .sch */
992 : /* see http://dds.cr.usgs.gov/srtm/version1/SRTM30/GTOPO30_Documentation */
993 0 : pszHeaderExt = "sch";
994 : }
995 :
996 1728 : if( poOpenInfo->papszSiblingFiles )
997 : {
998 : int iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
999 1720 : CPLFormFilename( NULL, osName, pszHeaderExt ) );
1000 1720 : if( iFile < 0 ) // return if there is no corresponding .hdr file
1001 1506 : return NULL;
1002 :
1003 : osHDRFilename =
1004 214 : CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1005 428 : NULL );
1006 : }
1007 : else
1008 : {
1009 8 : osHDRFilename = CPLFormCIFilename( osPath, osName, pszHeaderExt );
1010 : }
1011 :
1012 222 : bSelectedHDR = EQUAL( osHDRFilename, poOpenInfo->pszFilename );
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Do we have a .hdr file? */
1016 : /* -------------------------------------------------------------------- */
1017 : VSILFILE *fp;
1018 :
1019 222 : fp = VSIFOpenL( osHDRFilename, "r" );
1020 :
1021 222 : if( fp == NULL )
1022 : {
1023 8 : return NULL;
1024 : }
1025 :
1026 : /* -------------------------------------------------------------------- */
1027 : /* Is this file an ESRI header file? Read a few lines of text */
1028 : /* searching for something starting with nrows or ncols. */
1029 : /* -------------------------------------------------------------------- */
1030 : const char * pszLine;
1031 214 : int nRows = -1, nCols = -1, nBands = 1;
1032 214 : int nSkipBytes = 0;
1033 214 : double dfULXMap=0.5, dfULYMap = 0.5, dfYLLCorner = -123.456;
1034 214 : int bCenter = TRUE;
1035 214 : double dfXDim = 1.0, dfYDim = 1.0, dfNoData = 0.0;
1036 214 : int nLineCount = 0, bNoDataSet = FALSE;
1037 214 : GDALDataType eDataType = GDT_Byte;
1038 214 : int nBits = -1;
1039 214 : char chByteOrder = 'M';
1040 214 : char chPixelType = 'N'; // not defined
1041 214 : char szLayout[10] = "BIL";
1042 214 : char **papszHDR = NULL;
1043 214 : int bHasInternalProjection = FALSE;
1044 214 : int bHasMin = FALSE;
1045 214 : int bHasMax = FALSE;
1046 214 : double dfMin = 0, dfMax = 0;
1047 :
1048 2750 : while( (pszLine = CPLReadLineL( fp )) != NULL )
1049 : {
1050 : char **papszTokens;
1051 :
1052 2326 : nLineCount++;
1053 :
1054 2326 : if( nLineCount > 50 || strlen(pszLine) > 1000 )
1055 4 : break;
1056 :
1057 2322 : papszHDR = CSLAddString( papszHDR, pszLine );
1058 :
1059 2322 : papszTokens = CSLTokenizeStringComplex( pszLine, " \t", TRUE, FALSE );
1060 2322 : if( CSLCount( papszTokens ) < 2 )
1061 : {
1062 98 : CSLDestroy( papszTokens );
1063 98 : continue;
1064 : }
1065 :
1066 2224 : if( EQUAL(papszTokens[0],"ncols") )
1067 : {
1068 152 : nCols = atoi(papszTokens[1]);
1069 : }
1070 2072 : else if( EQUAL(papszTokens[0],"nrows") )
1071 : {
1072 152 : nRows = atoi(papszTokens[1]);
1073 : }
1074 1920 : else if( EQUAL(papszTokens[0],"skipbytes") )
1075 : {
1076 0 : nSkipBytes = atoi(papszTokens[1]);
1077 : }
1078 5688 : else if( EQUAL(papszTokens[0],"ulxmap")
1079 1848 : || EQUAL(papszTokens[0],"xllcorner")
1080 1844 : || EQUAL(papszTokens[0],"xllcenter") )
1081 : {
1082 76 : dfULXMap = atof(papszTokens[1]);
1083 76 : if( EQUAL(papszTokens[0],"xllcorner") )
1084 4 : bCenter = FALSE;
1085 : }
1086 1844 : else if( EQUAL(papszTokens[0],"ulymap") )
1087 : {
1088 72 : dfULYMap = atof(papszTokens[1]);
1089 : }
1090 3544 : else if( EQUAL(papszTokens[0],"yllcorner")
1091 1768 : || EQUAL(papszTokens[0],"yllcenter") )
1092 : {
1093 4 : dfYLLCorner = atof(papszTokens[1]);
1094 4 : if( EQUAL(papszTokens[0],"yllcorner") )
1095 4 : bCenter = FALSE;
1096 : }
1097 1768 : else if( EQUAL(papszTokens[0],"xdim") )
1098 : {
1099 72 : dfXDim = atof(papszTokens[1]);
1100 : }
1101 1696 : else if( EQUAL(papszTokens[0],"ydim") )
1102 : {
1103 72 : dfYDim = atof(papszTokens[1]);
1104 : }
1105 1624 : else if( EQUAL(papszTokens[0],"cellsize") )
1106 : {
1107 4 : dfXDim = dfYDim = atof(papszTokens[1]);
1108 : }
1109 1620 : else if( EQUAL(papszTokens[0],"nbands") )
1110 : {
1111 148 : nBands = atoi(papszTokens[1]);
1112 : }
1113 1472 : else if( EQUAL(papszTokens[0],"layout") )
1114 : {
1115 148 : strncpy( szLayout, papszTokens[1], sizeof(szLayout) );
1116 148 : szLayout[sizeof(szLayout)-1] = '\0';
1117 : }
1118 2656 : else if( EQUAL(papszTokens[0],"NODATA_value")
1119 1324 : || EQUAL(papszTokens[0],"NODATA") )
1120 : {
1121 8 : dfNoData = atof(papszTokens[1]);
1122 8 : bNoDataSet = TRUE;
1123 : }
1124 1316 : else if( EQUAL(papszTokens[0],"NBITS") )
1125 : {
1126 150 : nBits = atoi(papszTokens[1]);
1127 : }
1128 1166 : else if( EQUAL(papszTokens[0],"PIXELTYPE") )
1129 : {
1130 144 : chPixelType = (char) toupper(papszTokens[1][0]);
1131 : }
1132 1022 : else if( EQUAL(papszTokens[0],"byteorder") )
1133 : {
1134 152 : chByteOrder = (char) toupper(papszTokens[1][0]);
1135 : }
1136 :
1137 : /* http://www.worldclim.org/futdown.htm have the projection extensions */
1138 870 : else if( EQUAL(papszTokens[0],"Projection") )
1139 : {
1140 6 : bHasInternalProjection = TRUE;
1141 : }
1142 1728 : else if( EQUAL(papszTokens[0],"MinValue") ||
1143 862 : EQUAL(papszTokens[0],"MIN_VALUE") )
1144 : {
1145 2 : dfMin = atof(papszTokens[1]);
1146 2 : bHasMin = TRUE;
1147 : }
1148 1722 : else if( EQUAL(papszTokens[0],"MaxValue") ||
1149 860 : EQUAL(papszTokens[0],"MAX_VALUE") )
1150 : {
1151 2 : dfMax = atof(papszTokens[1]);
1152 2 : bHasMax = TRUE;
1153 : }
1154 :
1155 2224 : CSLDestroy( papszTokens );
1156 : }
1157 :
1158 214 : VSIFCloseL( fp );
1159 :
1160 : /* -------------------------------------------------------------------- */
1161 : /* Did we get the required keywords? If not we return with */
1162 : /* this never having been considered to be a match. This isn't */
1163 : /* an error! */
1164 : /* -------------------------------------------------------------------- */
1165 214 : if( nRows == -1 || nCols == -1 )
1166 : {
1167 62 : CSLDestroy( papszHDR );
1168 62 : return NULL;
1169 : }
1170 :
1171 152 : if (!GDALCheckDatasetDimensions(nCols, nRows) ||
1172 : !GDALCheckBandCount(nBands, FALSE))
1173 : {
1174 0 : CSLDestroy( papszHDR );
1175 0 : return NULL;
1176 : }
1177 :
1178 : /* -------------------------------------------------------------------- */
1179 : /* Has the user selected the .hdr file to open? */
1180 : /* -------------------------------------------------------------------- */
1181 152 : if( bSelectedHDR )
1182 : {
1183 : CPLError( CE_Failure, CPLE_AppDefined,
1184 : "The selected file is an ESRI BIL header file, but to\n"
1185 : "open ESRI BIL datasets, the data file should be selected\n"
1186 : "instead of the .hdr file. Please try again selecting\n"
1187 : "the data file (often with the extension .bil) corresponding\n"
1188 : "to the header file: %s\n",
1189 0 : poOpenInfo->pszFilename );
1190 0 : CSLDestroy( papszHDR );
1191 0 : return NULL;
1192 : }
1193 :
1194 : /* -------------------------------------------------------------------- */
1195 : /* If we aren't sure of the file type, check the data file */
1196 : /* size. If it is 4 bytes or more per pixel then we assume it */
1197 : /* is floating point data. */
1198 : /* -------------------------------------------------------------------- */
1199 152 : if( nBits == -1 && chPixelType == 'N' )
1200 : {
1201 : VSIStatBufL sStatBuf;
1202 2 : if( VSIStatL( poOpenInfo->pszFilename, &sStatBuf ) == 0 )
1203 : {
1204 2 : size_t nBytes = (size_t) (sStatBuf.st_size/nCols/nRows/nBands);
1205 2 : if( nBytes > 0 && nBytes != 3 )
1206 2 : nBits = nBytes*8;
1207 :
1208 2 : if( nBytes == 4 )
1209 2 : chPixelType = 'F';
1210 : }
1211 : }
1212 :
1213 : /* -------------------------------------------------------------------- */
1214 : /* If the extension is FLT it is likely a floating point file. */
1215 : /* -------------------------------------------------------------------- */
1216 152 : if( chPixelType == 'N' )
1217 : {
1218 6 : if( EQUAL( CPLGetExtension( poOpenInfo->pszFilename ), "FLT" ) )
1219 2 : chPixelType = 'F';
1220 : }
1221 :
1222 : /* -------------------------------------------------------------------- */
1223 : /* If we have a negative nodata value, let's assume that the */
1224 : /* pixel type is signed. This is necessary for datasets from */
1225 : /* http://www.worldclim.org/futdown.htm */
1226 : /* -------------------------------------------------------------------- */
1227 152 : if( bNoDataSet && dfNoData < 0 && chPixelType == 'N' )
1228 : {
1229 4 : chPixelType = 'S';
1230 : }
1231 :
1232 : /* -------------------------------------------------------------------- */
1233 : /* Create a corresponding GDALDataset. */
1234 : /* -------------------------------------------------------------------- */
1235 : EHdrDataset *poDS;
1236 :
1237 152 : poDS = new EHdrDataset();
1238 :
1239 304 : poDS->osHeaderExt = pszHeaderExt;
1240 :
1241 : /* -------------------------------------------------------------------- */
1242 : /* Capture some information from the file that is of interest. */
1243 : /* -------------------------------------------------------------------- */
1244 152 : poDS->nRasterXSize = nCols;
1245 152 : poDS->nRasterYSize = nRows;
1246 152 : poDS->papszHDR = papszHDR;
1247 :
1248 : /* -------------------------------------------------------------------- */
1249 : /* Open target binary file. */
1250 : /* -------------------------------------------------------------------- */
1251 152 : if( poOpenInfo->eAccess == GA_ReadOnly )
1252 84 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
1253 : else
1254 68 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
1255 :
1256 152 : if( poDS->fpImage == NULL )
1257 : {
1258 : CPLError( CE_Failure, CPLE_OpenFailed,
1259 : "Failed to open %s with write permission.\n%s",
1260 0 : osName.c_str(), VSIStrerror( errno ) );
1261 0 : delete poDS;
1262 0 : return NULL;
1263 : }
1264 :
1265 152 : poDS->eAccess = poOpenInfo->eAccess;
1266 :
1267 : /* -------------------------------------------------------------------- */
1268 : /* Figure out the data type. */
1269 : /* -------------------------------------------------------------------- */
1270 152 : if( nBits == 16 )
1271 : {
1272 30 : if ( chPixelType == 'S' )
1273 14 : eDataType = GDT_Int16;
1274 : else
1275 16 : eDataType = GDT_UInt16; // default
1276 : }
1277 122 : else if( nBits == 32 )
1278 : {
1279 62 : if( chPixelType == 'S' )
1280 16 : eDataType = GDT_Int32;
1281 46 : else if( chPixelType == 'F' )
1282 36 : eDataType = GDT_Float32;
1283 : else
1284 10 : eDataType = GDT_UInt32; // default
1285 : }
1286 60 : else if( nBits == 8 )
1287 : {
1288 54 : eDataType = GDT_Byte;
1289 54 : nBits = 8;
1290 : }
1291 12 : else if( nBits < 8 && nBits >= 1 )
1292 6 : eDataType = GDT_Byte;
1293 0 : else if( nBits == -1 )
1294 : {
1295 0 : if( chPixelType == 'F' )
1296 : {
1297 0 : eDataType = GDT_Float32;
1298 0 : nBits = 32;
1299 : }
1300 : else
1301 : {
1302 0 : eDataType = GDT_Byte;
1303 0 : nBits = 8;
1304 : }
1305 : }
1306 : else
1307 : {
1308 : CPLError( CE_Failure, CPLE_NotSupported,
1309 : "EHdr driver does not support %d NBITS value.",
1310 0 : nBits );
1311 0 : delete poDS;
1312 0 : return NULL;
1313 : }
1314 :
1315 : /* -------------------------------------------------------------------- */
1316 : /* Compute the line offset. */
1317 : /* -------------------------------------------------------------------- */
1318 152 : int nItemSize = GDALGetDataTypeSize(eDataType)/8;
1319 : int nPixelOffset, nLineOffset;
1320 : vsi_l_offset nBandOffset;
1321 :
1322 152 : if( EQUAL(szLayout,"BIP") )
1323 : {
1324 0 : nPixelOffset = nItemSize * nBands;
1325 0 : nLineOffset = nPixelOffset * nCols;
1326 0 : nBandOffset = (vsi_l_offset)nItemSize;
1327 : }
1328 152 : else if( EQUAL(szLayout,"BSQ") )
1329 : {
1330 0 : nPixelOffset = nItemSize;
1331 0 : nLineOffset = nPixelOffset * nCols;
1332 0 : nBandOffset = (vsi_l_offset)nLineOffset * nRows;
1333 : }
1334 : else /* assume BIL */
1335 : {
1336 152 : nPixelOffset = nItemSize;
1337 152 : nLineOffset = nItemSize * nBands * nCols;
1338 152 : nBandOffset = (vsi_l_offset)nItemSize * nCols;
1339 : }
1340 :
1341 152 : poDS->SetDescription( poOpenInfo->pszFilename );
1342 152 : poDS->PamInitialize();
1343 :
1344 : /* -------------------------------------------------------------------- */
1345 : /* Create band information objects. */
1346 : /* -------------------------------------------------------------------- */
1347 152 : poDS->nBands = nBands;
1348 404 : for( i = 0; i < poDS->nBands; i++ )
1349 : {
1350 : EHdrRasterBand *poBand;
1351 :
1352 : poBand =
1353 : new EHdrRasterBand( poDS, i+1, poDS->fpImage,
1354 : nSkipBytes + nBandOffset * i,
1355 : nPixelOffset, nLineOffset, eDataType,
1356 : #ifdef CPL_LSB
1357 : chByteOrder == 'I' || chByteOrder == 'L',
1358 : #else
1359 : chByteOrder == 'M',
1360 : #endif
1361 252 : nBits);
1362 :
1363 252 : if( bNoDataSet )
1364 8 : poBand->SetNoDataValue( dfNoData );
1365 :
1366 252 : if( bHasMin && bHasMax )
1367 : {
1368 2 : poBand->dfMin = dfMin;
1369 2 : poBand->dfMax = dfMax;
1370 2 : poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
1371 : }
1372 :
1373 252 : poDS->SetBand( i+1, poBand );
1374 : }
1375 :
1376 : /* -------------------------------------------------------------------- */
1377 : /* If we didn't get bounds in the .hdr, look for a worldfile. */
1378 : /* -------------------------------------------------------------------- */
1379 152 : if( dfYLLCorner != -123.456 )
1380 : {
1381 4 : if( bCenter )
1382 0 : dfULYMap = dfYLLCorner + (nRows-1) * dfYDim;
1383 : else
1384 4 : dfULYMap = dfYLLCorner + nRows * dfYDim;
1385 : }
1386 :
1387 152 : if( dfULYMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0 )
1388 : {
1389 76 : poDS->bGotTransform = TRUE;
1390 :
1391 76 : if( bCenter )
1392 : {
1393 72 : poDS->adfGeoTransform[0] = dfULXMap - dfXDim * 0.5;
1394 72 : poDS->adfGeoTransform[1] = dfXDim;
1395 72 : poDS->adfGeoTransform[2] = 0.0;
1396 72 : poDS->adfGeoTransform[3] = dfULYMap + dfYDim * 0.5;
1397 72 : poDS->adfGeoTransform[4] = 0.0;
1398 72 : poDS->adfGeoTransform[5] = - dfYDim;
1399 : }
1400 : else
1401 : {
1402 4 : poDS->adfGeoTransform[0] = dfULXMap;
1403 4 : poDS->adfGeoTransform[1] = dfXDim;
1404 4 : poDS->adfGeoTransform[2] = 0.0;
1405 4 : poDS->adfGeoTransform[3] = dfULYMap;
1406 4 : poDS->adfGeoTransform[4] = 0.0;
1407 4 : poDS->adfGeoTransform[5] = - dfYDim;
1408 : }
1409 : }
1410 :
1411 152 : if( !poDS->bGotTransform )
1412 : poDS->bGotTransform =
1413 : GDALReadWorldFile( poOpenInfo->pszFilename, 0,
1414 76 : poDS->adfGeoTransform );
1415 :
1416 152 : if( !poDS->bGotTransform )
1417 : poDS->bGotTransform =
1418 : GDALReadWorldFile( poOpenInfo->pszFilename, "wld",
1419 76 : poDS->adfGeoTransform );
1420 :
1421 : /* -------------------------------------------------------------------- */
1422 : /* Check for a .prj file. */
1423 : /* -------------------------------------------------------------------- */
1424 152 : const char *pszPrjFilename = CPLFormCIFilename( osPath, osName, "prj" );
1425 :
1426 152 : fp = VSIFOpenL( pszPrjFilename, "r" );
1427 :
1428 : /* .hdr files from http://www.worldclim.org/futdown.htm have the projection */
1429 : /* info in the .hdr file itself ! */
1430 152 : if (fp == NULL && bHasInternalProjection)
1431 : {
1432 2 : pszPrjFilename = osHDRFilename;
1433 2 : fp = VSIFOpenL( pszPrjFilename, "r" );
1434 : }
1435 :
1436 152 : if( fp != NULL )
1437 : {
1438 : char **papszLines;
1439 64 : OGRSpatialReference oSRS;
1440 :
1441 64 : VSIFCloseL( fp );
1442 :
1443 64 : papszLines = CSLLoad( pszPrjFilename );
1444 :
1445 64 : if( oSRS.importFromESRI( papszLines ) == OGRERR_NONE )
1446 : {
1447 : // If geographic values are in seconds, we must transform.
1448 : // Is there a code for minutes too?
1449 : char szResult[80];
1450 64 : if( oSRS.IsGeographic()
1451 : && EQUAL(OSR_GDS( szResult, sizeof(szResult),
1452 : papszLines, "Units", ""), "DS") )
1453 : {
1454 0 : poDS->adfGeoTransform[0] /= 3600.0;
1455 0 : poDS->adfGeoTransform[1] /= 3600.0;
1456 0 : poDS->adfGeoTransform[2] /= 3600.0;
1457 0 : poDS->adfGeoTransform[3] /= 3600.0;
1458 0 : poDS->adfGeoTransform[4] /= 3600.0;
1459 0 : poDS->adfGeoTransform[5] /= 3600.0;
1460 : }
1461 :
1462 64 : CPLFree( poDS->pszProjection );
1463 64 : oSRS.exportToWkt( &(poDS->pszProjection) );
1464 : }
1465 :
1466 64 : CSLDestroy( papszLines );
1467 : }
1468 : else
1469 : {
1470 : /* -------------------------------------------------------------------- */
1471 : /* Check for IMAGE.REP (Spatiocarte Defense 1.0) or name_of_image.rep */
1472 : /* if it's a GIS-GeoSPOT image */
1473 : /* For the specification of SPDF (in French), */
1474 : /* see http://eden.ign.fr/download/pub/doc/emabgi/spdf10.pdf/download */
1475 : /* -------------------------------------------------------------------- */
1476 88 : CPLString szImageRepFilename = GetImageRepFilename(poOpenInfo->pszFilename );
1477 88 : if (!szImageRepFilename.empty())
1478 : {
1479 0 : fp = VSIFOpenL( szImageRepFilename.c_str(), "r" );
1480 : }
1481 88 : if (fp != NULL)
1482 : {
1483 : const char *pszLine;
1484 0 : int bUTM = FALSE;
1485 0 : int bWGS84 = FALSE;
1486 0 : int bNorth = FALSE;
1487 0 : int bSouth = FALSE;
1488 0 : int utmZone = 0;
1489 :
1490 0 : while( (pszLine = CPLReadLineL( fp )) != NULL )
1491 : {
1492 0 : if (strncmp(pszLine, "PROJ_ID", strlen("PROJ_ID")) == 0 &&
1493 : strstr(pszLine, "UTM"))
1494 : {
1495 0 : bUTM = TRUE;
1496 : }
1497 0 : else if (strncmp(pszLine, "PROJ_ZONE", strlen("PROJ_ZONE")) == 0)
1498 : {
1499 0 : const char* c = strchr(pszLine, '"');
1500 0 : if (c)
1501 : {
1502 0 : c++;
1503 0 : if (*c >= '0' && *c <= '9')
1504 : {
1505 0 : utmZone = atoi(c);
1506 0 : if (utmZone >= 1 && utmZone <= 60)
1507 : {
1508 0 : if (strstr(pszLine, "Nord") || strstr(pszLine, "NORD"))
1509 : {
1510 0 : bNorth = TRUE;
1511 : }
1512 0 : else if (strstr(pszLine, "Sud") || strstr(pszLine, "SUD"))
1513 : {
1514 0 : bSouth = TRUE;
1515 : }
1516 : }
1517 : }
1518 : }
1519 : }
1520 0 : else if (strncmp(pszLine, "PROJ_CODE", strlen("PROJ_CODE")) == 0 &&
1521 : strstr(pszLine, "FR-MINDEF"))
1522 : {
1523 0 : const char* c = strchr(pszLine, 'A');
1524 0 : if (c)
1525 : {
1526 0 : c++;
1527 0 : if (*c >= '0' && *c <= '9')
1528 : {
1529 0 : utmZone = atoi(c);
1530 0 : if (utmZone >= 1 && utmZone <= 60)
1531 : {
1532 0 : if (c[1] == 'N' ||
1533 0 : (c[1] != '\0' && c[2] == 'N'))
1534 : {
1535 0 : bNorth = TRUE;
1536 : }
1537 0 : else if (c[1] == 'S' ||
1538 0 : (c[1] != '\0' && c[2] == 'S'))
1539 : {
1540 0 : bSouth = TRUE;
1541 : }
1542 : }
1543 : }
1544 : }
1545 : }
1546 0 : else if (strncmp(pszLine, "HORIZ_DATUM", strlen("HORIZ_DATUM")) == 0 &&
1547 : (strstr(pszLine, "WGS 84") || strstr(pszLine, "WGS84")))
1548 : {
1549 0 : bWGS84 = TRUE;
1550 : }
1551 0 : else if (strncmp(pszLine, "MAP_NUMBER", strlen("MAP_NUMBER")) == 0)
1552 : {
1553 0 : const char* c = strchr(pszLine, '"');
1554 0 : if (c)
1555 : {
1556 0 : char* pszMapNumber = CPLStrdup(c+1);
1557 0 : char* c2 = strchr(pszMapNumber, '"');
1558 0 : if (c2) *c2 = 0;
1559 0 : poDS->SetMetadataItem("SPDF_MAP_NUMBER", pszMapNumber);
1560 0 : CPLFree(pszMapNumber);
1561 : }
1562 : }
1563 0 : else if (strncmp(pszLine, "PRODUCTION_DATE", strlen("PRODUCTION_DATE")) == 0)
1564 : {
1565 0 : const char* c = pszLine + strlen("PRODUCTION_DATE");
1566 0 : while(*c == ' ')
1567 0 : c++;
1568 0 : if (*c)
1569 : {
1570 0 : poDS->SetMetadataItem("SPDF_PRODUCTION_DATE", c );
1571 : }
1572 : }
1573 : }
1574 :
1575 0 : VSIFCloseL( fp );
1576 :
1577 0 : if (utmZone != 0 && bUTM && bWGS84 && (bNorth || bSouth))
1578 : {
1579 : char projCSStr[64];
1580 0 : OGRSpatialReference oSRS;
1581 :
1582 : sprintf(projCSStr, "WGS 84 / UTM zone %d%c",
1583 0 : utmZone, (bNorth) ? 'N' : 'S');
1584 0 : oSRS.SetProjCS(projCSStr);
1585 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1586 0 : oSRS.SetUTM(utmZone, bNorth);
1587 0 : oSRS.SetAuthority("PROJCS", "EPSG", ((bNorth) ? 32600 : 32700) + utmZone);
1588 0 : oSRS.AutoIdentifyEPSG();
1589 :
1590 0 : CPLFree( poDS->pszProjection );
1591 0 : oSRS.exportToWkt( &(poDS->pszProjection) );
1592 : }
1593 : else
1594 : {
1595 0 : CPLError( CE_Warning, CPLE_NotSupported, "Cannot retrive projection from IMAGE.REP");
1596 : }
1597 88 : }
1598 : }
1599 :
1600 : /* -------------------------------------------------------------------- */
1601 : /* Check for a color table. */
1602 : /* -------------------------------------------------------------------- */
1603 152 : const char *pszCLRFilename = CPLFormCIFilename( osPath, osName, "clr" );
1604 :
1605 : /* Only read the .clr for byte, int16 or uint16 bands */
1606 152 : if (nItemSize <= 2)
1607 90 : fp = VSIFOpenL( pszCLRFilename, "r" );
1608 : else
1609 62 : fp = NULL;
1610 :
1611 152 : if( fp != NULL )
1612 : {
1613 8 : GDALColorTable oColorTable;
1614 8 : int bHasWarned = FALSE;
1615 :
1616 24 : while(TRUE)
1617 : {
1618 32 : const char *pszLine = CPLReadLineL(fp);
1619 32 : if ( !pszLine )
1620 : break;
1621 :
1622 24 : if( *pszLine == '#' || *pszLine == '!' )
1623 0 : continue;
1624 :
1625 : char **papszValues = CSLTokenizeString2(pszLine, "\t ",
1626 24 : CSLT_HONOURSTRINGS);
1627 : GDALColorEntry oEntry;
1628 :
1629 24 : if ( CSLCount(papszValues) >= 4 )
1630 : {
1631 24 : int nIndex = atoi( papszValues[0] ); // Index
1632 48 : if (nIndex >= 0 && nIndex < 65536)
1633 : {
1634 24 : oEntry.c1 = (short) atoi( papszValues[1] ); // Red
1635 24 : oEntry.c2 = (short) atoi( papszValues[2] ); // Green
1636 24 : oEntry.c3 = (short) atoi( papszValues[3] ); // Blue
1637 24 : oEntry.c4 = 255;
1638 :
1639 24 : oColorTable.SetColorEntry( nIndex, &oEntry );
1640 : }
1641 : else
1642 : {
1643 : /* Negative values are valid. At least we can find use of */
1644 : /* them here : http://www.ngdc.noaa.gov/mgg/topo/elev/esri/clr/ */
1645 : /* but there's no way of representing them with GDAL color */
1646 : /* table model */
1647 0 : if (!bHasWarned)
1648 0 : CPLDebug("EHdr", "Ignoring color index : %d", nIndex);
1649 0 : bHasWarned = TRUE;
1650 : }
1651 : }
1652 :
1653 24 : CSLDestroy( papszValues );
1654 : }
1655 :
1656 8 : VSIFCloseL( fp );
1657 :
1658 16 : for( i = 1; i <= poDS->nBands; i++ )
1659 : {
1660 8 : GDALRasterBand *poBand = poDS->GetRasterBand( i );
1661 8 : poBand->SetColorTable( &oColorTable );
1662 8 : poBand->SetColorInterpretation( GCI_PaletteIndex );
1663 : }
1664 :
1665 8 : poDS->bCLRDirty = FALSE;
1666 : }
1667 :
1668 : /* -------------------------------------------------------------------- */
1669 : /* Read statistics (.STX) */
1670 : /* -------------------------------------------------------------------- */
1671 152 : poDS->ReadSTX();
1672 :
1673 : /* -------------------------------------------------------------------- */
1674 : /* Initialize any PAM information. */
1675 : /* -------------------------------------------------------------------- */
1676 152 : poDS->TryLoadXML();
1677 :
1678 : /* -------------------------------------------------------------------- */
1679 : /* Check for overviews. */
1680 : /* -------------------------------------------------------------------- */
1681 152 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1682 :
1683 152 : return( poDS );
1684 : }
1685 :
1686 : /************************************************************************/
1687 : /* Create() */
1688 : /************************************************************************/
1689 :
1690 104 : GDALDataset *EHdrDataset::Create( const char * pszFilename,
1691 : int nXSize, int nYSize, int nBands,
1692 : GDALDataType eType,
1693 : char **papszParmList )
1694 :
1695 : {
1696 : /* -------------------------------------------------------------------- */
1697 : /* Verify input options. */
1698 : /* -------------------------------------------------------------------- */
1699 104 : if (nBands <= 0)
1700 : {
1701 : CPLError( CE_Failure, CPLE_NotSupported,
1702 2 : "EHdr driver does not support %d bands.\n", nBands);
1703 2 : return NULL;
1704 : }
1705 :
1706 102 : if( eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16
1707 : && eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_UInt32 )
1708 : {
1709 : CPLError( CE_Failure, CPLE_AppDefined,
1710 : "Attempt to create ESRI .hdr labelled dataset with an illegal\n"
1711 : "data type (%s).\n",
1712 30 : GDALGetDataTypeName(eType) );
1713 :
1714 30 : return NULL;
1715 : }
1716 :
1717 : /* -------------------------------------------------------------------- */
1718 : /* Try to create the file. */
1719 : /* -------------------------------------------------------------------- */
1720 : VSILFILE *fp;
1721 :
1722 72 : fp = VSIFOpenL( pszFilename, "wb" );
1723 :
1724 72 : if( fp == NULL )
1725 : {
1726 : CPLError( CE_Failure, CPLE_OpenFailed,
1727 : "Attempt to create file `%s' failed.\n",
1728 4 : pszFilename );
1729 4 : return NULL;
1730 : }
1731 :
1732 : /* -------------------------------------------------------------------- */
1733 : /* Just write out a couple of bytes to establish the binary */
1734 : /* file, and then close it. */
1735 : /* -------------------------------------------------------------------- */
1736 68 : VSIFWriteL( (void *) "\0\0", 2, 1, fp );
1737 68 : VSIFCloseL( fp );
1738 :
1739 : /* -------------------------------------------------------------------- */
1740 : /* Create the hdr filename. */
1741 : /* -------------------------------------------------------------------- */
1742 : char *pszHdrFilename;
1743 :
1744 : pszHdrFilename =
1745 68 : CPLStrdup( CPLResetExtension( pszFilename, "hdr" ) );
1746 :
1747 : /* -------------------------------------------------------------------- */
1748 : /* Open the file. */
1749 : /* -------------------------------------------------------------------- */
1750 68 : fp = VSIFOpenL( pszHdrFilename, "wt" );
1751 68 : if( fp == NULL )
1752 : {
1753 : CPLError( CE_Failure, CPLE_OpenFailed,
1754 : "Attempt to create file `%s' failed.\n",
1755 0 : pszHdrFilename );
1756 0 : CPLFree( pszHdrFilename );
1757 0 : return NULL;
1758 : }
1759 :
1760 : /* -------------------------------------------------------------------- */
1761 : /* Decide how many bits the file should have. */
1762 : /* -------------------------------------------------------------------- */
1763 : int nRowBytes;
1764 68 : int nBits = GDALGetDataTypeSize(eType);
1765 :
1766 68 : if( CSLFetchNameValue( papszParmList, "NBITS" ) != NULL )
1767 2 : nBits = atoi(CSLFetchNameValue( papszParmList, "NBITS" ));
1768 :
1769 68 : nRowBytes = (nBits * nXSize + 7) / 8;
1770 :
1771 : /* -------------------------------------------------------------------- */
1772 : /* Check for signed byte. */
1773 : /* -------------------------------------------------------------------- */
1774 68 : const char *pszPixelType = CSLFetchNameValue( papszParmList, "PIXELTYPE" );
1775 68 : if( pszPixelType == NULL )
1776 66 : pszPixelType = "";
1777 :
1778 : /* -------------------------------------------------------------------- */
1779 : /* Write out the raw definition for the dataset as a whole. */
1780 : /* -------------------------------------------------------------------- */
1781 68 : VSIFPrintfL( fp, "BYTEORDER I\n" );
1782 68 : VSIFPrintfL( fp, "LAYOUT BIL\n" );
1783 68 : VSIFPrintfL( fp, "NROWS %d\n", nYSize );
1784 68 : VSIFPrintfL( fp, "NCOLS %d\n", nXSize );
1785 68 : VSIFPrintfL( fp, "NBANDS %d\n", nBands );
1786 68 : VSIFPrintfL( fp, "NBITS %d\n", nBits );
1787 68 : VSIFPrintfL( fp, "BANDROWBYTES %d\n", nRowBytes );
1788 68 : VSIFPrintfL( fp, "TOTALROWBYTES %d\n", nRowBytes * nBands );
1789 :
1790 68 : if( eType == GDT_Float32 )
1791 10 : VSIFPrintfL( fp, "PIXELTYPE FLOAT\n");
1792 72 : else if( eType == GDT_Int16 || eType == GDT_Int32 )
1793 14 : VSIFPrintfL( fp, "PIXELTYPE SIGNEDINT\n");
1794 46 : else if( eType == GDT_Byte && EQUAL(pszPixelType,"SIGNEDBYTE") )
1795 2 : VSIFPrintfL( fp, "PIXELTYPE SIGNEDINT\n");
1796 : else
1797 42 : VSIFPrintfL( fp, "PIXELTYPE UNSIGNEDINT\n");
1798 :
1799 : /* -------------------------------------------------------------------- */
1800 : /* Cleanup */
1801 : /* -------------------------------------------------------------------- */
1802 68 : VSIFCloseL( fp );
1803 :
1804 68 : CPLFree( pszHdrFilename );
1805 :
1806 68 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
1807 : }
1808 :
1809 : /************************************************************************/
1810 : /* CreateCopy() */
1811 : /************************************************************************/
1812 :
1813 50 : GDALDataset *EHdrDataset::CreateCopy( const char * pszFilename,
1814 : GDALDataset * poSrcDS,
1815 : int bStrict, char ** papszOptions,
1816 : GDALProgressFunc pfnProgress,
1817 : void * pProgressData )
1818 :
1819 : {
1820 50 : char **papszAdjustedOptions = CSLDuplicate( papszOptions );
1821 : GDALDataset *poOutDS;
1822 :
1823 50 : int nBands = poSrcDS->GetRasterCount();
1824 50 : if (nBands == 0)
1825 : {
1826 : CPLError( CE_Failure, CPLE_NotSupported,
1827 2 : "EHdr driver does not support source dataset with zero band.\n");
1828 2 : return NULL;
1829 : }
1830 :
1831 : /* -------------------------------------------------------------------- */
1832 : /* Ensure we pass on NBITS and PIXELTYPE structure information. */
1833 : /* -------------------------------------------------------------------- */
1834 96 : if( poSrcDS->GetRasterBand(1)->GetMetadataItem( "NBITS",
1835 48 : "IMAGE_STRUCTURE" ) !=NULL
1836 : && CSLFetchNameValue( papszOptions, "NBITS" ) == NULL )
1837 : {
1838 : papszAdjustedOptions =
1839 : CSLSetNameValue( papszAdjustedOptions,
1840 : "NBITS",
1841 0 : poSrcDS->GetRasterBand(1)->GetMetadataItem("NBITS","IMAGE_STRUCTURE") );
1842 : }
1843 :
1844 96 : if( poSrcDS->GetRasterBand(1)->GetMetadataItem( "PIXELTYPE",
1845 48 : "IMAGE_STRUCTURE" ) !=NULL
1846 : && CSLFetchNameValue( papszOptions, "PIXELTYPE" ) == NULL )
1847 : {
1848 : papszAdjustedOptions =
1849 : CSLSetNameValue( papszAdjustedOptions,
1850 : "PIXELTYPE",
1851 2 : poSrcDS->GetRasterBand(1)->GetMetadataItem("PIXELTYPE","IMAGE_STRUCTURE") );
1852 : }
1853 :
1854 : /* -------------------------------------------------------------------- */
1855 : /* Proceed with normal copying using the default createcopy */
1856 : /* operators. */
1857 : /* -------------------------------------------------------------------- */
1858 48 : GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "EHdr" );
1859 :
1860 : poOutDS = poDriver->DefaultCreateCopy( pszFilename, poSrcDS, bStrict,
1861 : papszAdjustedOptions,
1862 48 : pfnProgress, pProgressData );
1863 48 : CSLDestroy( papszAdjustedOptions );
1864 :
1865 48 : return poOutDS;
1866 : }
1867 :
1868 : /************************************************************************/
1869 : /* GetMinimum() */
1870 : /************************************************************************/
1871 :
1872 6 : double EHdrRasterBand::GetMinimum( int *pbSuccess )
1873 : {
1874 6 : if( pbSuccess != NULL )
1875 6 : *pbSuccess = (minmaxmeanstddev & HAS_MIN_FLAG) != 0;
1876 :
1877 6 : if( minmaxmeanstddev & HAS_MIN_FLAG )
1878 4 : return dfMin;
1879 :
1880 2 : return RawRasterBand::GetMinimum( pbSuccess );
1881 : }
1882 :
1883 : /************************************************************************/
1884 : /* GetMaximum() */
1885 : /************************************************************************/
1886 :
1887 4 : double EHdrRasterBand::GetMaximum( int *pbSuccess )
1888 : {
1889 4 : if( pbSuccess != NULL )
1890 4 : *pbSuccess = (minmaxmeanstddev & HAS_MAX_FLAG) != 0;
1891 :
1892 4 : if( minmaxmeanstddev & HAS_MAX_FLAG )
1893 2 : return dfMax;
1894 :
1895 2 : return RawRasterBand::GetMaximum( pbSuccess );
1896 : }
1897 :
1898 : /************************************************************************/
1899 : /* GetStatistics() */
1900 : /************************************************************************/
1901 :
1902 4 : CPLErr EHdrRasterBand::GetStatistics( int bApproxOK, int bForce, double *pdfMin, double *pdfMax, double *pdfMean, double *pdfStdDev )
1903 : {
1904 4 : if( (minmaxmeanstddev & HAS_ALL_FLAGS) == HAS_ALL_FLAGS)
1905 : {
1906 2 : if ( pdfMin ) *pdfMin = dfMin;
1907 2 : if ( pdfMax ) *pdfMax = dfMax;
1908 2 : if ( pdfMean ) *pdfMean = dfMean;
1909 2 : if ( pdfStdDev ) *pdfStdDev = dfStdDev;
1910 2 : return CE_None;
1911 : }
1912 :
1913 : CPLErr eErr = RawRasterBand::GetStatistics( bApproxOK, bForce,
1914 : &dfMin, &dfMax,
1915 2 : &dfMean, &dfStdDev );
1916 :
1917 2 : if( eErr == CE_None )
1918 : {
1919 2 : EHdrDataset* poEDS = (EHdrDataset *) poDS;
1920 :
1921 2 : minmaxmeanstddev = HAS_ALL_FLAGS;
1922 :
1923 2 : if( poEDS->RewriteSTX() != CE_None )
1924 0 : RawRasterBand::SetStatistics( dfMin, dfMax, dfMean, dfStdDev );
1925 :
1926 2 : if( pdfMin )
1927 2 : *pdfMin = dfMin;
1928 2 : if( pdfMax )
1929 2 : *pdfMax = dfMax;
1930 2 : if( pdfMean )
1931 2 : *pdfMean = dfMean;
1932 2 : if( pdfStdDev )
1933 2 : *pdfStdDev = dfStdDev;
1934 : }
1935 :
1936 2 : return eErr;
1937 : }
1938 :
1939 : /************************************************************************/
1940 : /* SetStatistics() */
1941 : /************************************************************************/
1942 :
1943 2 : CPLErr EHdrRasterBand::SetStatistics( double dfMin, double dfMax, double dfMean, double dfStdDev )
1944 : {
1945 : // avoid churn if nothing is changing.
1946 2 : if( dfMin == this->dfMin
1947 : && dfMax == this->dfMax
1948 : && dfMean == this->dfMean
1949 : && dfStdDev == this->dfStdDev )
1950 0 : return CE_None;
1951 :
1952 2 : this->dfMin = dfMin;
1953 2 : this->dfMax = dfMax;
1954 2 : this->dfMean = dfMean;
1955 2 : this->dfStdDev = dfStdDev;
1956 :
1957 : // marks stats valid
1958 2 : minmaxmeanstddev = HAS_ALL_FLAGS;
1959 :
1960 2 : EHdrDataset* poEDS = (EHdrDataset *) poDS;
1961 :
1962 2 : if( poEDS->RewriteSTX() != CE_None )
1963 0 : return RawRasterBand::SetStatistics( dfMin, dfMax, dfMean, dfStdDev );
1964 : else
1965 2 : return CE_None;
1966 : }
1967 :
1968 : /************************************************************************/
1969 : /* SetColorTable() */
1970 : /************************************************************************/
1971 :
1972 12 : CPLErr EHdrRasterBand::SetColorTable( GDALColorTable *poNewCT )
1973 : {
1974 12 : CPLErr err = RawRasterBand::SetColorTable( poNewCT );
1975 12 : if( err != CE_None )
1976 0 : return err;
1977 :
1978 12 : ((EHdrDataset*)poDS)->bCLRDirty = TRUE;
1979 :
1980 12 : return CE_None;
1981 : }
1982 :
1983 : /************************************************************************/
1984 : /* GDALRegister_EHdr() */
1985 : /************************************************************************/
1986 :
1987 1135 : void GDALRegister_EHdr()
1988 :
1989 : {
1990 : GDALDriver *poDriver;
1991 :
1992 1135 : if( GDALGetDriverByName( "EHdr" ) == NULL )
1993 : {
1994 1093 : poDriver = new GDALDriver();
1995 :
1996 1093 : poDriver->SetDescription( "EHdr" );
1997 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1998 1093 : "ESRI .hdr Labelled" );
1999 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
2000 1093 : "frmt_various.html#EHdr" );
2001 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2002 1093 : "Byte Int16 UInt16 Int32 UInt32 Float32" );
2003 :
2004 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
2005 : "<CreationOptionList>"
2006 : " <Option name='NBITS' type='int' description='Special pixel bits (1-7)'/>"
2007 : " <Option name='PIXELTYPE' type='string' description='By setting this to SIGNEDBYTE, a new Byte file can be forced to be written as signed byte'/>"
2008 1093 : "</CreationOptionList>" );
2009 :
2010 1093 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2011 1093 : poDriver->pfnOpen = EHdrDataset::Open;
2012 1093 : poDriver->pfnCreate = EHdrDataset::Create;
2013 1093 : poDriver->pfnCreateCopy = EHdrDataset::CreateCopy;
2014 :
2015 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
2016 : }
2017 1135 : }
|
