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