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