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