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