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