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 50684 : int ZMapDataset::Identify(GDALOpenInfo *poOpenInfo)
233 : {
234 50684 : if (poOpenInfo->nHeaderBytes == 0)
235 47705 : return FALSE;
236 :
237 : /* -------------------------------------------------------------------- */
238 : /* Check that it looks roughly as a ZMap dataset */
239 : /* -------------------------------------------------------------------- */
240 2979 : const char *pszData =
241 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
242 :
243 : /* Skip comments line at the beginning */
244 2979 : int i = 0;
245 2979 : 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 2979 : if (pszData[i] != '@')
263 2953 : return FALSE;
264 26 : i++;
265 :
266 54 : 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 : CPLError(CE_Failure, CPLE_NotSupported,
295 : "The ZMAP driver does not support update access to existing"
296 : " datasets.");
297 0 : return nullptr;
298 : }
299 :
300 38 : auto poDS = std::make_unique<ZMapDataset>();
301 19 : poDS->m_fp.reset(poOpenInfo->fpL);
302 19 : poOpenInfo->fpL = nullptr;
303 :
304 : /* -------------------------------------------------------------------- */
305 : /* Find dataset characteristics */
306 : /* -------------------------------------------------------------------- */
307 :
308 : const char *pszLine;
309 19 : int nLine = 0;
310 19 : constexpr int MAX_HEADER_LINE = 1024;
311 76 : while ((pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE,
312 76 : nullptr)) != nullptr)
313 : {
314 76 : ++nLine;
315 76 : if (*pszLine == '!')
316 : {
317 57 : continue;
318 : }
319 : else
320 19 : break;
321 : }
322 : // cppcheck-suppress knownConditionTrueFalse
323 19 : if (pszLine == nullptr)
324 : {
325 0 : return nullptr;
326 : }
327 :
328 : /* Parse first header line */
329 38 : CPLStringList aosTokensFirstLine(CSLTokenizeString2(pszLine, ",", 0));
330 19 : if (aosTokensFirstLine.size() != 3)
331 : {
332 0 : return nullptr;
333 : }
334 :
335 19 : const int nValuesPerLine = atoi(aosTokensFirstLine[2]);
336 19 : if (nValuesPerLine <= 0)
337 : {
338 0 : CPLError(CE_Failure, CPLE_AppDefined,
339 : "Invalid/unsupported value for nValuesPerLine = %d",
340 : nValuesPerLine);
341 0 : return nullptr;
342 : }
343 :
344 : /* Parse second header line */
345 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
346 19 : ++nLine;
347 19 : if (pszLine == nullptr)
348 : {
349 0 : return nullptr;
350 : }
351 : const CPLStringList aosTokensSecondLine(
352 38 : CSLTokenizeString2(pszLine, ",", 0));
353 19 : if (aosTokensSecondLine.size() != 5)
354 : {
355 0 : return nullptr;
356 : }
357 :
358 19 : const int nFieldSize = atoi(aosTokensSecondLine[0]);
359 19 : const double dfNoDataValue = CPLAtofM(aosTokensSecondLine[1]);
360 19 : const int nDecimalCount = atoi(aosTokensSecondLine[3]);
361 19 : const int nColumnNumber = atoi(aosTokensSecondLine[4]);
362 :
363 19 : if (nFieldSize <= 0 || nFieldSize >= 40)
364 : {
365 0 : CPLError(CE_Failure, CPLE_AppDefined,
366 : "Invalid/unsupported value for nFieldSize = %d", nFieldSize);
367 0 : return nullptr;
368 : }
369 :
370 19 : if (nDecimalCount <= 0 || nDecimalCount >= nFieldSize)
371 : {
372 0 : CPLError(CE_Failure, CPLE_AppDefined,
373 : "Invalid/unsupported value for nDecimalCount = %d",
374 : nDecimalCount);
375 0 : return nullptr;
376 : }
377 :
378 19 : if (nColumnNumber != 1)
379 : {
380 0 : CPLError(CE_Failure, CPLE_AppDefined,
381 : "Invalid/unsupported value for nColumnNumber = %d",
382 : nColumnNumber);
383 0 : return nullptr;
384 : }
385 :
386 19 : if (nValuesPerLine <= 0 || nFieldSize > 1024 * 1024 / nValuesPerLine)
387 : {
388 0 : CPLError(CE_Failure, CPLE_AppDefined,
389 : "Invalid/unsupported value for nFieldSize = %d x "
390 : "nValuesPerLine = %d",
391 : nFieldSize, nValuesPerLine);
392 0 : return nullptr;
393 : }
394 :
395 : /* Parse third header line */
396 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
397 19 : ++nLine;
398 19 : if (pszLine == nullptr)
399 : {
400 0 : return nullptr;
401 : }
402 38 : const CPLStringList aosTokensThirdLine(CSLTokenizeString2(pszLine, ",", 0));
403 19 : if (aosTokensThirdLine.size() != 6)
404 : {
405 0 : return nullptr;
406 : }
407 :
408 19 : const int nRows = atoi(aosTokensThirdLine[0]);
409 19 : const int nCols = atoi(aosTokensThirdLine[1]);
410 19 : const double dfMinX = CPLAtofM(aosTokensThirdLine[2]);
411 19 : const double dfMaxX = CPLAtofM(aosTokensThirdLine[3]);
412 19 : const double dfMinY = CPLAtofM(aosTokensThirdLine[4]);
413 19 : const double dfMaxY = CPLAtofM(aosTokensThirdLine[5]);
414 :
415 19 : if (!GDALCheckDatasetDimensions(nCols, nRows) || nCols == 1 || nRows == 1)
416 : {
417 0 : return nullptr;
418 : }
419 :
420 : /* Ignore fourth header line */
421 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
422 19 : ++nLine;
423 19 : if (pszLine == nullptr)
424 : {
425 0 : return nullptr;
426 : }
427 :
428 : /* Check fifth header line */
429 19 : pszLine = CPLReadLine2L(poDS->m_fp.get(), MAX_HEADER_LINE, nullptr);
430 19 : ++nLine;
431 19 : if (pszLine == nullptr || pszLine[0] != '@')
432 : {
433 0 : return nullptr;
434 : }
435 :
436 : /* -------------------------------------------------------------------- */
437 : /* Create a corresponding GDALDataset. */
438 : /* -------------------------------------------------------------------- */
439 19 : poDS->m_nDataStartOff = VSIFTellL(poDS->m_fp.get());
440 19 : poDS->m_nValuesPerLine = nValuesPerLine;
441 19 : poDS->m_nFieldSize = nFieldSize;
442 19 : poDS->m_nDecimalCount = nDecimalCount;
443 19 : poDS->nRasterXSize = nCols;
444 19 : poDS->nRasterYSize = nRows;
445 19 : poDS->m_dfNoDataValue = dfNoDataValue;
446 19 : poDS->m_nFirstDataLine = nLine;
447 :
448 19 : if (CPLTestBool(CPLGetConfigOption("ZMAP_PIXEL_IS_POINT", "FALSE")))
449 : {
450 0 : const double dfStepX = (dfMaxX - dfMinX) / (nCols - 1);
451 0 : const double dfStepY = (dfMaxY - dfMinY) / (nRows - 1);
452 :
453 0 : poDS->m_adfGeoTransform[0] = dfMinX - dfStepX / 2;
454 0 : poDS->m_adfGeoTransform[1] = dfStepX;
455 0 : poDS->m_adfGeoTransform[3] = dfMaxY + dfStepY / 2;
456 0 : poDS->m_adfGeoTransform[5] = -dfStepY;
457 : }
458 : else
459 : {
460 19 : const double dfStepX = (dfMaxX - dfMinX) / nCols;
461 19 : const double dfStepY = (dfMaxY - dfMinY) / nRows;
462 :
463 19 : poDS->m_adfGeoTransform[0] = dfMinX;
464 19 : poDS->m_adfGeoTransform[1] = dfStepX;
465 19 : poDS->m_adfGeoTransform[3] = dfMaxY;
466 19 : poDS->m_adfGeoTransform[5] = -dfStepY;
467 : }
468 :
469 : /* -------------------------------------------------------------------- */
470 : /* Create band information objects. */
471 : /* -------------------------------------------------------------------- */
472 19 : poDS->nBands = 1;
473 19 : poDS->SetBand(1, std::make_unique<ZMapRasterBand>(poDS.get()));
474 :
475 : /* -------------------------------------------------------------------- */
476 : /* Initialize any PAM information. */
477 : /* -------------------------------------------------------------------- */
478 19 : poDS->SetDescription(poOpenInfo->pszFilename);
479 19 : poDS->TryLoadXML();
480 :
481 : /* -------------------------------------------------------------------- */
482 : /* Support overviews. */
483 : /* -------------------------------------------------------------------- */
484 19 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
485 19 : return poDS.release();
486 : }
487 :
488 : /************************************************************************/
489 : /* WriteRightJustified() */
490 : /************************************************************************/
491 :
492 1668 : static void WriteRightJustified(VSIVirtualHandleUniquePtr &fp,
493 : const char *pszValue, int nWidth)
494 : {
495 1668 : int nLen = (int)strlen(pszValue);
496 1668 : CPLAssert(nLen <= nWidth);
497 16550 : for (int i = 0; i < nWidth - nLen; i++)
498 14882 : fp->Write(" ", 1, 1);
499 1668 : fp->Write(pszValue, 1, nLen);
500 1668 : }
501 :
502 70 : static void WriteRightJustified(VSIVirtualHandleUniquePtr &fp, int nValue,
503 : int nWidth)
504 : {
505 140 : CPLString osValue(CPLSPrintf("%d", nValue));
506 70 : WriteRightJustified(fp, osValue.c_str(), nWidth);
507 70 : }
508 :
509 1584 : static void WriteRightJustified(VSIVirtualHandleUniquePtr &fp, double dfValue,
510 : int nWidth, int nDecimals = -1)
511 : {
512 : char szFormat[32];
513 1584 : if (nDecimals >= 0)
514 1584 : snprintf(szFormat, sizeof(szFormat), "%%.%df", nDecimals);
515 : else
516 0 : snprintf(szFormat, sizeof(szFormat), "%%g");
517 1584 : char *pszValue = const_cast<char *>(CPLSPrintf(szFormat, dfValue));
518 1584 : char *pszE = strchr(pszValue, 'e');
519 1584 : if (pszE)
520 0 : *pszE = 'E';
521 :
522 1584 : if (static_cast<int>(strlen(pszValue)) > nWidth)
523 : {
524 16 : CPLAssert(nDecimals >= 0);
525 16 : snprintf(szFormat, sizeof(szFormat), "%%.%dg", nDecimals);
526 16 : pszValue = const_cast<char *>(CPLSPrintf(szFormat, dfValue));
527 16 : pszE = strchr(pszValue, 'e');
528 16 : if (pszE)
529 14 : *pszE = 'E';
530 : }
531 :
532 3168 : CPLString osValue(pszValue);
533 1584 : WriteRightJustified(fp, osValue.c_str(), nWidth);
534 1584 : }
535 :
536 : /************************************************************************/
537 : /* CreateCopy() */
538 : /************************************************************************/
539 :
540 22 : GDALDataset *ZMapDataset::CreateCopy(const char *pszFilename,
541 : GDALDataset *poSrcDS, int bStrict,
542 : CPL_UNUSED char **papszOptions,
543 : GDALProgressFunc pfnProgress,
544 : void *pProgressData)
545 : {
546 : /* -------------------------------------------------------------------- */
547 : /* Some some rudimentary checks */
548 : /* -------------------------------------------------------------------- */
549 22 : int nBands = poSrcDS->GetRasterCount();
550 22 : if (nBands == 0)
551 : {
552 1 : CPLError(
553 : CE_Failure, CPLE_NotSupported,
554 : "ZMap driver does not support source dataset with zero band.\n");
555 1 : return nullptr;
556 : }
557 :
558 21 : if (nBands != 1)
559 : {
560 4 : CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
561 : "ZMap driver only uses the first band of the dataset.\n");
562 4 : if (bStrict)
563 4 : return nullptr;
564 : }
565 :
566 17 : if (pfnProgress && !pfnProgress(0.0, nullptr, pProgressData))
567 0 : return nullptr;
568 :
569 : /* -------------------------------------------------------------------- */
570 : /* Get source dataset info */
571 : /* -------------------------------------------------------------------- */
572 :
573 17 : const int nXSize = poSrcDS->GetRasterXSize();
574 17 : const int nYSize = poSrcDS->GetRasterYSize();
575 17 : if (nXSize == 1 || nYSize == 1)
576 : {
577 0 : return nullptr;
578 : }
579 :
580 : double adfGeoTransform[6];
581 17 : poSrcDS->GetGeoTransform(adfGeoTransform);
582 17 : if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
583 : {
584 0 : CPLError(CE_Failure, CPLE_NotSupported,
585 : "ZMap driver does not support CreateCopy() from skewed or "
586 : "rotated dataset.\n");
587 0 : return nullptr;
588 : }
589 :
590 : /* -------------------------------------------------------------------- */
591 : /* Create target file */
592 : /* -------------------------------------------------------------------- */
593 :
594 34 : auto fp = VSIVirtualHandleUniquePtr(VSIFOpenL(pszFilename, "wb"));
595 17 : if (fp == nullptr)
596 : {
597 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszFilename);
598 3 : return nullptr;
599 : }
600 :
601 14 : const int nFieldSize = 20;
602 14 : const int nValuesPerLine = 4;
603 14 : const int nDecimalCount = 7;
604 :
605 14 : int bHasNoDataValue = FALSE;
606 : double dfNoDataValue =
607 14 : poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoDataValue);
608 14 : if (!bHasNoDataValue)
609 13 : dfNoDataValue = 1.e30;
610 :
611 14 : fp->Printf("!\n");
612 14 : fp->Printf("! Created by GDAL.\n");
613 14 : fp->Printf("!\n");
614 14 : fp->Printf("@GRID FILE, GRID, %d\n", nValuesPerLine);
615 :
616 14 : WriteRightJustified(fp, nFieldSize, 10);
617 14 : fp->Printf(",");
618 14 : WriteRightJustified(fp, dfNoDataValue, nFieldSize, nDecimalCount);
619 14 : fp->Printf(",");
620 14 : WriteRightJustified(fp, "", 10);
621 14 : fp->Printf(",");
622 14 : WriteRightJustified(fp, nDecimalCount, 10);
623 14 : fp->Printf(",");
624 14 : WriteRightJustified(fp, 1, 10);
625 14 : fp->Printf("\n");
626 :
627 14 : WriteRightJustified(fp, nYSize, 10);
628 14 : fp->Printf(",");
629 14 : WriteRightJustified(fp, nXSize, 10);
630 14 : fp->Printf(",");
631 :
632 14 : if (CPLTestBool(CPLGetConfigOption("ZMAP_PIXEL_IS_POINT", "FALSE")))
633 : {
634 0 : WriteRightJustified(fp, adfGeoTransform[0] + adfGeoTransform[1] / 2, 14,
635 : 7);
636 0 : fp->Printf(",");
637 0 : WriteRightJustified(fp,
638 0 : adfGeoTransform[0] + adfGeoTransform[1] * nXSize -
639 0 : adfGeoTransform[1] / 2,
640 : 14, 7);
641 0 : fp->Printf(",");
642 0 : WriteRightJustified(fp,
643 0 : adfGeoTransform[3] + adfGeoTransform[5] * nYSize -
644 0 : adfGeoTransform[5] / 2,
645 : 14, 7);
646 0 : fp->Printf(",");
647 0 : WriteRightJustified(fp, adfGeoTransform[3] + adfGeoTransform[5] / 2, 14,
648 : 7);
649 : }
650 : else
651 : {
652 14 : WriteRightJustified(fp, adfGeoTransform[0], 14, 7);
653 14 : fp->Printf(",");
654 14 : WriteRightJustified(
655 14 : fp, adfGeoTransform[0] + adfGeoTransform[1] * nXSize, 14, 7);
656 14 : fp->Printf(",");
657 14 : WriteRightJustified(
658 14 : fp, adfGeoTransform[3] + adfGeoTransform[5] * nYSize, 14, 7);
659 14 : fp->Printf(",");
660 14 : WriteRightJustified(fp, adfGeoTransform[3], 14, 7);
661 : }
662 :
663 14 : fp->Printf("\n");
664 :
665 14 : fp->Printf("0.0, 0.0, 0.0\n");
666 14 : fp->Printf("@\n");
667 :
668 : /* -------------------------------------------------------------------- */
669 : /* Copy imagery */
670 : /* -------------------------------------------------------------------- */
671 28 : std::vector<double> adfLineBuffer(nYSize);
672 :
673 14 : CPLErr eErr = CE_None;
674 14 : const bool bEmitEOLAtEndOfColumn = CPLTestBool(
675 : CPLGetConfigOption("ZMAP_EMIT_EOL_AT_END_OF_COLUMN", "YES"));
676 14 : bool bEOLPrinted = false;
677 14 : int nValuesThisLine = 0;
678 148 : for (int i = 0; i < nXSize && eErr == CE_None; i++)
679 : {
680 268 : eErr = poSrcDS->GetRasterBand(1)->RasterIO(
681 134 : GF_Read, i, 0, 1, nYSize, adfLineBuffer.data(), 1, nYSize,
682 : GDT_Float64, 0, 0, nullptr);
683 134 : if (eErr != CE_None)
684 0 : break;
685 1648 : for (int j = 0; j < nYSize; j++)
686 : {
687 1514 : WriteRightJustified(fp, adfLineBuffer[j], nFieldSize,
688 : nDecimalCount);
689 1514 : ++nValuesThisLine;
690 1514 : if (nValuesThisLine == nValuesPerLine)
691 : {
692 322 : bEOLPrinted = true;
693 322 : nValuesThisLine = 0;
694 322 : fp->Printf("\n");
695 : }
696 : else
697 1192 : bEOLPrinted = false;
698 : }
699 134 : if (bEmitEOLAtEndOfColumn && !bEOLPrinted)
700 : {
701 112 : bEOLPrinted = true;
702 112 : nValuesThisLine = 0;
703 112 : fp->Printf("\n");
704 : }
705 :
706 268 : if (pfnProgress != nullptr &&
707 134 : !pfnProgress((i + 1) * 1.0 / nXSize, nullptr, pProgressData))
708 : {
709 0 : eErr = CE_Failure;
710 0 : break;
711 : }
712 : }
713 14 : if (!bEOLPrinted)
714 1 : fp->Printf("\n");
715 :
716 14 : if (eErr != CE_None || fp->Close() != 0)
717 0 : return nullptr;
718 :
719 14 : fp.reset();
720 28 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
721 14 : return ZMapDataset::Open(&oOpenInfo);
722 : }
723 :
724 : /************************************************************************/
725 : /* GetGeoTransform() */
726 : /************************************************************************/
727 :
728 1 : CPLErr ZMapDataset::GetGeoTransform(double *padfTransform)
729 :
730 : {
731 1 : memcpy(padfTransform, m_adfGeoTransform.data(), 6 * sizeof(double));
732 :
733 1 : return CE_None;
734 : }
735 :
736 : /************************************************************************/
737 : /* GDALRegister_ZMap() */
738 : /************************************************************************/
739 :
740 1595 : void GDALRegister_ZMap()
741 :
742 : {
743 1595 : if (GDALGetDriverByName("ZMap") != nullptr)
744 302 : return;
745 :
746 2586 : auto poDriver = std::make_unique<GDALDriver>();
747 :
748 1293 : poDriver->SetDescription("ZMap");
749 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
750 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ZMap Plus Grid");
751 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/zmap.html");
752 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dat");
753 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
754 :
755 1293 : poDriver->pfnOpen = ZMapDataset::Open;
756 1293 : poDriver->pfnIdentify = ZMapDataset::Identify;
757 1293 : poDriver->pfnCreateCopy = ZMapDataset::CreateCopy;
758 :
759 1293 : GetGDALDriverManager()->RegisterDriver(poDriver.release());
760 : }
|
