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