1 : /******************************************************************************
2 : * $Id: aaigriddataset.cpp 19956 2010-07-02 20:05:53Z rouault $
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 19956 2010-07-02 20:05:53Z rouault $");
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 : AAIGRasterBand::AAIGRasterBand( AAIGDataset *poDS, int nDataStart,
116 60 : GDALDataType eTypeIn )
117 :
118 : {
119 60 : this->poDS = poDS;
120 :
121 60 : nBand = 1;
122 60 : eDataType = eTypeIn;
123 :
124 60 : nBlockXSize = poDS->nRasterXSize;
125 60 : nBlockYSize = 1;
126 :
127 : panLineOffset =
128 60 : (GUIntBig *) VSICalloc( poDS->nRasterYSize, sizeof(GUIntBig) );
129 60 : 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 60 : panLineOffset[0] = nDataStart;
137 0 : }
138 :
139 : /************************************************************************/
140 : /* ~AAIGRasterBand() */
141 : /************************************************************************/
142 :
143 60 : AAIGRasterBand::~AAIGRasterBand()
144 :
145 : {
146 60 : CPLFree( panLineOffset );
147 60 : }
148 :
149 : /************************************************************************/
150 : /* IReadBlock() */
151 : /************************************************************************/
152 :
153 : CPLErr AAIGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
154 733 : void * pImage )
155 :
156 : {
157 733 : AAIGDataset *poODS = (AAIGDataset *) poDS;
158 : int iPixel;
159 :
160 733 : if( nBlockYOff < 0 || nBlockYOff > poODS->nRasterYSize - 1
161 : || nBlockXOff != 0 || panLineOffset == NULL)
162 0 : return CE_Failure;
163 :
164 733 : 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 733 : if( panLineOffset[nBlockYOff] == 0 )
174 0 : return CE_Failure;
175 :
176 :
177 733 : 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 31953 : for( iPixel = 0; iPixel < poODS->nRasterXSize; )
186 : {
187 : char szToken[500];
188 : char chNext;
189 30487 : int iTokenChar = 0;
190 :
191 : /* suck up any pre-white space. */
192 39874 : do {
193 39874 : chNext = poODS->Getc();
194 : } while( isspace( (unsigned char)chNext ) );
195 :
196 151258 : while( chNext != '\0' && !isspace((unsigned char)chNext) )
197 : {
198 90284 : 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 90284 : szToken[iTokenChar++] = chNext;
207 90284 : chNext = poODS->Getc();
208 : }
209 :
210 30487 : 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 30487 : szToken[iTokenChar] = '\0';
221 :
222 30487 : if( pImage != NULL )
223 : {
224 30412 : if( eDataType == GDT_Float32 )
225 1015 : ((float *) pImage)[iPixel] = (float) CPLAtofM(szToken);
226 : else
227 29397 : ((GInt32 *) pImage)[iPixel] = (GInt32) atoi(szToken);
228 : }
229 :
230 30487 : iPixel++;
231 : }
232 :
233 733 : if( nBlockYOff < poODS->nRasterYSize - 1 )
234 697 : panLineOffset[nBlockYOff + 1] = poODS->Tell();
235 :
236 733 : return CE_None;
237 : }
238 :
239 : /************************************************************************/
240 : /* GetNoDataValue() */
241 : /************************************************************************/
242 :
243 29 : double AAIGRasterBand::GetNoDataValue( int * pbSuccess )
244 :
245 : {
246 29 : AAIGDataset *poODS = (AAIGDataset *) poDS;
247 :
248 29 : if( pbSuccess )
249 26 : *pbSuccess = poODS->bNoDataSet;
250 :
251 29 : 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 60 : AAIGDataset::AAIGDataset()
281 :
282 : {
283 60 : papszPrj = NULL;
284 60 : pszProjection = CPLStrdup("");
285 60 : fp = NULL;
286 60 : bNoDataSet = FALSE;
287 60 : dfNoDataValue = -9999.0;
288 :
289 60 : adfGeoTransform[0] = 0.0;
290 60 : adfGeoTransform[1] = 1.0;
291 60 : adfGeoTransform[2] = 0.0;
292 60 : adfGeoTransform[3] = 0.0;
293 60 : adfGeoTransform[4] = 0.0;
294 60 : adfGeoTransform[5] = 1.0;
295 :
296 60 : nOffsetInBuffer = 256;
297 60 : nBufferOffset = 0;
298 60 : }
299 :
300 : /************************************************************************/
301 : /* ~AAIGDataset() */
302 : /************************************************************************/
303 :
304 60 : AAIGDataset::~AAIGDataset()
305 :
306 : {
307 60 : FlushCache();
308 :
309 60 : if( fp != NULL )
310 60 : VSIFCloseL( fp );
311 :
312 60 : CPLFree( pszProjection );
313 60 : CSLDestroy( papszPrj );
314 60 : }
315 :
316 : /************************************************************************/
317 : /* Tell() */
318 : /************************************************************************/
319 :
320 697 : GUIntBig AAIGDataset::Tell()
321 :
322 : {
323 697 : return nBufferOffset + nOffsetInBuffer;
324 : }
325 :
326 : /************************************************************************/
327 : /* Seek() */
328 : /************************************************************************/
329 :
330 733 : int AAIGDataset::Seek( GUIntBig nNewOffset )
331 :
332 : {
333 733 : nOffsetInBuffer = sizeof(achReadBuf);
334 733 : 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 130158 : char AAIGDataset::Getc()
345 :
346 : {
347 130158 : if( nOffsetInBuffer < (int) sizeof(achReadBuf) )
348 129225 : return achReadBuf[nOffsetInBuffer++];
349 :
350 933 : nBufferOffset = VSIFTellL( fp );
351 933 : int nRead = VSIFReadL( achReadBuf, 1, sizeof(achReadBuf), fp );
352 : unsigned int i;
353 20822 : for(i=nRead;i<sizeof(achReadBuf);i++)
354 19889 : achReadBuf[i] = '\0';
355 :
356 933 : nOffsetInBuffer = 0;
357 :
358 933 : 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 11334 : GDALDataset *AAIGDataset::Open( GDALOpenInfo * poOpenInfo )
381 :
382 : {
383 11334 : int i = 0;
384 11334 : int j = 0;
385 11334 : char **papszTokens = NULL;
386 :
387 : /* Default data type */
388 11334 : GDALDataType eDataType = GDT_Int32;
389 :
390 : /* -------------------------------------------------------------------- */
391 : /* Does this look like an AI grid file? */
392 : /* -------------------------------------------------------------------- */
393 11334 : 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 11274 : return NULL;
404 :
405 : papszTokens =
406 : CSLTokenizeString2( (const char *) poOpenInfo->pabyHeader,
407 60 : " \n\r\t", 0 );
408 60 : int nTokens = CSLCount(papszTokens);
409 :
410 : /* -------------------------------------------------------------------- */
411 : /* Create a corresponding GDALDataset. */
412 : /* -------------------------------------------------------------------- */
413 : AAIGDataset *poDS;
414 :
415 60 : poDS = new AAIGDataset();
416 :
417 : /* -------------------------------------------------------------------- */
418 : /* Parse the header. */
419 : /* -------------------------------------------------------------------- */
420 60 : double dfCellDX = 0;
421 60 : double dfCellDY = 0;
422 :
423 120 : 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 60 : poDS->nRasterXSize = atoi(papszTokens[i + 1]);
431 60 : 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 60 : poDS->nRasterYSize = atoi(papszTokens[i + 1]);
439 :
440 60 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
441 : {
442 0 : delete poDS;
443 0 : return NULL;
444 : }
445 :
446 60 : 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 = CPLAtofM( papszTokens[iDX+1] );
460 3 : dfCellDY = CPLAtofM( papszTokens[iDY+1] );
461 : }
462 : else
463 : {
464 57 : if (i + 1 >= nTokens)
465 : {
466 0 : CSLDestroy( papszTokens );
467 0 : delete poDS;
468 0 : return NULL;
469 : }
470 57 : dfCellDX = dfCellDY = CPLAtofM( papszTokens[i + 1] );
471 : }
472 :
473 60 : if ((i = CSLFindString( papszTokens, "xllcorner" )) >= 0 &&
474 : (j = CSLFindString( papszTokens, "yllcorner" )) >= 0 &&
475 : i + 1 < nTokens && j + 1 < nTokens)
476 : {
477 60 : poDS->adfGeoTransform[0] = CPLAtofM( papszTokens[i + 1] );
478 60 : poDS->adfGeoTransform[1] = dfCellDX;
479 60 : poDS->adfGeoTransform[2] = 0.0;
480 : poDS->adfGeoTransform[3] = CPLAtofM( papszTokens[j + 1] )
481 60 : + poDS->nRasterYSize * dfCellDY;
482 60 : poDS->adfGeoTransform[4] = 0.0;
483 60 : 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] = CPLAtofM(papszTokens[i + 1]) - 0.5 * dfCellDX;
492 0 : poDS->adfGeoTransform[1] = dfCellDX;
493 0 : poDS->adfGeoTransform[2] = 0.0;
494 : poDS->adfGeoTransform[3] = CPLAtofM( 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 60 : if( (i = CSLFindString( papszTokens, "NODATA_value" )) >= 0 &&
511 : i + 1 < nTokens)
512 : {
513 18 : const char* pszNoData = papszTokens[i + 1];
514 :
515 18 : poDS->bNoDataSet = TRUE;
516 18 : poDS->dfNoDataValue = CPLAtofM(pszNoData);
517 18 : if( strchr( pszNoData, '.' ) != NULL ||
518 : strchr( pszNoData, ',' ) != NULL ||
519 : INT_MIN > poDS->dfNoDataValue || poDS->dfNoDataValue > INT_MAX )
520 : {
521 1 : eDataType = GDT_Float32;
522 1 : poDS->dfNoDataValue = (double) (float) poDS->dfNoDataValue;
523 : }
524 : }
525 :
526 60 : CSLDestroy( papszTokens );
527 :
528 : /* -------------------------------------------------------------------- */
529 : /* Open file with large file API. */
530 : /* -------------------------------------------------------------------- */
531 :
532 60 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
533 60 : if( poDS->fp == NULL )
534 : {
535 : CPLError( CE_Failure, CPLE_OpenFailed,
536 : "VSIFOpenL(%s) failed unexpectedly.",
537 0 : poOpenInfo->pszFilename );
538 0 : delete poDS;
539 0 : return NULL;
540 : }
541 :
542 : /* -------------------------------------------------------------------- */
543 : /* Find the start of real data. */
544 : /* -------------------------------------------------------------------- */
545 : int nStartOfData;
546 :
547 7724 : for( i = 2; TRUE ; i++ )
548 : {
549 7724 : if( poOpenInfo->pabyHeader[i] == '\0' )
550 : {
551 : CPLError( CE_Failure, CPLE_AppDefined,
552 0 : "Couldn't find data values in ASCII Grid file.\n" );
553 0 : delete poDS;
554 0 : return NULL;
555 : }
556 :
557 7724 : if( poOpenInfo->pabyHeader[i-1] == '\n'
558 : || poOpenInfo->pabyHeader[i-2] == '\n'
559 : || poOpenInfo->pabyHeader[i-1] == '\r'
560 : || poOpenInfo->pabyHeader[i-2] == '\r' )
561 : {
562 612 : if( !isalpha(poOpenInfo->pabyHeader[i])
563 : && poOpenInfo->pabyHeader[i] != '\n'
564 : && poOpenInfo->pabyHeader[i] != '\r')
565 : {
566 60 : nStartOfData = i;
567 :
568 : /* Beginning of real data found. */
569 : break;
570 : }
571 : }
572 : }
573 :
574 : /* -------------------------------------------------------------------- */
575 : /* Recognize the type of data. */
576 : /* -------------------------------------------------------------------- */
577 60 : CPLAssert( NULL != poDS->fp );
578 :
579 60 : if( eDataType != GDT_Float32)
580 : {
581 : /* Allocate 100K chunk + 1 extra byte for NULL character. */
582 59 : const size_t nChunkSize = 1024 * 100;
583 59 : GByte* pabyChunk = (GByte *) CPLCalloc( nChunkSize + 1, sizeof(GByte) );
584 59 : pabyChunk[nChunkSize] = '\0';
585 :
586 59 : VSIFSeekL( poDS->fp, nStartOfData, SEEK_SET );
587 :
588 : /* Scan for dot in subsequent chunks of data. */
589 177 : while( !VSIFEofL( poDS->fp) )
590 : {
591 59 : VSIFReadL( pabyChunk, sizeof(GByte), nChunkSize, poDS->fp );
592 59 : CPLAssert( pabyChunk[nChunkSize] == '\0' );
593 :
594 5324885 : for(i = 0; i < (int)nChunkSize; i++)
595 : {
596 5324833 : GByte ch = pabyChunk[i];
597 5324833 : if (ch == '.' || ch == ',' || ch == 'e' || ch == 'E')
598 : {
599 7 : eDataType = GDT_Float32;
600 7 : break;
601 : }
602 : }
603 : }
604 :
605 : /* Deallocate chunk. */
606 59 : VSIFree( pabyChunk );
607 : }
608 :
609 : /* -------------------------------------------------------------------- */
610 : /* Create band information objects. */
611 : /* -------------------------------------------------------------------- */
612 60 : AAIGRasterBand* band = new AAIGRasterBand( poDS, nStartOfData, eDataType );
613 60 : poDS->SetBand( 1, band );
614 60 : if (band->panLineOffset == NULL)
615 : {
616 0 : delete poDS;
617 0 : return NULL;
618 : }
619 :
620 : /* -------------------------------------------------------------------- */
621 : /* Try to read projection file. */
622 : /* -------------------------------------------------------------------- */
623 : char *pszDirname, *pszBasename;
624 : VSIStatBufL sStatBuf;
625 :
626 60 : pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
627 60 : pszBasename = CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
628 :
629 60 : poDS->osPrjFilename = CPLFormFilename( pszDirname, pszBasename, "prj" );
630 60 : int nRet = VSIStatL( poDS->osPrjFilename, &sStatBuf );
631 :
632 : #ifndef WIN32
633 60 : if( nRet != 0 )
634 : {
635 23 : poDS->osPrjFilename = CPLFormFilename( pszDirname, pszBasename, "PRJ" );
636 23 : nRet = VSIStatL( poDS->osPrjFilename, &sStatBuf );
637 : }
638 : #endif
639 :
640 60 : if( nRet == 0 )
641 : {
642 38 : OGRSpatialReference oSRS;
643 :
644 38 : poDS->papszPrj = CSLLoad( poDS->osPrjFilename );
645 :
646 : CPLDebug( "AAIGrid", "Loaded SRS from %s",
647 38 : poDS->osPrjFilename.c_str() );
648 :
649 38 : if( oSRS.importFromESRI( poDS->papszPrj ) == OGRERR_NONE )
650 : {
651 : // If geographic values are in seconds, we must transform.
652 : // Is there a code for minutes too?
653 38 : if( oSRS.IsGeographic()
654 : && EQUAL(OSR_GDS( poDS->papszPrj, "Units", ""), "DS") )
655 : {
656 0 : poDS->adfGeoTransform[0] /= 3600.0;
657 0 : poDS->adfGeoTransform[1] /= 3600.0;
658 0 : poDS->adfGeoTransform[2] /= 3600.0;
659 0 : poDS->adfGeoTransform[3] /= 3600.0;
660 0 : poDS->adfGeoTransform[4] /= 3600.0;
661 0 : poDS->adfGeoTransform[5] /= 3600.0;
662 : }
663 :
664 38 : CPLFree( poDS->pszProjection );
665 38 : oSRS.exportToWkt( &(poDS->pszProjection) );
666 38 : }
667 : }
668 :
669 60 : CPLFree( pszDirname );
670 60 : CPLFree( pszBasename );
671 :
672 : /* -------------------------------------------------------------------- */
673 : /* Initialize any PAM information. */
674 : /* -------------------------------------------------------------------- */
675 60 : poDS->SetDescription( poOpenInfo->pszFilename );
676 60 : poDS->TryLoadXML();
677 :
678 60 : return( poDS );
679 : }
680 :
681 : /************************************************************************/
682 : /* GetGeoTransform() */
683 : /************************************************************************/
684 :
685 33 : CPLErr AAIGDataset::GetGeoTransform( double * padfTransform )
686 :
687 : {
688 33 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
689 33 : return( CE_None );
690 : }
691 :
692 : /************************************************************************/
693 : /* GetProjectionRef() */
694 : /************************************************************************/
695 :
696 45 : const char *AAIGDataset::GetProjectionRef()
697 :
698 : {
699 45 : return pszProjection;
700 : }
701 :
702 : /************************************************************************/
703 : /* AAIGCreateCopy() */
704 : /************************************************************************/
705 :
706 : static GDALDataset *
707 : AAIGCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
708 : int bStrict, char ** papszOptions,
709 25 : GDALProgressFunc pfnProgress, void * pProgressData )
710 :
711 : {
712 25 : int nBands = poSrcDS->GetRasterCount();
713 25 : int nXSize = poSrcDS->GetRasterXSize();
714 25 : int nYSize = poSrcDS->GetRasterYSize();
715 :
716 : /* -------------------------------------------------------------------- */
717 : /* Some rudimentary checks */
718 : /* -------------------------------------------------------------------- */
719 25 : if( nBands != 1 )
720 : {
721 : CPLError( CE_Failure, CPLE_NotSupported,
722 : "AAIG driver doesn't support %d bands. Must be 1 band.\n",
723 5 : nBands );
724 :
725 5 : return NULL;
726 : }
727 :
728 20 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
729 0 : return NULL;
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Create the dataset. */
733 : /* -------------------------------------------------------------------- */
734 : FILE *fpImage;
735 :
736 20 : fpImage = VSIFOpenL( pszFilename, "wt" );
737 20 : if( fpImage == NULL )
738 : {
739 : CPLError( CE_Failure, CPLE_OpenFailed,
740 : "Unable to create file %s.\n",
741 0 : pszFilename );
742 0 : return NULL;
743 : }
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Write ASCII Grid file header */
747 : /* -------------------------------------------------------------------- */
748 : double adfGeoTransform[6];
749 : char szHeader[2000];
750 : const char *pszForceCellsize =
751 20 : CSLFetchNameValue( papszOptions, "FORCE_CELLSIZE" );
752 :
753 20 : poSrcDS->GetGeoTransform( adfGeoTransform );
754 :
755 20 : if( ABS(adfGeoTransform[1]+adfGeoTransform[5]) < 0.0000001
756 : || ABS(adfGeoTransform[1]-adfGeoTransform[5]) < 0.0000001
757 : || (pszForceCellsize && CSLTestBoolean(pszForceCellsize)) )
758 : sprintf( szHeader,
759 : "ncols %d\n"
760 : "nrows %d\n"
761 : "xllcorner %.12f\n"
762 : "yllcorner %.12f\n"
763 : "cellsize %.12f\n",
764 : nXSize, nYSize,
765 : adfGeoTransform[0],
766 : adfGeoTransform[3] - nYSize * adfGeoTransform[1],
767 19 : adfGeoTransform[1] );
768 : else
769 : {
770 1 : if( pszForceCellsize == NULL )
771 : CPLError( CE_Warning, CPLE_AppDefined,
772 : "Producing a Golden Surfer style file with DX and DY instead\n"
773 : "of CELLSIZE since the input pixels are non-square. Use the\n"
774 : "FORCE_CELLSIZE=TRUE creation option to force use of DX for\n"
775 : "even though this will be distorted. Most ASCII Grid readers\n"
776 1 : "(ArcGIS included) do not support the DX and DY parameters.\n" );
777 : sprintf( szHeader,
778 : "ncols %d\n"
779 : "nrows %d\n"
780 : "xllcorner %.12f\n"
781 : "yllcorner %.12f\n"
782 : "dx %.12f\n"
783 : "dy %.12f\n",
784 : nXSize, nYSize,
785 : adfGeoTransform[0],
786 : adfGeoTransform[3] + nYSize * adfGeoTransform[5],
787 : adfGeoTransform[1],
788 1 : fabs(adfGeoTransform[5]) );
789 : }
790 :
791 : /* -------------------------------------------------------------------- */
792 : /* Handle nodata (optionally). */
793 : /* -------------------------------------------------------------------- */
794 20 : GDALRasterBand * poBand = poSrcDS->GetRasterBand( 1 );
795 : double dfNoData;
796 : int bSuccess;
797 :
798 : // Write `nodata' value to header if it is exists in source dataset
799 20 : dfNoData = poBand->GetNoDataValue( &bSuccess );
800 20 : if ( bSuccess )
801 : sprintf( szHeader+strlen(szHeader), "NODATA_value %6.20g\n",
802 0 : dfNoData );
803 :
804 20 : VSIFWriteL( szHeader, 1, strlen(szHeader), fpImage );
805 :
806 : /* -------------------------------------------------------------------- */
807 : /* Builds the format string used for printing float values. */
808 : /* -------------------------------------------------------------------- */
809 : char szFormatFloat[32];
810 20 : strcpy(szFormatFloat, " %6.20g");
811 : const char *pszDecimalPrecision =
812 20 : CSLFetchNameValue( papszOptions, "DECIMAL_PRECISION" );
813 20 : if (pszDecimalPrecision)
814 : {
815 1 : int nDecimal = atoi(pszDecimalPrecision);
816 1 : if (nDecimal >= 0)
817 1 : sprintf(szFormatFloat, " %%.%df", nDecimal);
818 : }
819 :
820 : /* -------------------------------------------------------------------- */
821 : /* Loop over image, copying image data. */
822 : /* -------------------------------------------------------------------- */
823 20 : int *panScanline = NULL;
824 20 : double *padfScanline = NULL;
825 : int bReadAsInt;
826 : int iLine, iPixel;
827 20 : CPLErr eErr = CE_None;
828 :
829 : bReadAsInt = ( poBand->GetRasterDataType() == GDT_Byte
830 : || poBand->GetRasterDataType() == GDT_Int16
831 : || poBand->GetRasterDataType() == GDT_UInt16
832 20 : || poBand->GetRasterDataType() == GDT_Int32 );
833 :
834 : // Write scanlines to output file
835 20 : if (bReadAsInt)
836 : panScanline = (int *) CPLMalloc( nXSize *
837 11 : GDALGetDataTypeSize(GDT_Int32) / 8 );
838 : else
839 : padfScanline = (double *) CPLMalloc( nXSize *
840 9 : GDALGetDataTypeSize(GDT_Float64) / 8 );
841 :
842 405 : for( iLine = 0; eErr == CE_None && iLine < nYSize; iLine++ )
843 : {
844 : eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
845 : (bReadAsInt) ? (void*)panScanline : (void*)padfScanline,
846 : nXSize, 1, (bReadAsInt) ? GDT_Int32 : GDT_Float64,
847 385 : 0, 0 );
848 :
849 385 : if( bReadAsInt )
850 : {
851 13760 : for ( iPixel = 0; iPixel < nXSize; iPixel++ )
852 : {
853 13485 : sprintf( szHeader, " %d", panScanline[iPixel] );
854 13485 : if( VSIFWriteL( szHeader, strlen(szHeader), 1, fpImage ) != 1 )
855 : {
856 0 : eErr = CE_Failure;
857 : CPLError( CE_Failure, CPLE_AppDefined,
858 0 : "Write failed, disk full?\n" );
859 0 : break;
860 : }
861 : }
862 : }
863 : else
864 : {
865 1610 : for ( iPixel = 0; iPixel < nXSize; iPixel++ )
866 : {
867 1500 : sprintf( szHeader, szFormatFloat, padfScanline[iPixel] );
868 1500 : if( VSIFWriteL( szHeader, strlen(szHeader), 1, fpImage ) != 1 )
869 : {
870 0 : eErr = CE_Failure;
871 : CPLError( CE_Failure, CPLE_AppDefined,
872 0 : "Write failed, disk full?\n" );
873 0 : break;
874 : }
875 : }
876 : }
877 385 : VSIFWriteL( (void *) "\n", 1, 1, fpImage );
878 :
879 385 : if( eErr == CE_None &&
880 : !pfnProgress((iLine + 1) / ((double) nYSize), NULL, pProgressData) )
881 : {
882 0 : eErr = CE_Failure;
883 : CPLError( CE_Failure, CPLE_UserInterrupt,
884 0 : "User terminated CreateCopy()" );
885 : }
886 : }
887 :
888 20 : CPLFree( panScanline );
889 20 : CPLFree( padfScanline );
890 20 : VSIFCloseL( fpImage );
891 :
892 : /* -------------------------------------------------------------------- */
893 : /* Try to write projection file. */
894 : /* -------------------------------------------------------------------- */
895 : const char *pszOriginalProjection;
896 :
897 20 : pszOriginalProjection = (char *)poSrcDS->GetProjectionRef();
898 20 : if( !EQUAL( pszOriginalProjection, "" ) )
899 : {
900 : char *pszDirname, *pszBasename;
901 : char *pszPrjFilename;
902 19 : char *pszESRIProjection = NULL;
903 : FILE *fp;
904 19 : OGRSpatialReference oSRS;
905 :
906 19 : pszDirname = CPLStrdup( CPLGetPath(pszFilename) );
907 19 : pszBasename = CPLStrdup( CPLGetBasename(pszFilename) );
908 :
909 19 : pszPrjFilename = CPLStrdup( CPLFormFilename( pszDirname, pszBasename, "prj" ) );
910 19 : fp = VSIFOpenL( pszPrjFilename, "wt" );
911 19 : if (fp != NULL)
912 : {
913 19 : oSRS.importFromWkt( (char **) &pszOriginalProjection );
914 19 : oSRS.morphToESRI();
915 19 : oSRS.exportToWkt( &pszESRIProjection );
916 19 : VSIFWriteL( pszESRIProjection, 1, strlen(pszESRIProjection), fp );
917 :
918 19 : VSIFCloseL( fp );
919 19 : CPLFree( pszESRIProjection );
920 : }
921 : else
922 : {
923 : CPLError( CE_Failure, CPLE_FileIO,
924 0 : "Unable to create file %s.\n", pszPrjFilename );
925 : }
926 19 : CPLFree( pszDirname );
927 19 : CPLFree( pszBasename );
928 19 : CPLFree( pszPrjFilename );
929 : }
930 :
931 : /* -------------------------------------------------------------------- */
932 : /* Re-open dataset, and copy any auxilary pam information. */
933 : /* -------------------------------------------------------------------- */
934 : GDALPamDataset *poDS = (GDALPamDataset *)
935 20 : GDALOpen( pszFilename, GA_ReadOnly );
936 :
937 20 : if( poDS )
938 20 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
939 :
940 20 : return poDS;
941 : }
942 :
943 : /************************************************************************/
944 : /* OSR_GDS() */
945 : /************************************************************************/
946 :
947 : static CPLString OSR_GDS( char **papszNV, const char * pszField,
948 11 : const char *pszDefaultValue )
949 :
950 : {
951 : int iLine;
952 :
953 11 : if( papszNV == NULL || papszNV[0] == NULL )
954 0 : return pszDefaultValue;
955 :
956 11 : for( iLine = 0;
957 : papszNV[iLine] != NULL &&
958 : !EQUALN(papszNV[iLine],pszField,strlen(pszField));
959 : iLine++ ) {}
960 :
961 11 : if( papszNV[iLine] == NULL )
962 11 : return pszDefaultValue;
963 : else
964 : {
965 0 : CPLString osResult;
966 : char **papszTokens;
967 :
968 0 : papszTokens = CSLTokenizeString(papszNV[iLine]);
969 :
970 0 : if( CSLCount(papszTokens) > 1 )
971 0 : osResult = papszTokens[1];
972 : else
973 0 : osResult = pszDefaultValue;
974 :
975 0 : CSLDestroy( papszTokens );
976 0 : return osResult;
977 : }
978 : }
979 :
980 : /************************************************************************/
981 : /* GDALRegister_AAIGrid() */
982 : /************************************************************************/
983 :
984 409 : void GDALRegister_AAIGrid()
985 :
986 : {
987 : GDALDriver *poDriver;
988 :
989 409 : if( GDALGetDriverByName( "AAIGrid" ) == NULL )
990 : {
991 392 : poDriver = new GDALDriver();
992 :
993 392 : poDriver->SetDescription( "AAIGrid" );
994 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
995 392 : "Arc/Info ASCII Grid" );
996 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
997 392 : "frmt_various.html#AAIGrid" );
998 392 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "asc" );
999 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1000 392 : "Byte UInt16 Int16 Int32 Float32" );
1001 :
1002 392 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1003 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
1004 : "<CreationOptionList>\n"
1005 : " <Option name='FORCE_CELLSIZE' type='boolean' description='Force use of CELLSIZE, default is FALSE.'/>\n"
1006 : " <Option name='DECIMAL_PRECISION' type='int' description='Number of decimal when writing floating-point numbers.'/>\n"
1007 392 : "</CreationOptionList>\n" );
1008 :
1009 392 : poDriver->pfnOpen = AAIGDataset::Open;
1010 392 : poDriver->pfnCreateCopy = AAIGCreateCopy;
1011 :
1012 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
1013 : }
1014 409 : }
1015 :
|
