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