1 : /******************************************************************************
2 : * $Id: xyzdataset.cpp 20000 2010-07-10 18:19:34Z rouault $
3 : *
4 : * Project: XYZ driver
5 : * Purpose: GDALDataset driver for XYZ dataset.
6 : * Author: Even Rouault, <even dot rouault at mines dash paris dot org>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2010, Even Rouault
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_vsi_virtual.h"
31 : #include "cpl_string.h"
32 : #include "gdal_pam.h"
33 :
34 : CPL_CVSID("$Id: xyzdataset.cpp 20000 2010-07-10 18:19:34Z rouault $");
35 :
36 : CPL_C_START
37 : void GDALRegister_XYZ(void);
38 : CPL_C_END
39 :
40 : /************************************************************************/
41 : /* ==================================================================== */
42 : /* XYZDataset */
43 : /* ==================================================================== */
44 : /************************************************************************/
45 :
46 : class XYZRasterBand;
47 :
48 : class XYZDataset : public GDALPamDataset
49 : {
50 : friend class XYZRasterBand;
51 :
52 : FILE *fp;
53 : int bHasHeaderLine;
54 : int nXIndex;
55 : int nYIndex;
56 : int nZIndex;
57 : int nMinTokens;
58 : int nLineNum; /* any line */
59 : int nDataLineNum; /* line with values (header line and empty lines ignored) */
60 : double adfGeoTransform[6];
61 :
62 : static int IdentifyEx( GDALOpenInfo *, int& );
63 :
64 : public:
65 : XYZDataset();
66 : virtual ~XYZDataset();
67 :
68 : virtual CPLErr GetGeoTransform( double * );
69 :
70 : static GDALDataset *Open( GDALOpenInfo * );
71 : static int Identify( GDALOpenInfo * );
72 : static GDALDataset *CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
73 : int bStrict, char ** papszOptions,
74 : GDALProgressFunc pfnProgress, void * pProgressData );
75 : };
76 :
77 : /************************************************************************/
78 : /* ==================================================================== */
79 : /* XYZRasterBand */
80 : /* ==================================================================== */
81 : /************************************************************************/
82 :
83 : class XYZRasterBand : public GDALPamRasterBand
84 18 : {
85 : friend class XYZDataset;
86 :
87 : public:
88 :
89 : XYZRasterBand( XYZDataset *, int, GDALDataType );
90 :
91 : virtual CPLErr IReadBlock( int, int, void * );
92 : };
93 :
94 :
95 : /************************************************************************/
96 : /* XYZRasterBand() */
97 : /************************************************************************/
98 :
99 18 : XYZRasterBand::XYZRasterBand( XYZDataset *poDS, int nBand, GDALDataType eDT )
100 :
101 : {
102 18 : this->poDS = poDS;
103 18 : this->nBand = nBand;
104 :
105 18 : eDataType = eDT;
106 :
107 18 : nBlockXSize = poDS->GetRasterXSize();
108 18 : nBlockYSize = 1;
109 18 : }
110 :
111 : /************************************************************************/
112 : /* IReadBlock() */
113 : /************************************************************************/
114 :
115 : CPLErr XYZRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
116 445 : void * pImage )
117 :
118 : {
119 445 : XYZDataset *poGDS = (XYZDataset *) poDS;
120 :
121 445 : int nLineInFile = nBlockYOff * nBlockXSize;
122 445 : if (poGDS->nDataLineNum > nLineInFile)
123 : {
124 7 : poGDS->nDataLineNum = 0;
125 7 : VSIFSeekL(poGDS->fp, 0, SEEK_SET);
126 7 : if (poGDS->bHasHeaderLine)
127 : {
128 5 : const char* pszLine = CPLReadLine2L(poGDS->fp, 100, 0);
129 5 : if (pszLine == NULL)
130 : {
131 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
132 0 : return CE_Failure;
133 : }
134 5 : poGDS->nLineNum ++;
135 : }
136 : }
137 :
138 902 : while(poGDS->nDataLineNum < nLineInFile)
139 : {
140 12 : const char* pszLine = CPLReadLine2L(poGDS->fp, 100, 0);
141 12 : if (pszLine == NULL)
142 : {
143 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
144 0 : return CE_Failure;
145 : }
146 12 : poGDS->nLineNum ++;
147 :
148 12 : const char* pszPtr = pszLine;
149 : char ch;
150 12 : int nCol = 0;
151 12 : int bLastWasSep = TRUE;
152 60 : while((ch = *pszPtr) != '\0')
153 : {
154 48 : if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
155 : {
156 12 : if (!bLastWasSep)
157 12 : nCol ++;
158 12 : bLastWasSep = TRUE;
159 : }
160 : else
161 : {
162 24 : bLastWasSep = FALSE;
163 : }
164 36 : pszPtr ++;
165 : }
166 :
167 : /* Skip empty line */
168 12 : if (nCol == 0 && bLastWasSep)
169 6 : continue;
170 :
171 6 : poGDS->nDataLineNum ++;
172 : }
173 :
174 : int i;
175 82053 : for(i=0;i<nBlockXSize;i++)
176 : {
177 : int nCol;
178 : int bLastWasSep;
179 81612 : do
180 : {
181 81612 : const char* pszLine = CPLReadLine2L(poGDS->fp, 100, 0);
182 81612 : if (pszLine == NULL)
183 : {
184 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
185 0 : return CE_Failure;
186 : }
187 81612 : poGDS->nLineNum ++;
188 :
189 81612 : const char* pszPtr = pszLine;
190 : char ch;
191 81612 : nCol = 0;
192 81612 : bLastWasSep = TRUE;
193 2799048 : while((ch = *pszPtr) != '\0')
194 : {
195 2799040 : if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
196 : {
197 163216 : if (!bLastWasSep)
198 163216 : nCol ++;
199 163216 : bLastWasSep = TRUE;
200 : }
201 : else
202 : {
203 2472608 : if (bLastWasSep && nCol == poGDS->nZIndex)
204 : {
205 81608 : double dfZ = CPLAtofM(pszPtr);
206 81608 : if (eDataType == GDT_Float32)
207 : {
208 80802 : ((float*)pImage)[i] = (float)dfZ;
209 : }
210 806 : else if (eDataType == GDT_Int32)
211 : {
212 0 : ((GInt32*)pImage)[i] = (GInt32)dfZ;
213 : }
214 806 : else if (eDataType == GDT_Int16)
215 : {
216 0 : ((GInt16*)pImage)[i] = (GInt16)dfZ;
217 : }
218 : else
219 : {
220 806 : ((GByte*)pImage)[i] = (GByte)dfZ;
221 : }
222 : }
223 2472608 : bLastWasSep = FALSE;
224 : }
225 2635824 : pszPtr ++;
226 : }
227 :
228 : /* Skip empty line */
229 : }
230 : while (nCol == 0 && bLastWasSep);
231 :
232 81608 : poGDS->nDataLineNum ++;
233 81608 : nCol ++;
234 81608 : if (nCol < poGDS->nMinTokens)
235 : {
236 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
237 0 : return CE_Failure;
238 : }
239 : }
240 445 : CPLAssert(poGDS->nDataLineNum == (nBlockYOff + 1) * nBlockXSize);
241 :
242 445 : return CE_None;
243 : }
244 :
245 : /************************************************************************/
246 : /* ~XYZDataset() */
247 : /************************************************************************/
248 :
249 18 : XYZDataset::XYZDataset()
250 : {
251 18 : fp = NULL;
252 18 : nDataLineNum = INT_MAX;
253 18 : nLineNum = 0;
254 18 : bHasHeaderLine = FALSE;
255 18 : nXIndex = -1;
256 18 : nYIndex = -1;
257 18 : nZIndex = -1;
258 18 : nMinTokens = 0;
259 18 : adfGeoTransform[0] = 0;
260 18 : adfGeoTransform[1] = 1;
261 18 : adfGeoTransform[2] = 0;
262 18 : adfGeoTransform[3] = 0;
263 18 : adfGeoTransform[4] = 0;
264 18 : adfGeoTransform[5] = 1;
265 18 : }
266 :
267 : /************************************************************************/
268 : /* ~XYZDataset() */
269 : /************************************************************************/
270 :
271 18 : XYZDataset::~XYZDataset()
272 :
273 : {
274 18 : FlushCache();
275 18 : if (fp)
276 18 : VSIFCloseL(fp);
277 18 : }
278 :
279 : /************************************************************************/
280 : /* Identify() */
281 : /************************************************************************/
282 :
283 8238 : int XYZDataset::Identify( GDALOpenInfo * poOpenInfo )
284 : {
285 : int bHasHeaderLine;
286 8238 : return IdentifyEx(poOpenInfo, bHasHeaderLine);
287 : }
288 :
289 : /************************************************************************/
290 : /* IdentifyEx() */
291 : /************************************************************************/
292 :
293 :
294 9411 : int XYZDataset::IdentifyEx( GDALOpenInfo * poOpenInfo, int& bHasHeaderLine )
295 :
296 : {
297 : int i;
298 :
299 9411 : bHasHeaderLine = FALSE;
300 :
301 9411 : CPLString osFilename(poOpenInfo->pszFilename);
302 :
303 9411 : GDALOpenInfo* poOpenInfoToDelete = NULL;
304 : /* GZipped .xyz files are common, so automagically open them */
305 : /* if the /vsigzip/ has not been explicitely passed */
306 9411 : if (strlen(poOpenInfo->pszFilename) > 6 &&
307 : EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "xyz.gz") &&
308 : !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
309 : {
310 0 : osFilename = "/vsigzip/";
311 0 : osFilename += poOpenInfo->pszFilename;
312 : poOpenInfo = poOpenInfoToDelete =
313 : new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
314 0 : poOpenInfo->papszSiblingFiles);
315 : }
316 :
317 9411 : if (poOpenInfo->nHeaderBytes == 0)
318 : {
319 9236 : delete poOpenInfoToDelete;
320 9236 : return FALSE;
321 : }
322 :
323 : /* -------------------------------------------------------------------- */
324 : /* Chech that it looks roughly as a XYZ dataset */
325 : /* -------------------------------------------------------------------- */
326 175 : const char* pszData = (const char*)poOpenInfo->pabyHeader;
327 13185 : for(i=0;i<poOpenInfo->nHeaderBytes;i++)
328 : {
329 13173 : char ch = pszData[i];
330 13173 : if (ch == 13 || ch == 10)
331 : {
332 18 : break;
333 : }
334 13155 : else if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
335 : ;
336 867 : else if ((ch >= '0' && ch <= '9') || ch == '.' || ch == '+' ||
337 : ch == '-' || ch == 'e' || ch == 'E')
338 : ;
339 461 : else if (ch == '"' || (ch >= 'a' && ch <= 'z') ||
340 : (ch >= 'A' && ch <= 'Z'))
341 158 : bHasHeaderLine = TRUE;
342 : else
343 : {
344 145 : delete poOpenInfoToDelete;
345 145 : return FALSE;
346 : }
347 : }
348 30 : int bHasFoundNewLine = FALSE;
349 30 : int bPrevWasSep = TRUE;
350 30 : int nCols = 0;
351 30 : int nMaxCols = 0;
352 16945 : for(;i<poOpenInfo->nHeaderBytes;i++)
353 : {
354 16915 : char ch = pszData[i];
355 17534 : if (ch == 13 || ch == 10)
356 : {
357 619 : bHasFoundNewLine = TRUE;
358 619 : if (!bPrevWasSep)
359 : {
360 595 : nCols ++;
361 595 : if (nCols > nMaxCols)
362 18 : nMaxCols = nCols;
363 : }
364 619 : bPrevWasSep = TRUE;
365 619 : nCols = 0;
366 : }
367 17500 : else if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
368 : {
369 1204 : if (!bPrevWasSep)
370 : {
371 1204 : nCols ++;
372 1204 : if (nCols > nMaxCols)
373 36 : nMaxCols = nCols;
374 : }
375 1204 : bPrevWasSep = TRUE;
376 : }
377 30184 : else if ((ch >= '0' && ch <= '9') || ch == '.' || ch == '+' ||
378 : ch == '-' || ch == 'e' || ch == 'E')
379 : {
380 15092 : bPrevWasSep = FALSE;
381 : }
382 : else
383 : {
384 0 : delete poOpenInfoToDelete;
385 0 : return FALSE;
386 : }
387 : }
388 :
389 30 : delete poOpenInfoToDelete;
390 30 : return bHasFoundNewLine && nMaxCols >= 3;
391 : }
392 :
393 : /************************************************************************/
394 : /* Open() */
395 : /************************************************************************/
396 :
397 1173 : GDALDataset *XYZDataset::Open( GDALOpenInfo * poOpenInfo )
398 :
399 : {
400 : int i;
401 : int bHasHeaderLine;
402 :
403 1173 : if (!IdentifyEx(poOpenInfo, bHasHeaderLine))
404 1155 : return NULL;
405 :
406 18 : CPLString osFilename(poOpenInfo->pszFilename);
407 :
408 : /* GZipped .xyz files are common, so automagically open them */
409 : /* if the /vsigzip/ has not been explicitely passed */
410 18 : if (strlen(poOpenInfo->pszFilename) > 6 &&
411 : EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "xyz.gz") &&
412 : !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
413 : {
414 0 : osFilename = "/vsigzip/";
415 0 : osFilename += poOpenInfo->pszFilename;
416 : }
417 :
418 : /* -------------------------------------------------------------------- */
419 : /* Find dataset characteristics */
420 : /* -------------------------------------------------------------------- */
421 18 : FILE* fp = VSIFOpenL(osFilename.c_str(), "rb");
422 18 : if (fp == NULL)
423 18 : return NULL;
424 :
425 : /* For better performance of CPLReadLine2L() we create a buffered reader */
426 : /* (except for /vsigzip/ since it has one internally) */
427 18 : if (!EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
428 18 : fp = (FILE*) VSICreateBufferedReaderHandle((VSIVirtualHandle*)fp);
429 :
430 : const char* pszLine;
431 18 : int nXIndex = -1, nYIndex = -1, nZIndex = -1;
432 18 : int nMinTokens = 0;
433 :
434 : /* -------------------------------------------------------------------- */
435 : /* Parse header line */
436 : /* -------------------------------------------------------------------- */
437 18 : if (bHasHeaderLine)
438 : {
439 4 : pszLine = CPLReadLine2L(fp, 100, NULL);
440 4 : if (pszLine == NULL)
441 : {
442 0 : VSIFCloseL(fp);
443 0 : return NULL;
444 : }
445 : char** papszTokens = CSLTokenizeString2( pszLine, " ,\t;",
446 4 : CSLT_HONOURSTRINGS );
447 4 : int nTokens = CSLCount(papszTokens);
448 4 : if (nTokens < 3)
449 : {
450 : CPLError(CE_Failure, CPLE_AppDefined,
451 : "At line %d, found %d tokens. Expected 3 at least",
452 0 : 1, nTokens);
453 0 : CSLDestroy(papszTokens);
454 0 : VSIFCloseL(fp);
455 0 : return NULL;
456 : }
457 : int i;
458 16 : for(i=0;i<nTokens;i++)
459 : {
460 16 : if (EQUAL(papszTokens[i], "x") ||
461 : EQUALN(papszTokens[i], "lon", 3) ||
462 : EQUALN(papszTokens[i], "east", 4))
463 4 : nXIndex = i;
464 12 : else if (EQUAL(papszTokens[i], "y") ||
465 : EQUALN(papszTokens[i], "lat", 3) ||
466 : EQUALN(papszTokens[i], "north", 5))
467 4 : nYIndex = i;
468 4 : else if (EQUAL(papszTokens[i], "z") ||
469 : EQUALN(papszTokens[i], "alt", 3) ||
470 : EQUAL(papszTokens[i], "height"))
471 4 : nZIndex = i;
472 : }
473 4 : CSLDestroy(papszTokens);
474 4 : papszTokens = NULL;
475 4 : if (nXIndex < 0 || nYIndex < 0 || nZIndex < 0)
476 : {
477 : CPLError(CE_Warning, CPLE_AppDefined,
478 0 : "Could not find one of the X, Y or Z column names in header line. Defaulting to the first 3 columns");
479 0 : nXIndex = 0;
480 0 : nYIndex = 1;
481 0 : nZIndex = 2;
482 : }
483 4 : nMinTokens = 1 + MAX(MAX(nXIndex, nYIndex), nZIndex);
484 : }
485 : else
486 : {
487 14 : nXIndex = 0;
488 14 : nYIndex = 1;
489 14 : nZIndex = 2;
490 14 : nMinTokens = 3;
491 : }
492 :
493 : /* -------------------------------------------------------------------- */
494 : /* Parse data lines */
495 : /* -------------------------------------------------------------------- */
496 :
497 18 : int nXSize = 0, nYSize = 0;
498 18 : int nLineNum = 0;
499 18 : int nDataLineNum = 0;
500 18 : double dfFirstX = 0;
501 18 : double dfX = 0, dfY = 0, dfZ = 0;
502 18 : double dfMinX = 0, dfMinY = 0, dfMaxX = 0, dfMaxY = 0;
503 18 : double dfLastX = 0, dfLastY = 0;
504 18 : double dfStepX = 0, dfStepY = 0;
505 18 : GDALDataType eDT = GDT_Byte;
506 123551 : while((pszLine = CPLReadLine2L(fp, 100, NULL)) != NULL)
507 : {
508 123515 : nLineNum ++;
509 :
510 123515 : const char* pszPtr = pszLine;
511 : char ch;
512 123515 : int nCol = 0;
513 123515 : int bLastWasSep = TRUE;
514 4241668 : while((ch = *pszPtr) != '\0')
515 : {
516 4241656 : if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
517 : {
518 247018 : if (!bLastWasSep)
519 247018 : nCol ++;
520 247018 : bLastWasSep = TRUE;
521 : }
522 : else
523 : {
524 3747620 : if (bLastWasSep)
525 : {
526 370527 : if (nCol == nXIndex)
527 123509 : dfX = CPLAtofM(pszPtr);
528 247018 : else if (nCol == nYIndex)
529 123509 : dfY = CPLAtofM(pszPtr);
530 123509 : else if (nCol == nZIndex && eDT != GDT_Float32)
531 : {
532 15131 : dfZ = CPLAtofM(pszPtr);
533 15131 : int nZ = (int)dfZ;
534 15131 : if ((double)nZ != dfZ)
535 : {
536 3 : eDT = GDT_Float32;
537 : }
538 15128 : else if ((eDT == GDT_Byte || eDT == GDT_Int16) && (nZ < 0 || nZ > 255))
539 : {
540 0 : if (nZ < -32768 || nZ > 32767)
541 0 : eDT = GDT_Int32;
542 : else
543 0 : eDT = GDT_Int16;
544 : }
545 : }
546 : }
547 3747620 : bLastWasSep = FALSE;
548 : }
549 3994638 : pszPtr ++;
550 : }
551 : /* skip empty lines */
552 123515 : if (bLastWasSep && nCol == 0)
553 : {
554 6 : continue;
555 : }
556 123509 : nDataLineNum ++;
557 123509 : nCol ++;
558 123509 : if (nCol < nMinTokens)
559 : {
560 : CPLError(CE_Failure, CPLE_AppDefined,
561 : "At line %d, found %d tokens. Expected %d at least",
562 0 : nLineNum, nCol, nMinTokens);
563 0 : VSIFCloseL(fp);
564 0 : return NULL;
565 : }
566 :
567 123509 : if (nDataLineNum == 1)
568 : {
569 18 : dfFirstX = dfMinX = dfMaxX = dfX;
570 18 : dfMinY = dfMaxY = dfY;
571 : }
572 : else
573 : {
574 123491 : if (dfX < dfMinX) dfMinX = dfX;
575 123491 : if (dfX > dfMaxX) dfMaxX = dfX;
576 123491 : if (dfY < dfMinY) dfMinY = dfY;
577 123491 : if (dfY > dfMaxY) dfMaxY = dfY;
578 : }
579 123509 : if (nDataLineNum == 2)
580 : {
581 18 : dfStepX = dfX - dfLastX;
582 18 : if (dfStepX <= 0)
583 : {
584 : CPLError(CE_Failure, CPLE_AppDefined,
585 : "Ungridded dataset: At line %d, X spacing was %f. Expected >0 value",
586 0 : nLineNum, dfStepX);
587 0 : VSIFCloseL(fp);
588 0 : return NULL;
589 : }
590 : }
591 123491 : else if (nDataLineNum > 2)
592 : {
593 123473 : double dfNewStepX = dfX - dfLastX;
594 123473 : double dfNewStepY = dfY - dfLastY;
595 123473 : if (dfNewStepY != 0)
596 : {
597 758 : nYSize ++;
598 758 : if (dfStepY == 0)
599 : {
600 18 : nXSize = nDataLineNum - 1;
601 18 : double dfAdjustedStepX = (dfMaxX - dfMinX) / (nXSize - 1);
602 18 : if (fabs(dfStepX - dfAdjustedStepX) > 1e-8)
603 : {
604 0 : CPLDebug("XYZ", "Adjusting stepx from %f to %f", dfStepX, dfAdjustedStepX);
605 : }
606 18 : dfStepX = dfAdjustedStepX;
607 : }
608 758 : if (dfStepY != 0 && fabs(dfX - dfFirstX) > 1e-8)
609 : {
610 : CPLError(CE_Failure, CPLE_AppDefined,
611 : "Ungridded dataset: At line %d, X is %f, where as %f was expected",
612 0 : nLineNum, dfX, dfFirstX);
613 0 : VSIFCloseL(fp);
614 0 : return NULL;
615 : }
616 758 : if (dfStepY != 0 && fabs(dfLastX - dfMaxX) > 1e-8)
617 : {
618 : CPLError(CE_Failure, CPLE_AppDefined,
619 : "Ungridded dataset: At line %d, X is %f, where as %f was expected",
620 0 : nLineNum - 1, dfLastX, dfMaxX);
621 0 : VSIFCloseL(fp);
622 0 : return NULL;
623 : }
624 : /*if (dfStepY != 0 && fabs(dfNewStepY - dfStepY) > 1e-8)
625 : {
626 : CPLError(CE_Failure, CPLE_AppDefined,
627 : "Ungridded dataset: At line %d, Y spacing was %f, whereas it was %f before",
628 : nLineNum, dfNewStepY, dfStepY);
629 : VSIFCloseL(fp);
630 : return NULL;
631 : }*/
632 758 : dfStepY = dfNewStepY;
633 : }
634 : else if (dfNewStepX != 0)
635 : {
636 : /*if (dfStepX != 0 && fabs(dfNewStepX - dfStepX) > 1e-8)
637 : {
638 : CPLError(CE_Failure, CPLE_AppDefined,
639 : "At line %d, X spacing was %f, whereas it was %f before",
640 : nLineNum, dfNewStepX, dfStepX);
641 : VSIFCloseL(fp);
642 : return NULL;
643 : }*/
644 : }
645 : }
646 123509 : dfLastX = dfX;
647 123509 : dfLastY = dfY;
648 : }
649 18 : nYSize ++;
650 :
651 18 : if (dfStepX == 0)
652 : {
653 0 : CPLError(CE_Failure, CPLE_AppDefined, "Couldn't determine X spacing");
654 0 : VSIFCloseL(fp);
655 0 : return NULL;
656 : }
657 :
658 18 : if (dfStepY == 0)
659 : {
660 0 : CPLError(CE_Failure, CPLE_AppDefined, "Couldn't determine Y spacing");
661 0 : VSIFCloseL(fp);
662 0 : return NULL;
663 : }
664 :
665 18 : double dfAdjustedStepY = ((dfStepY < 0) ? -1 : 1) * (dfMaxY - dfMinY) / (nYSize - 1);
666 18 : if (fabs(dfStepY - dfAdjustedStepY) > 1e-8)
667 : {
668 0 : CPLDebug("XYZ", "Adjusting stepy from %f to %f", dfStepY, dfAdjustedStepY);
669 : }
670 18 : dfStepY = dfAdjustedStepY;
671 :
672 : //CPLDebug("XYZ", "minx=%f maxx=%f stepx=%f", dfMinX, dfMaxX, dfStepX);
673 : //CPLDebug("XYZ", "miny=%f maxy=%f stepy=%f", dfMinY, dfMaxY, dfStepY);
674 :
675 18 : if (nDataLineNum != nXSize * nYSize)
676 : {
677 : CPLError(CE_Failure, CPLE_AppDefined, "Found %d lines. Expected %d",
678 0 : nDataLineNum,nXSize * nYSize);
679 0 : VSIFCloseL(fp);
680 0 : return NULL;
681 : }
682 :
683 18 : if (poOpenInfo->eAccess == GA_Update)
684 : {
685 : CPLError( CE_Failure, CPLE_NotSupported,
686 : "The XYZ driver does not support update access to existing"
687 0 : " datasets.\n" );
688 0 : VSIFCloseL(fp);
689 0 : return NULL;
690 : }
691 :
692 : /* -------------------------------------------------------------------- */
693 : /* Create a corresponding GDALDataset. */
694 : /* -------------------------------------------------------------------- */
695 : XYZDataset *poDS;
696 :
697 18 : poDS = new XYZDataset();
698 18 : poDS->fp = fp;
699 18 : poDS->bHasHeaderLine = bHasHeaderLine;
700 18 : poDS->nXIndex = nXIndex;
701 18 : poDS->nYIndex = nYIndex;
702 18 : poDS->nZIndex = nZIndex;
703 18 : poDS->nMinTokens = nMinTokens;
704 18 : poDS->nRasterXSize = nXSize;
705 18 : poDS->nRasterYSize = nYSize;
706 18 : poDS->adfGeoTransform[0] = dfMinX - dfStepX / 2;
707 18 : poDS->adfGeoTransform[1] = dfStepX;
708 : poDS->adfGeoTransform[3] = (dfStepY < 0) ? dfMaxY - dfStepY / 2 :
709 35 : dfMinY - dfStepY / 2;
710 18 : poDS->adfGeoTransform[5] = dfStepY;
711 :
712 18 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
713 : {
714 0 : delete poDS;
715 0 : return NULL;
716 : }
717 :
718 : /* -------------------------------------------------------------------- */
719 : /* Create band information objects. */
720 : /* -------------------------------------------------------------------- */
721 18 : poDS->nBands = 1;
722 36 : for( i = 0; i < poDS->nBands; i++ )
723 18 : poDS->SetBand( i+1, new XYZRasterBand( poDS, i+1, eDT ) );
724 :
725 : /* -------------------------------------------------------------------- */
726 : /* Initialize any PAM information. */
727 : /* -------------------------------------------------------------------- */
728 18 : poDS->SetDescription( poOpenInfo->pszFilename );
729 18 : poDS->TryLoadXML();
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Support overviews. */
733 : /* -------------------------------------------------------------------- */
734 18 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
735 18 : return( poDS );
736 : }
737 :
738 : /************************************************************************/
739 : /* CreateCopy() */
740 : /************************************************************************/
741 :
742 : GDALDataset* XYZDataset::CreateCopy( const char * pszFilename,
743 : GDALDataset *poSrcDS,
744 : int bStrict, char ** papszOptions,
745 : GDALProgressFunc pfnProgress,
746 18 : void * pProgressData )
747 : {
748 : /* -------------------------------------------------------------------- */
749 : /* Some some rudimentary checks */
750 : /* -------------------------------------------------------------------- */
751 18 : int nBands = poSrcDS->GetRasterCount();
752 18 : if (nBands == 0)
753 : {
754 : CPLError( CE_Failure, CPLE_NotSupported,
755 1 : "XYZ driver does not support source dataset with zero band.\n");
756 1 : return NULL;
757 : }
758 :
759 17 : if (nBands != 1)
760 : {
761 : CPLError( (bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
762 4 : "XYZ driver only uses the first band of the dataset.\n");
763 4 : if (bStrict)
764 4 : return NULL;
765 : }
766 :
767 13 : if( pfnProgress && !pfnProgress( 0.0, NULL, pProgressData ) )
768 0 : return NULL;
769 :
770 : /* -------------------------------------------------------------------- */
771 : /* Get source dataset info */
772 : /* -------------------------------------------------------------------- */
773 :
774 13 : int nXSize = poSrcDS->GetRasterXSize();
775 13 : int nYSize = poSrcDS->GetRasterYSize();
776 : double adfGeoTransform[6];
777 13 : poSrcDS->GetGeoTransform(adfGeoTransform);
778 13 : if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
779 : {
780 : CPLError( CE_Failure, CPLE_NotSupported,
781 0 : "XYZ driver does not support CreateCopy() from skewed or rotated dataset.\n");
782 0 : return NULL;
783 : }
784 :
785 13 : GDALDataType eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
786 : GDALDataType eReqDT;
787 18 : if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 ||
788 : eSrcDT == GDT_UInt16 || eSrcDT == GDT_Int32)
789 5 : eReqDT = GDT_Int32;
790 : else
791 8 : eReqDT = GDT_Float32;
792 :
793 : /* -------------------------------------------------------------------- */
794 : /* Create target file */
795 : /* -------------------------------------------------------------------- */
796 :
797 13 : FILE* fp = VSIFOpenL(pszFilename, "wb");
798 13 : if (fp == NULL)
799 : {
800 : CPLError( CE_Failure, CPLE_AppDefined,
801 0 : "Cannot create %s", pszFilename );
802 0 : return NULL;
803 : }
804 :
805 : /* -------------------------------------------------------------------- */
806 : /* Read creation options */
807 : /* -------------------------------------------------------------------- */
808 : const char* pszColSep =
809 13 : CSLFetchNameValue(papszOptions, "COLUMN_SEPARATOR");
810 13 : if (pszColSep == NULL)
811 12 : pszColSep = " ";
812 1 : else if (EQUAL(pszColSep, "COMMA"))
813 0 : pszColSep = ",";
814 1 : else if (EQUAL(pszColSep, "SPACE"))
815 0 : pszColSep = " ";
816 1 : else if (EQUAL(pszColSep, "SEMICOLON"))
817 0 : pszColSep = ";";
818 1 : else if (EQUAL(pszColSep, "\\t") || EQUAL(pszColSep, "TAB"))
819 0 : pszColSep = "\t";
820 :
821 : const char* pszAddHeaderLine =
822 13 : CSLFetchNameValue(papszOptions, "ADD_HEADER_LINE");
823 13 : if (pszAddHeaderLine != NULL && CSLTestBoolean(pszAddHeaderLine))
824 : {
825 1 : VSIFPrintfL(fp, "X%sY%sZ\n", pszColSep, pszColSep);
826 : }
827 :
828 : /* -------------------------------------------------------------------- */
829 : /* Copy imagery */
830 : /* -------------------------------------------------------------------- */
831 13 : void* pLineBuffer = (void*) CPLMalloc(nXSize * sizeof(int));
832 : int i, j;
833 344 : for(j=0;j<nYSize;j++)
834 : {
835 : CPLErr eErr = poSrcDS->GetRasterBand(1)->RasterIO(
836 : GF_Read, 0, j, nXSize, 1,
837 : pLineBuffer, nXSize, 1,
838 331 : eReqDT, 0, 0);
839 331 : if (eErr != CE_None)
840 0 : break;
841 331 : double dfY = adfGeoTransform[3] + (j + 0.5) * adfGeoTransform[5];
842 42232 : for(i=0;i<nXSize;i++)
843 : {
844 41901 : double dfX = adfGeoTransform[0] + (i + 0.5) * adfGeoTransform[1];
845 41901 : VSIFPrintfL(fp, "%.18g%s%.18g%s", dfX, pszColSep, dfY, pszColSep);
846 41901 : if (eReqDT == GDT_Int32)
847 800 : VSIFPrintfL(fp, "%d\n", ((int*)pLineBuffer)[i]);
848 : else
849 41101 : VSIFPrintfL(fp, "%.18g\n", ((float*)pLineBuffer)[i]);
850 : }
851 331 : if (!pfnProgress( (j+1) * 1.0 / nYSize, NULL, pProgressData))
852 0 : break;
853 : }
854 13 : CPLFree(pLineBuffer);
855 13 : VSIFCloseL(fp);
856 :
857 13 : return (GDALDataset*) GDALOpen(pszFilename, GA_ReadOnly);
858 : }
859 :
860 : /************************************************************************/
861 : /* GetGeoTransform() */
862 : /************************************************************************/
863 :
864 1 : CPLErr XYZDataset::GetGeoTransform( double * padfTransform )
865 :
866 : {
867 1 : memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
868 :
869 1 : return( CE_None );
870 : }
871 :
872 : /************************************************************************/
873 : /* GDALRegister_XYZ() */
874 : /************************************************************************/
875 :
876 409 : void GDALRegister_XYZ()
877 :
878 : {
879 : GDALDriver *poDriver;
880 :
881 409 : if( GDALGetDriverByName( "XYZ" ) == NULL )
882 : {
883 392 : poDriver = new GDALDriver();
884 :
885 392 : poDriver->SetDescription( "XYZ" );
886 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
887 392 : "ASCII Gridded XYZ" );
888 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
889 392 : "frmt_xyz.html" );
890 392 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "xyz" );
891 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
892 : "<CreationOptionList>"
893 : " <Option name='COLUMN_SEPARATOR' type='string' default=' ' description='Separator between fields.'/>"
894 : " <Option name='ADD_HEADER_LINE' type='boolean' default='false' description='Add an header line with column names.'/>"
895 392 : "</CreationOptionList>");
896 :
897 392 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
898 :
899 392 : poDriver->pfnOpen = XYZDataset::Open;
900 392 : poDriver->pfnIdentify = XYZDataset::Identify;
901 392 : poDriver->pfnCreateCopy = XYZDataset::CreateCopy;
902 :
903 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
904 : }
905 409 : }
906 :
|
