1 : /******************************************************************************
2 : * $Id: gsagdataset.cpp 25215 2012-11-08 08:25:05Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Implements the Golden Software ASCII Grid Format.
6 : * Author: Kevin Locke, kwl7@cornell.edu
7 : * (Based largely on aaigriddataset.cpp by Frank Warmerdam)
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2006, Kevin Locke <kwl7@cornell.edu>
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "cpl_conv.h"
32 :
33 : #include <sstream>
34 : #include <float.h>
35 : #include <limits.h>
36 : #include <assert.h>
37 :
38 : #include "gdal_pam.h"
39 :
40 : #ifndef DBL_MAX
41 : # ifdef __DBL_MAX__
42 : # define DBL_MAX __DBL_MAX__
43 : # else
44 : # define DBL_MAX 1.7976931348623157E+308
45 : # endif /* __DBL_MAX__ */
46 : #endif /* DBL_MAX */
47 :
48 : #ifndef INT_MAX
49 : # define INT_MAX 2147483647
50 : #endif /* INT_MAX */
51 :
52 : CPL_CVSID("$Id: gsagdataset.cpp 25215 2012-11-08 08:25:05Z rouault $");
53 :
54 : CPL_C_START
55 : void GDALRegister_GSAG(void);
56 : CPL_C_END
57 :
58 : /************************************************************************/
59 : /* ==================================================================== */
60 : /* GSAGDataset */
61 : /* ==================================================================== */
62 : /************************************************************************/
63 :
64 : class GSAGRasterBand;
65 :
66 : class GSAGDataset : public GDALPamDataset
67 : {
68 : friend class GSAGRasterBand;
69 :
70 : static const double dfNODATA_VALUE;
71 : static const int nFIELD_PRECISION;
72 : static const size_t nMAX_HEADER_SIZE;
73 :
74 : static CPLErr ShiftFileContents( VSILFILE *, vsi_l_offset, int, const char * );
75 :
76 : VSILFILE *fp;
77 : size_t nMinMaxZOffset;
78 : char szEOL[3];
79 :
80 : CPLErr UpdateHeader();
81 :
82 : public:
83 : GSAGDataset( const char *pszEOL = "\x0D\x0A" );
84 : ~GSAGDataset();
85 :
86 : static int Identify( GDALOpenInfo * );
87 : static GDALDataset *Open( GDALOpenInfo * );
88 : static GDALDataset *CreateCopy( const char *pszFilename,
89 : GDALDataset *poSrcDS,
90 : int bStrict, char **papszOptions,
91 : GDALProgressFunc pfnProgress,
92 : void *pProgressData );
93 :
94 : CPLErr GetGeoTransform( double *padfGeoTransform );
95 : CPLErr SetGeoTransform( double *padfGeoTransform );
96 : };
97 :
98 : /* NOTE: This is not mentioned in the spec, but Surfer 8 uses this value */
99 : const double GSAGDataset::dfNODATA_VALUE = 1.70141E+38;
100 :
101 : const int GSAGDataset::nFIELD_PRECISION = 14;
102 : const size_t GSAGDataset::nMAX_HEADER_SIZE = 200;
103 :
104 : /************************************************************************/
105 : /* ==================================================================== */
106 : /* GSAGRasterBand */
107 : /* ==================================================================== */
108 : /************************************************************************/
109 :
110 : class GSAGRasterBand : public GDALPamRasterBand
111 : {
112 : friend class GSAGDataset;
113 :
114 : double dfMinX;
115 : double dfMaxX;
116 : double dfMinY;
117 : double dfMaxY;
118 : double dfMinZ;
119 : double dfMaxZ;
120 :
121 : vsi_l_offset *panLineOffset;
122 : int nLastReadLine;
123 :
124 : double *padfRowMinZ;
125 : double *padfRowMaxZ;
126 : int nMinZRow;
127 : int nMaxZRow;
128 :
129 : CPLErr ScanForMinMaxZ();
130 :
131 : public:
132 :
133 : GSAGRasterBand( GSAGDataset *, int, vsi_l_offset );
134 : ~GSAGRasterBand();
135 :
136 : CPLErr IReadBlock( int, int, void * );
137 : CPLErr IWriteBlock( int, int, void * );
138 :
139 : double GetNoDataValue( int *pbSuccess = NULL );
140 : double GetMinimum( int *pbSuccess = NULL );
141 : double GetMaximum( int *pbSuccess = NULL );
142 : };
143 :
144 : /************************************************************************/
145 : /* AlmostEqual() */
146 : /* This function is needed because in release mode "1.70141E+38" is not */
147 : /* parsed as 1.70141E+38 in the last bit of the mantissa. */
148 : /* See http://gcc.gnu.org/ml/gcc/2003-08/msg01195.html for some */
149 : /* explanation. */
150 : /************************************************************************/
151 :
152 400 : bool AlmostEqual( double dfVal1, double dfVal2 )
153 :
154 : {
155 400 : const double dfTOLERANCE = 0.0000000001;
156 400 : if( dfVal1 == 0.0 || dfVal2 == 0.0 )
157 0 : return fabs(dfVal1 - dfVal2) < dfTOLERANCE;
158 400 : return fabs((dfVal1 - dfVal2)/dfVal1) < dfTOLERANCE;
159 : }
160 :
161 : /************************************************************************/
162 : /* GSAGRasterBand() */
163 : /************************************************************************/
164 :
165 16 : GSAGRasterBand::GSAGRasterBand( GSAGDataset *poDS, int nBand,
166 : vsi_l_offset nDataStart ) :
167 : padfRowMinZ(NULL),
168 : padfRowMaxZ(NULL),
169 : nMinZRow(-1),
170 16 : nMaxZRow(-1)
171 :
172 : {
173 16 : this->poDS = poDS;
174 16 : this->nBand = nBand;
175 :
176 16 : eDataType = GDT_Float64;
177 :
178 16 : nBlockXSize = poDS->GetRasterXSize();
179 16 : nBlockYSize = 1;
180 :
181 : panLineOffset =
182 16 : (vsi_l_offset *)VSICalloc( poDS->nRasterYSize+1, sizeof(vsi_l_offset) );
183 16 : if( panLineOffset == NULL )
184 : {
185 : CPLError(CE_Failure, CPLE_OutOfMemory,
186 : "GSAGRasterBand::GSAGRasterBand : Out of memory allocating %d * %d bytes",
187 0 : (int) poDS->nRasterYSize+1, (int) sizeof(vsi_l_offset) );
188 0 : return;
189 : }
190 :
191 16 : panLineOffset[poDS->nRasterYSize-1] = nDataStart;
192 16 : nLastReadLine = poDS->nRasterYSize;
193 0 : }
194 :
195 : /************************************************************************/
196 : /* ~GSAGRasterBand() */
197 : /************************************************************************/
198 :
199 16 : GSAGRasterBand::~GSAGRasterBand()
200 : {
201 16 : CPLFree( panLineOffset );
202 16 : if( padfRowMinZ != NULL )
203 0 : CPLFree( padfRowMinZ );
204 16 : if( padfRowMaxZ != NULL )
205 0 : CPLFree( padfRowMaxZ );
206 16 : }
207 :
208 : /************************************************************************/
209 : /* ScanForMinMaxZ() */
210 : /************************************************************************/
211 :
212 0 : CPLErr GSAGRasterBand::ScanForMinMaxZ()
213 :
214 : {
215 0 : double *padfRowValues = (double *)VSIMalloc2( nBlockXSize, sizeof(double) );
216 0 : if( padfRowValues == NULL )
217 : {
218 : CPLError( CE_Failure, CPLE_OutOfMemory,
219 0 : "Unable to allocate memory for grid row values.\n" );
220 0 : return CE_Failure;
221 : }
222 :
223 0 : double dfNewMinZ = DBL_MAX;
224 0 : double dfNewMaxZ = -DBL_MAX;
225 0 : int nNewMinZRow = 0;
226 0 : int nNewMaxZRow = 0;
227 :
228 : /* Since we have to scan, lets calc. statistics too */
229 0 : double dfSum = 0.0;
230 0 : double dfSum2 = 0.0;
231 0 : unsigned long nValuesRead = 0;
232 0 : for( int iRow=0; iRow<nRasterYSize; iRow++ )
233 : {
234 0 : CPLErr eErr = IReadBlock( 0, iRow, padfRowValues );
235 0 : if( eErr != CE_None )
236 : {
237 0 : VSIFree( padfRowValues );
238 0 : return eErr;
239 : }
240 :
241 0 : padfRowMinZ[iRow] = DBL_MAX;
242 0 : padfRowMaxZ[iRow] = -DBL_MAX;
243 0 : for( int iCell=0; iCell<nRasterXSize; iCell++ )
244 : {
245 0 : if( AlmostEqual(padfRowValues[iCell], GSAGDataset::dfNODATA_VALUE) )
246 0 : continue;
247 :
248 0 : if( padfRowValues[iCell] < padfRowMinZ[iRow] )
249 0 : padfRowMinZ[iRow] = padfRowValues[iCell];
250 :
251 0 : if( padfRowValues[iCell] > padfRowMaxZ[iRow] )
252 0 : padfRowMaxZ[iRow] = padfRowValues[iCell];
253 :
254 0 : dfSum += padfRowValues[iCell];
255 0 : dfSum2 += padfRowValues[iCell] * padfRowValues[iCell];
256 0 : nValuesRead++;
257 : }
258 :
259 0 : if( padfRowMinZ[iRow] < dfNewMinZ )
260 : {
261 0 : dfNewMinZ = padfRowMinZ[iRow];
262 0 : nNewMinZRow = iRow;
263 : }
264 :
265 0 : if( padfRowMaxZ[iRow] > dfNewMaxZ )
266 : {
267 0 : dfNewMaxZ = padfRowMaxZ[iRow];
268 0 : nNewMaxZRow = iRow;
269 : }
270 : }
271 :
272 0 : VSIFree( padfRowValues );
273 :
274 0 : if( nValuesRead == 0 )
275 : {
276 0 : dfMinZ = 0.0;
277 0 : dfMaxZ = 0.0;
278 0 : nMinZRow = 0;
279 0 : nMaxZRow = 0;
280 0 : return CE_None;
281 : }
282 :
283 0 : dfMinZ = dfNewMinZ;
284 0 : dfMaxZ = dfNewMaxZ;
285 0 : nMinZRow = nNewMinZRow;
286 0 : nMaxZRow = nNewMaxZRow;
287 :
288 0 : double dfMean = dfSum / nValuesRead;
289 0 : double dfStdDev = sqrt((dfSum2 / nValuesRead) - (dfMean * dfMean));
290 0 : SetStatistics( dfMinZ, dfMaxZ, dfMean, dfStdDev );
291 :
292 0 : return CE_None;
293 : }
294 :
295 : /************************************************************************/
296 : /* IReadBlock() */
297 : /************************************************************************/
298 :
299 156 : CPLErr GSAGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
300 : void * pImage )
301 : {
302 : static size_t nMaxLineSize = 128;
303 156 : double *pdfImage = (double *)pImage;
304 156 : GSAGDataset *poGDS = (GSAGDataset *)poDS;
305 :
306 156 : assert( poGDS != NULL );
307 :
308 156 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
309 0 : return CE_Failure;
310 :
311 156 : if( panLineOffset[nBlockYOff] == 0 )
312 : {
313 : // Discover the last read block
314 80 : for ( int iFoundLine = nLastReadLine - 1; iFoundLine > nBlockYOff; iFoundLine--)
315 : {
316 76 : if( IReadBlock( nBlockXOff, iFoundLine, NULL) != CE_None )
317 0 : return CE_Failure;
318 : }
319 : }
320 :
321 156 : if( panLineOffset[nBlockYOff] == 0 )
322 0 : return CE_Failure;
323 156 : if( VSIFSeekL( poGDS->fp, panLineOffset[nBlockYOff], SEEK_SET ) != 0 )
324 : {
325 : CPLError( CE_Failure, CPLE_FileIO,
326 : "Can't seek to offset %ld to read grid row %d.",
327 0 : (long) panLineOffset[nBlockYOff], nBlockYOff );
328 0 : return CE_Failure;
329 : }
330 :
331 : size_t nLineBufSize;
332 156 : char *szLineBuf = NULL;
333 : size_t nCharsRead;
334 156 : size_t nCharsExamined = 0;
335 : /* If we know the offsets, we can just read line directly */
336 232 : if( (nBlockYOff > 0) && ( panLineOffset[nBlockYOff-1] != 0 ) )
337 : {
338 76 : assert(panLineOffset[nBlockYOff-1] > panLineOffset[nBlockYOff]);
339 76 : nLineBufSize = (size_t) (panLineOffset[nBlockYOff-1]
340 76 : - panLineOffset[nBlockYOff] + 1);
341 : }
342 : else
343 : {
344 80 : nLineBufSize = nMaxLineSize;
345 : }
346 :
347 156 : szLineBuf = (char *)VSIMalloc( nLineBufSize );
348 156 : if( szLineBuf == NULL )
349 : {
350 : CPLError( CE_Failure, CPLE_OutOfMemory,
351 0 : "Unable to read block, unable to allocate line buffer.\n" );
352 0 : return CE_Failure;
353 : }
354 :
355 156 : nCharsRead = VSIFReadL( szLineBuf, 1, nLineBufSize - 1, poGDS->fp );
356 156 : if( nCharsRead == 0 )
357 : {
358 0 : VSIFree( szLineBuf );
359 : CPLError( CE_Failure, CPLE_FileIO,
360 : "Can't read grid row %d at offset %ld.\n",
361 0 : nBlockYOff, (long) panLineOffset[nBlockYOff] );
362 0 : return CE_Failure;
363 : }
364 156 : szLineBuf[nCharsRead] = '\0';
365 :
366 156 : char *szStart = szLineBuf;
367 156 : char *szEnd = szStart;
368 3276 : for( int iCell=0; iCell<nBlockXSize; szStart = szEnd )
369 : {
370 3120 : double dfValue = CPLStrtod( szStart, &szEnd );
371 3120 : if( szStart == szEnd )
372 : {
373 : /* No number found */
374 :
375 : /* Check if this was an expected failure */
376 0 : while( isspace( (unsigned char)*szStart ) )
377 0 : szStart++;
378 :
379 : /* Found sign at end of input, seek back to re-read it */
380 0 : if ( (*szStart == '-' || *szStart == '+') && *(szStart+1) == '\0' )
381 : {
382 0 : if( VSIFSeekL( poGDS->fp,
383 : VSIFTellL( poGDS->fp)-1,
384 : SEEK_SET ) != 0 )
385 : {
386 0 : VSIFree( szLineBuf );
387 : CPLError( CE_Failure, CPLE_FileIO,
388 : "Unable to seek in grid row %d "
389 : "(offset %ld, seek %d).\n",
390 : nBlockYOff,
391 : (long) VSIFTellL(poGDS->fp),
392 0 : -1 );
393 :
394 0 : return CE_Failure;
395 : }
396 : }
397 0 : else if( *szStart != '\0' )
398 : {
399 0 : szEnd = szStart;
400 0 : while( !isspace( (unsigned char)*szEnd ) && *szEnd != '\0' )
401 0 : szEnd++;
402 0 : char cOldEnd = *szEnd;
403 0 : *szEnd = '\0';
404 :
405 : CPLError( CE_Warning, CPLE_FileIO,
406 : "Unexpected value in grid row %d (expected floating "
407 : "point value, found \"%s\").\n",
408 0 : nBlockYOff, szStart );
409 :
410 0 : *szEnd = cOldEnd;
411 :
412 0 : szEnd = szStart;
413 0 : while( !isdigit( *szEnd ) && *szEnd != '.' && *szEnd != '\0' )
414 0 : szEnd++;
415 :
416 0 : continue;
417 : }
418 0 : else if( static_cast<size_t>(szStart - szLineBuf) != nCharsRead )
419 : {
420 : CPLError( CE_Warning, CPLE_FileIO,
421 : "Unexpected ASCII null-character in grid row %d at "
422 : "offset %ld.\n",
423 : nBlockYOff,
424 0 : (long) (szStart - szLineBuf) );
425 :
426 0 : while( *szStart == '\0' &&
427 : static_cast<size_t>(szStart - szLineBuf) < nCharsRead )
428 0 : szStart++;
429 :
430 0 : szEnd = szStart;
431 0 : continue;
432 : }
433 :
434 0 : nCharsExamined += szStart - szLineBuf;
435 0 : nCharsRead = VSIFReadL( szLineBuf, 1, nLineBufSize - 1, poGDS->fp );
436 0 : if( nCharsRead == 0 )
437 : {
438 0 : VSIFree( szLineBuf );
439 : CPLError( CE_Failure, CPLE_FileIO,
440 : "Can't read portion of grid row %d at offset %ld.",
441 0 : nBlockYOff, (long) panLineOffset[nBlockYOff] );
442 0 : return CE_Failure;
443 : }
444 0 : szLineBuf[nCharsRead] = '\0';
445 0 : szStart = szEnd = szLineBuf;
446 0 : continue;
447 : }
448 3120 : else if( *szEnd == '\0'
449 : || (*szEnd == '.' && *(szEnd+1) == '\0')
450 : || (*szEnd == '-' && *(szEnd+1) == '\0')
451 : || (*szEnd == '+' && *(szEnd+1) == '\0')
452 : || (*szEnd == 'E' && *(szEnd+1) == '\0')
453 : || (*szEnd == 'E' && *(szEnd+1) == '-' && *(szEnd+2) == '\0')
454 : || (*szEnd == 'E' && *(szEnd+1) == '+' && *(szEnd+2) == '\0')
455 : || (*szEnd == 'e' && *(szEnd+1) == '\0')
456 : || (*szEnd == 'e' && *(szEnd+1) == '-' && *(szEnd+2) == '\0')
457 : || (*szEnd == 'e' && *(szEnd+1) == '+' && *(szEnd+2) == '\0'))
458 : {
459 : /* Number was interrupted by a nul character */
460 0 : while( *szEnd != '\0' )
461 0 : szEnd++;
462 :
463 0 : if( static_cast<size_t>(szEnd - szLineBuf) != nCharsRead )
464 : {
465 : CPLError( CE_Warning, CPLE_FileIO,
466 : "Unexpected ASCII null-character in grid row %d at "
467 : "offset %ld.\n",
468 : nBlockYOff,
469 0 : (long) (szStart - szLineBuf) );
470 :
471 0 : while( *szEnd == '\0' &&
472 : static_cast<size_t>(szStart - szLineBuf) < nCharsRead )
473 0 : szEnd++;
474 :
475 0 : continue;
476 : }
477 :
478 : /* End of buffer, could be interrupting a number */
479 0 : if( VSIFSeekL( poGDS->fp, szStart - szEnd, SEEK_CUR ) != 0 )
480 : {
481 0 : VSIFree( szLineBuf );
482 : CPLError( CE_Failure, CPLE_FileIO,
483 : "Unable to seek in grid row %d (offset %ld, seek %d)"
484 : ".\n", nBlockYOff,
485 : (long) VSIFTellL(poGDS->fp),
486 0 : (int) (szStart - szEnd) );
487 :
488 0 : return CE_Failure;
489 : }
490 0 : nCharsExamined += szStart - szLineBuf;
491 0 : nCharsRead = VSIFReadL( szLineBuf, 1, nLineBufSize - 1, poGDS->fp );
492 0 : szLineBuf[nCharsRead] = '\0';
493 :
494 0 : if( nCharsRead == 0 )
495 : {
496 0 : VSIFree( szLineBuf );
497 : CPLError( CE_Failure, CPLE_FileIO,
498 : "Can't read portion of grid row %d at offset %ld.",
499 0 : nBlockYOff, (long) panLineOffset[nBlockYOff] );
500 0 : return CE_Failure;
501 : }
502 0 : else if( nCharsRead > static_cast<size_t>(szEnd - szStart) )
503 : {
504 : /* Read new data, this was not really the end */
505 0 : szEnd = szStart = szLineBuf;
506 0 : continue;
507 : }
508 :
509 : /* This is really the last value and has no tailing newline */
510 0 : szStart = szLineBuf;
511 0 : szEnd = szLineBuf + nCharsRead;
512 : }
513 :
514 3120 : if( pdfImage != NULL )
515 : {
516 1600 : *(pdfImage+iCell) = dfValue;
517 : }
518 :
519 3120 : iCell++;
520 : }
521 :
522 468 : while( *szEnd == ' ' )
523 156 : szEnd++;
524 :
525 156 : if( *szEnd != '\0' && *szEnd != poGDS->szEOL[0] )
526 : CPLDebug( "GSAG", "Grid row %d does not end with a newline. "
527 0 : "Possible skew.\n", nBlockYOff );
528 :
529 986 : while( isspace( (unsigned char)*szEnd ) )
530 674 : szEnd++;
531 :
532 156 : nCharsExamined += szEnd - szLineBuf;
533 :
534 156 : if( nCharsExamined >= nMaxLineSize )
535 0 : nMaxLineSize = nCharsExamined + 1;
536 :
537 156 : if( nBlockYOff > 0 )
538 152 : panLineOffset[nBlockYOff - 1] =
539 152 : panLineOffset[nBlockYOff] + nCharsExamined;
540 :
541 156 : nLastReadLine = nBlockYOff;
542 :
543 156 : VSIFree( szLineBuf );
544 :
545 156 : return CE_None;
546 : }
547 :
548 : /************************************************************************/
549 : /* IWriteBlock() */
550 : /************************************************************************/
551 :
552 0 : CPLErr GSAGRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
553 : void * pImage )
554 :
555 : {
556 0 : if( eAccess == GA_ReadOnly )
557 : {
558 : CPLError( CE_Failure, CPLE_NoWriteAccess,
559 0 : "Unable to write block, dataset opened read only.\n" );
560 0 : return CE_Failure;
561 : }
562 :
563 0 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
564 0 : return CE_Failure;
565 :
566 0 : GSAGDataset *poGDS = (GSAGDataset *)poDS;
567 0 : assert( poGDS != NULL );
568 :
569 0 : if( padfRowMinZ == NULL || padfRowMaxZ == NULL
570 : || nMinZRow < 0 || nMaxZRow < 0 )
571 : {
572 0 : padfRowMinZ = (double *)VSIMalloc2( nRasterYSize,sizeof(double) );
573 0 : if( padfRowMinZ == NULL )
574 : {
575 : CPLError( CE_Failure, CPLE_OutOfMemory,
576 0 : "Unable to allocate space for row minimums array.\n" );
577 0 : return CE_Failure;
578 : }
579 :
580 0 : padfRowMaxZ = (double *)VSIMalloc2( nRasterYSize,sizeof(double) );
581 0 : if( padfRowMaxZ == NULL )
582 : {
583 0 : VSIFree( padfRowMinZ );
584 0 : padfRowMinZ = NULL;
585 : CPLError( CE_Failure, CPLE_OutOfMemory,
586 0 : "Unable to allocate space for row maximums array.\n" );
587 0 : return CE_Failure;
588 : }
589 :
590 0 : CPLErr eErr = ScanForMinMaxZ();
591 0 : if( eErr != CE_None )
592 0 : return eErr;
593 : }
594 :
595 0 : if( panLineOffset[nBlockYOff+1] == 0 )
596 0 : IReadBlock( nBlockXOff, nBlockYOff, NULL );
597 :
598 0 : if( panLineOffset[nBlockYOff+1] == 0 || panLineOffset[nBlockYOff] == 0 )
599 0 : return CE_Failure;
600 :
601 0 : std::ostringstream ssOutBuf;
602 0 : ssOutBuf.precision( GSAGDataset::nFIELD_PRECISION );
603 0 : ssOutBuf.setf( std::ios::uppercase );
604 :
605 0 : double *pdfImage = (double *)pImage;
606 0 : padfRowMinZ[nBlockYOff] = DBL_MAX;
607 0 : padfRowMaxZ[nBlockYOff] = -DBL_MAX;
608 0 : for( int iCell=0; iCell<nBlockXSize; )
609 : {
610 0 : for( int iCol=0; iCol<10 && iCell<nBlockXSize; iCol++, iCell++ )
611 : {
612 0 : if( AlmostEqual( pdfImage[iCell], GSAGDataset::dfNODATA_VALUE ) )
613 : {
614 0 : if( pdfImage[iCell] < padfRowMinZ[nBlockYOff] )
615 0 : padfRowMinZ[nBlockYOff] = pdfImage[iCell];
616 :
617 0 : if( pdfImage[iCell] > padfRowMaxZ[nBlockYOff] )
618 0 : padfRowMaxZ[nBlockYOff] = pdfImage[iCell];
619 : }
620 :
621 0 : ssOutBuf << pdfImage[iCell] << " ";
622 : }
623 0 : ssOutBuf << poGDS->szEOL;
624 : }
625 0 : ssOutBuf << poGDS->szEOL;
626 :
627 0 : CPLString sOut = ssOutBuf.str();
628 0 : if( sOut.length() != panLineOffset[nBlockYOff+1]-panLineOffset[nBlockYOff] )
629 : {
630 0 : int nShiftSize = (int) (sOut.length() - (panLineOffset[nBlockYOff+1]
631 0 : - panLineOffset[nBlockYOff]));
632 0 : if( nBlockYOff != poGDS->nRasterYSize
633 : && GSAGDataset::ShiftFileContents( poGDS->fp,
634 : panLineOffset[nBlockYOff+1],
635 : nShiftSize,
636 : poGDS->szEOL ) != CE_None )
637 : {
638 : CPLError( CE_Failure, CPLE_FileIO,
639 : "Failure writing block, "
640 0 : "unable to shift file contents.\n" );
641 0 : return CE_Failure;
642 : }
643 :
644 0 : for( size_t iLine=nBlockYOff+1;
645 : iLine < static_cast<unsigned>(poGDS->nRasterYSize+1)
646 0 : && panLineOffset[iLine] != 0; iLine++ )
647 0 : panLineOffset[iLine] += nShiftSize;
648 : }
649 :
650 0 : if( VSIFSeekL( poGDS->fp, panLineOffset[nBlockYOff], SEEK_SET ) != 0 )
651 : {
652 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to seek to grid line.\n" );
653 0 : return CE_Failure;
654 : }
655 :
656 0 : if( VSIFWriteL( sOut.c_str(), 1, sOut.length(),
657 : poGDS->fp ) != sOut.length() )
658 : {
659 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to write grid block.\n" );
660 0 : return CE_Failure;
661 : }
662 :
663 : /* Update header as needed */
664 0 : bool bHeaderNeedsUpdate = false;
665 0 : if( nMinZRow == nBlockYOff && padfRowMinZ[nBlockYOff] > dfMinZ )
666 : {
667 0 : double dfNewMinZ = -DBL_MAX;
668 0 : for( int iRow=0; iRow<nRasterYSize; iRow++ )
669 : {
670 0 : if( padfRowMinZ[iRow] < dfNewMinZ )
671 : {
672 0 : dfNewMinZ = padfRowMinZ[iRow];
673 0 : nMinZRow = iRow;
674 : }
675 : }
676 :
677 0 : if( dfNewMinZ != dfMinZ )
678 : {
679 0 : dfMinZ = dfNewMinZ;
680 0 : bHeaderNeedsUpdate = true;
681 : }
682 : }
683 :
684 0 : if( nMaxZRow == nBlockYOff && padfRowMaxZ[nBlockYOff] < dfMaxZ )
685 : {
686 0 : double dfNewMaxZ = -DBL_MAX;
687 0 : for( int iRow=0; iRow<nRasterYSize; iRow++ )
688 : {
689 0 : if( padfRowMaxZ[iRow] > dfNewMaxZ )
690 : {
691 0 : dfNewMaxZ = padfRowMaxZ[iRow];
692 0 : nMaxZRow = iRow;
693 : }
694 : }
695 :
696 0 : if( dfNewMaxZ != dfMaxZ )
697 : {
698 0 : dfMaxZ = dfNewMaxZ;
699 0 : bHeaderNeedsUpdate = true;
700 : }
701 : }
702 :
703 0 : if( padfRowMinZ[nBlockYOff] < dfMinZ || padfRowMaxZ[nBlockYOff] > dfMaxZ )
704 : {
705 0 : if( padfRowMinZ[nBlockYOff] < dfMinZ )
706 : {
707 0 : dfMinZ = padfRowMinZ[nBlockYOff];
708 0 : nMinZRow = nBlockYOff;
709 : }
710 :
711 0 : if( padfRowMaxZ[nBlockYOff] > dfMaxZ )
712 : {
713 0 : dfMaxZ = padfRowMaxZ[nBlockYOff];
714 0 : nMaxZRow = nBlockYOff;
715 : }
716 :
717 0 : bHeaderNeedsUpdate = true;
718 : }
719 :
720 0 : if( bHeaderNeedsUpdate && dfMaxZ > dfMinZ )
721 : {
722 0 : CPLErr eErr = poGDS->UpdateHeader();
723 0 : return eErr;
724 : }
725 :
726 0 : return CE_None;
727 : }
728 :
729 : /************************************************************************/
730 : /* GetNoDataValue() */
731 : /************************************************************************/
732 :
733 8 : double GSAGRasterBand::GetNoDataValue( int * pbSuccess )
734 : {
735 8 : if( pbSuccess )
736 7 : *pbSuccess = TRUE;
737 :
738 8 : return GSAGDataset::dfNODATA_VALUE;
739 : }
740 :
741 : /************************************************************************/
742 : /* GetMinimum() */
743 : /************************************************************************/
744 :
745 0 : double GSAGRasterBand::GetMinimum( int *pbSuccess )
746 : {
747 0 : if( pbSuccess )
748 0 : *pbSuccess = TRUE;
749 :
750 0 : return dfMinZ;
751 : }
752 :
753 : /************************************************************************/
754 : /* GetMaximum() */
755 : /************************************************************************/
756 :
757 0 : double GSAGRasterBand::GetMaximum( int *pbSuccess )
758 : {
759 0 : if( pbSuccess )
760 0 : *pbSuccess = TRUE;
761 :
762 0 : return dfMaxZ;
763 : }
764 :
765 : /************************************************************************/
766 : /* ==================================================================== */
767 : /* GSAGDataset */
768 : /* ==================================================================== */
769 : /************************************************************************/
770 :
771 : /************************************************************************/
772 : /* GSAGDataset() */
773 : /************************************************************************/
774 :
775 16 : GSAGDataset::GSAGDataset( const char *pszEOL )
776 :
777 : {
778 16 : if( pszEOL == NULL || EQUAL(pszEOL, "") )
779 : {
780 0 : CPLDebug( "GSAG", "GSAGDataset() created with invalid EOL string.\n" );
781 0 : this->szEOL[0] = '\x0D';
782 0 : this->szEOL[1] = '\x0A';
783 0 : this->szEOL[2] = '\0';
784 : }
785 : else
786 : {
787 16 : strncpy(this->szEOL, pszEOL, sizeof(this->szEOL));
788 16 : this->szEOL[sizeof(this->szEOL) - 1] = '\0';
789 : }
790 16 : }
791 :
792 : /************************************************************************/
793 : /* ~GSAGDataset() */
794 : /************************************************************************/
795 :
796 16 : GSAGDataset::~GSAGDataset()
797 :
798 : {
799 16 : FlushCache();
800 16 : if( fp != NULL )
801 16 : VSIFCloseL( fp );
802 16 : }
803 :
804 : /************************************************************************/
805 : /* Identify() */
806 : /************************************************************************/
807 :
808 12490 : int GSAGDataset::Identify( GDALOpenInfo * poOpenInfo )
809 :
810 : {
811 : /* Check for signature */
812 12506 : if( poOpenInfo->nHeaderBytes < 5
813 : || !EQUALN((const char *) poOpenInfo->pabyHeader,"DSAA",4)
814 16 : || ( poOpenInfo->pabyHeader[4] != '\x0D'
815 0 : && poOpenInfo->pabyHeader[4] != '\x0A' ))
816 : {
817 12474 : return FALSE;
818 : }
819 16 : return TRUE;
820 : }
821 :
822 : /************************************************************************/
823 : /* Open() */
824 : /************************************************************************/
825 :
826 2397 : GDALDataset *GSAGDataset::Open( GDALOpenInfo * poOpenInfo )
827 :
828 : {
829 2397 : if( !Identify(poOpenInfo) )
830 : {
831 2381 : return NULL;
832 : }
833 :
834 : /* Identify the end of line marker (should be \x0D\x0A, but try some others)
835 : * (note that '\x0D' == '\r' and '\x0A' == '\n' on most systems) */
836 : char szEOL[3];
837 16 : szEOL[0] = poOpenInfo->pabyHeader[4];
838 16 : szEOL[1] = poOpenInfo->pabyHeader[5];
839 16 : szEOL[2] = '\0';
840 16 : if( szEOL[1] != '\x0D' && szEOL[1] != '\x0A' )
841 0 : szEOL[1] = '\0';
842 :
843 : /* -------------------------------------------------------------------- */
844 : /* Create a corresponding GDALDataset. */
845 : /* -------------------------------------------------------------------- */
846 16 : GSAGDataset *poDS = new GSAGDataset( szEOL );
847 :
848 : /* -------------------------------------------------------------------- */
849 : /* Open file with large file API. */
850 : /* -------------------------------------------------------------------- */
851 :
852 16 : poDS->eAccess = poOpenInfo->eAccess;
853 16 : if( poOpenInfo->eAccess == GA_ReadOnly )
854 4 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
855 : else
856 12 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
857 :
858 16 : if( poDS->fp == NULL )
859 : {
860 : CPLError( CE_Failure, CPLE_OpenFailed,
861 : "VSIFOpenL(%s) failed unexpectedly.",
862 0 : poOpenInfo->pszFilename );
863 0 : delete poDS;
864 0 : return NULL;
865 : }
866 :
867 : /* -------------------------------------------------------------------- */
868 : /* Read the header. */
869 : /* -------------------------------------------------------------------- */
870 : char *pabyHeader;
871 16 : bool bMustFreeHeader = false;
872 16 : if( poOpenInfo->nHeaderBytes >= static_cast<int>(nMAX_HEADER_SIZE) )
873 : {
874 16 : pabyHeader = (char *)poOpenInfo->pabyHeader;
875 : }
876 : else
877 : {
878 0 : bMustFreeHeader = true;
879 0 : pabyHeader = (char *)VSIMalloc( nMAX_HEADER_SIZE );
880 0 : if( pabyHeader == NULL )
881 : {
882 : CPLError( CE_Failure, CPLE_OutOfMemory,
883 0 : "Unable to open dataset, unable to header buffer.\n" );
884 0 : return NULL;
885 : }
886 :
887 0 : size_t nRead = VSIFReadL( pabyHeader, 1, nMAX_HEADER_SIZE-1, poDS->fp );
888 0 : pabyHeader[nRead] = '\0';
889 : }
890 :
891 16 : const char *szErrorMsg = NULL;
892 16 : const char *szStart = pabyHeader + 5;
893 : char *szEnd;
894 : double dfTemp;
895 : double dfMinX;
896 : double dfMaxX;
897 : double dfMinY;
898 : double dfMaxY;
899 : double dfMinZ;
900 : double dfMaxZ;
901 :
902 : /* Parse number of X axis grid rows */
903 16 : long nTemp = strtol( szStart, &szEnd, 10 );
904 16 : if( szStart == szEnd || nTemp < 0l )
905 : {
906 0 : szErrorMsg = "Unable to parse the number of X axis grid columns.\n";
907 0 : goto error;
908 : }
909 16 : else if( nTemp > INT_MAX )
910 : {
911 : CPLError( CE_Warning, CPLE_AppDefined,
912 0 : "Number of X axis grid columns not representable.\n" );
913 0 : poDS->nRasterXSize = INT_MAX;
914 : }
915 16 : else if ( nTemp == 0 )
916 : {
917 0 : szErrorMsg = "Number of X axis grid columns is zero, which is invalid.\n";
918 0 : goto error;
919 : }
920 : else
921 : {
922 16 : poDS->nRasterXSize = static_cast<int>(nTemp);
923 : }
924 16 : szStart = szEnd;
925 :
926 : /* Parse number of Y axis grid rows */
927 16 : nTemp = strtol( szStart, &szEnd, 10 );
928 16 : if( szStart == szEnd || nTemp < 0l )
929 : {
930 0 : szErrorMsg = "Unable to parse the number of Y axis grid rows.\n";
931 0 : goto error;
932 : }
933 16 : else if( nTemp > INT_MAX )
934 : {
935 : CPLError( CE_Warning, CPLE_AppDefined,
936 0 : "Number of Y axis grid rows not representable.\n" );
937 0 : poDS->nRasterYSize = INT_MAX;
938 : }
939 16 : else if ( nTemp == 0)
940 : {
941 0 : szErrorMsg = "Number of Y axis grid rows is zero, which is invalid.\n";
942 0 : goto error;
943 : }
944 : else
945 : {
946 16 : poDS->nRasterYSize = static_cast<int>(nTemp);
947 : }
948 16 : szStart = szEnd;
949 :
950 : /* Parse the minimum X value of the grid */
951 16 : dfTemp = CPLStrtod( szStart, &szEnd );
952 16 : if( szStart == szEnd )
953 : {
954 0 : szErrorMsg = "Unable to parse the minimum X value.\n";
955 0 : goto error;
956 : }
957 : else
958 : {
959 16 : dfMinX = dfTemp;
960 : }
961 16 : szStart = szEnd;
962 :
963 : /* Parse the maximum X value of the grid */
964 16 : dfTemp = CPLStrtod( szStart, &szEnd );
965 16 : if( szStart == szEnd )
966 : {
967 0 : szErrorMsg = "Unable to parse the maximum X value.\n";
968 0 : goto error;
969 : }
970 : else
971 : {
972 16 : dfMaxX = dfTemp;
973 : }
974 16 : szStart = szEnd;
975 :
976 : /* Parse the minimum Y value of the grid */
977 16 : dfTemp = CPLStrtod( szStart, &szEnd );
978 16 : if( szStart == szEnd )
979 : {
980 0 : szErrorMsg = "Unable to parse the minimum Y value.\n";
981 0 : goto error;
982 : }
983 : else
984 : {
985 16 : dfMinY = dfTemp;
986 : }
987 16 : szStart = szEnd;
988 :
989 : /* Parse the maximum Y value of the grid */
990 16 : dfTemp = CPLStrtod( szStart, &szEnd );
991 16 : if( szStart == szEnd )
992 : {
993 0 : szErrorMsg = "Unable to parse the maximum Y value.\n";
994 0 : goto error;
995 : }
996 : else
997 : {
998 16 : dfMaxY = dfTemp;
999 : }
1000 16 : szStart = szEnd;
1001 :
1002 : /* Parse the minimum Z value of the grid */
1003 64 : while( isspace( (unsigned char)*szStart ) )
1004 32 : szStart++;
1005 16 : poDS->nMinMaxZOffset = szStart - pabyHeader;
1006 :
1007 16 : dfTemp = CPLStrtod( szStart, &szEnd );
1008 16 : if( szStart == szEnd )
1009 : {
1010 0 : szErrorMsg = "Unable to parse the minimum Z value.\n";
1011 0 : goto error;
1012 : }
1013 : else
1014 : {
1015 16 : dfMinZ = dfTemp;
1016 : }
1017 16 : szStart = szEnd;
1018 :
1019 : /* Parse the maximum Z value of the grid */
1020 16 : dfTemp = CPLStrtod( szStart, &szEnd );
1021 16 : if( szStart == szEnd )
1022 : {
1023 0 : szErrorMsg = "Unable to parse the maximum Z value.\n";
1024 0 : goto error;
1025 : }
1026 : else
1027 : {
1028 16 : dfMaxZ = dfTemp;
1029 : }
1030 :
1031 64 : while( isspace((unsigned char)*szEnd) )
1032 32 : szEnd++;
1033 :
1034 : /* -------------------------------------------------------------------- */
1035 : /* Create band information objects. */
1036 : /* -------------------------------------------------------------------- */
1037 : {
1038 16 : GSAGRasterBand *poBand = new GSAGRasterBand( poDS, 1, szEnd-pabyHeader );
1039 16 : if( poBand->panLineOffset == NULL )
1040 : {
1041 0 : delete poBand;
1042 0 : goto error;
1043 : }
1044 :
1045 16 : poBand->dfMinX = dfMinX;
1046 16 : poBand->dfMaxX = dfMaxX;
1047 16 : poBand->dfMinY = dfMinY;
1048 16 : poBand->dfMaxY = dfMaxY;
1049 16 : poBand->dfMinZ = dfMinZ;
1050 16 : poBand->dfMaxZ = dfMaxZ;
1051 :
1052 16 : poDS->SetBand( 1, poBand );
1053 : }
1054 :
1055 16 : if( bMustFreeHeader )
1056 : {
1057 0 : CPLFree( pabyHeader );
1058 : }
1059 :
1060 : /* -------------------------------------------------------------------- */
1061 : /* Initialize any PAM information. */
1062 : /* -------------------------------------------------------------------- */
1063 16 : poDS->SetDescription( poOpenInfo->pszFilename );
1064 16 : poDS->TryLoadXML();
1065 :
1066 : /* -------------------------------------------------------------------- */
1067 : /* Check for external overviews. */
1068 : /* -------------------------------------------------------------------- */
1069 16 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename, poOpenInfo->papszSiblingFiles );
1070 :
1071 16 : return( poDS );
1072 :
1073 : error:
1074 0 : if ( bMustFreeHeader )
1075 : {
1076 0 : CPLFree( pabyHeader );
1077 : }
1078 :
1079 0 : delete poDS;
1080 :
1081 0 : if (szErrorMsg)
1082 0 : CPLError( CE_Failure, CPLE_AppDefined, "%s", szErrorMsg );
1083 0 : return NULL;
1084 : }
1085 :
1086 : /************************************************************************/
1087 : /* GetGeoTransform() */
1088 : /************************************************************************/
1089 :
1090 17 : CPLErr GSAGDataset::GetGeoTransform( double *padfGeoTransform )
1091 : {
1092 17 : if( padfGeoTransform == NULL )
1093 0 : return CE_Failure;
1094 :
1095 17 : GSAGRasterBand *poGRB = (GSAGRasterBand *)GetRasterBand( 1 );
1096 :
1097 17 : if( poGRB == NULL )
1098 : {
1099 0 : padfGeoTransform[0] = 0;
1100 0 : padfGeoTransform[1] = 1;
1101 0 : padfGeoTransform[2] = 0;
1102 0 : padfGeoTransform[3] = 0;
1103 0 : padfGeoTransform[4] = 0;
1104 0 : padfGeoTransform[5] = 1;
1105 0 : return CE_Failure;
1106 : }
1107 :
1108 : /* check if we have a PAM GeoTransform stored */
1109 17 : CPLPushErrorHandler( CPLQuietErrorHandler );
1110 17 : CPLErr eErr = GDALPamDataset::GetGeoTransform( padfGeoTransform );
1111 17 : CPLPopErrorHandler();
1112 :
1113 17 : if( eErr == CE_None )
1114 0 : return CE_None;
1115 :
1116 : /* calculate pixel size first */
1117 17 : padfGeoTransform[1] = (poGRB->dfMaxX - poGRB->dfMinX)/(nRasterXSize - 1);
1118 17 : padfGeoTransform[5] = (poGRB->dfMinY - poGRB->dfMaxY)/(nRasterYSize - 1);
1119 :
1120 : /* then calculate image origin */
1121 17 : padfGeoTransform[0] = poGRB->dfMinX - padfGeoTransform[1] / 2;
1122 17 : padfGeoTransform[3] = poGRB->dfMaxY - padfGeoTransform[5] / 2;
1123 :
1124 : /* tilt/rotation does not supported by the GS grids */
1125 17 : padfGeoTransform[4] = 0.0;
1126 17 : padfGeoTransform[2] = 0.0;
1127 :
1128 17 : return CE_None;
1129 : }
1130 :
1131 : /************************************************************************/
1132 : /* SetGeoTransform() */
1133 : /************************************************************************/
1134 :
1135 0 : CPLErr GSAGDataset::SetGeoTransform( double *padfGeoTransform )
1136 : {
1137 0 : if( eAccess == GA_ReadOnly )
1138 : {
1139 : CPLError( CE_Failure, CPLE_NoWriteAccess,
1140 0 : "Unable to set GeoTransform, dataset opened read only.\n" );
1141 0 : return CE_Failure;
1142 : }
1143 :
1144 0 : GSAGRasterBand *poGRB = (GSAGRasterBand *)GetRasterBand( 1 );
1145 :
1146 0 : if( poGRB == NULL || padfGeoTransform == NULL)
1147 0 : return CE_Failure;
1148 :
1149 : /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
1150 0 : CPLErr eErr = CE_None;
1151 : /*if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0
1152 : || padfGeoTransform[1] < 0.0 || padfGeoTransform[5] < 0.0 )
1153 : eErr = GDALPamDataset::SetGeoTransform( padfGeoTransform );*/
1154 :
1155 0 : if( eErr != CE_None )
1156 0 : return eErr;
1157 :
1158 0 : double dfOldMinX = poGRB->dfMinX;
1159 0 : double dfOldMaxX = poGRB->dfMaxX;
1160 0 : double dfOldMinY = poGRB->dfMinY;
1161 0 : double dfOldMaxY = poGRB->dfMaxY;
1162 :
1163 0 : poGRB->dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
1164 : poGRB->dfMaxX =
1165 0 : padfGeoTransform[1] * (nRasterXSize - 0.5) + padfGeoTransform[0];
1166 : poGRB->dfMinY =
1167 0 : padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
1168 0 : poGRB->dfMaxY = padfGeoTransform[3] + padfGeoTransform[5] / 2;
1169 :
1170 0 : eErr = UpdateHeader();
1171 :
1172 0 : if( eErr != CE_None )
1173 : {
1174 0 : poGRB->dfMinX = dfOldMinX;
1175 0 : poGRB->dfMaxX = dfOldMaxX;
1176 0 : poGRB->dfMinY = dfOldMinY;
1177 0 : poGRB->dfMaxY = dfOldMaxY;
1178 : }
1179 :
1180 0 : return eErr;
1181 : }
1182 :
1183 : /************************************************************************/
1184 : /* ShiftFileContents() */
1185 : /************************************************************************/
1186 12 : CPLErr GSAGDataset::ShiftFileContents( VSILFILE *fp, vsi_l_offset nShiftStart,
1187 : int nShiftSize, const char *pszEOL )
1188 : {
1189 : /* nothing to do for zero-shift */
1190 12 : if( nShiftSize == 0 )
1191 0 : return CE_None;
1192 :
1193 : /* make sure start location is sane */
1194 12 : if( nShiftStart < 0
1195 : || (nShiftSize < 0
1196 : && nShiftStart < static_cast<vsi_l_offset>(-nShiftSize)) )
1197 0 : nShiftStart = (nShiftSize > 0) ? 0 : -nShiftSize;
1198 :
1199 : /* get offset at end of file */
1200 12 : if( VSIFSeekL( fp, 0, SEEK_END ) != 0 )
1201 : {
1202 : CPLError( CE_Failure, CPLE_FileIO,
1203 0 : "Unable to seek to end of grid file.\n" );
1204 0 : return CE_Failure;
1205 : }
1206 :
1207 12 : vsi_l_offset nOldEnd = VSIFTellL( fp );
1208 :
1209 : /* If shifting past end, just zero-pad as necessary */
1210 12 : if( nShiftStart >= nOldEnd )
1211 : {
1212 0 : if( nShiftSize < 0 )
1213 : {
1214 0 : if( nShiftStart + nShiftSize >= nOldEnd )
1215 0 : return CE_None;
1216 :
1217 0 : if( VSIFSeekL( fp, nShiftStart + nShiftSize, SEEK_SET ) != 0 )
1218 : {
1219 : CPLError( CE_Failure, CPLE_FileIO,
1220 0 : "Unable to seek near end of file.\n" );
1221 0 : return CE_Failure;
1222 : }
1223 :
1224 : /* ftruncate()? */
1225 0 : for( vsi_l_offset nPos = nShiftStart + nShiftSize;
1226 : nPos > nOldEnd; nPos++ )
1227 : {
1228 0 : if( VSIFWriteL( (void *)" ", 1, 1, fp ) != 1 )
1229 : {
1230 : CPLError( CE_Failure, CPLE_FileIO,
1231 : "Unable to write padding to grid file "
1232 0 : "(Out of space?).\n" );
1233 0 : return CE_Failure;
1234 : }
1235 : }
1236 :
1237 0 : return CE_None;
1238 : }
1239 : else
1240 : {
1241 0 : for( vsi_l_offset nPos = nOldEnd;
1242 : nPos < nShiftStart + nShiftSize; nPos++ )
1243 : {
1244 0 : if( VSIFWriteL( (void *)" ", 1, 1, fp ) != 1 )
1245 : {
1246 : CPLError( CE_Failure, CPLE_FileIO,
1247 : "Unable to write padding to grid file "
1248 0 : "(Out of space?).\n" );
1249 0 : return CE_Failure;
1250 : }
1251 : }
1252 0 : return CE_None;
1253 : }
1254 : }
1255 :
1256 : /* prepare buffer for real shifting */
1257 12 : size_t nBufferSize = (1024 >= abs(nShiftSize)*2) ? 1024 : abs(nShiftSize)*2;
1258 12 : char *pabyBuffer = (char *)VSIMalloc( nBufferSize );
1259 12 : if( pabyBuffer == NULL)
1260 : {
1261 : CPLError( CE_Failure, CPLE_OutOfMemory,
1262 0 : "Unable to allocate space for shift buffer.\n" );
1263 0 : return CE_Failure;
1264 : }
1265 :
1266 12 : if( VSIFSeekL( fp, nShiftStart, SEEK_SET ) != 0 )
1267 : {
1268 0 : VSIFree( pabyBuffer );
1269 : CPLError( CE_Failure, CPLE_FileIO,
1270 0 : "Unable to seek to start of shift in grid file.\n" );
1271 0 : return CE_Failure;
1272 : }
1273 :
1274 : size_t nRead;
1275 12 : size_t nOverlap = (nShiftSize > 0) ? nShiftSize : 0;
1276 : /* If there is overlap, fill buffer with the overlap to start */
1277 12 : if( nOverlap > 0)
1278 : {
1279 0 : nRead = VSIFReadL( (void *)pabyBuffer, 1, nOverlap, fp );
1280 0 : if( nRead < nOverlap && !VSIFEofL( fp ) )
1281 : {
1282 0 : VSIFree( pabyBuffer );
1283 : CPLError( CE_Failure, CPLE_FileIO,
1284 0 : "Error reading grid file.\n" );
1285 0 : return CE_Failure;
1286 : }
1287 :
1288 : /* overwrite the new space with ' ' */
1289 0 : if( VSIFSeekL( fp, nShiftStart, SEEK_SET ) != 0 )
1290 : {
1291 0 : VSIFree( pabyBuffer );
1292 : CPLError( CE_Failure, CPLE_FileIO,
1293 0 : "Unable to seek to start of shift in grid file.\n" );
1294 0 : return CE_Failure;
1295 : }
1296 :
1297 0 : for( int iFill=0; iFill<nShiftSize; iFill++ )
1298 : {
1299 0 : if( VSIFWriteL( (void *)" ", 1, 1, fp ) != 1 )
1300 : {
1301 0 : VSIFree( pabyBuffer );
1302 : CPLError( CE_Failure, CPLE_FileIO,
1303 : "Unable to write padding to grid file "
1304 0 : "(Out of space?).\n" );
1305 0 : return CE_Failure;
1306 : }
1307 : }
1308 :
1309 : /* if we have already read the entire file, finish it off */
1310 0 : if( VSIFTellL( fp ) >= nOldEnd )
1311 : {
1312 0 : if( VSIFWriteL( (void *)pabyBuffer, 1, nRead, fp ) != nRead )
1313 : {
1314 0 : VSIFree( pabyBuffer );
1315 : CPLError( CE_Failure, CPLE_FileIO,
1316 0 : "Unable to write to grid file (Out of space?).\n" );
1317 0 : return CE_Failure;
1318 : }
1319 :
1320 0 : VSIFree( pabyBuffer );
1321 0 : return CE_None;
1322 : }
1323 : }
1324 :
1325 : /* iterate over the remainder of the file and shift as requested */
1326 12 : bool bEOF = false;
1327 37 : while( !bEOF )
1328 : {
1329 : nRead = VSIFReadL( (void *)(pabyBuffer+nOverlap), 1,
1330 13 : nBufferSize - nOverlap, fp );
1331 :
1332 13 : if( VSIFEofL( fp ) )
1333 12 : bEOF = true;
1334 : else
1335 1 : bEOF = false;
1336 :
1337 13 : if( nRead == 0 && !bEOF )
1338 : {
1339 0 : VSIFree( pabyBuffer );
1340 : CPLError( CE_Failure, CPLE_FileIO,
1341 0 : "Unable to read from grid file (possible corruption).\n");
1342 0 : return CE_Failure;
1343 : }
1344 :
1345 : /* FIXME: Should use SEEK_CUR, review integer promotions... */
1346 13 : if( VSIFSeekL( fp, VSIFTellL(fp)-nRead+nShiftSize-nOverlap,
1347 : SEEK_SET ) != 0 )
1348 : {
1349 0 : VSIFree( pabyBuffer );
1350 : CPLError( CE_Failure, CPLE_FileIO,
1351 0 : "Unable to seek in grid file (possible corruption).\n" );
1352 0 : return CE_Failure;
1353 : }
1354 :
1355 13 : size_t nWritten = VSIFWriteL( (void *)pabyBuffer, 1, nRead, fp );
1356 13 : if( nWritten != nRead )
1357 : {
1358 0 : VSIFree( pabyBuffer );
1359 : CPLError( CE_Failure, CPLE_FileIO,
1360 0 : "Unable to write to grid file (out of space?).\n" );
1361 0 : return CE_Failure;
1362 : }
1363 :
1364 : /* shift overlapped contents to the front of the buffer if necessary */
1365 13 : if( nOverlap > 0)
1366 0 : memmove(pabyBuffer, pabyBuffer+nRead, nOverlap);
1367 : }
1368 :
1369 : /* write the remainder of the buffer or overwrite leftovers and finish */
1370 12 : if( nShiftSize > 0 )
1371 : {
1372 0 : size_t nTailSize = nOverlap;
1373 0 : while( nTailSize > 0 && isspace( (unsigned char)pabyBuffer[nTailSize-1] ) )
1374 0 : nTailSize--;
1375 :
1376 0 : if( VSIFWriteL( (void *)pabyBuffer, 1, nTailSize, fp ) != nTailSize )
1377 : {
1378 0 : VSIFree( pabyBuffer );
1379 : CPLError( CE_Failure, CPLE_FileIO,
1380 0 : "Unable to write to grid file (out of space?).\n" );
1381 0 : return CE_Failure;
1382 : }
1383 :
1384 0 : if( VSIFWriteL( (void *)pszEOL, 1, strlen(pszEOL), fp )
1385 : != strlen(pszEOL) )
1386 : {
1387 0 : VSIFree( pabyBuffer );
1388 : CPLError( CE_Failure, CPLE_FileIO,
1389 0 : "Unable to write to grid file (out of space?).\n" );
1390 0 : return CE_Failure;
1391 : }
1392 : }
1393 : else
1394 : {
1395 : /* FIXME: ftruncate()? */
1396 : /* FIXME: Should use SEEK_CUR, review integer promotions... */
1397 12 : if( VSIFSeekL( fp, VSIFTellL(fp)-strlen(pszEOL), SEEK_SET ) != 0 )
1398 : {
1399 0 : VSIFree( pabyBuffer );
1400 : CPLError( CE_Failure, CPLE_FileIO,
1401 0 : "Unable to seek in grid file.\n" );
1402 0 : return CE_Failure;
1403 : }
1404 :
1405 301 : for( int iPadding=0; iPadding<-nShiftSize; iPadding++ )
1406 : {
1407 289 : if( VSIFWriteL( (void *)" ", 1, 1, fp ) != 1 )
1408 : {
1409 0 : VSIFree( pabyBuffer );
1410 : CPLError( CE_Failure, CPLE_FileIO,
1411 0 : "Error writing to grid file.\n" );
1412 0 : return CE_Failure;
1413 : }
1414 : }
1415 :
1416 12 : if( VSIFWriteL( (void *)pszEOL, 1, strlen(pszEOL), fp )
1417 : != strlen(pszEOL) )
1418 : {
1419 0 : VSIFree( pabyBuffer );
1420 : CPLError( CE_Failure, CPLE_FileIO,
1421 0 : "Unable to write to grid file (out of space?).\n" );
1422 0 : return CE_Failure;
1423 : }
1424 : }
1425 :
1426 12 : VSIFree( pabyBuffer );
1427 12 : return CE_None;
1428 : }
1429 :
1430 : /************************************************************************/
1431 : /* UpdateHeader() */
1432 : /************************************************************************/
1433 :
1434 0 : CPLErr GSAGDataset::UpdateHeader()
1435 :
1436 : {
1437 0 : GSAGRasterBand *poBand = (GSAGRasterBand *)GetRasterBand( 1 );
1438 0 : if( poBand == NULL )
1439 : {
1440 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to open raster band.\n" );
1441 0 : return CE_Failure;
1442 : }
1443 :
1444 0 : std::ostringstream ssOutBuf;
1445 0 : ssOutBuf.precision( nFIELD_PRECISION );
1446 0 : ssOutBuf.setf( std::ios::uppercase );
1447 :
1448 : /* signature */
1449 0 : ssOutBuf << "DSAA" << szEOL;
1450 :
1451 : /* columns rows */
1452 0 : ssOutBuf << nRasterXSize << " " << nRasterYSize << szEOL;
1453 :
1454 : /* x range */
1455 0 : ssOutBuf << poBand->dfMinX << " " << poBand->dfMaxX << szEOL;
1456 :
1457 : /* y range */
1458 0 : ssOutBuf << poBand->dfMinY << " " << poBand->dfMaxY << szEOL;
1459 :
1460 : /* z range */
1461 0 : ssOutBuf << poBand->dfMinZ << " " << poBand->dfMaxZ << szEOL;
1462 :
1463 0 : CPLString sOut = ssOutBuf.str();
1464 0 : if( sOut.length() != poBand->panLineOffset[0] )
1465 : {
1466 0 : int nShiftSize = (int) (sOut.length() - poBand->panLineOffset[0]);
1467 0 : if( ShiftFileContents( fp, poBand->panLineOffset[0], nShiftSize,
1468 : szEOL ) != CE_None )
1469 : {
1470 : CPLError( CE_Failure, CPLE_FileIO,
1471 : "Unable to update grid header, "
1472 0 : "failure shifting file contents.\n" );
1473 0 : return CE_Failure;
1474 : }
1475 :
1476 0 : for( size_t iLine=0;
1477 : iLine < static_cast<unsigned>(nRasterYSize+1)
1478 0 : && poBand->panLineOffset[iLine] != 0;
1479 : iLine++ )
1480 0 : poBand->panLineOffset[iLine] += nShiftSize;
1481 : }
1482 :
1483 0 : if( VSIFSeekL( fp, 0, SEEK_SET ) != 0 )
1484 : {
1485 : CPLError( CE_Failure, CPLE_FileIO,
1486 0 : "Unable to seek to start of grid file.\n" );
1487 0 : return CE_Failure;
1488 : }
1489 :
1490 0 : if( VSIFWriteL( sOut.c_str(), 1, sOut.length(), fp ) != sOut.length() )
1491 : {
1492 : CPLError( CE_Failure, CPLE_FileIO,
1493 0 : "Unable to update file header. Disk full?\n" );
1494 0 : return CE_Failure;
1495 : }
1496 :
1497 0 : return CE_None;
1498 : }
1499 :
1500 : /************************************************************************/
1501 : /* CreateCopy() */
1502 : /************************************************************************/
1503 :
1504 30 : GDALDataset *GSAGDataset::CreateCopy( const char *pszFilename,
1505 : GDALDataset *poSrcDS,
1506 : int bStrict, char **papszOptions,
1507 : GDALProgressFunc pfnProgress,
1508 : void *pProgressData )
1509 : {
1510 30 : if( pfnProgress == NULL )
1511 0 : pfnProgress = GDALDummyProgress;
1512 :
1513 30 : int nBands = poSrcDS->GetRasterCount();
1514 30 : if (nBands == 0)
1515 : {
1516 : CPLError( CE_Failure, CPLE_NotSupported,
1517 1 : "GSAG driver does not support source dataset with zero band.\n");
1518 1 : return NULL;
1519 : }
1520 29 : else if (nBands > 1)
1521 : {
1522 4 : if( bStrict )
1523 : {
1524 : CPLError( CE_Failure, CPLE_NotSupported,
1525 : "Unable to create copy, Golden Software ASCII Grid "
1526 4 : "format only supports one raster band.\n" );
1527 4 : return NULL;
1528 : }
1529 : else
1530 : CPLError( CE_Warning, CPLE_NotSupported,
1531 : "Golden Software ASCII Grid format only supports one "
1532 0 : "raster band, first band will be copied.\n" );
1533 : }
1534 :
1535 25 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1536 : {
1537 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated\n" );
1538 0 : return NULL;
1539 : }
1540 :
1541 25 : VSILFILE *fp = VSIFOpenL( pszFilename, "w+b" );
1542 :
1543 25 : if( fp == NULL )
1544 : {
1545 : CPLError( CE_Failure, CPLE_OpenFailed,
1546 : "Attempt to create file '%s' failed.\n",
1547 13 : pszFilename );
1548 13 : return NULL;
1549 : }
1550 :
1551 12 : int nXSize = poSrcDS->GetRasterXSize();
1552 12 : int nYSize = poSrcDS->GetRasterYSize();
1553 : double adfGeoTransform[6];
1554 :
1555 12 : poSrcDS->GetGeoTransform( adfGeoTransform );
1556 :
1557 12 : std::ostringstream ssHeader;
1558 12 : ssHeader.precision( nFIELD_PRECISION );
1559 12 : ssHeader.setf( std::ios::uppercase );
1560 :
1561 12 : ssHeader << "DSAA\x0D\x0A";
1562 :
1563 12 : ssHeader << nXSize << " " << nYSize << "\x0D\x0A";
1564 :
1565 12 : ssHeader << adfGeoTransform[0] + adfGeoTransform[1] / 2 << " "
1566 24 : << adfGeoTransform[1] * (nXSize - 0.5) + adfGeoTransform[0]
1567 36 : << "\x0D\x0A";
1568 :
1569 24 : ssHeader << adfGeoTransform[5] * (nYSize - 0.5) + adfGeoTransform[3] << " "
1570 12 : << adfGeoTransform[3] + adfGeoTransform[5] / 2
1571 36 : << "\x0D\x0A";
1572 :
1573 12 : if( VSIFWriteL( (void *)ssHeader.str().c_str(), 1, ssHeader.str().length(),
1574 : fp ) != ssHeader.str().length() )
1575 : {
1576 0 : VSIFCloseL( fp );
1577 : CPLError( CE_Failure, CPLE_FileIO,
1578 0 : "Unable to create copy, writing header failed.\n" );
1579 0 : return NULL;
1580 : }
1581 :
1582 : /* Save the location and write placeholders for the min/max Z value */
1583 12 : vsi_l_offset nRangeStart = VSIFTellL( fp );
1584 12 : const char *szDummyRange = "0.0000000000001 0.0000000000001\x0D\x0A";
1585 12 : size_t nDummyRangeLen = strlen( szDummyRange );
1586 12 : if( VSIFWriteL( (void *)szDummyRange, 1, nDummyRangeLen,
1587 : fp ) != nDummyRangeLen )
1588 : {
1589 0 : VSIFCloseL( fp );
1590 : CPLError( CE_Failure, CPLE_FileIO,
1591 0 : "Unable to create copy, writing header failed.\n" );
1592 0 : return NULL;
1593 : }
1594 :
1595 : /* -------------------------------------------------------------------- */
1596 : /* Copy band data. */
1597 : /* -------------------------------------------------------------------- */
1598 12 : double *pdfData = (double *)VSIMalloc2( nXSize, sizeof( double ) );
1599 12 : if( pdfData == NULL )
1600 : {
1601 0 : VSIFCloseL( fp );
1602 : CPLError( CE_Failure, CPLE_OutOfMemory,
1603 0 : "Unable to create copy, unable to allocate line buffer.\n" );
1604 0 : return NULL;
1605 : }
1606 :
1607 12 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
1608 : int bSrcHasNDValue;
1609 12 : double dfSrcNoDataValue = poSrcBand->GetNoDataValue( &bSrcHasNDValue );
1610 12 : double dfMin = DBL_MAX;
1611 12 : double dfMax = -DBL_MAX;
1612 142 : for( int iRow=0; iRow<nYSize; iRow++ )
1613 : {
1614 : CPLErr eErr = poSrcBand->RasterIO( GF_Read, 0, nYSize-iRow-1,
1615 : nXSize, 1, pdfData,
1616 130 : nXSize, 1, GDT_Float64, 0, 0 );
1617 :
1618 130 : if( eErr != CE_None )
1619 : {
1620 0 : VSIFCloseL( fp );
1621 0 : VSIFree( pdfData );
1622 0 : return NULL;
1623 : }
1624 :
1625 410 : for( int iCol=0; iCol<nXSize; )
1626 : {
1627 150 : for( int iCount=0;
1628 : iCount<10 && iCol<nXSize;
1629 : iCount++, iCol++ )
1630 : {
1631 1500 : double dfValue = pdfData[iCol];
1632 :
1633 1500 : if( bSrcHasNDValue && AlmostEqual( dfValue, dfSrcNoDataValue ) )
1634 : {
1635 0 : dfValue = dfNODATA_VALUE;
1636 : }
1637 : else
1638 : {
1639 1500 : if( dfValue > dfMax )
1640 14 : dfMax = dfValue;
1641 :
1642 1500 : if( dfValue < dfMin )
1643 19 : dfMin = dfValue;
1644 : }
1645 :
1646 1500 : std::ostringstream ssOut;
1647 1500 : ssOut.precision(nFIELD_PRECISION);
1648 1500 : ssOut.setf( std::ios::uppercase );
1649 1500 : ssOut << dfValue << " ";
1650 1500 : CPLString sOut = ssOut.str();
1651 :
1652 1500 : if( VSIFWriteL( sOut.c_str(), 1, sOut.length(), fp )
1653 : != sOut.length() )
1654 : {
1655 0 : VSIFCloseL( fp );
1656 0 : VSIFree( pdfData );
1657 : CPLError( CE_Failure, CPLE_FileIO,
1658 0 : "Unable to write grid cell. Disk full?\n" );
1659 0 : return NULL;
1660 : }
1661 : }
1662 :
1663 150 : if( VSIFWriteL( (void *)"\x0D\x0A", 1, 2, fp ) != 2 )
1664 : {
1665 0 : VSIFCloseL( fp );
1666 0 : VSIFree( pdfData );
1667 : CPLError( CE_Failure, CPLE_FileIO,
1668 0 : "Unable to finish write of grid line. Disk full?\n" );
1669 0 : return NULL;
1670 : }
1671 : }
1672 :
1673 130 : if( VSIFWriteL( (void *)"\x0D\x0A", 1, 2, fp ) != 2 )
1674 : {
1675 0 : VSIFCloseL( fp );
1676 0 : VSIFree( pdfData );
1677 : CPLError( CE_Failure, CPLE_FileIO,
1678 0 : "Unable to finish write of grid row. Disk full?\n" );
1679 0 : return NULL;
1680 : }
1681 :
1682 130 : if( !pfnProgress( static_cast<double>(iRow + 1)/nYSize,
1683 : NULL, pProgressData ) )
1684 : {
1685 0 : VSIFCloseL( fp );
1686 0 : VSIFree( pdfData );
1687 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1688 0 : return NULL;
1689 : }
1690 : }
1691 :
1692 12 : VSIFree( pdfData );
1693 :
1694 : /* write out the min and max values */
1695 12 : std::ostringstream ssRange;
1696 12 : ssRange.precision( nFIELD_PRECISION );
1697 12 : ssRange.setf( std::ios::uppercase );
1698 12 : ssRange << dfMin << " " << dfMax << "\x0D\x0A";
1699 12 : if( ssRange.str().length() != nDummyRangeLen )
1700 : {
1701 12 : int nShiftSize = ssRange.str().length() - nDummyRangeLen;
1702 12 : if( ShiftFileContents( fp, nRangeStart + nDummyRangeLen,
1703 : nShiftSize, "\x0D\x0A" ) != CE_None )
1704 : {
1705 0 : VSIFCloseL( fp );
1706 : CPLError( CE_Failure, CPLE_FileIO,
1707 0 : "Unable to shift file contents.\n" );
1708 0 : return NULL;
1709 : }
1710 : }
1711 :
1712 12 : if( VSIFSeekL( fp, nRangeStart, SEEK_SET ) != 0 )
1713 : {
1714 0 : VSIFCloseL( fp );
1715 : CPLError( CE_Failure, CPLE_FileIO,
1716 0 : "Unable to seek to start of grid file copy.\n" );
1717 0 : return NULL;
1718 : }
1719 :
1720 12 : if( VSIFWriteL( (void *)ssRange.str().c_str(), 1, ssRange.str().length(),
1721 : fp ) != ssRange.str().length() )
1722 : {
1723 0 : VSIFCloseL( fp );
1724 : CPLError( CE_Failure, CPLE_FileIO,
1725 0 : "Unable to write range information.\n" );
1726 0 : return NULL;
1727 : }
1728 :
1729 12 : VSIFCloseL( fp );
1730 :
1731 : GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen( pszFilename,
1732 12 : GA_Update );
1733 12 : if (poDS)
1734 : {
1735 12 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1736 : }
1737 12 : return poDS;
1738 : }
1739 :
1740 : /************************************************************************/
1741 : /* GDALRegister_GSAG() */
1742 : /************************************************************************/
1743 :
1744 582 : void GDALRegister_GSAG()
1745 :
1746 : {
1747 : GDALDriver *poDriver;
1748 :
1749 582 : if( GDALGetDriverByName( "GSAG" ) == NULL )
1750 : {
1751 561 : poDriver = new GDALDriver();
1752 :
1753 561 : poDriver->SetDescription( "GSAG" );
1754 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1755 561 : "Golden Software ASCII Grid (.grd)" );
1756 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1757 561 : "frmt_various.html#GSAG" );
1758 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "grd" );
1759 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1760 : "Byte Int16 UInt16 Int32 UInt32 "
1761 561 : "Float32 Float64" );
1762 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1763 :
1764 561 : poDriver->pfnIdentify = GSAGDataset::Identify;
1765 561 : poDriver->pfnOpen = GSAGDataset::Open;
1766 561 : poDriver->pfnCreateCopy = GSAGDataset::CreateCopy;
1767 :
1768 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1769 : }
1770 582 : }
|
