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