Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: ZMap driver
4 : * Purpose: GDALDataset driver for ZMap dataset.
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2011-2012, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_string.h"
14 : #include "cpl_vsi_virtual.h"
15 : #include "gdal_frmts.h"
16 : #include "gdal_pam.h"
17 :
18 : #include <array>
19 : #include <cmath>
20 : #include <deque>
21 :
22 : /************************************************************************/
23 : /* ==================================================================== */
24 : /* ZMapDataset */
25 : /* ==================================================================== */
26 : /************************************************************************/
27 :
28 : class ZMapRasterBand;
29 :
30 19 : class ZMapDataset final : public GDALPamDataset
31 : {
32 : friend class ZMapRasterBand;
33 :
34 : VSIVirtualHandleUniquePtr m_fp{};
35 : int m_nValuesPerLine = 0;
36 : int m_nFieldSize = 0;
37 : int m_nDecimalCount = 0;
38 : int m_nColNum = -1;
39 : double m_dfNoDataValue = 0;
40 : vsi_l_offset m_nDataStartOff = 0;
41 : std::array<double, 6> m_adfGeoTransform = {{0, 1, 0, 0, 0, 1}};
42 : int m_nFirstDataLine = 0;
43 : int m_nCurLine = 0;
44 : std::deque<double> m_odfQueue{};
45 :
46 : public:
47 : ZMapDataset();
48 : virtual ~ZMapDataset();
49 :
50 : virtual CPLErr GetGeoTransform(double *) override;
51 :
52 : static GDALDataset *Open(GDALOpenInfo *);
53 : static int Identify(GDALOpenInfo *);
54 : static GDALDataset *CreateCopy(const char *pszFilename,
55 : GDALDataset *poSrcDS, int bStrict,
56 : char **papszOptions,
57 : GDALProgressFunc pfnProgress,
58 : void *pProgressData);
59 : };
60 :
61 : /************************************************************************/
62 : /* ==================================================================== */
63 : /* ZMapRasterBand */
64 : /* ==================================================================== */
65 : /************************************************************************/
66 :
67 : class ZMapRasterBand final : public GDALPamRasterBand
68 : {
69 : friend class ZMapDataset;
70 :
71 : public:
72 : explicit ZMapRasterBand(ZMapDataset *);
73 :
74 : virtual CPLErr IReadBlock(int, int, void *) override;
75 : virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
76 : };
77 :
78 : /************************************************************************/
79 : /* ZMapRasterBand() */
80 : /************************************************************************/
81 :
82 19 : ZMapRasterBand::ZMapRasterBand(ZMapDataset *poDSIn)
83 :
84 : {
85 19 : poDS = poDSIn;
86 19 : nBand = 1;
87 :
88 19 : eDataType = GDT_Float64;
89 :
90 : // The format is column oriented! That is we have first the value of
91 : // pixel (col=0, line=0), then the one of (col=0, line=1), etc.
92 19 : nBlockXSize = 1;
93 19 : nBlockYSize = poDSIn->GetRasterYSize();
94 19 : }
95 :
96 : /************************************************************************/
97 : /* IReadBlock() */
98 : /************************************************************************/
99 :
100 42 : CPLErr ZMapRasterBand::IReadBlock(int nBlockXOff, CPL_UNUSED int nBlockYOff,
101 : void *pImage)
102 : {
103 42 : ZMapDataset *poGDS = cpl::down_cast<ZMapDataset *>(poDS);
104 :
105 : // If seeking backwards in term of columns, reset reading to the first
106 : // column
107 42 : if (nBlockXOff < poGDS->m_nColNum + 1)
108 : {
109 0 : poGDS->m_fp->Seek(poGDS->m_nDataStartOff, SEEK_SET);
110 0 : poGDS->m_nColNum = -1;
111 0 : poGDS->m_nCurLine = poGDS->m_nFirstDataLine;
112 0 : poGDS->m_odfQueue.clear();
113 : }
114 :
115 42 : if (nBlockXOff > poGDS->m_nColNum + 1)
116 : {
117 0 : for (int i = poGDS->m_nColNum + 1; i < nBlockXOff; i++)
118 : {
119 0 : if (IReadBlock(i, 0, nullptr) != CE_None)
120 0 : return CE_Failure;
121 : }
122 : }
123 :
124 42 : int iRow = 0;
125 42 : const double dfExp = std::pow(10.0, poGDS->m_nDecimalCount);
126 42 : double *padfImage = reinterpret_cast<double *>(pImage);
127 :
128 : // If we have previously read too many values, start by consuming the
129 : // queue
130 45 : while (iRow < nRasterYSize && !poGDS->m_odfQueue.empty())
131 : {
132 3 : if (padfImage)
133 3 : padfImage[iRow] = poGDS->m_odfQueue.front();
134 3 : ++iRow;
135 3 : poGDS->m_odfQueue.pop_front();
136 : }
137 :
138 : // Now read as many lines as needed to finish filling the column buffer
139 245 : while (iRow < nRasterYSize)
140 : {
141 203 : constexpr int MARGIN = 16; // Should be at least 2 for \r\n
142 203 : char *pszLine = const_cast<char *>(CPLReadLine2L(
143 : poGDS->m_fp.get(),
144 203 : poGDS->m_nValuesPerLine * poGDS->m_nFieldSize + MARGIN, nullptr));
145 203 : ++poGDS->m_nCurLine;
146 203 : if (pszLine == nullptr)
147 0 : return CE_Failure;
148 :
149 : // Each line should have at most m_nValuesPerLine values of size
150 : // m_nFieldSize
151 203 : const int nLineLen = static_cast<int>(strlen(pszLine));
152 203 : if ((nLineLen % poGDS->m_nFieldSize) != 0)
153 : {
154 0 : CPLError(CE_Failure, CPLE_AppDefined,
155 : "Line %d has length %d, which is not a multiple of %d",
156 : poGDS->m_nCurLine, nLineLen, poGDS->m_nFieldSize);
157 0 : return CE_Failure;
158 : }
159 :
160 203 : const int nValuesThisLine = nLineLen / poGDS->m_nFieldSize;
161 203 : if (nValuesThisLine > poGDS->m_nValuesPerLine)
162 : {
163 0 : CPLError(
164 : CE_Failure, CPLE_AppDefined,
165 : "Line %d has %d values, whereas the maximum expected is %d",
166 : poGDS->m_nCurLine, nValuesThisLine, poGDS->m_nValuesPerLine);
167 0 : return CE_Failure;
168 : }
169 :
170 1013 : for (int iValueThisLine = 0; iValueThisLine < nValuesThisLine;
171 : iValueThisLine++)
172 : {
173 810 : char *pszValue = pszLine + iValueThisLine * poGDS->m_nFieldSize;
174 810 : const char chSaved = pszValue[poGDS->m_nFieldSize];
175 810 : pszValue[poGDS->m_nFieldSize] = 0;
176 810 : const double dfVal = strchr(pszValue, '.') != nullptr
177 810 : ? CPLAtofM(pszValue)
178 0 : : atoi(pszValue) * dfExp;
179 810 : pszValue[poGDS->m_nFieldSize] = chSaved;
180 810 : if (iRow < nRasterYSize)
181 : {
182 807 : if (padfImage)
183 807 : padfImage[iRow] = dfVal;
184 807 : ++iRow;
185 : }
186 : else
187 : {
188 3 : poGDS->m_odfQueue.push_back(dfVal);
189 : }
190 : }
191 : }
192 :
193 42 : poGDS->m_nColNum++;
194 :
195 42 : return CE_None;
196 : }
197 :
198 : /************************************************************************/
199 : /* GetNoDataValue() */
200 : /************************************************************************/
201 :
202 3 : double ZMapRasterBand::GetNoDataValue(int *pbSuccess)
203 : {
204 3 : ZMapDataset *poGDS = cpl::down_cast<ZMapDataset *>(poDS);
205 :
206 3 : if (pbSuccess)
207 3 : *pbSuccess = TRUE;
208 :
209 3 : return poGDS->m_dfNoDataValue;
210 : }
211 :
212 : /************************************************************************/
213 : /* ~ZMapDataset() */
214 : /************************************************************************/
215 :
216 : ZMapDataset::ZMapDataset() = default;
217 :
218 : /************************************************************************/
219 : /* ~ZMapDataset() */
220 : /************************************************************************/
221 :
222 38 : ZMapDataset::~ZMapDataset()
223 :
224 : {
225 19 : FlushCache(true);
226 38 : }
227 :
228 : /************************************************************************/
229 : /* Identify() */
230 : /************************************************************************/
231 :
232 51697 : int ZMapDataset::Identify(GDALOpenInfo *poOpenInfo)
233 : {
234 51697 : if (poOpenInfo->nHeaderBytes == 0)
235 48639 : return FALSE;
236 :
237 : /* -------------------------------------------------------------------- */
238 : /* Check that it looks roughly as a ZMap dataset */
239 : /* -------------------------------------------------------------------- */
240 3058 : const char *pszData =
241 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
242 :
243 : /* Skip comments line at the beginning */
244 3058 : int i = 0;
245 3058 : if (pszData[i] == '!')
246 : {
247 33 : i++;
248 706 : for (; i < poOpenInfo->nHeaderBytes; i++)
249 : {
250 706 : char ch = pszData[i];
251 706 : if (ch == 13 || ch == 10)
252 : {
253 89 : i++;
254 89 : if (ch == 13 && pszData[i] == 10)
255 0 : i++;
256 89 : if (pszData[i] != '!')
257 33 : break;
258 : }
259 : }
260 : }
261 :
262 3058 : if (pszData[i] != '@')
263 3030 : return FALSE;
264 28 : i++;
265 :
266 56 : const CPLStringList aosTokens(CSLTokenizeString2(pszData + i, ",", 0));
267 28 : if (aosTokens.size() < 3)
268 : {
269 4 : return FALSE;
270 : }
271 :
272 24 : const char *pszToken = aosTokens[1];
273 48 : while (*pszToken == ' ')
274 24 : pszToken++;
275 :
276 24 : return STARTS_WITH(pszToken, "GRID");
277 : }
278 :
279 : /************************************************************************/
280 : /* Open() */
281 : /************************************************************************/
282 :
283 19 : GDALDataset *ZMapDataset::Open(GDALOpenInfo *poOpenInfo)
284 :
285 : {
286 19 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
287 0 : return nullptr;
288 :
289 : /* -------------------------------------------------------------------- */
290 : /* Confirm the requested access is supported. */
291 : /* -------------------------------------------------------------------- */
292 19 : if (poOpenInfo->eAccess == GA_Update)
293 : {
294 0 : ReportUpdateNotSupportedByDriver("ZMAP");
295 0 : return nullptr;
296 : }
297 :
298 38 : auto poDS = std::make_unique<ZMapDataset>();
299 19 : poDS->m_fp.reset(poOpenInfo->fpL);
300 19 : poOpenInfo->fpL = nullptr;
301 :
302 : /* -------------------------------------------------------------------- */
303 : /* Find dataset characteristics */
304 : /* -------------------------------------------------------------------- */
305 :
306 : const char *pszLine;
307 19 : int nLine = 0;
308 19 : constexpr int MAX_HEADER_LINE = 1024;
309 76 : while ((pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE,
310 76 : nullptr)) != nullptr)
311 : {
312 76 : ++nLine;
313 76 : if (*pszLine == '!')
314 : {
315 57 : continue;
316 : }
317 : else
318 19 : break;
319 : }
320 : // cppcheck-suppress knownConditionTrueFalse
321 19 : if (pszLine == nullptr)
322 : {
323 0 : return nullptr;
324 : }
325 :
326 : /* Parse first header line */
327 38 : CPLStringList aosTokensFirstLine(CSLTokenizeString2(pszLine, ",", 0));
328 19 : if (aosTokensFirstLine.size() != 3)
329 : {
330 0 : return nullptr;
331 : }
332 :
333 19 : const int nValuesPerLine = atoi(aosTokensFirstLine[2]);
334 19 : if (nValuesPerLine <= 0)
335 : {
336 0 : CPLError(CE_Failure, CPLE_AppDefined,
337 : "Invalid/unsupported value for nValuesPerLine = %d",
338 : nValuesPerLine);
339 0 : return nullptr;
340 : }
341 :
342 : /* Parse second header line */
343 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
344 19 : ++nLine;
345 19 : if (pszLine == nullptr)
346 : {
347 0 : return nullptr;
348 : }
349 : const CPLStringList aosTokensSecondLine(
350 38 : CSLTokenizeString2(pszLine, ",", 0));
351 19 : if (aosTokensSecondLine.size() != 5)
352 : {
353 0 : return nullptr;
354 : }
355 :
356 19 : const int nFieldSize = atoi(aosTokensSecondLine[0]);
357 19 : const double dfNoDataValue = CPLAtofM(aosTokensSecondLine[1]);
358 19 : const int nDecimalCount = atoi(aosTokensSecondLine[3]);
359 19 : const int nColumnNumber = atoi(aosTokensSecondLine[4]);
360 :
361 19 : if (nFieldSize <= 0 || nFieldSize >= 40)
362 : {
363 0 : CPLError(CE_Failure, CPLE_AppDefined,
364 : "Invalid/unsupported value for nFieldSize = %d", nFieldSize);
365 0 : return nullptr;
366 : }
367 :
368 19 : if (nDecimalCount <= 0 || nDecimalCount >= nFieldSize)
369 : {
370 0 : CPLError(CE_Failure, CPLE_AppDefined,
371 : "Invalid/unsupported value for nDecimalCount = %d",
372 : nDecimalCount);
373 0 : return nullptr;
374 : }
375 :
376 19 : if (nColumnNumber != 1)
377 : {
378 0 : CPLError(CE_Failure, CPLE_AppDefined,
379 : "Invalid/unsupported value for nColumnNumber = %d",
380 : nColumnNumber);
381 0 : return nullptr;
382 : }
383 :
384 19 : if (nValuesPerLine <= 0 || nFieldSize > 1024 * 1024 / nValuesPerLine)
385 : {
386 0 : CPLError(CE_Failure, CPLE_AppDefined,
387 : "Invalid/unsupported value for nFieldSize = %d x "
388 : "nValuesPerLine = %d",
389 : nFieldSize, nValuesPerLine);
390 0 : return nullptr;
391 : }
392 :
393 : /* Parse third header line */
394 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
395 19 : ++nLine;
396 19 : if (pszLine == nullptr)
397 : {
398 0 : return nullptr;
399 : }
400 38 : const CPLStringList aosTokensThirdLine(CSLTokenizeString2(pszLine, ",", 0));
401 19 : if (aosTokensThirdLine.size() != 6)
402 : {
403 0 : return nullptr;
404 : }
405 :
406 19 : const int nRows = atoi(aosTokensThirdLine[0]);
407 19 : const int nCols = atoi(aosTokensThirdLine[1]);
408 19 : const double dfMinX = CPLAtofM(aosTokensThirdLine[2]);
409 19 : const double dfMaxX = CPLAtofM(aosTokensThirdLine[3]);
410 19 : const double dfMinY = CPLAtofM(aosTokensThirdLine[4]);
411 19 : const double dfMaxY = CPLAtofM(aosTokensThirdLine[5]);
412 :
413 19 : if (!GDALCheckDatasetDimensions(nCols, nRows) || nCols == 1 || nRows == 1)
414 : {
415 0 : return nullptr;
416 : }
417 :
418 : /* Ignore fourth header line */
419 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
420 19 : ++nLine;
421 19 : if (pszLine == nullptr)
422 : {
423 0 : return nullptr;
424 : }
425 :
426 : /* Check fifth header line */
427 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
428 19 : ++nLine;
429 19 : if (pszLine == nullptr || pszLine[0] != '@')
430 : {
431 0 : return nullptr;
432 : }
433 :
434 : /* -------------------------------------------------------------------- */
435 : /* Create a corresponding GDALDataset. */
436 : /* -------------------------------------------------------------------- */
437 19 : poDS->m_nDataStartOff = VSIFTellL(poDS->m_fp.get());
438 19 : poDS->m_nValuesPerLine = nValuesPerLine;
439 19 : poDS->m_nFieldSize = nFieldSize;
440 19 : poDS->m_nDecimalCount = nDecimalCount;
441 19 : poDS->nRasterXSize = nCols;
442 19 : poDS->nRasterYSize = nRows;
443 19 : poDS->m_dfNoDataValue = dfNoDataValue;
444 19 : poDS->m_nFirstDataLine = nLine;
445 :
446 19 : if (CPLTestBool(CPLGetConfigOption("ZMAP_PIXEL_IS_POINT", "FALSE")))
447 : {
448 0 : const double dfStepX = (dfMaxX - dfMinX) / (nCols - 1);
449 0 : const double dfStepY = (dfMaxY - dfMinY) / (nRows - 1);
450 :
451 0 : poDS->m_adfGeoTransform[0] = dfMinX - dfStepX / 2;
452 0 : poDS->m_adfGeoTransform[1] = dfStepX;
453 0 : poDS->m_adfGeoTransform[3] = dfMaxY + dfStepY / 2;
454 0 : poDS->m_adfGeoTransform[5] = -dfStepY;
455 : }
456 : else
457 : {
458 19 : const double dfStepX = (dfMaxX - dfMinX) / nCols;
459 19 : const double dfStepY = (dfMaxY - dfMinY) / nRows;
460 :
461 19 : poDS->m_adfGeoTransform[0] = dfMinX;
462 19 : poDS->m_adfGeoTransform[1] = dfStepX;
463 19 : poDS->m_adfGeoTransform[3] = dfMaxY;
464 19 : poDS->m_adfGeoTransform[5] = -dfStepY;
465 : }
466 :
467 : /* -------------------------------------------------------------------- */
468 : /* Create band information objects. */
469 : /* -------------------------------------------------------------------- */
470 19 : poDS->nBands = 1;
471 19 : poDS->SetBand(1, std::make_unique<ZMapRasterBand>(poDS.get()));
472 :
473 : /* -------------------------------------------------------------------- */
474 : /* Initialize any PAM information. */
475 : /* -------------------------------------------------------------------- */
476 19 : poDS->SetDescription(poOpenInfo->pszFilename);
477 19 : poDS->TryLoadXML();
478 :
479 : /* -------------------------------------------------------------------- */
480 : /* Support overviews. */
481 : /* -------------------------------------------------------------------- */
482 19 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
483 19 : return poDS.release();
484 : }
485 :
486 : /************************************************************************/
487 : /* WriteRightJustified() */
488 : /************************************************************************/
489 :
490 1668 : static void WriteRightJustified(VSIVirtualHandleUniquePtr &fp,
491 : const char *pszValue, int nWidth)
492 : {
493 1668 : int nLen = (int)strlen(pszValue);
494 1668 : CPLAssert(nLen <= nWidth);
495 16550 : for (int i = 0; i < nWidth - nLen; i++)
496 14882 : fp->Write(" ", 1, 1);
497 1668 : fp->Write(pszValue, 1, nLen);
498 1668 : }
499 :
500 70 : static void WriteRightJustified(VSIVirtualHandleUniquePtr &fp, int nValue,
501 : int nWidth)
502 : {
503 140 : CPLString osValue(CPLSPrintf("%d", nValue));
504 70 : WriteRightJustified(fp, osValue.c_str(), nWidth);
505 70 : }
506 :
507 1584 : static void WriteRightJustified(VSIVirtualHandleUniquePtr &fp, double dfValue,
508 : int nWidth, int nDecimals = -1)
509 : {
510 : char szFormat[32];
511 1584 : if (nDecimals >= 0)
512 1584 : snprintf(szFormat, sizeof(szFormat), "%%.%df", nDecimals);
513 : else
514 0 : snprintf(szFormat, sizeof(szFormat), "%%g");
515 1584 : char *pszValue = const_cast<char *>(CPLSPrintf(szFormat, dfValue));
516 1584 : char *pszE = strchr(pszValue, 'e');
517 1584 : if (pszE)
518 0 : *pszE = 'E';
519 :
520 1584 : if (static_cast<int>(strlen(pszValue)) > nWidth)
521 : {
522 16 : CPLAssert(nDecimals >= 0);
523 16 : snprintf(szFormat, sizeof(szFormat), "%%.%dg", nDecimals);
524 16 : pszValue = const_cast<char *>(CPLSPrintf(szFormat, dfValue));
525 16 : pszE = strchr(pszValue, 'e');
526 16 : if (pszE)
527 14 : *pszE = 'E';
528 : }
529 :
530 3168 : CPLString osValue(pszValue);
531 1584 : WriteRightJustified(fp, osValue.c_str(), nWidth);
532 1584 : }
533 :
534 : /************************************************************************/
535 : /* CreateCopy() */
536 : /************************************************************************/
537 :
538 22 : GDALDataset *ZMapDataset::CreateCopy(const char *pszFilename,
539 : GDALDataset *poSrcDS, int bStrict,
540 : CPL_UNUSED char **papszOptions,
541 : GDALProgressFunc pfnProgress,
542 : void *pProgressData)
543 : {
544 : /* -------------------------------------------------------------------- */
545 : /* Some some rudimentary checks */
546 : /* -------------------------------------------------------------------- */
547 22 : int nBands = poSrcDS->GetRasterCount();
548 22 : if (nBands == 0)
549 : {
550 1 : CPLError(
551 : CE_Failure, CPLE_NotSupported,
552 : "ZMap driver does not support source dataset with zero band.\n");
553 1 : return nullptr;
554 : }
555 :
556 21 : if (nBands != 1)
557 : {
558 4 : CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
559 : "ZMap driver only uses the first band of the dataset.\n");
560 4 : if (bStrict)
561 4 : return nullptr;
562 : }
563 :
564 17 : if (pfnProgress && !pfnProgress(0.0, nullptr, pProgressData))
565 0 : return nullptr;
566 :
567 : /* -------------------------------------------------------------------- */
568 : /* Get source dataset info */
569 : /* -------------------------------------------------------------------- */
570 :
571 17 : const int nXSize = poSrcDS->GetRasterXSize();
572 17 : const int nYSize = poSrcDS->GetRasterYSize();
573 17 : if (nXSize == 1 || nYSize == 1)
574 : {
575 0 : return nullptr;
576 : }
577 :
578 : double adfGeoTransform[6];
579 17 : poSrcDS->GetGeoTransform(adfGeoTransform);
580 17 : if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
581 : {
582 0 : CPLError(CE_Failure, CPLE_NotSupported,
583 : "ZMap driver does not support CreateCopy() from skewed or "
584 : "rotated dataset.\n");
585 0 : return nullptr;
586 : }
587 :
588 : /* -------------------------------------------------------------------- */
589 : /* Create target file */
590 : /* -------------------------------------------------------------------- */
591 :
592 34 : auto fp = VSIVirtualHandleUniquePtr(VSIFOpenL(pszFilename, "wb"));
593 17 : if (fp == nullptr)
594 : {
595 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszFilename);
596 3 : return nullptr;
597 : }
598 :
599 14 : const int nFieldSize = 20;
600 14 : const int nValuesPerLine = 4;
601 14 : const int nDecimalCount = 7;
602 :
603 14 : int bHasNoDataValue = FALSE;
604 : double dfNoDataValue =
605 14 : poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoDataValue);
606 14 : if (!bHasNoDataValue)
607 13 : dfNoDataValue = 1.e30;
608 :
609 14 : fp->Printf("!\n");
610 14 : fp->Printf("! Created by GDAL.\n");
611 14 : fp->Printf("!\n");
612 14 : fp->Printf("@GRID FILE, GRID, %d\n", nValuesPerLine);
613 :
614 14 : WriteRightJustified(fp, nFieldSize, 10);
615 14 : fp->Printf(",");
616 14 : WriteRightJustified(fp, dfNoDataValue, nFieldSize, nDecimalCount);
617 14 : fp->Printf(",");
618 14 : WriteRightJustified(fp, "", 10);
619 14 : fp->Printf(",");
620 14 : WriteRightJustified(fp, nDecimalCount, 10);
621 14 : fp->Printf(",");
622 14 : WriteRightJustified(fp, 1, 10);
623 14 : fp->Printf("\n");
624 :
625 14 : WriteRightJustified(fp, nYSize, 10);
626 14 : fp->Printf(",");
627 14 : WriteRightJustified(fp, nXSize, 10);
628 14 : fp->Printf(",");
629 :
630 14 : if (CPLTestBool(CPLGetConfigOption("ZMAP_PIXEL_IS_POINT", "FALSE")))
631 : {
632 0 : WriteRightJustified(fp, adfGeoTransform[0] + adfGeoTransform[1] / 2, 14,
633 : 7);
634 0 : fp->Printf(",");
635 0 : WriteRightJustified(fp,
636 0 : adfGeoTransform[0] + adfGeoTransform[1] * nXSize -
637 0 : adfGeoTransform[1] / 2,
638 : 14, 7);
639 0 : fp->Printf(",");
640 0 : WriteRightJustified(fp,
641 0 : adfGeoTransform[3] + adfGeoTransform[5] * nYSize -
642 0 : adfGeoTransform[5] / 2,
643 : 14, 7);
644 0 : fp->Printf(",");
645 0 : WriteRightJustified(fp, adfGeoTransform[3] + adfGeoTransform[5] / 2, 14,
646 : 7);
647 : }
648 : else
649 : {
650 14 : WriteRightJustified(fp, adfGeoTransform[0], 14, 7);
651 14 : fp->Printf(",");
652 14 : WriteRightJustified(
653 14 : fp, adfGeoTransform[0] + adfGeoTransform[1] * nXSize, 14, 7);
654 14 : fp->Printf(",");
655 14 : WriteRightJustified(
656 14 : fp, adfGeoTransform[3] + adfGeoTransform[5] * nYSize, 14, 7);
657 14 : fp->Printf(",");
658 14 : WriteRightJustified(fp, adfGeoTransform[3], 14, 7);
659 : }
660 :
661 14 : fp->Printf("\n");
662 :
663 14 : fp->Printf("0.0, 0.0, 0.0\n");
664 14 : fp->Printf("@\n");
665 :
666 : /* -------------------------------------------------------------------- */
667 : /* Copy imagery */
668 : /* -------------------------------------------------------------------- */
669 28 : std::vector<double> adfLineBuffer(nYSize);
670 :
671 14 : CPLErr eErr = CE_None;
672 14 : const bool bEmitEOLAtEndOfColumn = CPLTestBool(
673 : CPLGetConfigOption("ZMAP_EMIT_EOL_AT_END_OF_COLUMN", "YES"));
674 14 : bool bEOLPrinted = false;
675 14 : int nValuesThisLine = 0;
676 148 : for (int i = 0; i < nXSize && eErr == CE_None; i++)
677 : {
678 268 : eErr = poSrcDS->GetRasterBand(1)->RasterIO(
679 134 : GF_Read, i, 0, 1, nYSize, adfLineBuffer.data(), 1, nYSize,
680 : GDT_Float64, 0, 0, nullptr);
681 134 : if (eErr != CE_None)
682 0 : break;
683 1648 : for (int j = 0; j < nYSize; j++)
684 : {
685 1514 : WriteRightJustified(fp, adfLineBuffer[j], nFieldSize,
686 : nDecimalCount);
687 1514 : ++nValuesThisLine;
688 1514 : if (nValuesThisLine == nValuesPerLine)
689 : {
690 322 : bEOLPrinted = true;
691 322 : nValuesThisLine = 0;
692 322 : fp->Printf("\n");
693 : }
694 : else
695 1192 : bEOLPrinted = false;
696 : }
697 134 : if (bEmitEOLAtEndOfColumn && !bEOLPrinted)
698 : {
699 112 : bEOLPrinted = true;
700 112 : nValuesThisLine = 0;
701 112 : fp->Printf("\n");
702 : }
703 :
704 268 : if (pfnProgress != nullptr &&
705 134 : !pfnProgress((i + 1) * 1.0 / nXSize, nullptr, pProgressData))
706 : {
707 0 : eErr = CE_Failure;
708 0 : break;
709 : }
710 : }
711 14 : if (!bEOLPrinted)
712 1 : fp->Printf("\n");
713 :
714 14 : if (eErr != CE_None || fp->Close() != 0)
715 0 : return nullptr;
716 :
717 14 : fp.reset();
718 28 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
719 14 : return ZMapDataset::Open(&oOpenInfo);
720 : }
721 :
722 : /************************************************************************/
723 : /* GetGeoTransform() */
724 : /************************************************************************/
725 :
726 1 : CPLErr ZMapDataset::GetGeoTransform(double *padfTransform)
727 :
728 : {
729 1 : memcpy(padfTransform, m_adfGeoTransform.data(), 6 * sizeof(double));
730 :
731 1 : return CE_None;
732 : }
733 :
734 : /************************************************************************/
735 : /* GDALRegister_ZMap() */
736 : /************************************************************************/
737 :
738 1686 : void GDALRegister_ZMap()
739 :
740 : {
741 1686 : if (GDALGetDriverByName("ZMap") != nullptr)
742 302 : return;
743 :
744 2768 : auto poDriver = std::make_unique<GDALDriver>();
745 :
746 1384 : poDriver->SetDescription("ZMap");
747 1384 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
748 1384 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ZMap Plus Grid");
749 1384 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/zmap.html");
750 1384 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dat");
751 1384 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
752 :
753 1384 : poDriver->pfnOpen = ZMapDataset::Open;
754 1384 : poDriver->pfnIdentify = ZMapDataset::Identify;
755 1384 : poDriver->pfnCreateCopy = ZMapDataset::CreateCopy;
756 :
757 1384 : GetGDALDriverManager()->RegisterDriver(poDriver.release());
758 : }
|
