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, VSIFTellL(poGDS->fp) + szStart - szEnd,
458 0 : SEEK_SET) != 0)
459 : {
460 0 : VSIFree(szLineBuf);
461 0 : CPLError(CE_Failure, CPLE_FileIO,
462 : "Unable to seek in grid row %d (offset %ld, seek %d)"
463 : ".\n",
464 0 : nBlockYOff, (long)VSIFTellL(poGDS->fp),
465 0 : (int)(szStart - szEnd));
466 :
467 0 : return CE_Failure;
468 : }
469 0 : nCharsExamined += szStart - szLineBuf;
470 0 : nCharsRead = VSIFReadL(szLineBuf, 1, nLineBufSize - 1, poGDS->fp);
471 0 : szLineBuf[nCharsRead] = '\0';
472 :
473 0 : if (nCharsRead == 0)
474 : {
475 0 : VSIFree(szLineBuf);
476 0 : CPLError(CE_Failure, CPLE_FileIO,
477 : "Can't read portion of grid row %d at offset %ld.",
478 0 : nBlockYOff, (long)panLineOffset[nBlockYOff]);
479 0 : return CE_Failure;
480 : }
481 0 : else if (nCharsRead > static_cast<size_t>(szEnd - szStart))
482 : {
483 : /* Read new data, this was not really the end */
484 0 : szEnd = szLineBuf;
485 0 : continue;
486 : }
487 :
488 : /* This is really the last value and has no tailing newline */
489 0 : szEnd = szLineBuf + nCharsRead;
490 : }
491 :
492 3120 : if (pdfImage != nullptr)
493 : {
494 1600 : *(pdfImage + iCell) = dfValue;
495 : }
496 :
497 3120 : iCell++;
498 : }
499 :
500 312 : while (*szEnd == ' ')
501 156 : szEnd++;
502 :
503 156 : if (*szEnd != '\0' && *szEnd != poGDS->szEOL[0])
504 0 : CPLDebug("GSAG",
505 : "Grid row %d does not end with a newline. "
506 : "Possible skew.\n",
507 : nBlockYOff);
508 :
509 830 : while (isspace((unsigned char)*szEnd))
510 674 : szEnd++;
511 :
512 156 : nCharsExamined += szEnd - szLineBuf;
513 :
514 156 : if (nCharsExamined >= nMaxLineSize)
515 0 : nMaxLineSize = nCharsExamined + 1;
516 :
517 156 : if (nBlockYOff > 0)
518 : {
519 152 : vsi_l_offset nNewOffset = panLineOffset[nBlockYOff] + nCharsExamined;
520 152 : if (panLineOffset[nBlockYOff - 1] == 0)
521 : {
522 76 : panLineOffset[nBlockYOff - 1] = nNewOffset;
523 : }
524 76 : else if (panLineOffset[nBlockYOff - 1] != nNewOffset)
525 : {
526 0 : CPLError(
527 : CE_Failure, CPLE_AppDefined,
528 : "Coding error: previous offset for line %d was " CPL_FRMT_GUIB
529 : ", new offset would be " CPL_FRMT_GUIB,
530 0 : nBlockYOff - 1, panLineOffset[nBlockYOff - 1], nNewOffset);
531 : }
532 : }
533 :
534 156 : nLastReadLine = nBlockYOff;
535 :
536 156 : VSIFree(szLineBuf);
537 :
538 156 : return CE_None;
539 : }
540 :
541 : /************************************************************************/
542 : /* IWriteBlock() */
543 : /************************************************************************/
544 :
545 0 : CPLErr GSAGRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)
546 :
547 : {
548 0 : if (eAccess == GA_ReadOnly)
549 : {
550 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
551 : "Unable to write block, dataset opened read only.\n");
552 0 : return CE_Failure;
553 : }
554 :
555 0 : if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
556 0 : return CE_Failure;
557 :
558 0 : GSAGDataset *poGDS = cpl::down_cast<GSAGDataset *>(poDS);
559 :
560 0 : if (padfRowMinZ == nullptr || padfRowMaxZ == nullptr || nMinZRow < 0 ||
561 0 : nMaxZRow < 0)
562 : {
563 0 : padfRowMinZ =
564 0 : (double *)VSI_MALLOC2_VERBOSE(nRasterYSize, sizeof(double));
565 0 : if (padfRowMinZ == nullptr)
566 : {
567 0 : return CE_Failure;
568 : }
569 :
570 0 : padfRowMaxZ =
571 0 : (double *)VSI_MALLOC2_VERBOSE(nRasterYSize, sizeof(double));
572 0 : if (padfRowMaxZ == nullptr)
573 : {
574 0 : VSIFree(padfRowMinZ);
575 0 : padfRowMinZ = nullptr;
576 0 : return CE_Failure;
577 : }
578 :
579 0 : CPLErr eErr = ScanForMinMaxZ();
580 0 : if (eErr != CE_None)
581 0 : return eErr;
582 : }
583 :
584 0 : if (panLineOffset[nBlockYOff + 1] == 0)
585 0 : IReadBlock(nBlockXOff, nBlockYOff, nullptr);
586 :
587 0 : if (panLineOffset[nBlockYOff + 1] == 0 || panLineOffset[nBlockYOff] == 0)
588 0 : return CE_Failure;
589 :
590 0 : std::ostringstream ssOutBuf;
591 0 : ssOutBuf.precision(GSAGDataset::nFIELD_PRECISION);
592 0 : ssOutBuf.setf(std::ios::uppercase);
593 :
594 0 : double *pdfImage = (double *)pImage;
595 0 : padfRowMinZ[nBlockYOff] = std::numeric_limits<double>::max();
596 0 : padfRowMaxZ[nBlockYOff] = std::numeric_limits<double>::lowest();
597 0 : for (int iCell = 0; iCell < nBlockXSize;)
598 : {
599 0 : for (int iCol = 0; iCol < 10 && iCell < nBlockXSize; iCol++, iCell++)
600 : {
601 0 : if (AlmostEqual(pdfImage[iCell], GSAGDataset::dfNODATA_VALUE))
602 : {
603 0 : if (pdfImage[iCell] < padfRowMinZ[nBlockYOff])
604 0 : padfRowMinZ[nBlockYOff] = pdfImage[iCell];
605 :
606 0 : if (pdfImage[iCell] > padfRowMaxZ[nBlockYOff])
607 0 : padfRowMaxZ[nBlockYOff] = pdfImage[iCell];
608 : }
609 :
610 0 : ssOutBuf << pdfImage[iCell] << " ";
611 : }
612 0 : ssOutBuf << poGDS->szEOL;
613 : }
614 0 : ssOutBuf << poGDS->szEOL;
615 :
616 0 : CPLString sOut = ssOutBuf.str();
617 0 : if (sOut.length() !=
618 0 : panLineOffset[nBlockYOff + 1] - panLineOffset[nBlockYOff])
619 : {
620 0 : int nShiftSize = (int)(sOut.length() - (panLineOffset[nBlockYOff + 1] -
621 0 : panLineOffset[nBlockYOff]));
622 0 : if (nBlockYOff != poGDS->nRasterYSize &&
623 0 : GSAGDataset::ShiftFileContents(poGDS->fp,
624 0 : panLineOffset[nBlockYOff + 1],
625 0 : nShiftSize, poGDS->szEOL) != CE_None)
626 : {
627 0 : CPLError(CE_Failure, CPLE_FileIO,
628 : "Failure writing block, "
629 : "unable to shift file contents.\n");
630 0 : return CE_Failure;
631 : }
632 :
633 0 : for (size_t iLine = nBlockYOff + 1;
634 0 : iLine < static_cast<unsigned>(poGDS->nRasterYSize + 1) &&
635 0 : panLineOffset[iLine] != 0;
636 : iLine++)
637 0 : panLineOffset[iLine] += nShiftSize;
638 : }
639 :
640 0 : if (VSIFSeekL(poGDS->fp, panLineOffset[nBlockYOff], SEEK_SET) != 0)
641 : {
642 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to seek to grid line.\n");
643 0 : return CE_Failure;
644 : }
645 :
646 0 : if (VSIFWriteL(sOut.c_str(), 1, sOut.length(), poGDS->fp) != sOut.length())
647 : {
648 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to write grid block.\n");
649 0 : return CE_Failure;
650 : }
651 :
652 : /* Update header as needed */
653 0 : bool bHeaderNeedsUpdate = false;
654 0 : if (nMinZRow == nBlockYOff && padfRowMinZ[nBlockYOff] > dfMinZ)
655 : {
656 0 : double dfNewMinZ = std::numeric_limits<double>::lowest();
657 0 : for (int iRow = 0; iRow < nRasterYSize; iRow++)
658 : {
659 0 : if (padfRowMinZ[iRow] < dfNewMinZ)
660 : {
661 0 : dfNewMinZ = padfRowMinZ[iRow];
662 0 : nMinZRow = iRow;
663 : }
664 : }
665 :
666 0 : if (dfNewMinZ != dfMinZ)
667 : {
668 0 : dfMinZ = dfNewMinZ;
669 0 : bHeaderNeedsUpdate = true;
670 : }
671 : }
672 :
673 0 : if (nMaxZRow == nBlockYOff && padfRowMaxZ[nBlockYOff] < dfMaxZ)
674 : {
675 0 : double dfNewMaxZ = std::numeric_limits<double>::lowest();
676 0 : for (int iRow = 0; iRow < nRasterYSize; iRow++)
677 : {
678 0 : if (padfRowMaxZ[iRow] > dfNewMaxZ)
679 : {
680 0 : dfNewMaxZ = padfRowMaxZ[iRow];
681 0 : nMaxZRow = iRow;
682 : }
683 : }
684 :
685 0 : if (dfNewMaxZ != dfMaxZ)
686 : {
687 0 : dfMaxZ = dfNewMaxZ;
688 0 : bHeaderNeedsUpdate = true;
689 : }
690 : }
691 :
692 0 : if (padfRowMinZ[nBlockYOff] < dfMinZ || padfRowMaxZ[nBlockYOff] > dfMaxZ)
693 : {
694 0 : if (padfRowMinZ[nBlockYOff] < dfMinZ)
695 : {
696 0 : dfMinZ = padfRowMinZ[nBlockYOff];
697 0 : nMinZRow = nBlockYOff;
698 : }
699 :
700 0 : if (padfRowMaxZ[nBlockYOff] > dfMaxZ)
701 : {
702 0 : dfMaxZ = padfRowMaxZ[nBlockYOff];
703 0 : nMaxZRow = nBlockYOff;
704 : }
705 :
706 0 : bHeaderNeedsUpdate = true;
707 : }
708 :
709 0 : if (bHeaderNeedsUpdate && dfMaxZ > dfMinZ)
710 : {
711 0 : CPLErr eErr = poGDS->UpdateHeader();
712 0 : return eErr;
713 : }
714 :
715 0 : return CE_None;
716 : }
717 :
718 : /************************************************************************/
719 : /* GetNoDataValue() */
720 : /************************************************************************/
721 :
722 8 : double GSAGRasterBand::GetNoDataValue(int *pbSuccess)
723 : {
724 8 : if (pbSuccess)
725 7 : *pbSuccess = TRUE;
726 :
727 8 : return GSAGDataset::dfNODATA_VALUE;
728 : }
729 :
730 : /************************************************************************/
731 : /* GetMinimum() */
732 : /************************************************************************/
733 :
734 0 : double GSAGRasterBand::GetMinimum(int *pbSuccess)
735 : {
736 0 : if (pbSuccess)
737 0 : *pbSuccess = TRUE;
738 :
739 0 : return dfMinZ;
740 : }
741 :
742 : /************************************************************************/
743 : /* GetMaximum() */
744 : /************************************************************************/
745 :
746 0 : double GSAGRasterBand::GetMaximum(int *pbSuccess)
747 : {
748 0 : if (pbSuccess)
749 0 : *pbSuccess = TRUE;
750 :
751 0 : return dfMaxZ;
752 : }
753 :
754 : /************************************************************************/
755 : /* ==================================================================== */
756 : /* GSAGDataset */
757 : /* ==================================================================== */
758 : /************************************************************************/
759 :
760 : /************************************************************************/
761 : /* GSAGDataset() */
762 : /************************************************************************/
763 :
764 17 : GSAGDataset::GSAGDataset(const char *pszEOL) : fp(nullptr), nMinMaxZOffset(0)
765 : {
766 17 : if (pszEOL == nullptr || EQUAL(pszEOL, ""))
767 : {
768 0 : CPLDebug("GSAG", "GSAGDataset() created with invalid EOL string.\n");
769 0 : szEOL[0] = '\x0D';
770 0 : szEOL[1] = '\x0A';
771 0 : szEOL[2] = '\0';
772 0 : return;
773 : }
774 :
775 17 : snprintf(szEOL, sizeof(szEOL), "%s", pszEOL);
776 : }
777 :
778 : /************************************************************************/
779 : /* ~GSAGDataset() */
780 : /************************************************************************/
781 :
782 34 : GSAGDataset::~GSAGDataset()
783 :
784 : {
785 17 : FlushCache(true);
786 17 : if (fp != nullptr)
787 17 : VSIFCloseL(fp);
788 34 : }
789 :
790 : /************************************************************************/
791 : /* Identify() */
792 : /************************************************************************/
793 :
794 60281 : int GSAGDataset::Identify(GDALOpenInfo *poOpenInfo)
795 :
796 : {
797 : /* Check for signature */
798 60281 : if (poOpenInfo->nHeaderBytes < 5 ||
799 4651 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "DSAA") ||
800 34 : (poOpenInfo->pabyHeader[4] != '\x0D' &&
801 0 : poOpenInfo->pabyHeader[4] != '\x0A'))
802 : {
803 60247 : return FALSE;
804 : }
805 34 : return TRUE;
806 : }
807 :
808 : /************************************************************************/
809 : /* Open() */
810 : /************************************************************************/
811 :
812 17 : GDALDataset *GSAGDataset::Open(GDALOpenInfo *poOpenInfo)
813 :
814 : {
815 17 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
816 : {
817 0 : return nullptr;
818 : }
819 :
820 : /* Identify the end of line marker (should be \x0D\x0A, but try some others)
821 : * (note that '\x0D' == '\r' and '\x0A' == '\n' on most systems) */
822 : char szEOL[3];
823 17 : szEOL[0] = poOpenInfo->pabyHeader[4];
824 17 : szEOL[1] = poOpenInfo->pabyHeader[5];
825 17 : szEOL[2] = '\0';
826 17 : if (szEOL[1] != '\x0D' && szEOL[1] != '\x0A')
827 0 : szEOL[1] = '\0';
828 :
829 : /* -------------------------------------------------------------------- */
830 : /* Create a corresponding GDALDataset. */
831 : /* -------------------------------------------------------------------- */
832 17 : GSAGDataset *poDS = new GSAGDataset(szEOL);
833 17 : poDS->eAccess = poOpenInfo->eAccess;
834 17 : poDS->fp = poOpenInfo->fpL;
835 17 : poOpenInfo->fpL = nullptr;
836 :
837 : /* -------------------------------------------------------------------- */
838 : /* Read the header. */
839 : /* -------------------------------------------------------------------- */
840 17 : char *pabyHeader = nullptr;
841 17 : bool bMustFreeHeader = false;
842 17 : if (poOpenInfo->nHeaderBytes >= static_cast<int>(nMAX_HEADER_SIZE))
843 : {
844 17 : pabyHeader = (char *)poOpenInfo->pabyHeader;
845 : }
846 : else
847 : {
848 0 : bMustFreeHeader = true;
849 0 : pabyHeader = (char *)VSI_MALLOC_VERBOSE(nMAX_HEADER_SIZE);
850 0 : if (pabyHeader == nullptr)
851 : {
852 0 : delete poDS;
853 0 : return nullptr;
854 : }
855 :
856 0 : size_t nRead = VSIFReadL(pabyHeader, 1, nMAX_HEADER_SIZE - 1, poDS->fp);
857 0 : pabyHeader[nRead] = '\0';
858 : }
859 :
860 17 : const char *szErrorMsg = nullptr;
861 17 : const char *szStart = pabyHeader + 5;
862 17 : char *szEnd = nullptr;
863 : double dfTemp;
864 :
865 : /* Parse number of X axis grid rows */
866 17 : long nTemp = strtol(szStart, &szEnd, 10);
867 17 : if (szStart == szEnd || nTemp < 0l)
868 : {
869 0 : szErrorMsg = "Unable to parse the number of X axis grid columns.\n";
870 0 : goto error;
871 : }
872 17 : else if (nTemp > std::numeric_limits<int>::max())
873 : {
874 0 : CPLError(CE_Warning, CPLE_AppDefined,
875 : "Number of X axis grid columns not representable.\n");
876 0 : poDS->nRasterXSize = std::numeric_limits<int>::max();
877 : }
878 17 : else if (nTemp == 0)
879 : {
880 0 : szErrorMsg =
881 : "Number of X axis grid columns is zero, which is invalid.\n";
882 0 : goto error;
883 : }
884 : else
885 : {
886 17 : poDS->nRasterXSize = static_cast<int>(nTemp);
887 : }
888 17 : szStart = szEnd;
889 :
890 : /* Parse number of Y axis grid rows */
891 17 : nTemp = strtol(szStart, &szEnd, 10);
892 17 : if (szStart == szEnd || nTemp < 0l)
893 : {
894 0 : szErrorMsg = "Unable to parse the number of Y axis grid rows.\n";
895 0 : goto error;
896 : }
897 17 : else if (nTemp > std::numeric_limits<int>::max() - 1)
898 : {
899 0 : CPLError(CE_Warning, CPLE_AppDefined,
900 : "Number of Y axis grid rows not representable.\n");
901 0 : poDS->nRasterYSize = std::numeric_limits<int>::max() - 1;
902 : }
903 17 : else if (nTemp == 0)
904 : {
905 0 : szErrorMsg = "Number of Y axis grid rows is zero, which is invalid.\n";
906 0 : goto error;
907 : }
908 : else
909 : {
910 17 : poDS->nRasterYSize = static_cast<int>(nTemp);
911 : }
912 17 : szStart = szEnd;
913 :
914 : /* Parse the minimum X value of the grid */
915 : double dfMinX;
916 17 : dfTemp = CPLStrtod(szStart, &szEnd);
917 17 : if (szStart == szEnd)
918 : {
919 0 : szErrorMsg = "Unable to parse the minimum X value.\n";
920 0 : goto error;
921 : }
922 : else
923 : {
924 17 : dfMinX = dfTemp;
925 : }
926 17 : szStart = szEnd;
927 :
928 : /* Parse the maximum X value of the grid */
929 : double dfMaxX;
930 17 : dfTemp = CPLStrtod(szStart, &szEnd);
931 17 : if (szStart == szEnd)
932 : {
933 0 : szErrorMsg = "Unable to parse the maximum X value.\n";
934 0 : goto error;
935 : }
936 : else
937 : {
938 17 : dfMaxX = dfTemp;
939 : }
940 17 : szStart = szEnd;
941 :
942 : /* Parse the minimum Y value of the grid */
943 : double dfMinY;
944 17 : dfTemp = CPLStrtod(szStart, &szEnd);
945 17 : if (szStart == szEnd)
946 : {
947 0 : szErrorMsg = "Unable to parse the minimum Y value.\n";
948 0 : goto error;
949 : }
950 : else
951 : {
952 17 : dfMinY = dfTemp;
953 : }
954 17 : szStart = szEnd;
955 :
956 : /* Parse the maximum Y value of the grid */
957 : double dfMaxY;
958 17 : dfTemp = CPLStrtod(szStart, &szEnd);
959 17 : if (szStart == szEnd)
960 : {
961 0 : szErrorMsg = "Unable to parse the maximum Y value.\n";
962 0 : goto error;
963 : }
964 : else
965 : {
966 17 : dfMaxY = dfTemp;
967 : }
968 17 : szStart = szEnd;
969 :
970 : /* Parse the minimum Z value of the grid */
971 51 : while (isspace((unsigned char)*szStart))
972 34 : szStart++;
973 17 : poDS->nMinMaxZOffset = szStart - pabyHeader;
974 :
975 : double dfMinZ;
976 17 : dfTemp = CPLStrtod(szStart, &szEnd);
977 17 : if (szStart == szEnd)
978 : {
979 0 : szErrorMsg = "Unable to parse the minimum Z value.\n";
980 0 : goto error;
981 : }
982 : else
983 : {
984 17 : dfMinZ = dfTemp;
985 : }
986 17 : szStart = szEnd;
987 :
988 : /* Parse the maximum Z value of the grid */
989 : double dfMaxZ;
990 17 : dfTemp = CPLStrtod(szStart, &szEnd);
991 17 : if (szStart == szEnd)
992 : {
993 0 : szErrorMsg = "Unable to parse the maximum Z value.\n";
994 0 : goto error;
995 : }
996 : else
997 : {
998 17 : dfMaxZ = dfTemp;
999 : }
1000 :
1001 51 : while (isspace((unsigned char)*szEnd))
1002 34 : szEnd++;
1003 :
1004 : /* -------------------------------------------------------------------- */
1005 : /* Create band information objects. */
1006 : /* -------------------------------------------------------------------- */
1007 : {
1008 : GSAGRasterBand *poBand =
1009 17 : new GSAGRasterBand(poDS, 1, szEnd - pabyHeader);
1010 17 : if (poBand->panLineOffset == nullptr)
1011 : {
1012 0 : delete poBand;
1013 0 : goto error;
1014 : }
1015 :
1016 17 : poBand->dfMinX = dfMinX;
1017 17 : poBand->dfMaxX = dfMaxX;
1018 17 : poBand->dfMinY = dfMinY;
1019 17 : poBand->dfMaxY = dfMaxY;
1020 17 : poBand->dfMinZ = dfMinZ;
1021 17 : poBand->dfMaxZ = dfMaxZ;
1022 :
1023 17 : poDS->SetBand(1, poBand);
1024 : }
1025 :
1026 17 : if (bMustFreeHeader)
1027 : {
1028 0 : CPLFree(pabyHeader);
1029 : }
1030 :
1031 : /* -------------------------------------------------------------------- */
1032 : /* Initialize any PAM information. */
1033 : /* -------------------------------------------------------------------- */
1034 17 : poDS->SetDescription(poOpenInfo->pszFilename);
1035 17 : poDS->TryLoadXML();
1036 :
1037 : /* -------------------------------------------------------------------- */
1038 : /* Check for external overviews. */
1039 : /* -------------------------------------------------------------------- */
1040 34 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
1041 17 : poOpenInfo->GetSiblingFiles());
1042 :
1043 17 : return poDS;
1044 :
1045 0 : error:
1046 0 : if (bMustFreeHeader)
1047 : {
1048 0 : CPLFree(pabyHeader);
1049 : }
1050 :
1051 0 : delete poDS;
1052 :
1053 0 : if (szErrorMsg)
1054 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", szErrorMsg);
1055 0 : return nullptr;
1056 : }
1057 :
1058 : /************************************************************************/
1059 : /* GetGeoTransform() */
1060 : /************************************************************************/
1061 :
1062 17 : CPLErr GSAGDataset::GetGeoTransform(GDALGeoTransform >) const
1063 : {
1064 17 : gt = GDALGeoTransform();
1065 :
1066 : const GSAGRasterBand *poGRB =
1067 17 : cpl::down_cast<const GSAGRasterBand *>(GetRasterBand(1));
1068 :
1069 : /* check if we have a PAM GeoTransform stored */
1070 17 : CPLPushErrorHandler(CPLQuietErrorHandler);
1071 17 : CPLErr eErr = GDALPamDataset::GetGeoTransform(gt);
1072 17 : CPLPopErrorHandler();
1073 :
1074 17 : if (eErr == CE_None)
1075 0 : return CE_None;
1076 :
1077 17 : if (nRasterXSize == 1 || nRasterYSize == 1)
1078 0 : return CE_Failure;
1079 :
1080 : /* calculate pixel size first */
1081 17 : gt[1] = (poGRB->dfMaxX - poGRB->dfMinX) / (nRasterXSize - 1);
1082 17 : gt[5] = (poGRB->dfMinY - poGRB->dfMaxY) / (nRasterYSize - 1);
1083 :
1084 : /* then calculate image origin */
1085 17 : gt[0] = poGRB->dfMinX - gt[1] / 2;
1086 17 : gt[3] = poGRB->dfMaxY - gt[5] / 2;
1087 :
1088 : /* tilt/rotation does not supported by the GS grids */
1089 17 : gt[4] = 0.0;
1090 17 : gt[2] = 0.0;
1091 :
1092 17 : return CE_None;
1093 : }
1094 :
1095 : /************************************************************************/
1096 : /* SetGeoTransform() */
1097 : /************************************************************************/
1098 :
1099 0 : CPLErr GSAGDataset::SetGeoTransform(const GDALGeoTransform >)
1100 : {
1101 0 : if (eAccess == GA_ReadOnly)
1102 : {
1103 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
1104 : "Unable to set GeoTransform, dataset opened read only.\n");
1105 0 : return CE_Failure;
1106 : }
1107 :
1108 0 : GSAGRasterBand *poGRB = cpl::down_cast<GSAGRasterBand *>(GetRasterBand(1));
1109 :
1110 : /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
1111 : /*if( gt[2] != 0.0 || gt[4] != 0.0
1112 : || gt[1] < 0.0 || gt[5] < 0.0 )
1113 : eErr = GDALPamDataset::SetGeoTransform( gt );*/
1114 : // if( eErr != CE_None )
1115 : // return eErr;
1116 :
1117 0 : const double dfOldMinX = poGRB->dfMinX;
1118 0 : const double dfOldMaxX = poGRB->dfMaxX;
1119 0 : const double dfOldMinY = poGRB->dfMinY;
1120 0 : const double dfOldMaxY = poGRB->dfMaxY;
1121 :
1122 0 : poGRB->dfMinX = gt[0] + gt[1] / 2;
1123 0 : poGRB->dfMaxX = gt[1] * (nRasterXSize - 0.5) + gt[0];
1124 0 : poGRB->dfMinY = gt[5] * (nRasterYSize - 0.5) + gt[3];
1125 0 : poGRB->dfMaxY = gt[3] + gt[5] / 2;
1126 :
1127 0 : CPLErr eErr = UpdateHeader();
1128 :
1129 0 : if (eErr != CE_None)
1130 : {
1131 0 : poGRB->dfMinX = dfOldMinX;
1132 0 : poGRB->dfMaxX = dfOldMaxX;
1133 0 : poGRB->dfMinY = dfOldMinY;
1134 0 : poGRB->dfMaxY = dfOldMaxY;
1135 : }
1136 :
1137 0 : return eErr;
1138 : }
1139 :
1140 : /************************************************************************/
1141 : /* ShiftFileContents() */
1142 : /************************************************************************/
1143 12 : CPLErr GSAGDataset::ShiftFileContents(VSILFILE *fp, vsi_l_offset nShiftStart,
1144 : int nShiftSize, const char *pszEOL)
1145 : {
1146 : /* nothing to do for zero-shift */
1147 12 : if (nShiftSize == 0)
1148 0 : return CE_None;
1149 :
1150 : /* make sure start location is sane */
1151 : /* Tautology is always false. nShiftStart is unsigned. */
1152 12 : if (/* nShiftStart < 0
1153 : || */
1154 12 : (nShiftSize < 0 &&
1155 12 : nShiftStart < static_cast<vsi_l_offset>(-nShiftSize)))
1156 0 : nShiftStart = /*(nShiftSize > 0) ? 0 :*/ -nShiftSize;
1157 :
1158 : /* get offset at end of file */
1159 12 : if (VSIFSeekL(fp, 0, SEEK_END) != 0)
1160 : {
1161 0 : CPLError(CE_Failure, CPLE_FileIO,
1162 : "Unable to seek to end of grid file.\n");
1163 0 : return CE_Failure;
1164 : }
1165 :
1166 12 : vsi_l_offset nOldEnd = VSIFTellL(fp);
1167 :
1168 : /* If shifting past end, just zero-pad as necessary */
1169 12 : if (nShiftStart >= nOldEnd)
1170 : {
1171 0 : if (nShiftSize < 0)
1172 : {
1173 0 : if (nShiftStart + nShiftSize >= nOldEnd)
1174 0 : return CE_None;
1175 :
1176 0 : VSIFTruncateL(fp, nShiftStart + nShiftSize);
1177 :
1178 0 : return CE_None;
1179 : }
1180 : else
1181 : {
1182 0 : for (vsi_l_offset nPos = nOldEnd; nPos < nShiftStart + nShiftSize;
1183 : nPos++)
1184 : {
1185 0 : if (VSIFWriteL((void *)" ", 1, 1, fp) != 1)
1186 : {
1187 0 : CPLError(CE_Failure, CPLE_FileIO,
1188 : "Unable to write padding to grid file "
1189 : "(Out of space?).\n");
1190 0 : return CE_Failure;
1191 : }
1192 : }
1193 0 : return CE_None;
1194 : }
1195 : }
1196 :
1197 : /* prepare buffer for real shifting */
1198 12 : size_t nBufferSize =
1199 12 : (1024 >= abs(nShiftSize) * 2) ? 1024 : abs(nShiftSize) * 2;
1200 12 : char *pabyBuffer = (char *)VSI_MALLOC_VERBOSE(nBufferSize);
1201 12 : if (pabyBuffer == nullptr)
1202 : {
1203 0 : return CE_Failure;
1204 : }
1205 :
1206 12 : if (VSIFSeekL(fp, nShiftStart, SEEK_SET) != 0)
1207 : {
1208 0 : VSIFree(pabyBuffer);
1209 0 : CPLError(CE_Failure, CPLE_FileIO,
1210 : "Unable to seek to start of shift in grid file.\n");
1211 0 : return CE_Failure;
1212 : }
1213 :
1214 : size_t nRead;
1215 12 : size_t nOverlap = (nShiftSize > 0) ? nShiftSize : 0;
1216 : /* If there is overlap, fill buffer with the overlap to start */
1217 12 : if (nOverlap > 0)
1218 : {
1219 0 : nRead = VSIFReadL((void *)pabyBuffer, 1, nOverlap, fp);
1220 0 : if (nRead < nOverlap && !VSIFEofL(fp))
1221 : {
1222 0 : VSIFree(pabyBuffer);
1223 0 : CPLError(CE_Failure, CPLE_FileIO, "Error reading grid file.\n");
1224 0 : return CE_Failure;
1225 : }
1226 :
1227 : /* overwrite the new space with ' ' */
1228 0 : if (VSIFSeekL(fp, nShiftStart, SEEK_SET) != 0)
1229 : {
1230 0 : VSIFree(pabyBuffer);
1231 0 : CPLError(CE_Failure, CPLE_FileIO,
1232 : "Unable to seek to start of shift in grid file.\n");
1233 0 : return CE_Failure;
1234 : }
1235 :
1236 0 : for (int iFill = 0; iFill < nShiftSize; iFill++)
1237 : {
1238 0 : if (VSIFWriteL((void *)" ", 1, 1, fp) != 1)
1239 : {
1240 0 : VSIFree(pabyBuffer);
1241 0 : CPLError(CE_Failure, CPLE_FileIO,
1242 : "Unable to write padding to grid file "
1243 : "(Out of space?).\n");
1244 0 : return CE_Failure;
1245 : }
1246 : }
1247 :
1248 : /* if we have already read the entire file, finish it off */
1249 0 : if (VSIFTellL(fp) >= nOldEnd)
1250 : {
1251 0 : if (VSIFWriteL((void *)pabyBuffer, 1, nRead, fp) != nRead)
1252 : {
1253 0 : VSIFree(pabyBuffer);
1254 0 : CPLError(CE_Failure, CPLE_FileIO,
1255 : "Unable to write to grid file (Out of space?).\n");
1256 0 : return CE_Failure;
1257 : }
1258 :
1259 0 : VSIFree(pabyBuffer);
1260 0 : return CE_None;
1261 : }
1262 : }
1263 :
1264 : /* iterate over the remainder of the file and shift as requested */
1265 12 : bool bEOF = false;
1266 25 : while (!bEOF)
1267 : {
1268 13 : nRead = VSIFReadL((void *)(pabyBuffer + nOverlap), 1,
1269 : nBufferSize - nOverlap, fp);
1270 :
1271 13 : if (VSIFEofL(fp))
1272 12 : bEOF = true;
1273 : else
1274 1 : bEOF = false;
1275 :
1276 13 : if (nRead == 0 && !bEOF)
1277 : {
1278 0 : VSIFree(pabyBuffer);
1279 0 : CPLError(CE_Failure, CPLE_FileIO,
1280 : "Unable to read from grid file (possible corruption).\n");
1281 0 : return CE_Failure;
1282 : }
1283 :
1284 : /* FIXME: Should use SEEK_CUR, review integer promotions... */
1285 : vsi_l_offset nNewPos =
1286 : (nShiftSize >= 0)
1287 13 : ? VSIFTellL(fp) + nShiftSize - nRead - nOverlap
1288 13 : : VSIFTellL(fp) - (-nShiftSize) - nRead - nOverlap;
1289 13 : if (VSIFSeekL(fp, nNewPos, SEEK_SET) != 0)
1290 : {
1291 0 : VSIFree(pabyBuffer);
1292 0 : CPLError(CE_Failure, CPLE_FileIO,
1293 : "Unable to seek in grid file (possible corruption).\n");
1294 0 : return CE_Failure;
1295 : }
1296 :
1297 13 : size_t nWritten = VSIFWriteL((void *)pabyBuffer, 1, nRead, fp);
1298 13 : if (nWritten != nRead)
1299 : {
1300 0 : VSIFree(pabyBuffer);
1301 0 : CPLError(CE_Failure, CPLE_FileIO,
1302 : "Unable to write to grid file (out of space?).\n");
1303 0 : return CE_Failure;
1304 : }
1305 :
1306 : /* shift overlapped contents to the front of the buffer if necessary */
1307 13 : if (nOverlap > 0)
1308 0 : memmove(pabyBuffer, pabyBuffer + nRead, nOverlap);
1309 : }
1310 :
1311 : /* write the remainder of the buffer or overwrite leftovers and finish */
1312 12 : if (nShiftSize > 0)
1313 : {
1314 0 : size_t nTailSize = nOverlap;
1315 0 : while (nTailSize > 0 &&
1316 0 : isspace((unsigned char)pabyBuffer[nTailSize - 1]))
1317 0 : nTailSize--;
1318 :
1319 0 : if (VSIFWriteL((void *)pabyBuffer, 1, nTailSize, fp) != nTailSize)
1320 : {
1321 0 : VSIFree(pabyBuffer);
1322 0 : CPLError(CE_Failure, CPLE_FileIO,
1323 : "Unable to write to grid file (out of space?).\n");
1324 0 : return CE_Failure;
1325 : }
1326 :
1327 0 : if (VSIFWriteL((void *)pszEOL, 1, strlen(pszEOL), fp) != strlen(pszEOL))
1328 : {
1329 0 : VSIFree(pabyBuffer);
1330 0 : CPLError(CE_Failure, CPLE_FileIO,
1331 : "Unable to write to grid file (out of space?).\n");
1332 0 : return CE_Failure;
1333 : }
1334 : }
1335 : else
1336 : {
1337 : /* FIXME: ftruncate()? */
1338 : /* FIXME: Should use SEEK_CUR, review integer promotions... */
1339 12 : if (VSIFSeekL(fp, VSIFTellL(fp) - strlen(pszEOL), SEEK_SET) != 0)
1340 : {
1341 0 : VSIFree(pabyBuffer);
1342 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to seek in grid file.\n");
1343 0 : return CE_Failure;
1344 : }
1345 :
1346 301 : for (int iPadding = 0; iPadding < -nShiftSize; iPadding++)
1347 : {
1348 289 : if (VSIFWriteL((void *)" ", 1, 1, fp) != 1)
1349 : {
1350 0 : VSIFree(pabyBuffer);
1351 0 : CPLError(CE_Failure, CPLE_FileIO,
1352 : "Error writing to grid file.\n");
1353 0 : return CE_Failure;
1354 : }
1355 : }
1356 :
1357 12 : if (VSIFWriteL((void *)pszEOL, 1, strlen(pszEOL), fp) != strlen(pszEOL))
1358 : {
1359 0 : VSIFree(pabyBuffer);
1360 0 : CPLError(CE_Failure, CPLE_FileIO,
1361 : "Unable to write to grid file (out of space?).\n");
1362 0 : return CE_Failure;
1363 : }
1364 : }
1365 :
1366 12 : VSIFree(pabyBuffer);
1367 12 : return CE_None;
1368 : }
1369 :
1370 : /************************************************************************/
1371 : /* UpdateHeader() */
1372 : /************************************************************************/
1373 :
1374 0 : CPLErr GSAGDataset::UpdateHeader()
1375 :
1376 : {
1377 0 : GSAGRasterBand *poBand = cpl::down_cast<GSAGRasterBand *>(GetRasterBand(1));
1378 :
1379 0 : std::ostringstream ssOutBuf;
1380 0 : ssOutBuf.precision(nFIELD_PRECISION);
1381 0 : ssOutBuf.setf(std::ios::uppercase);
1382 :
1383 : /* signature */
1384 0 : ssOutBuf << "DSAA" << szEOL;
1385 :
1386 : /* columns rows */
1387 0 : ssOutBuf << nRasterXSize << " " << nRasterYSize << szEOL;
1388 :
1389 : /* x range */
1390 0 : ssOutBuf << poBand->dfMinX << " " << poBand->dfMaxX << szEOL;
1391 :
1392 : /* y range */
1393 0 : ssOutBuf << poBand->dfMinY << " " << poBand->dfMaxY << szEOL;
1394 :
1395 : /* z range */
1396 0 : ssOutBuf << poBand->dfMinZ << " " << poBand->dfMaxZ << szEOL;
1397 :
1398 0 : CPLString sOut = ssOutBuf.str();
1399 0 : if (sOut.length() != poBand->panLineOffset[0])
1400 : {
1401 0 : int nShiftSize = (int)(sOut.length() - poBand->panLineOffset[0]);
1402 0 : if (ShiftFileContents(fp, poBand->panLineOffset[0], nShiftSize,
1403 0 : szEOL) != CE_None)
1404 : {
1405 0 : CPLError(CE_Failure, CPLE_FileIO,
1406 : "Unable to update grid header, "
1407 : "failure shifting file contents.\n");
1408 0 : return CE_Failure;
1409 : }
1410 :
1411 0 : for (size_t iLine = 0;
1412 0 : iLine < static_cast<unsigned>(nRasterYSize + 1) &&
1413 0 : poBand->panLineOffset[iLine] != 0;
1414 : iLine++)
1415 0 : poBand->panLineOffset[iLine] += nShiftSize;
1416 : }
1417 :
1418 0 : if (VSIFSeekL(fp, 0, SEEK_SET) != 0)
1419 : {
1420 0 : CPLError(CE_Failure, CPLE_FileIO,
1421 : "Unable to seek to start of grid file.\n");
1422 0 : return CE_Failure;
1423 : }
1424 :
1425 0 : if (VSIFWriteL(sOut.c_str(), 1, sOut.length(), fp) != sOut.length())
1426 : {
1427 0 : CPLError(CE_Failure, CPLE_FileIO,
1428 : "Unable to update file header. Disk full?\n");
1429 0 : return CE_Failure;
1430 : }
1431 :
1432 0 : return CE_None;
1433 : }
1434 :
1435 : /************************************************************************/
1436 : /* CreateCopy() */
1437 : /************************************************************************/
1438 :
1439 30 : GDALDataset *GSAGDataset::CreateCopy(const char *pszFilename,
1440 : GDALDataset *poSrcDS, int bStrict,
1441 : CPL_UNUSED char **papszOptions,
1442 : GDALProgressFunc pfnProgress,
1443 : void *pProgressData)
1444 : {
1445 30 : if (pfnProgress == nullptr)
1446 0 : pfnProgress = GDALDummyProgress;
1447 :
1448 30 : int nBands = poSrcDS->GetRasterCount();
1449 30 : if (nBands == 0)
1450 : {
1451 1 : CPLError(
1452 : CE_Failure, CPLE_NotSupported,
1453 : "GSAG driver does not support source dataset with zero band.\n");
1454 1 : return nullptr;
1455 : }
1456 29 : else if (nBands > 1)
1457 : {
1458 4 : if (bStrict)
1459 : {
1460 4 : CPLError(CE_Failure, CPLE_NotSupported,
1461 : "Unable to create copy, Golden Software ASCII Grid "
1462 : "format only supports one raster band.\n");
1463 4 : return nullptr;
1464 : }
1465 : else
1466 0 : CPLError(CE_Warning, CPLE_NotSupported,
1467 : "Golden Software ASCII Grid format only supports one "
1468 : "raster band, first band will be copied.\n");
1469 : }
1470 :
1471 25 : if (!pfnProgress(0.0, nullptr, pProgressData))
1472 : {
1473 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated\n");
1474 0 : return nullptr;
1475 : }
1476 :
1477 25 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
1478 :
1479 25 : if (fp == nullptr)
1480 : {
1481 3 : CPLError(CE_Failure, CPLE_OpenFailed,
1482 : "Attempt to create file '%s' failed.\n", pszFilename);
1483 3 : return nullptr;
1484 : }
1485 :
1486 22 : const int nXSize = poSrcDS->GetRasterXSize();
1487 22 : const int nYSize = poSrcDS->GetRasterYSize();
1488 22 : GDALGeoTransform gt;
1489 :
1490 22 : poSrcDS->GetGeoTransform(gt);
1491 :
1492 44 : std::ostringstream ssHeader;
1493 22 : ssHeader.precision(nFIELD_PRECISION);
1494 22 : ssHeader.setf(std::ios::uppercase);
1495 :
1496 22 : ssHeader << "DSAA\x0D\x0A";
1497 :
1498 22 : ssHeader << nXSize << " " << nYSize << "\x0D\x0A";
1499 :
1500 22 : ssHeader << gt[0] + gt[1] / 2 << " " << gt[1] * (nXSize - 0.5) + gt[0]
1501 22 : << "\x0D\x0A";
1502 :
1503 22 : ssHeader << gt[5] * (nYSize - 0.5) + gt[3] << " " << gt[3] + gt[5] / 2
1504 22 : << "\x0D\x0A";
1505 :
1506 44 : if (VSIFWriteL((void *)ssHeader.str().c_str(), 1, ssHeader.str().length(),
1507 22 : fp) != ssHeader.str().length())
1508 : {
1509 1 : VSIFCloseL(fp);
1510 1 : CPLError(CE_Failure, CPLE_FileIO,
1511 : "Unable to create copy, writing header failed.\n");
1512 1 : return nullptr;
1513 : }
1514 :
1515 : /* Save the location and write placeholders for the min/max Z value */
1516 21 : vsi_l_offset nRangeStart = VSIFTellL(fp);
1517 21 : const char *szDummyRange = "0.0000000000001 0.0000000000001\x0D\x0A";
1518 21 : size_t nDummyRangeLen = strlen(szDummyRange);
1519 21 : if (VSIFWriteL((void *)szDummyRange, 1, nDummyRangeLen, fp) !=
1520 : nDummyRangeLen)
1521 : {
1522 1 : VSIFCloseL(fp);
1523 1 : CPLError(CE_Failure, CPLE_FileIO,
1524 : "Unable to create copy, writing header failed.\n");
1525 1 : return nullptr;
1526 : }
1527 :
1528 : /* -------------------------------------------------------------------- */
1529 : /* Copy band data. */
1530 : /* -------------------------------------------------------------------- */
1531 20 : double *pdfData = (double *)VSI_MALLOC2_VERBOSE(nXSize, sizeof(double));
1532 20 : if (pdfData == nullptr)
1533 : {
1534 0 : VSIFCloseL(fp);
1535 0 : return nullptr;
1536 : }
1537 :
1538 20 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
1539 : int bSrcHasNDValue;
1540 20 : double dfSrcNoDataValue = poSrcBand->GetNoDataValue(&bSrcHasNDValue);
1541 20 : double dfMin = std::numeric_limits<double>::max();
1542 20 : double dfMax = std::numeric_limits<double>::lowest();
1543 184 : for (int iRow = 0; iRow < nYSize; iRow++)
1544 : {
1545 : CPLErr eErr =
1546 172 : poSrcBand->RasterIO(GF_Read, 0, nYSize - iRow - 1, nXSize, 1,
1547 : pdfData, nXSize, 1, GDT_Float64, 0, 0, nullptr);
1548 :
1549 172 : if (eErr != CE_None)
1550 : {
1551 0 : VSIFCloseL(fp);
1552 0 : VSIFree(pdfData);
1553 0 : return nullptr;
1554 : }
1555 :
1556 356 : for (int iCol = 0; iCol < nXSize;)
1557 : {
1558 2076 : for (int iCount = 0; iCount < 10 && iCol < nXSize; iCount++, iCol++)
1559 : {
1560 1892 : double dfValue = pdfData[iCol];
1561 :
1562 1892 : if (bSrcHasNDValue && AlmostEqual(dfValue, dfSrcNoDataValue))
1563 : {
1564 0 : dfValue = dfNODATA_VALUE;
1565 : }
1566 : else
1567 : {
1568 1892 : if (dfValue > dfMax)
1569 22 : dfMax = dfValue;
1570 :
1571 1892 : if (dfValue < dfMin)
1572 27 : dfMin = dfValue;
1573 : }
1574 :
1575 1892 : std::ostringstream ssOut;
1576 1892 : ssOut.precision(nFIELD_PRECISION);
1577 1892 : ssOut.setf(std::ios::uppercase);
1578 1892 : ssOut << dfValue << " ";
1579 1892 : CPLString sOut = ssOut.str();
1580 :
1581 3784 : if (VSIFWriteL(sOut.c_str(), 1, sOut.length(), fp) !=
1582 1892 : sOut.length())
1583 : {
1584 8 : VSIFCloseL(fp);
1585 8 : VSIFree(pdfData);
1586 8 : CPLError(CE_Failure, CPLE_FileIO,
1587 : "Unable to write grid cell. Disk full?\n");
1588 8 : return nullptr;
1589 : }
1590 : }
1591 :
1592 184 : if (VSIFWriteL((void *)"\x0D\x0A", 1, 2, fp) != 2)
1593 : {
1594 0 : VSIFCloseL(fp);
1595 0 : VSIFree(pdfData);
1596 0 : CPLError(CE_Failure, CPLE_FileIO,
1597 : "Unable to finish write of grid line. Disk full?\n");
1598 0 : return nullptr;
1599 : }
1600 : }
1601 :
1602 164 : if (VSIFWriteL((void *)"\x0D\x0A", 1, 2, fp) != 2)
1603 : {
1604 0 : VSIFCloseL(fp);
1605 0 : VSIFree(pdfData);
1606 0 : CPLError(CE_Failure, CPLE_FileIO,
1607 : "Unable to finish write of grid row. Disk full?\n");
1608 0 : return nullptr;
1609 : }
1610 :
1611 164 : if (!pfnProgress(static_cast<double>(iRow + 1) / nYSize, nullptr,
1612 : pProgressData))
1613 : {
1614 0 : VSIFCloseL(fp);
1615 0 : VSIFree(pdfData);
1616 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1617 0 : return nullptr;
1618 : }
1619 : }
1620 :
1621 12 : VSIFree(pdfData);
1622 :
1623 : /* write out the min and max values */
1624 24 : std::ostringstream ssRange;
1625 12 : ssRange.precision(nFIELD_PRECISION);
1626 12 : ssRange.setf(std::ios::uppercase);
1627 12 : ssRange << dfMin << " " << dfMax << "\x0D\x0A";
1628 12 : if (ssRange.str().length() != nDummyRangeLen)
1629 : {
1630 12 : int nShiftSize = static_cast<int>(ssRange.str().length()) -
1631 12 : static_cast<int>(nDummyRangeLen);
1632 12 : if (ShiftFileContents(fp, nRangeStart + nDummyRangeLen, nShiftSize,
1633 12 : "\x0D\x0A") != CE_None)
1634 : {
1635 0 : VSIFCloseL(fp);
1636 0 : CPLError(CE_Failure, CPLE_FileIO,
1637 : "Unable to shift file contents.\n");
1638 0 : return nullptr;
1639 : }
1640 : }
1641 :
1642 12 : if (VSIFSeekL(fp, nRangeStart, SEEK_SET) != 0)
1643 : {
1644 0 : VSIFCloseL(fp);
1645 0 : CPLError(CE_Failure, CPLE_FileIO,
1646 : "Unable to seek to start of grid file copy.\n");
1647 0 : return nullptr;
1648 : }
1649 :
1650 24 : if (VSIFWriteL((void *)ssRange.str().c_str(), 1, ssRange.str().length(),
1651 12 : fp) != ssRange.str().length())
1652 : {
1653 0 : VSIFCloseL(fp);
1654 0 : CPLError(CE_Failure, CPLE_FileIO,
1655 : "Unable to write range information.\n");
1656 0 : return nullptr;
1657 : }
1658 :
1659 12 : VSIFCloseL(fp);
1660 :
1661 12 : GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen(pszFilename, GA_Update);
1662 12 : if (poDS)
1663 : {
1664 12 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1665 : }
1666 12 : return poDS;
1667 : }
1668 :
1669 : /************************************************************************/
1670 : /* GDALRegister_GSAG() */
1671 : /************************************************************************/
1672 :
1673 2038 : void GDALRegister_GSAG()
1674 :
1675 : {
1676 2038 : if (GDALGetDriverByName("GSAG") != nullptr)
1677 283 : return;
1678 :
1679 1755 : GDALDriver *poDriver = new GDALDriver();
1680 :
1681 1755 : poDriver->SetDescription("GSAG");
1682 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1683 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1684 1755 : "Golden Software ASCII Grid (.grd)");
1685 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gsag.html");
1686 1755 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1687 1755 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1688 : "Byte Int16 UInt16 Int32 UInt32 "
1689 1755 : "Float32 Float64");
1690 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1691 :
1692 1755 : poDriver->pfnIdentify = GSAGDataset::Identify;
1693 1755 : poDriver->pfnOpen = GSAGDataset::Open;
1694 1755 : poDriver->pfnCreateCopy = GSAGDataset::CreateCopy;
1695 :
1696 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
1697 : }
|
