1 : /******************************************************************************
2 : * $Id: ehdrdataset.cpp 24766 2012-08-11 18:50:50Z 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 24766 2012-08-11 18:50:50Z 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 126 : {
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 126 : 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 126 : minmaxmeanstddev(0)
164 : {
165 126 : EHdrDataset* poEDS = (EHdrDataset*)poDS;
166 :
167 126 : if (nBits < 8)
168 : {
169 3 : nStartBit = atoi(poEDS->GetKeyValue("SKIPBYTES")) * 8;
170 3 : 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 3 : nPixelOffsetBits = nBits;
180 3 : nLineOffsetBits = atoi(poEDS->GetKeyValue("TOTALROWBYTES")) * 8;
181 3 : if( nLineOffsetBits == 0 )
182 0 : nLineOffsetBits = nPixelOffsetBits * poDS->GetRasterXSize();
183 :
184 3 : nBlockXSize = poDS->GetRasterXSize();
185 3 : nBlockYSize = 1;
186 :
187 : SetMetadataItem( "NBITS",
188 : CPLString().Printf( "%d", nBits ),
189 3 : "IMAGE_STRUCTURE" );
190 : }
191 :
192 126 : if( eDataType == GDT_Byte
193 : && EQUAL(poEDS->GetKeyValue("PIXELTYPE",""),"SIGNEDINT") )
194 : SetMetadataItem( "PIXELTYPE", "SIGNEDBYTE",
195 2 : "IMAGE_STRUCTURE" );
196 126 : }
197 :
198 : /************************************************************************/
199 : /* IReadBlock() */
200 : /************************************************************************/
201 :
202 256 : CPLErr EHdrRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
203 : void * pImage )
204 :
205 : {
206 256 : if (nBits >= 8)
207 224 : 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 32 : nLineBytes = (nPixelOffsetBits*nBlockXSize + 7)/8;
218 32 : nLineStart = (nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) / 8;
219 : iBitOffset = (int)
220 32 : ((nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) % 8);
221 :
222 : /* -------------------------------------------------------------------- */
223 : /* Read data into buffer. */
224 : /* -------------------------------------------------------------------- */
225 32 : pabyBuffer = (GByte *) CPLCalloc(nLineBytes,1);
226 :
227 32 : 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 32 : int iPixel = 0, iX;
241 :
242 1056 : for( iX = 0; iX < nBlockXSize; iX++ )
243 : {
244 1024 : int nOutWord = 0, iBit;
245 :
246 2048 : for( iBit = 0; iBit < nBits; iBit++ )
247 : {
248 1024 : if( pabyBuffer[iBitOffset>>3] & (0x80 >>(iBitOffset & 7)) )
249 200 : nOutWord |= (1 << (nBits - 1 - iBit));
250 1024 : iBitOffset++;
251 : }
252 :
253 1024 : iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
254 :
255 1024 : ((GByte *) pImage)[iPixel++] = (GByte) nOutWord;
256 : }
257 :
258 32 : CPLFree( pabyBuffer );
259 :
260 32 : return CE_None;
261 : }
262 :
263 : /************************************************************************/
264 : /* IWriteBlock() */
265 : /************************************************************************/
266 :
267 472 : CPLErr EHdrRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
268 : void * pImage )
269 :
270 : {
271 472 : if (nBits >= 8)
272 440 : 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 32 : nLineBytes = (nPixelOffsetBits*nBlockXSize + 7)/8;
283 32 : nLineStart = (nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) / 8;
284 : iBitOffset = (int)
285 32 : ((nStartBit + ((vsi_l_offset)nLineOffsetBits) * nBlockYOff) % 8);
286 :
287 : /* -------------------------------------------------------------------- */
288 : /* Read data into buffer. */
289 : /* -------------------------------------------------------------------- */
290 32 : pabyBuffer = (GByte *) CPLCalloc(nLineBytes,1);
291 :
292 32 : 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 32 : VSIFReadL( pabyBuffer, 1, nLineBytes, GetFPL() );
302 :
303 : /* -------------------------------------------------------------------- */
304 : /* Copy data, promoting to 8bit. */
305 : /* -------------------------------------------------------------------- */
306 32 : int iPixel = 0, iX;
307 :
308 1056 : for( iX = 0; iX < nBlockXSize; iX++ )
309 : {
310 : int iBit;
311 1024 : int nOutWord = ((GByte *) pImage)[iPixel++];
312 :
313 2048 : for( iBit = 0; iBit < nBits; iBit++ )
314 : {
315 1024 : if( nOutWord & (1 << (nBits - 1 - iBit)) )
316 200 : pabyBuffer[iBitOffset>>3] |= (0x80 >>(iBitOffset & 7));
317 : else
318 824 : pabyBuffer[iBitOffset>>3] &= ~((0x80 >>(iBitOffset & 7)));
319 :
320 1024 : iBitOffset++;
321 : }
322 :
323 1024 : iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
324 : }
325 :
326 : /* -------------------------------------------------------------------- */
327 : /* Write the data back out. */
328 : /* -------------------------------------------------------------------- */
329 32 : 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 32 : CPLFree( pabyBuffer );
340 :
341 32 : return CE_None;
342 : }
343 :
344 : /************************************************************************/
345 : /* IRasterIO() */
346 : /************************************************************************/
347 :
348 408 : 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 408 : if (nBits >= 8)
357 : return RawRasterBand::IRasterIO( eRWFlag,
358 : nXOff, nYOff, nXSize, nYSize,
359 : pData, nBufXSize, nBufYSize,
360 375 : 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 33 : eBufType, nPixelSpace, nLineSpace );
368 : }
369 :
370 : /************************************************************************/
371 : /* OSR_GDS() */
372 : /************************************************************************/
373 :
374 16 : 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 16 : if( papszNV == NULL || papszNV[0] == NULL )
382 0 : return pszDefaultValue;
383 :
384 125 : for( iLine = 0;
385 62 : papszNV[iLine] != NULL &&
386 47 : !EQUALN(papszNV[iLine],pszField,strlen(pszField));
387 : iLine++ ) {}
388 :
389 16 : if( papszNV[iLine] == NULL )
390 15 : return pszDefaultValue;
391 : else
392 : {
393 : char **papszTokens;
394 :
395 1 : papszTokens = CSLTokenizeString(papszNV[iLine]);
396 :
397 1 : if( CSLCount(papszTokens) > 1 )
398 1 : strncpy( pszResult, papszTokens[1], nResultLen);
399 : else
400 0 : strncpy( pszResult, pszDefaultValue, nResultLen);
401 1 : pszResult[nResultLen-1] = '\0';
402 :
403 1 : CSLDestroy( papszTokens );
404 1 : return pszResult;
405 : }
406 : }
407 :
408 :
409 : /************************************************************************/
410 : /* ==================================================================== */
411 : /* EHdrDataset */
412 : /* ==================================================================== */
413 : /************************************************************************/
414 :
415 : /************************************************************************/
416 : /* EHdrDataset() */
417 : /************************************************************************/
418 :
419 76 : EHdrDataset::EHdrDataset()
420 : {
421 76 : fpImage = NULL;
422 76 : pszProjection = CPLStrdup("");
423 76 : bGotTransform = FALSE;
424 76 : adfGeoTransform[0] = 0.0;
425 76 : adfGeoTransform[1] = 1.0;
426 76 : adfGeoTransform[2] = 0.0;
427 76 : adfGeoTransform[3] = 0.0;
428 76 : adfGeoTransform[4] = 0.0;
429 76 : adfGeoTransform[5] = 1.0;
430 76 : papszHDR = NULL;
431 76 : bHDRDirty = FALSE;
432 76 : bCLRDirty = FALSE;
433 76 : osHeaderExt = "hdr";
434 76 : }
435 :
436 : /************************************************************************/
437 : /* ~EHdrDataset() */
438 : /************************************************************************/
439 :
440 76 : EHdrDataset::~EHdrDataset()
441 :
442 : {
443 76 : FlushCache();
444 :
445 76 : if( nBands > 0 && GetAccess() == GA_Update )
446 : {
447 : int bNoDataSet;
448 : double dfNoData;
449 34 : RawRasterBand *poBand = (RawRasterBand *) GetRasterBand( 1 );
450 :
451 34 : dfNoData = poBand->GetNoDataValue(&bNoDataSet);
452 34 : if( bNoDataSet )
453 : {
454 : ResetKeyValue( "NODATA",
455 1 : CPLString().Printf( "%.8g", dfNoData ) );
456 : }
457 :
458 34 : if( bCLRDirty )
459 2 : RewriteColorTable( poBand->GetColorTable() );
460 :
461 34 : if( bHDRDirty )
462 33 : RewriteHDR();
463 : }
464 :
465 76 : if( fpImage != NULL )
466 76 : VSIFCloseL( fpImage );
467 :
468 76 : CPLFree( pszProjection );
469 76 : CSLDestroy( papszHDR );
470 76 : }
471 :
472 : /************************************************************************/
473 : /* GetKeyValue() */
474 : /************************************************************************/
475 :
476 66 : const char *EHdrDataset::GetKeyValue( const char *pszKey,
477 : const char *pszDefault )
478 :
479 : {
480 : int i;
481 :
482 602 : for( i = 0; papszHDR[i] != NULL; i++ )
483 : {
484 662 : if( EQUALN(pszKey,papszHDR[i],strlen(pszKey))
485 63 : && isspace((unsigned char)papszHDR[i][strlen(pszKey)]) )
486 : {
487 63 : const char *pszValue = papszHDR[i] + strlen(pszKey);
488 492 : while( isspace((unsigned char)*pszValue) )
489 366 : pszValue++;
490 :
491 63 : return pszValue;
492 : }
493 : }
494 :
495 3 : 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 129 : void EHdrDataset::ResetKeyValue( const char *pszKey, const char *pszValue )
506 :
507 : {
508 : int i;
509 : char szNewLine[82];
510 :
511 129 : if( strlen(pszValue) > 65 )
512 : {
513 0 : CPLAssert( strlen(pszValue) <= 65 );
514 0 : return;
515 : }
516 :
517 129 : sprintf( szNewLine, "%-15s%s", pszKey, pszValue );
518 :
519 1482 : for( i = CSLCount(papszHDR)-1; i >= 0; i-- )
520 : {
521 1353 : 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 129 : bHDRDirty = TRUE;
534 129 : papszHDR = CSLAddString( papszHDR, szNewLine );
535 : }
536 :
537 : /************************************************************************/
538 : /* RewriteColorTable() */
539 : /************************************************************************/
540 :
541 2 : void EHdrDataset::RewriteColorTable( GDALColorTable *poTable )
542 :
543 : {
544 2 : CPLString osCLRFilename = CPLResetExtension( GetDescription(), "clr" );
545 2 : if( poTable )
546 : {
547 2 : VSILFILE *fp = VSIFOpenL( osCLRFilename, "wt" );
548 2 : if( fp != NULL )
549 : {
550 8 : for( int iColor = 0; iColor < poTable->GetColorEntryCount(); iColor++ )
551 : {
552 6 : CPLString oLine;
553 : GDALColorEntry sEntry;
554 :
555 6 : poTable->GetColorEntryAsRGB( iColor, &sEntry );
556 :
557 : // I wish we had a way to mark transparency.
558 : oLine.Printf( "%3d %3d %3d %3d\n",
559 6 : iColor, sEntry.c1, sEntry.c2, sEntry.c3 );
560 6 : VSIFWriteL( (void *) oLine.c_str(), 1, strlen(oLine), fp );
561 : }
562 2 : 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 2 : }
574 :
575 : /************************************************************************/
576 : /* GetProjectionRef() */
577 : /************************************************************************/
578 :
579 16 : const char *EHdrDataset::GetProjectionRef()
580 :
581 : {
582 16 : if (pszProjection && strlen(pszProjection) > 0)
583 16 : return pszProjection;
584 :
585 0 : return GDALPamDataset::GetProjectionRef();
586 : }
587 :
588 : /************************************************************************/
589 : /* SetProjection() */
590 : /************************************************************************/
591 :
592 30 : CPLErr EHdrDataset::SetProjection( const char *pszSRS )
593 :
594 : {
595 : /* -------------------------------------------------------------------- */
596 : /* Reset coordinate system on the dataset. */
597 : /* -------------------------------------------------------------------- */
598 30 : CPLFree( pszProjection );
599 30 : pszProjection = CPLStrdup( pszSRS );
600 :
601 30 : if( strlen(pszSRS) == 0 )
602 0 : return CE_None;
603 :
604 : /* -------------------------------------------------------------------- */
605 : /* Convert to ESRI WKT. */
606 : /* -------------------------------------------------------------------- */
607 30 : OGRSpatialReference oSRS( pszSRS );
608 30 : char *pszESRI_SRS = NULL;
609 :
610 30 : oSRS.morphToESRI();
611 30 : oSRS.exportToWkt( &pszESRI_SRS );
612 :
613 : /* -------------------------------------------------------------------- */
614 : /* Write to .prj file. */
615 : /* -------------------------------------------------------------------- */
616 30 : CPLString osPrjFilename = CPLResetExtension( GetDescription(), "prj" );
617 : VSILFILE *fp;
618 :
619 30 : fp = VSIFOpenL( osPrjFilename.c_str(), "wt" );
620 30 : if( fp != NULL )
621 : {
622 30 : VSIFWriteL( pszESRI_SRS, 1, strlen(pszESRI_SRS), fp );
623 30 : VSIFWriteL( (void *) "\n", 1, 1, fp );
624 30 : VSIFCloseL( fp );
625 : }
626 :
627 30 : CPLFree( pszESRI_SRS );
628 :
629 30 : return CE_None;
630 : }
631 :
632 : /************************************************************************/
633 : /* GetGeoTransform() */
634 : /************************************************************************/
635 :
636 24 : CPLErr EHdrDataset::GetGeoTransform( double * padfTransform )
637 :
638 : {
639 24 : if( bGotTransform )
640 : {
641 24 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
642 24 : return CE_None;
643 : }
644 : else
645 : {
646 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
647 : }
648 : }
649 :
650 : /************************************************************************/
651 : /* SetGeoTransform() */
652 : /************************************************************************/
653 :
654 32 : CPLErr EHdrDataset::SetGeoTransform( double *padfGeoTransform )
655 :
656 : {
657 : /* -------------------------------------------------------------------- */
658 : /* We only support non-rotated images with info in the .HDR file. */
659 : /* -------------------------------------------------------------------- */
660 64 : if( padfGeoTransform[2] != 0.0
661 32 : || padfGeoTransform[4] != 0.0 )
662 : {
663 0 : return GDALPamDataset::SetGeoTransform( padfGeoTransform );
664 : }
665 :
666 : /* -------------------------------------------------------------------- */
667 : /* Record new geotransform. */
668 : /* -------------------------------------------------------------------- */
669 32 : bGotTransform = TRUE;
670 32 : memcpy( adfGeoTransform, padfGeoTransform, sizeof(double) * 6 );
671 :
672 : /* -------------------------------------------------------------------- */
673 : /* Strip out all old geotransform keywords from HDR records. */
674 : /* -------------------------------------------------------------------- */
675 : int i;
676 320 : for( i = CSLCount(papszHDR)-1; i >= 0; i-- )
677 : {
678 1152 : if( EQUALN(papszHDR[i],"ul",2)
679 288 : || EQUALN(papszHDR[i]+1,"ll",2)
680 288 : || EQUALN(papszHDR[i],"cell",4)
681 288 : || 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 32 : CPLString oValue;
691 :
692 32 : oValue.Printf( "%.15g", adfGeoTransform[0] + adfGeoTransform[1] * 0.5 );
693 32 : ResetKeyValue( "ULXMAP", oValue );
694 :
695 32 : oValue.Printf( "%.15g", adfGeoTransform[3] + adfGeoTransform[5] * 0.5 );
696 32 : ResetKeyValue( "ULYMAP", oValue );
697 :
698 32 : oValue.Printf( "%.15g", adfGeoTransform[1] );
699 32 : ResetKeyValue( "XDIM", oValue );
700 :
701 32 : oValue.Printf( "%.15g", fabs(adfGeoTransform[5]) );
702 32 : ResetKeyValue( "YDIM", oValue );
703 :
704 32 : return CE_None;
705 : }
706 :
707 : /************************************************************************/
708 : /* RewriteHDR() */
709 : /************************************************************************/
710 :
711 33 : CPLErr EHdrDataset::RewriteHDR()
712 :
713 : {
714 33 : CPLString osPath = CPLGetPath( GetDescription() );
715 33 : CPLString osName = CPLGetBasename( GetDescription() );
716 33 : CPLString osHDRFilename = CPLFormCIFilename( osPath, osName, osHeaderExt );
717 :
718 : /* -------------------------------------------------------------------- */
719 : /* Write .hdr file. */
720 : /* -------------------------------------------------------------------- */
721 : VSILFILE *fp;
722 : int i;
723 :
724 33 : fp = VSIFOpenL( osHDRFilename, "wt" );
725 :
726 33 : 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 459 : for( i = 0; papszHDR[i] != NULL; i++ )
735 : {
736 426 : VSIFWriteL( papszHDR[i], 1, strlen(papszHDR[i]), fp );
737 426 : VSIFWriteL( (void *) "\n", 1, 1, fp );
738 : }
739 :
740 33 : VSIFCloseL( fp );
741 :
742 33 : bHDRDirty = FALSE;
743 :
744 33 : return CE_None;
745 : }
746 :
747 : /************************************************************************/
748 : /* RewriteSTX() */
749 : /************************************************************************/
750 :
751 2 : CPLErr EHdrDataset::RewriteSTX()
752 : {
753 2 : CPLString osPath = CPLGetPath( GetDescription() );
754 2 : CPLString osName = CPLGetBasename( GetDescription() );
755 2 : CPLString osSTXFilename = CPLFormCIFilename( osPath, osName, "stx" );
756 :
757 : /* -------------------------------------------------------------------- */
758 : /* Write .stx file. */
759 : /* -------------------------------------------------------------------- */
760 : VSILFILE *fp;
761 2 : fp = VSIFOpenL( osSTXFilename, "wt" );
762 2 : 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 4 : for (int i = 0; i < nBands; ++i)
770 : {
771 2 : EHdrRasterBand* poBand = (EHdrRasterBand*)papoBands[i];
772 2 : VSIFPrintfL( fp, "%d %.10f %.10f ", i+1, poBand->dfMin, poBand->dfMax );
773 2 : if ( poBand->minmaxmeanstddev & HAS_MEAN_FLAG )
774 2 : VSIFPrintfL( fp, "%.10f ", poBand->dfMean);
775 : else
776 0 : VSIFPrintfL( fp, "# ");
777 :
778 2 : if ( poBand->minmaxmeanstddev & HAS_STDDEV_FLAG )
779 2 : VSIFPrintfL( fp, "%.10f\n", poBand->dfStdDev);
780 : else
781 0 : VSIFPrintfL( fp, "#\n");
782 : }
783 :
784 2 : VSIFCloseL( fp );
785 :
786 2 : return CE_None;
787 : }
788 :
789 : /************************************************************************/
790 : /* ReadSTX() */
791 : /************************************************************************/
792 :
793 76 : CPLErr EHdrDataset::ReadSTX()
794 : {
795 76 : CPLString osPath = CPLGetPath( GetDescription() );
796 76 : CPLString osName = CPLGetBasename( GetDescription() );
797 76 : CPLString osSTXFilename = CPLFormCIFilename( osPath, osName, "stx" );
798 :
799 : /* -------------------------------------------------------------------- */
800 : /* Read .stx file. */
801 : /* -------------------------------------------------------------------- */
802 : VSILFILE *fp;
803 76 : if ((fp = VSIFOpenL( osSTXFilename, "rt" )) != NULL)
804 : {
805 : const char * pszLine;
806 6 : while( (pszLine = CPLReadLineL( fp )) != NULL )
807 : {
808 : char **papszTokens;
809 2 : papszTokens = CSLTokenizeStringComplex( pszLine, " \t", TRUE, FALSE );
810 2 : int nTokens = CSLCount( papszTokens );
811 2 : if( nTokens >= 5 )
812 : {
813 2 : int i = atoi(papszTokens[0]);
814 2 : if (i > 0 && i <= nBands)
815 : {
816 2 : EHdrRasterBand* poBand = (EHdrRasterBand*)papoBands[i-1];
817 2 : poBand->dfMin = atof(papszTokens[1]);
818 2 : poBand->dfMax = atof(papszTokens[2]);
819 :
820 2 : int bNoDataSet = FALSE;
821 2 : double dfNoData = poBand->GetNoDataValue(&bNoDataSet);
822 2 : 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 2 : poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
834 : // reads optional mean and stddev
835 2 : if ( !EQUAL(papszTokens[3], "#") )
836 : {
837 2 : poBand->dfMean = atof(papszTokens[3]);
838 2 : poBand->minmaxmeanstddev |= HAS_MEAN_FLAG;
839 : }
840 2 : if ( !EQUAL(papszTokens[4], "#") )
841 : {
842 2 : poBand->dfStdDev = atof(papszTokens[4]);
843 2 : poBand->minmaxmeanstddev |= HAS_STDDEV_FLAG;
844 : }
845 :
846 2 : if( nTokens >= 6 && !EQUAL(papszTokens[5], "#") )
847 0 : poBand->SetMetadataItem( "STRETCHMIN", papszTokens[5], "RENDERING_HINTS" );
848 :
849 2 : if( nTokens >= 7 && !EQUAL(papszTokens[6], "#") )
850 0 : poBand->SetMetadataItem( "STRETCHMAX", papszTokens[6], "RENDERING_HINTS" );
851 : }
852 : }
853 :
854 2 : CSLDestroy( papszTokens );
855 : }
856 :
857 2 : VSIFCloseL( fp );
858 : }
859 :
860 76 : 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 53 : CPLString EHdrDataset::GetImageRepFilename(const char* pszFilename)
876 : {
877 : VSIStatBufL sStatBuf;
878 :
879 53 : CPLString osPath = CPLGetPath( pszFilename );
880 53 : CPLString osName = CPLGetBasename( pszFilename );
881 : CPLString osREPFilename =
882 53 : CPLFormCIFilename( osPath, osName, "rep" );
883 53 : if( VSIStatExL( osREPFilename.c_str(), &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
884 0 : return osREPFilename;
885 :
886 53 : 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 53 : return CPLString();
919 : }
920 :
921 : /************************************************************************/
922 : /* GetFileList() */
923 : /************************************************************************/
924 :
925 9 : char **EHdrDataset::GetFileList()
926 :
927 : {
928 : VSIStatBufL sStatBuf;
929 9 : CPLString osPath = CPLGetPath( GetDescription() );
930 9 : CPLString osName = CPLGetBasename( GetDescription() );
931 9 : char **papszFileList = NULL;
932 :
933 : // Main data file, etc.
934 9 : papszFileList = GDALPamDataset::GetFileList();
935 :
936 : // Header file.
937 9 : CPLString osFilename = CPLFormCIFilename( osPath, osName, osHeaderExt );
938 9 : papszFileList = CSLAddString( papszFileList, osFilename );
939 :
940 : // Statistics file
941 9 : osFilename = CPLFormCIFilename( osPath, osName, "stx" );
942 9 : if( VSIStatExL( osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
943 1 : papszFileList = CSLAddString( papszFileList, osFilename );
944 :
945 : // color table file.
946 9 : osFilename = CPLFormCIFilename( osPath, osName, "clr" );
947 9 : if( VSIStatExL( osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
948 2 : papszFileList = CSLAddString( papszFileList, osFilename );
949 :
950 : // projections file.
951 9 : osFilename = CPLFormCIFilename( osPath, osName, "prj" );
952 9 : if( VSIStatExL( osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
953 5 : papszFileList = CSLAddString( papszFileList, osFilename );
954 :
955 9 : CPLString imageRepFilename = GetImageRepFilename( GetDescription() );
956 9 : if (!imageRepFilename.empty())
957 0 : papszFileList = CSLAddString( papszFileList, imageRepFilename.c_str() );
958 :
959 9 : return papszFileList;
960 : }
961 :
962 : /************************************************************************/
963 : /* Open() */
964 : /************************************************************************/
965 :
966 12225 : 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 12225 : if( poOpenInfo->nHeaderBytes < 2 )
975 11319 : return NULL;
976 :
977 : /* -------------------------------------------------------------------- */
978 : /* Now we need to tear apart the filename to form a .HDR */
979 : /* filename. */
980 : /* -------------------------------------------------------------------- */
981 906 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
982 906 : CPLString osName = CPLGetBasename( poOpenInfo->pszFilename );
983 906 : CPLString osHDRFilename;
984 :
985 906 : const char* pszHeaderExt = "hdr";
986 906 : 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 906 : if( poOpenInfo->papszSiblingFiles )
997 : {
998 : int iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
999 902 : CPLFormFilename( NULL, osName, pszHeaderExt ) );
1000 902 : if( iFile < 0 ) // return if there is no corresponding .hdr file
1001 795 : return NULL;
1002 :
1003 : osHDRFilename =
1004 107 : CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1005 214 : NULL );
1006 : }
1007 : else
1008 : {
1009 4 : osHDRFilename = CPLFormCIFilename( osPath, osName, pszHeaderExt );
1010 : }
1011 :
1012 111 : bSelectedHDR = EQUAL( osHDRFilename, poOpenInfo->pszFilename );
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Do we have a .hdr file? */
1016 : /* -------------------------------------------------------------------- */
1017 : VSILFILE *fp;
1018 :
1019 111 : fp = VSIFOpenL( osHDRFilename, "r" );
1020 :
1021 111 : if( fp == NULL )
1022 : {
1023 4 : 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 107 : int nRows = -1, nCols = -1, nBands = 1;
1032 107 : int nSkipBytes = 0;
1033 107 : double dfULXMap=0.5, dfULYMap = 0.5, dfYLLCorner = -123.456;
1034 107 : int bCenter = TRUE;
1035 107 : double dfXDim = 1.0, dfYDim = 1.0, dfNoData = 0.0;
1036 107 : int nLineCount = 0, bNoDataSet = FALSE;
1037 107 : GDALDataType eDataType = GDT_Byte;
1038 107 : int nBits = -1;
1039 107 : char chByteOrder = 'M';
1040 107 : char chPixelType = 'N'; // not defined
1041 107 : char szLayout[10] = "BIL";
1042 107 : char **papszHDR = NULL;
1043 107 : int bHasInternalProjection = FALSE;
1044 107 : int bHasMin = FALSE;
1045 107 : int bHasMax = FALSE;
1046 107 : double dfMin = 0, dfMax = 0;
1047 :
1048 1375 : while( (pszLine = CPLReadLineL( fp )) != NULL )
1049 : {
1050 : char **papszTokens;
1051 :
1052 1163 : nLineCount++;
1053 :
1054 1163 : if( nLineCount > 50 || strlen(pszLine) > 1000 )
1055 2 : break;
1056 :
1057 1161 : papszHDR = CSLAddString( papszHDR, pszLine );
1058 :
1059 1161 : papszTokens = CSLTokenizeStringComplex( pszLine, " \t", TRUE, FALSE );
1060 1161 : if( CSLCount( papszTokens ) < 2 )
1061 : {
1062 49 : CSLDestroy( papszTokens );
1063 49 : continue;
1064 : }
1065 :
1066 1112 : if( EQUAL(papszTokens[0],"ncols") )
1067 : {
1068 76 : nCols = atoi(papszTokens[1]);
1069 : }
1070 1036 : else if( EQUAL(papszTokens[0],"nrows") )
1071 : {
1072 76 : nRows = atoi(papszTokens[1]);
1073 : }
1074 960 : else if( EQUAL(papszTokens[0],"skipbytes") )
1075 : {
1076 0 : nSkipBytes = atoi(papszTokens[1]);
1077 : }
1078 2844 : else if( EQUAL(papszTokens[0],"ulxmap")
1079 924 : || EQUAL(papszTokens[0],"xllcorner")
1080 922 : || EQUAL(papszTokens[0],"xllcenter") )
1081 : {
1082 38 : dfULXMap = CPLAtofM(papszTokens[1]);
1083 38 : if( EQUAL(papszTokens[0],"xllcorner") )
1084 2 : bCenter = FALSE;
1085 : }
1086 922 : else if( EQUAL(papszTokens[0],"ulymap") )
1087 : {
1088 36 : dfULYMap = CPLAtofM(papszTokens[1]);
1089 : }
1090 1772 : else if( EQUAL(papszTokens[0],"yllcorner")
1091 884 : || EQUAL(papszTokens[0],"yllcenter") )
1092 : {
1093 2 : dfYLLCorner = CPLAtofM(papszTokens[1]);
1094 2 : if( EQUAL(papszTokens[0],"yllcorner") )
1095 2 : bCenter = FALSE;
1096 : }
1097 884 : else if( EQUAL(papszTokens[0],"xdim") )
1098 : {
1099 36 : dfXDim = CPLAtofM(papszTokens[1]);
1100 : }
1101 848 : else if( EQUAL(papszTokens[0],"ydim") )
1102 : {
1103 36 : dfYDim = CPLAtofM(papszTokens[1]);
1104 : }
1105 812 : else if( EQUAL(papszTokens[0],"cellsize") )
1106 : {
1107 2 : dfXDim = dfYDim = CPLAtofM(papszTokens[1]);
1108 : }
1109 810 : else if( EQUAL(papszTokens[0],"nbands") )
1110 : {
1111 74 : nBands = atoi(papszTokens[1]);
1112 : }
1113 736 : else if( EQUAL(papszTokens[0],"layout") )
1114 : {
1115 74 : strncpy( szLayout, papszTokens[1], sizeof(szLayout) );
1116 74 : szLayout[sizeof(szLayout)-1] = '\0';
1117 : }
1118 1328 : else if( EQUAL(papszTokens[0],"NODATA_value")
1119 662 : || EQUAL(papszTokens[0],"NODATA") )
1120 : {
1121 4 : dfNoData = CPLAtofM(papszTokens[1]);
1122 4 : bNoDataSet = TRUE;
1123 : }
1124 658 : else if( EQUAL(papszTokens[0],"NBITS") )
1125 : {
1126 75 : nBits = atoi(papszTokens[1]);
1127 : }
1128 583 : else if( EQUAL(papszTokens[0],"PIXELTYPE") )
1129 : {
1130 72 : chPixelType = (char) toupper(papszTokens[1][0]);
1131 : }
1132 511 : else if( EQUAL(papszTokens[0],"byteorder") )
1133 : {
1134 76 : chByteOrder = (char) toupper(papszTokens[1][0]);
1135 : }
1136 :
1137 : /* http://www.worldclim.org/futdown.htm have the projection extensions */
1138 435 : else if( EQUAL(papszTokens[0],"Projection") )
1139 : {
1140 3 : bHasInternalProjection = TRUE;
1141 : }
1142 864 : else if( EQUAL(papszTokens[0],"MinValue") ||
1143 431 : EQUAL(papszTokens[0],"MIN_VALUE") )
1144 : {
1145 1 : dfMin = CPLAtofM(papszTokens[1]);
1146 1 : bHasMin = TRUE;
1147 : }
1148 861 : else if( EQUAL(papszTokens[0],"MaxValue") ||
1149 430 : EQUAL(papszTokens[0],"MAX_VALUE") )
1150 : {
1151 1 : dfMax = CPLAtofM(papszTokens[1]);
1152 1 : bHasMax = TRUE;
1153 : }
1154 :
1155 1112 : CSLDestroy( papszTokens );
1156 : }
1157 :
1158 107 : 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 107 : if( nRows == -1 || nCols == -1 )
1166 : {
1167 31 : CSLDestroy( papszHDR );
1168 31 : return NULL;
1169 : }
1170 :
1171 76 : 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 76 : 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 76 : if( nBits == -1 && chPixelType == 'N' )
1200 : {
1201 : VSIStatBufL sStatBuf;
1202 1 : if( VSIStatL( poOpenInfo->pszFilename, &sStatBuf ) == 0 )
1203 : {
1204 1 : size_t nBytes = (size_t) (sStatBuf.st_size/nCols/nRows/nBands);
1205 1 : if( nBytes > 0 && nBytes != 3 )
1206 1 : nBits = nBytes*8;
1207 :
1208 1 : if( nBytes == 4 )
1209 1 : chPixelType = 'F';
1210 : }
1211 : }
1212 :
1213 : /* -------------------------------------------------------------------- */
1214 : /* If the extension is FLT it is likely a floating point file. */
1215 : /* -------------------------------------------------------------------- */
1216 76 : if( chPixelType == 'N' )
1217 : {
1218 3 : if( EQUAL( CPLGetExtension( poOpenInfo->pszFilename ), "FLT" ) )
1219 1 : 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 76 : if( bNoDataSet && dfNoData < 0 && chPixelType == 'N' )
1228 : {
1229 2 : chPixelType = 'S';
1230 : }
1231 :
1232 : /* -------------------------------------------------------------------- */
1233 : /* Create a corresponding GDALDataset. */
1234 : /* -------------------------------------------------------------------- */
1235 : EHdrDataset *poDS;
1236 :
1237 76 : poDS = new EHdrDataset();
1238 :
1239 152 : poDS->osHeaderExt = pszHeaderExt;
1240 :
1241 : /* -------------------------------------------------------------------- */
1242 : /* Capture some information from the file that is of interest. */
1243 : /* -------------------------------------------------------------------- */
1244 76 : poDS->nRasterXSize = nCols;
1245 76 : poDS->nRasterYSize = nRows;
1246 76 : poDS->papszHDR = papszHDR;
1247 :
1248 : /* -------------------------------------------------------------------- */
1249 : /* Open target binary file. */
1250 : /* -------------------------------------------------------------------- */
1251 76 : if( poOpenInfo->eAccess == GA_ReadOnly )
1252 42 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
1253 : else
1254 34 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
1255 :
1256 76 : 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 76 : poDS->eAccess = poOpenInfo->eAccess;
1266 :
1267 : /* -------------------------------------------------------------------- */
1268 : /* Figure out the data type. */
1269 : /* -------------------------------------------------------------------- */
1270 76 : if( nBits == 16 )
1271 : {
1272 15 : if ( chPixelType == 'S' )
1273 7 : eDataType = GDT_Int16;
1274 : else
1275 8 : eDataType = GDT_UInt16; // default
1276 : }
1277 61 : else if( nBits == 32 )
1278 : {
1279 31 : if( chPixelType == 'S' )
1280 8 : eDataType = GDT_Int32;
1281 23 : else if( chPixelType == 'F' )
1282 18 : eDataType = GDT_Float32;
1283 : else
1284 5 : eDataType = GDT_UInt32; // default
1285 : }
1286 30 : else if( nBits == 8 )
1287 : {
1288 27 : eDataType = GDT_Byte;
1289 27 : nBits = 8;
1290 : }
1291 6 : else if( nBits < 8 && nBits >= 1 )
1292 3 : 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 76 : int nItemSize = GDALGetDataTypeSize(eDataType)/8;
1319 : int nPixelOffset, nLineOffset;
1320 : vsi_l_offset nBandOffset;
1321 :
1322 76 : if( EQUAL(szLayout,"BIP") )
1323 : {
1324 0 : nPixelOffset = nItemSize * nBands;
1325 0 : nLineOffset = nPixelOffset * nCols;
1326 0 : nBandOffset = (vsi_l_offset)nItemSize;
1327 : }
1328 76 : 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 76 : nPixelOffset = nItemSize;
1337 76 : nLineOffset = nItemSize * nBands * nCols;
1338 76 : nBandOffset = (vsi_l_offset)nItemSize * nCols;
1339 : }
1340 :
1341 76 : poDS->SetDescription( poOpenInfo->pszFilename );
1342 76 : poDS->PamInitialize();
1343 :
1344 : /* -------------------------------------------------------------------- */
1345 : /* Create band information objects. */
1346 : /* -------------------------------------------------------------------- */
1347 76 : poDS->nBands = nBands;
1348 202 : 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 126 : nBits);
1362 :
1363 126 : if( bNoDataSet )
1364 4 : poBand->SetNoDataValue( dfNoData );
1365 :
1366 126 : if( bHasMin && bHasMax )
1367 : {
1368 1 : poBand->dfMin = dfMin;
1369 1 : poBand->dfMax = dfMax;
1370 1 : poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
1371 : }
1372 :
1373 126 : poDS->SetBand( i+1, poBand );
1374 : }
1375 :
1376 : /* -------------------------------------------------------------------- */
1377 : /* If we didn't get bounds in the .hdr, look for a worldfile. */
1378 : /* -------------------------------------------------------------------- */
1379 76 : if( dfYLLCorner != -123.456 )
1380 : {
1381 2 : if( bCenter )
1382 0 : dfULYMap = dfYLLCorner + (nRows-1) * dfYDim;
1383 : else
1384 2 : dfULYMap = dfYLLCorner + nRows * dfYDim;
1385 : }
1386 :
1387 76 : if( dfULYMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0 )
1388 : {
1389 38 : poDS->bGotTransform = TRUE;
1390 :
1391 38 : if( bCenter )
1392 : {
1393 36 : poDS->adfGeoTransform[0] = dfULXMap - dfXDim * 0.5;
1394 36 : poDS->adfGeoTransform[1] = dfXDim;
1395 36 : poDS->adfGeoTransform[2] = 0.0;
1396 36 : poDS->adfGeoTransform[3] = dfULYMap + dfYDim * 0.5;
1397 36 : poDS->adfGeoTransform[4] = 0.0;
1398 36 : poDS->adfGeoTransform[5] = - dfYDim;
1399 : }
1400 : else
1401 : {
1402 2 : poDS->adfGeoTransform[0] = dfULXMap;
1403 2 : poDS->adfGeoTransform[1] = dfXDim;
1404 2 : poDS->adfGeoTransform[2] = 0.0;
1405 2 : poDS->adfGeoTransform[3] = dfULYMap;
1406 2 : poDS->adfGeoTransform[4] = 0.0;
1407 2 : poDS->adfGeoTransform[5] = - dfYDim;
1408 : }
1409 : }
1410 :
1411 76 : if( !poDS->bGotTransform )
1412 : poDS->bGotTransform =
1413 : GDALReadWorldFile( poOpenInfo->pszFilename, 0,
1414 38 : poDS->adfGeoTransform );
1415 :
1416 76 : if( !poDS->bGotTransform )
1417 : poDS->bGotTransform =
1418 : GDALReadWorldFile( poOpenInfo->pszFilename, "wld",
1419 38 : poDS->adfGeoTransform );
1420 :
1421 : /* -------------------------------------------------------------------- */
1422 : /* Check for a .prj file. */
1423 : /* -------------------------------------------------------------------- */
1424 76 : const char *pszPrjFilename = CPLFormCIFilename( osPath, osName, "prj" );
1425 :
1426 76 : 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 76 : if (fp == NULL && bHasInternalProjection)
1431 : {
1432 1 : pszPrjFilename = osHDRFilename;
1433 1 : fp = VSIFOpenL( pszPrjFilename, "r" );
1434 : }
1435 :
1436 76 : if( fp != NULL )
1437 : {
1438 : char **papszLines;
1439 32 : OGRSpatialReference oSRS;
1440 :
1441 32 : VSIFCloseL( fp );
1442 :
1443 32 : papszLines = CSLLoad( pszPrjFilename );
1444 :
1445 32 : 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 32 : 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 32 : CPLFree( poDS->pszProjection );
1463 32 : oSRS.exportToWkt( &(poDS->pszProjection) );
1464 : }
1465 :
1466 32 : 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 44 : CPLString szImageRepFilename = GetImageRepFilename(poOpenInfo->pszFilename );
1477 44 : if (!szImageRepFilename.empty())
1478 : {
1479 0 : fp = VSIFOpenL( szImageRepFilename.c_str(), "r" );
1480 : }
1481 44 : 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 44 : }
1598 : }
1599 :
1600 : /* -------------------------------------------------------------------- */
1601 : /* Check for a color table. */
1602 : /* -------------------------------------------------------------------- */
1603 76 : const char *pszCLRFilename = CPLFormCIFilename( osPath, osName, "clr" );
1604 :
1605 : /* Only read the .clr for byte, int16 or uint16 bands */
1606 76 : if (nItemSize <= 2)
1607 45 : fp = VSIFOpenL( pszCLRFilename, "r" );
1608 : else
1609 31 : fp = NULL;
1610 :
1611 76 : if( fp != NULL )
1612 : {
1613 4 : GDALColorTable oColorTable;
1614 4 : int bHasWarned = FALSE;
1615 :
1616 12 : while(TRUE)
1617 : {
1618 16 : const char *pszLine = CPLReadLineL(fp);
1619 16 : if ( !pszLine )
1620 : break;
1621 :
1622 12 : if( *pszLine == '#' || *pszLine == '!' )
1623 0 : continue;
1624 :
1625 : char **papszValues = CSLTokenizeString2(pszLine, "\t ",
1626 12 : CSLT_HONOURSTRINGS);
1627 : GDALColorEntry oEntry;
1628 :
1629 12 : if ( CSLCount(papszValues) >= 4 )
1630 : {
1631 12 : int nIndex = atoi( papszValues[0] ); // Index
1632 24 : if (nIndex >= 0 && nIndex < 65536)
1633 : {
1634 12 : oEntry.c1 = (short) atoi( papszValues[1] ); // Red
1635 12 : oEntry.c2 = (short) atoi( papszValues[2] ); // Green
1636 12 : oEntry.c3 = (short) atoi( papszValues[3] ); // Blue
1637 12 : oEntry.c4 = 255;
1638 :
1639 12 : 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 12 : CSLDestroy( papszValues );
1654 : }
1655 :
1656 4 : VSIFCloseL( fp );
1657 :
1658 8 : for( i = 1; i <= poDS->nBands; i++ )
1659 : {
1660 4 : GDALRasterBand *poBand = poDS->GetRasterBand( i );
1661 4 : poBand->SetColorTable( &oColorTable );
1662 4 : poBand->SetColorInterpretation( GCI_PaletteIndex );
1663 : }
1664 :
1665 4 : poDS->bCLRDirty = FALSE;
1666 : }
1667 :
1668 : /* -------------------------------------------------------------------- */
1669 : /* Read statistics (.STX) */
1670 : /* -------------------------------------------------------------------- */
1671 76 : poDS->ReadSTX();
1672 :
1673 : /* -------------------------------------------------------------------- */
1674 : /* Initialize any PAM information. */
1675 : /* -------------------------------------------------------------------- */
1676 76 : poDS->TryLoadXML();
1677 :
1678 : /* -------------------------------------------------------------------- */
1679 : /* Check for overviews. */
1680 : /* -------------------------------------------------------------------- */
1681 76 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1682 :
1683 76 : return( poDS );
1684 : }
1685 :
1686 : /************************************************************************/
1687 : /* Create() */
1688 : /************************************************************************/
1689 :
1690 62 : 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 62 : if (nBands <= 0)
1700 : {
1701 : CPLError( CE_Failure, CPLE_NotSupported,
1702 1 : "EHdr driver does not support %d bands.\n", nBands);
1703 1 : return NULL;
1704 : }
1705 :
1706 61 : 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 15 : GDALGetDataTypeName(eType) );
1713 :
1714 15 : return NULL;
1715 : }
1716 :
1717 : /* -------------------------------------------------------------------- */
1718 : /* Try to create the file. */
1719 : /* -------------------------------------------------------------------- */
1720 : VSILFILE *fp;
1721 :
1722 46 : fp = VSIFOpenL( pszFilename, "wb" );
1723 :
1724 46 : if( fp == NULL )
1725 : {
1726 : CPLError( CE_Failure, CPLE_OpenFailed,
1727 : "Attempt to create file `%s' failed.\n",
1728 12 : pszFilename );
1729 12 : 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 34 : VSIFWriteL( (void *) "\0\0", 2, 1, fp );
1737 34 : VSIFCloseL( fp );
1738 :
1739 : /* -------------------------------------------------------------------- */
1740 : /* Create the hdr filename. */
1741 : /* -------------------------------------------------------------------- */
1742 : char *pszHdrFilename;
1743 :
1744 : pszHdrFilename =
1745 34 : CPLStrdup( CPLResetExtension( pszFilename, "hdr" ) );
1746 :
1747 : /* -------------------------------------------------------------------- */
1748 : /* Open the file. */
1749 : /* -------------------------------------------------------------------- */
1750 34 : fp = VSIFOpenL( pszHdrFilename, "wt" );
1751 34 : 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 34 : int nBits = GDALGetDataTypeSize(eType);
1765 :
1766 34 : if( CSLFetchNameValue( papszParmList, "NBITS" ) != NULL )
1767 1 : nBits = atoi(CSLFetchNameValue( papszParmList, "NBITS" ));
1768 :
1769 34 : nRowBytes = (nBits * nXSize + 7) / 8;
1770 :
1771 : /* -------------------------------------------------------------------- */
1772 : /* Check for signed byte. */
1773 : /* -------------------------------------------------------------------- */
1774 34 : const char *pszPixelType = CSLFetchNameValue( papszParmList, "PIXELTYPE" );
1775 34 : if( pszPixelType == NULL )
1776 33 : pszPixelType = "";
1777 :
1778 : /* -------------------------------------------------------------------- */
1779 : /* Write out the raw definition for the dataset as a whole. */
1780 : /* -------------------------------------------------------------------- */
1781 34 : VSIFPrintfL( fp, "BYTEORDER I\n" );
1782 34 : VSIFPrintfL( fp, "LAYOUT BIL\n" );
1783 34 : VSIFPrintfL( fp, "NROWS %d\n", nYSize );
1784 34 : VSIFPrintfL( fp, "NCOLS %d\n", nXSize );
1785 34 : VSIFPrintfL( fp, "NBANDS %d\n", nBands );
1786 34 : VSIFPrintfL( fp, "NBITS %d\n", nBits );
1787 34 : VSIFPrintfL( fp, "BANDROWBYTES %d\n", nRowBytes );
1788 34 : VSIFPrintfL( fp, "TOTALROWBYTES %d\n", nRowBytes * nBands );
1789 :
1790 34 : if( eType == GDT_Float32 )
1791 5 : VSIFPrintfL( fp, "PIXELTYPE FLOAT\n");
1792 36 : else if( eType == GDT_Int16 || eType == GDT_Int32 )
1793 7 : VSIFPrintfL( fp, "PIXELTYPE SIGNEDINT\n");
1794 23 : else if( eType == GDT_Byte && EQUAL(pszPixelType,"SIGNEDBYTE") )
1795 1 : VSIFPrintfL( fp, "PIXELTYPE SIGNEDINT\n");
1796 : else
1797 21 : VSIFPrintfL( fp, "PIXELTYPE UNSIGNEDINT\n");
1798 :
1799 : /* -------------------------------------------------------------------- */
1800 : /* Cleanup */
1801 : /* -------------------------------------------------------------------- */
1802 34 : VSIFCloseL( fp );
1803 :
1804 34 : CPLFree( pszHdrFilename );
1805 :
1806 34 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
1807 : }
1808 :
1809 : /************************************************************************/
1810 : /* CreateCopy() */
1811 : /************************************************************************/
1812 :
1813 35 : GDALDataset *EHdrDataset::CreateCopy( const char * pszFilename,
1814 : GDALDataset * poSrcDS,
1815 : int bStrict, char ** papszOptions,
1816 : GDALProgressFunc pfnProgress,
1817 : void * pProgressData )
1818 :
1819 : {
1820 35 : char **papszAdjustedOptions = CSLDuplicate( papszOptions );
1821 : GDALDataset *poOutDS;
1822 :
1823 35 : int nBands = poSrcDS->GetRasterCount();
1824 35 : if (nBands == 0)
1825 : {
1826 : CPLError( CE_Failure, CPLE_NotSupported,
1827 1 : "EHdr driver does not support source dataset with zero band.\n");
1828 1 : return NULL;
1829 : }
1830 :
1831 : /* -------------------------------------------------------------------- */
1832 : /* Ensure we pass on NBITS and PIXELTYPE structure information. */
1833 : /* -------------------------------------------------------------------- */
1834 68 : if( poSrcDS->GetRasterBand(1)->GetMetadataItem( "NBITS",
1835 34 : "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 68 : if( poSrcDS->GetRasterBand(1)->GetMetadataItem( "PIXELTYPE",
1845 34 : "IMAGE_STRUCTURE" ) !=NULL
1846 : && CSLFetchNameValue( papszOptions, "PIXELTYPE" ) == NULL )
1847 : {
1848 : papszAdjustedOptions =
1849 : CSLSetNameValue( papszAdjustedOptions,
1850 : "PIXELTYPE",
1851 1 : poSrcDS->GetRasterBand(1)->GetMetadataItem("PIXELTYPE","IMAGE_STRUCTURE") );
1852 : }
1853 :
1854 : /* -------------------------------------------------------------------- */
1855 : /* Proceed with normal copying using the default createcopy */
1856 : /* operators. */
1857 : /* -------------------------------------------------------------------- */
1858 34 : GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "EHdr" );
1859 :
1860 : poOutDS = poDriver->DefaultCreateCopy( pszFilename, poSrcDS, bStrict,
1861 : papszAdjustedOptions,
1862 34 : pfnProgress, pProgressData );
1863 34 : CSLDestroy( papszAdjustedOptions );
1864 :
1865 34 : return poOutDS;
1866 : }
1867 :
1868 : /************************************************************************/
1869 : /* GetMinimum() */
1870 : /************************************************************************/
1871 :
1872 3 : double EHdrRasterBand::GetMinimum( int *pbSuccess )
1873 : {
1874 3 : if( pbSuccess != NULL )
1875 3 : *pbSuccess = (minmaxmeanstddev & HAS_MIN_FLAG) != 0;
1876 :
1877 3 : if( minmaxmeanstddev & HAS_MIN_FLAG )
1878 2 : return dfMin;
1879 :
1880 1 : return RawRasterBand::GetMinimum( pbSuccess );
1881 : }
1882 :
1883 : /************************************************************************/
1884 : /* GetMaximum() */
1885 : /************************************************************************/
1886 :
1887 2 : double EHdrRasterBand::GetMaximum( int *pbSuccess )
1888 : {
1889 2 : if( pbSuccess != NULL )
1890 2 : *pbSuccess = (minmaxmeanstddev & HAS_MAX_FLAG) != 0;
1891 :
1892 2 : if( minmaxmeanstddev & HAS_MAX_FLAG )
1893 1 : return dfMax;
1894 :
1895 1 : return RawRasterBand::GetMaximum( pbSuccess );
1896 : }
1897 :
1898 : /************************************************************************/
1899 : /* GetStatistics() */
1900 : /************************************************************************/
1901 :
1902 2 : CPLErr EHdrRasterBand::GetStatistics( int bApproxOK, int bForce, double *pdfMin, double *pdfMax, double *pdfMean, double *pdfStdDev )
1903 : {
1904 2 : if( (minmaxmeanstddev & HAS_ALL_FLAGS) == HAS_ALL_FLAGS)
1905 : {
1906 1 : if ( pdfMin ) *pdfMin = dfMin;
1907 1 : if ( pdfMax ) *pdfMax = dfMax;
1908 1 : if ( pdfMean ) *pdfMean = dfMean;
1909 1 : if ( pdfStdDev ) *pdfStdDev = dfStdDev;
1910 1 : return CE_None;
1911 : }
1912 :
1913 : CPLErr eErr = RawRasterBand::GetStatistics( bApproxOK, bForce,
1914 : &dfMin, &dfMax,
1915 1 : &dfMean, &dfStdDev );
1916 :
1917 1 : if( eErr == CE_None )
1918 : {
1919 1 : EHdrDataset* poEDS = (EHdrDataset *) poDS;
1920 :
1921 1 : minmaxmeanstddev = HAS_ALL_FLAGS;
1922 :
1923 1 : if( poEDS->RewriteSTX() != CE_None )
1924 0 : RawRasterBand::SetStatistics( dfMin, dfMax, dfMean, dfStdDev );
1925 :
1926 1 : if( pdfMin )
1927 1 : *pdfMin = dfMin;
1928 1 : if( pdfMax )
1929 1 : *pdfMax = dfMax;
1930 1 : if( pdfMean )
1931 1 : *pdfMean = dfMean;
1932 1 : if( pdfStdDev )
1933 1 : *pdfStdDev = dfStdDev;
1934 : }
1935 :
1936 1 : return eErr;
1937 : }
1938 :
1939 : /************************************************************************/
1940 : /* SetStatistics() */
1941 : /************************************************************************/
1942 :
1943 1 : CPLErr EHdrRasterBand::SetStatistics( double dfMin, double dfMax, double dfMean, double dfStdDev )
1944 : {
1945 : // avoid churn if nothing is changing.
1946 1 : if( dfMin == this->dfMin
1947 : && dfMax == this->dfMax
1948 : && dfMean == this->dfMean
1949 : && dfStdDev == this->dfStdDev )
1950 0 : return CE_None;
1951 :
1952 1 : this->dfMin = dfMin;
1953 1 : this->dfMax = dfMax;
1954 1 : this->dfMean = dfMean;
1955 1 : this->dfStdDev = dfStdDev;
1956 :
1957 : // marks stats valid
1958 1 : minmaxmeanstddev = HAS_ALL_FLAGS;
1959 :
1960 1 : EHdrDataset* poEDS = (EHdrDataset *) poDS;
1961 :
1962 1 : if( poEDS->RewriteSTX() != CE_None )
1963 0 : return RawRasterBand::SetStatistics( dfMin, dfMax, dfMean, dfStdDev );
1964 : else
1965 1 : return CE_None;
1966 : }
1967 :
1968 : /************************************************************************/
1969 : /* SetColorTable() */
1970 : /************************************************************************/
1971 :
1972 6 : CPLErr EHdrRasterBand::SetColorTable( GDALColorTable *poNewCT )
1973 : {
1974 6 : CPLErr err = RawRasterBand::SetColorTable( poNewCT );
1975 6 : if( err != CE_None )
1976 0 : return err;
1977 :
1978 6 : ((EHdrDataset*)poDS)->bCLRDirty = TRUE;
1979 :
1980 6 : return CE_None;
1981 : }
1982 :
1983 : /************************************************************************/
1984 : /* GDALRegister_EHdr() */
1985 : /************************************************************************/
1986 :
1987 582 : void GDALRegister_EHdr()
1988 :
1989 : {
1990 : GDALDriver *poDriver;
1991 :
1992 582 : if( GDALGetDriverByName( "EHdr" ) == NULL )
1993 : {
1994 561 : poDriver = new GDALDriver();
1995 :
1996 561 : poDriver->SetDescription( "EHdr" );
1997 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1998 561 : "ESRI .hdr Labelled" );
1999 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
2000 561 : "frmt_various.html#EHdr" );
2001 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2002 561 : "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 561 : "</CreationOptionList>" );
2009 :
2010 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2011 561 : poDriver->pfnOpen = EHdrDataset::Open;
2012 561 : poDriver->pfnCreate = EHdrDataset::Create;
2013 561 : poDriver->pfnCreateCopy = EHdrDataset::CreateCopy;
2014 :
2015 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
2016 : }
2017 582 : }
|
