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