1 : /******************************************************************************
2 : * $Id: aaigriddataset.cpp 17974 2009-11-07 16:46:40Z chaitanya $
3 : *
4 : * Project: GDAL
5 : * Purpose: Implements Arc/Info ASCII Grid Format.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2001, 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 "gdal_pam.h"
31 : #include <ctype.h>
32 : #include <limits.h>
33 : #include "cpl_string.h"
34 : #include "ogr_spatialref.h"
35 :
36 : CPL_CVSID("$Id: aaigriddataset.cpp 17974 2009-11-07 16:46:40Z chaitanya $");
37 :
38 : CPL_C_START
39 : void GDALRegister_AAIGrid(void);
40 : CPL_C_END
41 :
42 : static CPLString OSR_GDS( char **papszNV, const char * pszField,
43 : const char *pszDefaultValue );
44 :
45 : /************************************************************************/
46 : /* ==================================================================== */
47 : /* AAIGDataset */
48 : /* ==================================================================== */
49 : /************************************************************************/
50 :
51 : class AAIGRasterBand;
52 :
53 : class CPL_DLL AAIGDataset : public GDALPamDataset
54 : {
55 : friend class AAIGRasterBand;
56 :
57 : FILE *fp;
58 :
59 : double adfGeoTransform[6];
60 : char **papszPrj;
61 : CPLString osPrjFilename;
62 : char *pszProjection;
63 :
64 : int bNoDataSet;
65 : double dfNoDataValue;
66 :
67 : unsigned char achReadBuf[256];
68 : GUIntBig nBufferOffset;
69 : int nOffsetInBuffer;
70 :
71 : char Getc();
72 : GUIntBig Tell();
73 : int Seek( GUIntBig nOffset );
74 :
75 : public:
76 : AAIGDataset();
77 : ~AAIGDataset();
78 :
79 : virtual char **GetFileList(void);
80 :
81 : static GDALDataset *Open( GDALOpenInfo * );
82 : static CPLErr Delete( const char *pszFilename );
83 : static CPLErr Remove( const char *pszFilename, int bRepError );
84 :
85 : virtual CPLErr GetGeoTransform( double * );
86 : virtual const char *GetProjectionRef(void);
87 : };
88 :
89 : /************************************************************************/
90 : /* ==================================================================== */
91 : /* AAIGRasterBand */
92 : /* ==================================================================== */
93 : /************************************************************************/
94 :
95 : class AAIGRasterBand : public GDALPamRasterBand
96 : {
97 : friend class AAIGDataset;
98 :
99 : GUIntBig *panLineOffset;
100 :
101 : public:
102 :
103 : AAIGRasterBand( AAIGDataset *, int, GDALDataType );
104 : virtual ~AAIGRasterBand();
105 :
106 : virtual double GetNoDataValue( int * );
107 : virtual CPLErr SetNoDataValue( double );
108 : virtual CPLErr IReadBlock( int, int, void * );
109 : };
110 :
111 : /************************************************************************/
112 : /* AAIGRasterBand() */
113 : /************************************************************************/
114 :
115 54 : AAIGRasterBand::AAIGRasterBand( AAIGDataset *poDS, int nDataStart,
116 54 : GDALDataType eTypeIn )
117 :
118 : {
119 54 : this->poDS = poDS;
120 :
121 54 : nBand = 1;
122 54 : eDataType = eTypeIn;
123 :
124 54 : nBlockXSize = poDS->nRasterXSize;
125 54 : nBlockYSize = 1;
126 :
127 : panLineOffset =
128 54 : (GUIntBig *) VSICalloc( poDS->nRasterYSize, sizeof(GUIntBig) );
129 54 : if (panLineOffset == NULL)
130 : {
131 : CPLError(CE_Failure, CPLE_OutOfMemory,
132 : "AAIGRasterBand::AAIGRasterBand : Out of memory (nRasterYSize = %d)",
133 0 : poDS->nRasterYSize);
134 0 : return;
135 : }
136 54 : panLineOffset[0] = nDataStart;
137 0 : }
138 :
139 : /************************************************************************/
140 : /* ~AAIGRasterBand() */
141 : /************************************************************************/
142 :
143 108 : AAIGRasterBand::~AAIGRasterBand()
144 :
145 : {
146 54 : CPLFree( panLineOffset );
147 108 : }
148 :
149 : /************************************************************************/
150 : /* IReadBlock() */
151 : /************************************************************************/
152 :
153 688 : CPLErr AAIGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
154 : void * pImage )
155 :
156 : {
157 688 : AAIGDataset *poODS = (AAIGDataset *) poDS;
158 : int iPixel;
159 :
160 688 : if( nBlockYOff < 0 || nBlockYOff > poODS->nRasterYSize - 1
161 : || nBlockXOff != 0 || panLineOffset == NULL)
162 0 : return CE_Failure;
163 :
164 688 : if( panLineOffset[nBlockYOff] == 0 )
165 : {
166 : int iPrevLine;
167 :
168 6 : for( iPrevLine = 1; iPrevLine <= nBlockYOff; iPrevLine++ )
169 5 : if( panLineOffset[iPrevLine] == 0 )
170 5 : IReadBlock( nBlockXOff, iPrevLine-1, NULL );
171 : }
172 :
173 688 : if( panLineOffset[nBlockYOff] == 0 )
174 0 : return CE_Failure;
175 :
176 :
177 688 : if( poODS->Seek( panLineOffset[nBlockYOff] ) != 0 )
178 : {
179 : CPLError( CE_Failure, CPLE_FileIO,
180 : "Can't seek to offset %lu in input file to read data.",
181 0 : (long unsigned int)panLineOffset[nBlockYOff] );
182 0 : return CE_Failure;
183 : }
184 :
185 31518 : for( iPixel = 0; iPixel < poODS->nRasterXSize; )
186 : {
187 : char szToken[500];
188 : char chNext;
189 30142 : int iTokenChar = 0;
190 :
191 : /* suck up any pre-white space. */
192 38200 : do {
193 38200 : chNext = poODS->Getc();
194 : } while( isspace( (unsigned char)chNext ) );
195 :
196 148633 : while( chNext != '\0' && !isspace((unsigned char)chNext) )
197 : {
198 88349 : if( iTokenChar == sizeof(szToken)-2 )
199 : {
200 : CPLError( CE_Failure, CPLE_FileIO,
201 : "Token too long at scanline %d.",
202 0 : nBlockYOff );
203 0 : return CE_Failure;
204 : }
205 :
206 88349 : szToken[iTokenChar++] = chNext;
207 88349 : chNext = poODS->Getc();
208 : }
209 :
210 30142 : if( chNext == '\0' &&
211 : (iPixel != poODS->nRasterXSize - 1 ||
212 : nBlockYOff != poODS->nRasterYSize - 1) )
213 : {
214 : CPLError( CE_Failure, CPLE_FileIO,
215 : "File short, can't read line %d.",
216 0 : nBlockYOff );
217 0 : return CE_Failure;
218 : }
219 :
220 30142 : szToken[iTokenChar] = '\0';
221 :
222 30142 : if( pImage != NULL )
223 : {
224 30067 : if( eDataType == GDT_Float32 )
225 835 : ((float *) pImage)[iPixel] = (float) atof(szToken);
226 : else
227 29232 : ((GInt32 *) pImage)[iPixel] = (GInt32) atoi(szToken);
228 : }
229 :
230 30142 : iPixel++;
231 : }
232 :
233 688 : if( nBlockYOff < poODS->nRasterYSize - 1 )
234 658 : panLineOffset[nBlockYOff + 1] = poODS->Tell();
235 :
236 688 : return CE_None;
237 : }
238 :
239 : /************************************************************************/
240 : /* GetNoDataValue() */
241 : /************************************************************************/
242 :
243 22 : double AAIGRasterBand::GetNoDataValue( int * pbSuccess )
244 :
245 : {
246 22 : AAIGDataset *poODS = (AAIGDataset *) poDS;
247 :
248 22 : if( pbSuccess )
249 21 : *pbSuccess = poODS->bNoDataSet;
250 :
251 22 : return poODS->dfNoDataValue;
252 : }
253 :
254 : /************************************************************************/
255 : /* SetNoDataValue() */
256 : /************************************************************************/
257 :
258 0 : CPLErr AAIGRasterBand::SetNoDataValue( double dfNoData )
259 :
260 : {
261 0 : AAIGDataset *poODS = (AAIGDataset *) poDS;
262 :
263 0 : poODS->bNoDataSet = TRUE;
264 0 : poODS->dfNoDataValue = dfNoData;
265 :
266 0 : return CE_None;
267 : }
268 :
269 : /************************************************************************/
270 : /* ==================================================================== */
271 : /* AAIGDataset */
272 : /* ==================================================================== */
273 : /************************************************************************/
274 :
275 :
276 : /************************************************************************/
277 : /* AAIGDataset() */
278 : /************************************************************************/
279 :
280 54 : AAIGDataset::AAIGDataset()
281 :
282 : {
283 54 : papszPrj = NULL;
284 54 : pszProjection = CPLStrdup("");
285 54 : fp = NULL;
286 54 : bNoDataSet = FALSE;
287 54 : dfNoDataValue = -9999.0;
288 :
289 54 : adfGeoTransform[0] = 0.0;
290 54 : adfGeoTransform[1] = 1.0;
291 54 : adfGeoTransform[2] = 0.0;
292 54 : adfGeoTransform[3] = 0.0;
293 54 : adfGeoTransform[4] = 0.0;
294 54 : adfGeoTransform[5] = 1.0;
295 :
296 54 : nOffsetInBuffer = 256;
297 54 : nBufferOffset = 0;
298 54 : }
299 :
300 : /************************************************************************/
301 : /* ~AAIGDataset() */
302 : /************************************************************************/
303 :
304 108 : AAIGDataset::~AAIGDataset()
305 :
306 : {
307 54 : FlushCache();
308 :
309 54 : if( fp != NULL )
310 54 : VSIFCloseL( fp );
311 :
312 54 : CPLFree( pszProjection );
313 54 : CSLDestroy( papszPrj );
314 108 : }
315 :
316 : /************************************************************************/
317 : /* Tell() */
318 : /************************************************************************/
319 :
320 658 : GUIntBig AAIGDataset::Tell()
321 :
322 : {
323 658 : return nBufferOffset + nOffsetInBuffer;
324 : }
325 :
326 : /************************************************************************/
327 : /* Seek() */
328 : /************************************************************************/
329 :
330 688 : int AAIGDataset::Seek( GUIntBig nNewOffset )
331 :
332 : {
333 688 : nOffsetInBuffer = sizeof(achReadBuf);
334 688 : return VSIFSeekL( fp, nNewOffset, SEEK_SET );
335 : }
336 :
337 : /************************************************************************/
338 : /* Getc() */
339 : /* */
340 : /* Read a single character from the input file (efficiently we */
341 : /* hope). */
342 : /************************************************************************/
343 :
344 126549 : char AAIGDataset::Getc()
345 :
346 : {
347 126549 : if( nOffsetInBuffer < (int) sizeof(achReadBuf) )
348 125661 : return achReadBuf[nOffsetInBuffer++];
349 :
350 888 : nBufferOffset = VSIFTellL( fp );
351 888 : int nRead = VSIFReadL( achReadBuf, 1, sizeof(achReadBuf), fp );
352 : unsigned int i;
353 15206 : for(i=nRead;i<sizeof(achReadBuf);i++)
354 14318 : achReadBuf[i] = '\0';
355 :
356 888 : nOffsetInBuffer = 0;
357 :
358 888 : return achReadBuf[nOffsetInBuffer++];
359 : }
360 :
361 : /************************************************************************/
362 : /* GetFileList() */
363 : /************************************************************************/
364 :
365 9 : char **AAIGDataset::GetFileList()
366 :
367 : {
368 9 : char **papszFileList = GDALPamDataset::GetFileList();
369 :
370 9 : if( papszPrj != NULL )
371 8 : papszFileList = CSLAddString( papszFileList, osPrjFilename );
372 :
373 9 : return papszFileList;
374 : }
375 :
376 : /************************************************************************/
377 : /* Open() */
378 : /************************************************************************/
379 :
380 9915 : GDALDataset *AAIGDataset::Open( GDALOpenInfo * poOpenInfo )
381 :
382 : {
383 9915 : int i = 0;
384 9915 : int j = 0;
385 9915 : char **papszTokens = NULL;
386 :
387 : /* Default data type */
388 9915 : GDALDataType eDataType = GDT_Int32;
389 :
390 : /* -------------------------------------------------------------------- */
391 : /* Does this look like an AI grid file? */
392 : /* -------------------------------------------------------------------- */
393 9915 : if( poOpenInfo->nHeaderBytes < 100
394 : || !( EQUALN((const char *) poOpenInfo->pabyHeader,"ncols",5) ||
395 : EQUALN((const char *) poOpenInfo->pabyHeader,"nrows",5) ||
396 : EQUALN((const char *) poOpenInfo->pabyHeader,"xllcorner",9)||
397 : EQUALN((const char *) poOpenInfo->pabyHeader,"yllcorner",9)||
398 : EQUALN((const char *) poOpenInfo->pabyHeader,"xllcenter",9)||
399 : EQUALN((const char *) poOpenInfo->pabyHeader,"yllcenter",9)||
400 : EQUALN((const char *) poOpenInfo->pabyHeader,"dx",2)||
401 : EQUALN((const char *) poOpenInfo->pabyHeader,"dy",2)||
402 : EQUALN((const char *) poOpenInfo->pabyHeader,"cellsize",8)) )
403 9861 : return NULL;
404 :
405 : papszTokens =
406 : CSLTokenizeString2( (const char *) poOpenInfo->pabyHeader,
407 54 : " \n\r\t", 0 );
408 54 : int nTokens = CSLCount(papszTokens);
409 :
410 : /* -------------------------------------------------------------------- */
411 : /* Create a corresponding GDALDataset. */
412 : /* -------------------------------------------------------------------- */
413 : AAIGDataset *poDS;
414 :
415 54 : poDS = new AAIGDataset();
416 :
417 : /* -------------------------------------------------------------------- */
418 : /* Parse the header. */
419 : /* -------------------------------------------------------------------- */
420 54 : double dfCellDX = 0;
421 54 : double dfCellDY = 0;
422 :
423 108 : if ( (i = CSLFindString( papszTokens, "ncols" )) < 0 ||
424 : i + 1 >= nTokens)
425 : {
426 0 : CSLDestroy( papszTokens );
427 0 : delete poDS;
428 0 : return NULL;
429 : }
430 54 : poDS->nRasterXSize = atoi(papszTokens[i + 1]);
431 54 : if ( (i = CSLFindString( papszTokens, "nrows" )) < 0 ||
432 : i + 1 >= nTokens)
433 : {
434 0 : CSLDestroy( papszTokens );
435 0 : delete poDS;
436 0 : return NULL;
437 : }
438 54 : poDS->nRasterYSize = atoi(papszTokens[i + 1]);
439 :
440 54 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
441 : {
442 0 : delete poDS;
443 0 : return NULL;
444 : }
445 :
446 54 : if ( (i = CSLFindString( papszTokens, "cellsize" )) < 0 )
447 : {
448 : int iDX, iDY;
449 3 : if( (iDX = CSLFindString(papszTokens,"dx")) < 0
450 : || (iDY = CSLFindString(papszTokens,"dy")) < 0
451 : || iDX+1 >= nTokens
452 : || iDY+1 >= nTokens)
453 : {
454 0 : CSLDestroy( papszTokens );
455 0 : delete poDS;
456 0 : return NULL;
457 : }
458 :
459 3 : dfCellDX = atof( papszTokens[iDX+1] );
460 3 : dfCellDY = atof( papszTokens[iDY+1] );
461 : }
462 : else
463 : {
464 51 : if (i + 1 >= nTokens)
465 : {
466 0 : CSLDestroy( papszTokens );
467 0 : delete poDS;
468 0 : return NULL;
469 : }
470 51 : dfCellDX = dfCellDY = atof( papszTokens[i + 1] );
471 : }
472 :
473 54 : if ((i = CSLFindString( papszTokens, "xllcorner" )) >= 0 &&
474 : (j = CSLFindString( papszTokens, "yllcorner" )) >= 0 &&
475 : i + 1 < nTokens && j + 1 < nTokens)
476 : {
477 54 : poDS->adfGeoTransform[0] = atof( papszTokens[i + 1] );
478 54 : poDS->adfGeoTransform[1] = dfCellDX;
479 54 : poDS->adfGeoTransform[2] = 0.0;
480 54 : poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] )
481 108 : + poDS->nRasterYSize * dfCellDY;
482 54 : poDS->adfGeoTransform[4] = 0.0;
483 54 : poDS->adfGeoTransform[5] = - dfCellDY;
484 : }
485 0 : else if ((i = CSLFindString( papszTokens, "xllcenter" )) >= 0 &&
486 : (j = CSLFindString( papszTokens, "yllcenter" )) >= 0 &&
487 : i + 1 < nTokens && j + 1 < nTokens)
488 : {
489 0 : poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
490 :
491 0 : poDS->adfGeoTransform[0] = atof(papszTokens[i + 1]) - 0.5 * dfCellDX;
492 0 : poDS->adfGeoTransform[1] = dfCellDX;
493 0 : poDS->adfGeoTransform[2] = 0.0;
494 0 : poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] )
495 : - 0.5 * dfCellDY
496 0 : + poDS->nRasterYSize * dfCellDY;
497 0 : poDS->adfGeoTransform[4] = 0.0;
498 0 : poDS->adfGeoTransform[5] = - dfCellDY;
499 : }
500 : else
501 : {
502 0 : poDS->adfGeoTransform[0] = 0.0;
503 0 : poDS->adfGeoTransform[1] = dfCellDX;
504 0 : poDS->adfGeoTransform[2] = 0.0;
505 0 : poDS->adfGeoTransform[3] = 0.0;
506 0 : poDS->adfGeoTransform[4] = 0.0;
507 0 : poDS->adfGeoTransform[5] = - dfCellDY;
508 : }
509 :
510 54 : if( (i = CSLFindString( papszTokens, "NODATA_value" )) >= 0 &&
511 : i + 1 < nTokens)
512 : {
513 12 : const char* pszNoData = papszTokens[i + 1];
514 :
515 12 : poDS->bNoDataSet = TRUE;
516 12 : poDS->dfNoDataValue = atof(pszNoData);
517 12 : if( strchr( pszNoData, '.' ) != NULL )
518 : {
519 1 : eDataType = GDT_Float32;
520 : }
521 : }
522 :
523 54 : CSLDestroy( papszTokens );
524 :
525 : /* -------------------------------------------------------------------- */
526 : /* Open file with large file API. */
527 : /* -------------------------------------------------------------------- */
528 :
529 54 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
530 54 : if( poDS->fp == NULL )
531 : {
532 : CPLError( CE_Failure, CPLE_OpenFailed,
533 : "VSIFOpenL(%s) failed unexpectedly.",
534 0 : poOpenInfo->pszFilename );
535 0 : delete poDS;
536 0 : return NULL;
537 : }
538 :
539 : /* -------------------------------------------------------------------- */
540 : /* Find the start of real data. */
541 : /* -------------------------------------------------------------------- */
542 : int nStartOfData;
543 :
544 6891 : for( i = 2; TRUE ; i++ )
545 : {
546 6891 : if( poOpenInfo->pabyHeader[i] == '\0' )
547 : {
548 : CPLError( CE_Failure, CPLE_AppDefined,
549 0 : "Couldn't find data values in ASCII Grid file.\n" );
550 0 : delete poDS;
551 0 : return NULL;
552 : }
553 :
554 26223 : if( poOpenInfo->pabyHeader[i-1] == '\n'
555 6606 : || poOpenInfo->pabyHeader[i-2] == '\n'
556 6375 : || poOpenInfo->pabyHeader[i-1] == '\r'
557 6351 : || poOpenInfo->pabyHeader[i-2] == '\r' )
558 : {
559 672 : if( !isalpha(poOpenInfo->pabyHeader[i])
560 78 : && poOpenInfo->pabyHeader[i] != '\n'
561 54 : && poOpenInfo->pabyHeader[i] != '\r')
562 : {
563 54 : nStartOfData = i;
564 :
565 : /* Beginning of real data found. */
566 : break;
567 : }
568 : }
569 : }
570 :
571 : /* -------------------------------------------------------------------- */
572 : /* Recognize the type of data. */
573 : /* -------------------------------------------------------------------- */
574 : CPLAssert( NULL != poDS->fp );
575 :
576 : /* Use bigger data type. */
577 66 : if( poDS->bNoDataSet
578 : && ( INT_MIN > poDS->dfNoDataValue || poDS->dfNoDataValue > INT_MAX) )
579 : {
580 0 : eDataType = GDT_Float32;
581 : }
582 : else
583 : {
584 : /* Allocate 100K chunk + 1 extra byte for NULL character. */
585 54 : const size_t nChunkSize = 1024 * 100;
586 54 : GByte* pabyChunk = (GByte *) CPLCalloc( nChunkSize + 1, sizeof(GByte) );
587 54 : pabyChunk[nChunkSize] = '\0';
588 :
589 54 : VSIFSeekL( poDS->fp, nStartOfData, SEEK_SET );
590 :
591 : /* Scan for dot in subsequent chunks of data. */
592 54 : while( !VSIFEofL( poDS->fp) )
593 : {
594 54 : VSIFReadL( pabyChunk, sizeof(GByte), nChunkSize, poDS->fp );
595 : CPLAssert( pabyChunk[nChunkSize] == '\0' );
596 :
597 54 : if( strchr( (const char *)pabyChunk, '.' ) != NULL )
598 : {
599 6 : eDataType = GDT_Float32;
600 6 : break;
601 : }
602 : }
603 :
604 : /* Deallocate chunk. */
605 54 : VSIFree( pabyChunk );
606 : }
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* Create band information objects. */
610 : /* -------------------------------------------------------------------- */
611 54 : AAIGRasterBand* band = new AAIGRasterBand( poDS, nStartOfData, eDataType );
612 54 : poDS->SetBand( 1, band );
613 54 : if (band->panLineOffset == NULL)
614 : {
615 0 : delete poDS;
616 0 : return NULL;
617 : }
618 :
619 : /* -------------------------------------------------------------------- */
620 : /* Try to read projection file. */
621 : /* -------------------------------------------------------------------- */
622 : char *pszDirname, *pszBasename;
623 : VSIStatBufL sStatBuf;
624 :
625 54 : pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
626 54 : pszBasename = CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
627 :
628 54 : poDS->osPrjFilename = CPLFormFilename( pszDirname, pszBasename, "prj" );
629 54 : int nRet = VSIStatL( poDS->osPrjFilename, &sStatBuf );
630 :
631 : #ifndef WIN32
632 54 : if( nRet != 0 )
633 : {
634 17 : poDS->osPrjFilename = CPLFormFilename( pszDirname, pszBasename, "PRJ" );
635 17 : nRet = VSIStatL( poDS->osPrjFilename, &sStatBuf );
636 : }
637 : #endif
638 :
639 54 : if( nRet == 0 )
640 : {
641 38 : OGRSpatialReference oSRS;
642 :
643 38 : poDS->papszPrj = CSLLoad( poDS->osPrjFilename );
644 :
645 : CPLDebug( "AAIGrid", "Loaded SRS from %s",
646 38 : poDS->osPrjFilename.c_str() );
647 :
648 38 : if( oSRS.importFromESRI( poDS->papszPrj ) == OGRERR_NONE )
649 : {
650 : // If geographic values are in seconds, we must transform.
651 : // Is there a code for minutes too?
652 38 : if( oSRS.IsGeographic()
653 : && EQUAL(OSR_GDS( poDS->papszPrj, "Units", ""), "DS") )
654 : {
655 0 : poDS->adfGeoTransform[0] /= 3600.0;
656 0 : poDS->adfGeoTransform[1] /= 3600.0;
657 0 : poDS->adfGeoTransform[2] /= 3600.0;
658 0 : poDS->adfGeoTransform[3] /= 3600.0;
659 0 : poDS->adfGeoTransform[4] /= 3600.0;
660 0 : poDS->adfGeoTransform[5] /= 3600.0;
661 : }
662 :
663 38 : CPLFree( poDS->pszProjection );
664 38 : oSRS.exportToWkt( &(poDS->pszProjection) );
665 38 : }
666 : }
667 :
668 54 : CPLFree( pszDirname );
669 54 : CPLFree( pszBasename );
670 :
671 : /* -------------------------------------------------------------------- */
672 : /* Initialize any PAM information. */
673 : /* -------------------------------------------------------------------- */
674 54 : poDS->SetDescription( poOpenInfo->pszFilename );
675 54 : poDS->TryLoadXML();
676 :
677 54 : return( poDS );
678 : }
679 :
680 : /************************************************************************/
681 : /* GetGeoTransform() */
682 : /************************************************************************/
683 :
684 26 : CPLErr AAIGDataset::GetGeoTransform( double * padfTransform )
685 :
686 : {
687 26 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
688 26 : return( CE_None );
689 : }
690 :
691 : /************************************************************************/
692 : /* GetProjectionRef() */
693 : /************************************************************************/
694 :
695 41 : const char *AAIGDataset::GetProjectionRef()
696 :
697 : {
698 41 : return pszProjection;
699 : }
700 :
701 : /************************************************************************/
702 : /* AAIGCreateCopy() */
703 : /************************************************************************/
704 :
705 : static GDALDataset *
706 25 : AAIGCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
707 : int bStrict, char ** papszOptions,
708 : GDALProgressFunc pfnProgress, void * pProgressData )
709 :
710 : {
711 25 : int nBands = poSrcDS->GetRasterCount();
712 25 : int nXSize = poSrcDS->GetRasterXSize();
713 25 : int nYSize = poSrcDS->GetRasterYSize();
714 :
715 : /* -------------------------------------------------------------------- */
716 : /* Some rudimentary checks */
717 : /* -------------------------------------------------------------------- */
718 25 : if( nBands != 1 )
719 : {
720 : CPLError( CE_Failure, CPLE_NotSupported,
721 : "AAIG driver doesn't support %d bands. Must be 1 band.\n",
722 5 : nBands );
723 :
724 5 : return NULL;
725 : }
726 :
727 20 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
728 0 : return NULL;
729 :
730 : /* -------------------------------------------------------------------- */
731 : /* Create the dataset. */
732 : /* -------------------------------------------------------------------- */
733 : FILE *fpImage;
734 :
735 20 : fpImage = VSIFOpenL( pszFilename, "wt" );
736 20 : if( fpImage == NULL )
737 : {
738 : CPLError( CE_Failure, CPLE_OpenFailed,
739 : "Unable to create file %s.\n",
740 0 : pszFilename );
741 0 : return NULL;
742 : }
743 :
744 : /* -------------------------------------------------------------------- */
745 : /* Write ASCII Grid file header */
746 : /* -------------------------------------------------------------------- */
747 : double adfGeoTransform[6];
748 : char szHeader[2000];
749 : const char *pszForceCellsize =
750 20 : CSLFetchNameValue( papszOptions, "FORCE_CELLSIZE" );
751 :
752 20 : poSrcDS->GetGeoTransform( adfGeoTransform );
753 :
754 22 : if( ABS(adfGeoTransform[1]+adfGeoTransform[5]) < 0.0000001
755 2 : || ABS(adfGeoTransform[1]-adfGeoTransform[5]) < 0.0000001
756 : || (pszForceCellsize && CSLTestBoolean(pszForceCellsize)) )
757 : sprintf( szHeader,
758 : "ncols %d\n"
759 : "nrows %d\n"
760 : "xllcorner %.12f\n"
761 : "yllcorner %.12f\n"
762 : "cellsize %.12f\n",
763 : nXSize, nYSize,
764 : adfGeoTransform[0],
765 38 : adfGeoTransform[3] - nYSize * adfGeoTransform[1],
766 57 : adfGeoTransform[1] );
767 : else
768 : {
769 1 : if( pszForceCellsize == NULL )
770 : CPLError( CE_Warning, CPLE_AppDefined,
771 : "Producing a Golden Surfer style file with DX and DY instead\n"
772 : "of CELLSIZE since the input pixels are non-square. Use the\n"
773 : "FORCE_CELLSIZE=TRUE creation option to force use of DX for\n"
774 : "even though this will be distorted. Most ASCII Grid readers\n"
775 1 : "(ArcGIS included) do not support the DX and DY parameters.\n" );
776 : sprintf( szHeader,
777 : "ncols %d\n"
778 : "nrows %d\n"
779 : "xllcorner %.12f\n"
780 : "yllcorner %.12f\n"
781 : "dx %.12f\n"
782 : "dy %.12f\n",
783 : nXSize, nYSize,
784 : adfGeoTransform[0],
785 2 : adfGeoTransform[3] + nYSize * adfGeoTransform[5],
786 : adfGeoTransform[1],
787 3 : fabs(adfGeoTransform[5]) );
788 : }
789 :
790 : /* -------------------------------------------------------------------- */
791 : /* Handle nodata (optionally). */
792 : /* -------------------------------------------------------------------- */
793 20 : GDALRasterBand * poBand = poSrcDS->GetRasterBand( 1 );
794 : double dfNoData;
795 : int bSuccess;
796 :
797 : // Write `nodata' value to header if it is exists in source dataset
798 20 : dfNoData = poBand->GetNoDataValue( &bSuccess );
799 20 : if ( bSuccess )
800 : sprintf( szHeader+strlen(szHeader), "NODATA_value %6.20g\n",
801 0 : dfNoData );
802 :
803 20 : VSIFWriteL( szHeader, 1, strlen(szHeader), fpImage );
804 :
805 : /* -------------------------------------------------------------------- */
806 : /* Builds the format string used for printing float values. */
807 : /* -------------------------------------------------------------------- */
808 : char szFormatFloat[32];
809 20 : strcpy(szFormatFloat, " %6.20g");
810 : const char *pszDecimalPrecision =
811 20 : CSLFetchNameValue( papszOptions, "DECIMAL_PRECISION" );
812 20 : if (pszDecimalPrecision)
813 : {
814 1 : int nDecimal = atoi(pszDecimalPrecision);
815 1 : if (nDecimal >= 0)
816 1 : sprintf(szFormatFloat, " %%.%df", nDecimal);
817 : }
818 :
819 : /* -------------------------------------------------------------------- */
820 : /* Loop over image, copying image data. */
821 : /* -------------------------------------------------------------------- */
822 20 : int *panScanline = NULL;
823 20 : double *padfScanline = NULL;
824 : int bReadAsInt;
825 : int iLine, iPixel;
826 20 : CPLErr eErr = CE_None;
827 :
828 : bReadAsInt = ( poBand->GetRasterDataType() == GDT_Byte
829 : || poBand->GetRasterDataType() == GDT_Int16
830 : || poBand->GetRasterDataType() == GDT_UInt16
831 20 : || poBand->GetRasterDataType() == GDT_Int32 );
832 :
833 : // Write scanlines to output file
834 20 : if (bReadAsInt)
835 : panScanline = (int *) CPLMalloc( nXSize *
836 11 : GDALGetDataTypeSize(GDT_Int32) / 8 );
837 : else
838 : padfScanline = (double *) CPLMalloc( nXSize *
839 9 : GDALGetDataTypeSize(GDT_Float64) / 8 );
840 :
841 405 : for( iLine = 0; eErr == CE_None && iLine < nYSize; iLine++ )
842 : {
843 : eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
844 : (bReadAsInt) ? (void*)panScanline : (void*)padfScanline,
845 : nXSize, 1, (bReadAsInt) ? GDT_Int32 : GDT_Float64,
846 385 : 0, 0 );
847 :
848 385 : if( bReadAsInt )
849 : {
850 13760 : for ( iPixel = 0; iPixel < nXSize; iPixel++ )
851 : {
852 13485 : sprintf( szHeader, " %d", panScanline[iPixel] );
853 13485 : if( VSIFWriteL( szHeader, strlen(szHeader), 1, fpImage ) != 1 )
854 : {
855 0 : eErr = CE_Failure;
856 : CPLError( CE_Failure, CPLE_AppDefined,
857 0 : "Write failed, disk full?\n" );
858 0 : break;
859 : }
860 : }
861 : }
862 : else
863 : {
864 1610 : for ( iPixel = 0; iPixel < nXSize; iPixel++ )
865 : {
866 1500 : sprintf( szHeader, szFormatFloat, padfScanline[iPixel] );
867 1500 : if( VSIFWriteL( szHeader, strlen(szHeader), 1, fpImage ) != 1 )
868 : {
869 0 : eErr = CE_Failure;
870 : CPLError( CE_Failure, CPLE_AppDefined,
871 0 : "Write failed, disk full?\n" );
872 0 : break;
873 : }
874 : }
875 : }
876 385 : VSIFWriteL( (void *) "\n", 1, 1, fpImage );
877 :
878 385 : if( eErr == CE_None &&
879 : !pfnProgress((iLine + 1) / ((double) nYSize), NULL, pProgressData) )
880 : {
881 0 : eErr = CE_Failure;
882 : CPLError( CE_Failure, CPLE_UserInterrupt,
883 0 : "User terminated CreateCopy()" );
884 : }
885 : }
886 :
887 20 : CPLFree( panScanline );
888 20 : CPLFree( padfScanline );
889 20 : VSIFCloseL( fpImage );
890 :
891 : /* -------------------------------------------------------------------- */
892 : /* Try to write projection file. */
893 : /* -------------------------------------------------------------------- */
894 : const char *pszOriginalProjection;
895 :
896 20 : pszOriginalProjection = (char *)poSrcDS->GetProjectionRef();
897 20 : if( !EQUAL( pszOriginalProjection, "" ) )
898 : {
899 : char *pszDirname, *pszBasename;
900 : char *pszPrjFilename;
901 19 : char *pszESRIProjection = NULL;
902 : FILE *fp;
903 19 : OGRSpatialReference oSRS;
904 :
905 19 : pszDirname = CPLStrdup( CPLGetPath(pszFilename) );
906 19 : pszBasename = CPLStrdup( CPLGetBasename(pszFilename) );
907 :
908 19 : pszPrjFilename = CPLStrdup( CPLFormFilename( pszDirname, pszBasename, "prj" ) );
909 19 : fp = VSIFOpenL( pszPrjFilename, "wt" );
910 19 : if (fp != NULL)
911 : {
912 19 : oSRS.importFromWkt( (char **) &pszOriginalProjection );
913 19 : oSRS.morphToESRI();
914 19 : oSRS.exportToWkt( &pszESRIProjection );
915 19 : VSIFWriteL( pszESRIProjection, 1, strlen(pszESRIProjection), fp );
916 :
917 19 : VSIFCloseL( fp );
918 19 : CPLFree( pszESRIProjection );
919 : }
920 : else
921 : {
922 : CPLError( CE_Failure, CPLE_FileIO,
923 0 : "Unable to create file %s.\n", pszPrjFilename );
924 : }
925 19 : CPLFree( pszDirname );
926 19 : CPLFree( pszBasename );
927 19 : CPLFree( pszPrjFilename );
928 : }
929 :
930 : /* -------------------------------------------------------------------- */
931 : /* Re-open dataset, and copy any auxilary pam information. */
932 : /* -------------------------------------------------------------------- */
933 : GDALPamDataset *poDS = (GDALPamDataset *)
934 20 : GDALOpen( pszFilename, GA_ReadOnly );
935 :
936 20 : if( poDS )
937 20 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
938 :
939 20 : return poDS;
940 : }
941 :
942 : /************************************************************************/
943 : /* OSR_GDS() */
944 : /************************************************************************/
945 :
946 11 : static CPLString OSR_GDS( char **papszNV, const char * pszField,
947 : const char *pszDefaultValue )
948 :
949 : {
950 : int iLine;
951 :
952 11 : if( papszNV == NULL || papszNV[0] == NULL )
953 0 : return pszDefaultValue;
954 :
955 44 : for( iLine = 0;
956 22 : papszNV[iLine] != NULL &&
957 11 : !EQUALN(papszNV[iLine],pszField,strlen(pszField));
958 : iLine++ ) {}
959 :
960 11 : if( papszNV[iLine] == NULL )
961 11 : return pszDefaultValue;
962 : else
963 : {
964 0 : CPLString osResult;
965 : char **papszTokens;
966 :
967 0 : papszTokens = CSLTokenizeString(papszNV[iLine]);
968 :
969 0 : if( CSLCount(papszTokens) > 1 )
970 0 : osResult = papszTokens[1];
971 : else
972 0 : osResult = pszDefaultValue;
973 :
974 0 : CSLDestroy( papszTokens );
975 0 : return osResult;
976 : }
977 : }
978 :
979 : /************************************************************************/
980 : /* GDALRegister_AAIGrid() */
981 : /************************************************************************/
982 :
983 338 : void GDALRegister_AAIGrid()
984 :
985 : {
986 : GDALDriver *poDriver;
987 :
988 338 : if( GDALGetDriverByName( "AAIGrid" ) == NULL )
989 : {
990 336 : poDriver = new GDALDriver();
991 :
992 336 : poDriver->SetDescription( "AAIGrid" );
993 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
994 336 : "Arc/Info ASCII Grid" );
995 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
996 336 : "frmt_various.html#AAIGrid" );
997 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "asc" );
998 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
999 336 : "Byte UInt16 Int16 Int32 Float32" );
1000 :
1001 336 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1002 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
1003 : "<CreationOptionList>\n"
1004 : " <Option name='FORCE_CELLSIZE' type='boolean' description='Force use of CELLSIZE, default is FALSE.'/>\n"
1005 : " <Option name='DECIMAL_PRECISION' type='int' description='Number of decimal when writing floating-point numbers.'/>\n"
1006 336 : "</CreationOptionList>\n" );
1007 :
1008 336 : poDriver->pfnOpen = AAIGDataset::Open;
1009 336 : poDriver->pfnCreateCopy = AAIGCreateCopy;
1010 :
1011 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
1012 : }
1013 338 : }
1014 :
|
