1 : /******************************************************************************
2 : * $Id: xyzdataset.cpp 22596 2011-06-27 18:53:18Z 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 22596 2011-06-27 18:53:18Z 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 : VSILFILE *fp;
53 : int bHasHeaderLine;
54 : int nCommentLineCount;
55 : int nXIndex;
56 : int nYIndex;
57 : int nZIndex;
58 : int nMinTokens;
59 : int nLineNum; /* any line */
60 : int nDataLineNum; /* line with values (header line and empty lines ignored) */
61 : double adfGeoTransform[6];
62 :
63 : static int IdentifyEx( GDALOpenInfo *, int&, int& nCommentLineCount );
64 :
65 : public:
66 : XYZDataset();
67 : virtual ~XYZDataset();
68 :
69 : virtual CPLErr GetGeoTransform( double * );
70 :
71 : static GDALDataset *Open( GDALOpenInfo * );
72 : static int Identify( GDALOpenInfo * );
73 : static GDALDataset *CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
74 : int bStrict, char ** papszOptions,
75 : GDALProgressFunc pfnProgress, void * pProgressData );
76 : };
77 :
78 : /************************************************************************/
79 : /* ==================================================================== */
80 : /* XYZRasterBand */
81 : /* ==================================================================== */
82 : /************************************************************************/
83 :
84 : class XYZRasterBand : public GDALPamRasterBand
85 36 : {
86 : friend class XYZDataset;
87 :
88 : public:
89 :
90 : XYZRasterBand( XYZDataset *, int, GDALDataType );
91 :
92 : virtual CPLErr IReadBlock( int, int, void * );
93 : };
94 :
95 :
96 : /************************************************************************/
97 : /* XYZRasterBand() */
98 : /************************************************************************/
99 :
100 36 : XYZRasterBand::XYZRasterBand( XYZDataset *poDS, int nBand, GDALDataType eDT )
101 :
102 : {
103 36 : this->poDS = poDS;
104 36 : this->nBand = nBand;
105 :
106 36 : eDataType = eDT;
107 :
108 36 : nBlockXSize = poDS->GetRasterXSize();
109 36 : nBlockYSize = 1;
110 36 : }
111 :
112 : /************************************************************************/
113 : /* IReadBlock() */
114 : /************************************************************************/
115 :
116 890 : CPLErr XYZRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
117 : void * pImage )
118 :
119 : {
120 890 : XYZDataset *poGDS = (XYZDataset *) poDS;
121 :
122 890 : if (poGDS->fp == NULL)
123 0 : return CE_Failure;
124 :
125 890 : int nLineInFile = nBlockYOff * nBlockXSize;
126 890 : if (poGDS->nDataLineNum > nLineInFile)
127 : {
128 14 : poGDS->nDataLineNum = 0;
129 14 : VSIFSeekL(poGDS->fp, 0, SEEK_SET);
130 :
131 14 : for(int i=0;i<poGDS->nCommentLineCount;i++)
132 0 : CPLReadLine2L(poGDS->fp, 100, NULL);
133 :
134 14 : if (poGDS->bHasHeaderLine)
135 : {
136 10 : const char* pszLine = CPLReadLine2L(poGDS->fp, 100, 0);
137 10 : if (pszLine == NULL)
138 : {
139 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
140 0 : return CE_Failure;
141 : }
142 10 : poGDS->nLineNum ++;
143 : }
144 : }
145 :
146 1804 : while(poGDS->nDataLineNum < nLineInFile)
147 : {
148 24 : const char* pszLine = CPLReadLine2L(poGDS->fp, 100, 0);
149 24 : if (pszLine == NULL)
150 : {
151 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
152 0 : return CE_Failure;
153 : }
154 24 : poGDS->nLineNum ++;
155 :
156 24 : const char* pszPtr = pszLine;
157 : char ch;
158 24 : int nCol = 0;
159 24 : int bLastWasSep = TRUE;
160 120 : while((ch = *pszPtr) != '\0')
161 : {
162 96 : if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
163 : {
164 24 : if (!bLastWasSep)
165 24 : nCol ++;
166 24 : bLastWasSep = TRUE;
167 : }
168 : else
169 : {
170 48 : bLastWasSep = FALSE;
171 : }
172 72 : pszPtr ++;
173 : }
174 :
175 : /* Skip empty line */
176 24 : if (nCol == 0 && bLastWasSep)
177 12 : continue;
178 :
179 12 : poGDS->nDataLineNum ++;
180 : }
181 :
182 : int i;
183 164106 : for(i=0;i<nBlockXSize;i++)
184 : {
185 : int nCol;
186 : int bLastWasSep;
187 163224 : do
188 : {
189 163224 : const char* pszLine = CPLReadLine2L(poGDS->fp, 100, 0);
190 163224 : if (pszLine == NULL)
191 : {
192 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
193 0 : return CE_Failure;
194 : }
195 163224 : poGDS->nLineNum ++;
196 :
197 163224 : const char* pszPtr = pszLine;
198 : char ch;
199 163224 : nCol = 0;
200 163224 : bLastWasSep = TRUE;
201 5598096 : while((ch = *pszPtr) != '\0')
202 : {
203 5598080 : if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
204 : {
205 326432 : if (!bLastWasSep)
206 326432 : nCol ++;
207 326432 : bLastWasSep = TRUE;
208 : }
209 : else
210 : {
211 4945216 : if (bLastWasSep && nCol == poGDS->nZIndex)
212 : {
213 163216 : double dfZ = CPLAtofM(pszPtr);
214 163216 : if (eDataType == GDT_Float32)
215 : {
216 161604 : ((float*)pImage)[i] = (float)dfZ;
217 : }
218 1612 : else if (eDataType == GDT_Int32)
219 : {
220 0 : ((GInt32*)pImage)[i] = (GInt32)dfZ;
221 : }
222 1612 : else if (eDataType == GDT_Int16)
223 : {
224 0 : ((GInt16*)pImage)[i] = (GInt16)dfZ;
225 : }
226 : else
227 : {
228 1612 : ((GByte*)pImage)[i] = (GByte)dfZ;
229 : }
230 : }
231 4945216 : bLastWasSep = FALSE;
232 : }
233 5271648 : pszPtr ++;
234 : }
235 :
236 : /* Skip empty line */
237 : }
238 : while (nCol == 0 && bLastWasSep);
239 :
240 163216 : poGDS->nDataLineNum ++;
241 163216 : nCol ++;
242 163216 : if (nCol < poGDS->nMinTokens)
243 : {
244 0 : memset(pImage, 0, nBlockXSize * (GDALGetDataTypeSize(eDataType) / 8));
245 0 : return CE_Failure;
246 : }
247 : }
248 890 : CPLAssert(poGDS->nDataLineNum == (nBlockYOff + 1) * nBlockXSize);
249 :
250 890 : return CE_None;
251 : }
252 :
253 : /************************************************************************/
254 : /* ~XYZDataset() */
255 : /************************************************************************/
256 :
257 36 : XYZDataset::XYZDataset()
258 : {
259 36 : fp = NULL;
260 36 : nDataLineNum = INT_MAX;
261 36 : nLineNum = 0;
262 36 : bHasHeaderLine = FALSE;
263 36 : nXIndex = -1;
264 36 : nYIndex = -1;
265 36 : nZIndex = -1;
266 36 : nMinTokens = 0;
267 36 : adfGeoTransform[0] = 0;
268 36 : adfGeoTransform[1] = 1;
269 36 : adfGeoTransform[2] = 0;
270 36 : adfGeoTransform[3] = 0;
271 36 : adfGeoTransform[4] = 0;
272 36 : adfGeoTransform[5] = 1;
273 36 : }
274 :
275 : /************************************************************************/
276 : /* ~XYZDataset() */
277 : /************************************************************************/
278 :
279 36 : XYZDataset::~XYZDataset()
280 :
281 : {
282 36 : FlushCache();
283 36 : if (fp)
284 36 : VSIFCloseL(fp);
285 36 : }
286 :
287 : /************************************************************************/
288 : /* Identify() */
289 : /************************************************************************/
290 :
291 18894 : int XYZDataset::Identify( GDALOpenInfo * poOpenInfo )
292 : {
293 : int bHasHeaderLine, nCommentLineCount;
294 18894 : return IdentifyEx(poOpenInfo, bHasHeaderLine, nCommentLineCount);
295 : }
296 :
297 : /************************************************************************/
298 : /* IdentifyEx() */
299 : /************************************************************************/
300 :
301 :
302 21904 : int XYZDataset::IdentifyEx( GDALOpenInfo * poOpenInfo,
303 : int& bHasHeaderLine,
304 : int& nCommentLineCount)
305 :
306 : {
307 : int i;
308 :
309 21904 : bHasHeaderLine = FALSE;
310 21904 : nCommentLineCount = 0;
311 :
312 21904 : CPLString osFilename(poOpenInfo->pszFilename);
313 :
314 21904 : GDALOpenInfo* poOpenInfoToDelete = NULL;
315 : /* GZipped .xyz files are common, so automagically open them */
316 : /* if the /vsigzip/ has not been explicitely passed */
317 21904 : if (strlen(poOpenInfo->pszFilename) > 6 &&
318 : EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "xyz.gz") &&
319 : !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
320 : {
321 0 : osFilename = "/vsigzip/";
322 0 : osFilename += poOpenInfo->pszFilename;
323 : poOpenInfo = poOpenInfoToDelete =
324 : new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
325 0 : poOpenInfo->papszSiblingFiles);
326 : }
327 :
328 21904 : if (poOpenInfo->nHeaderBytes == 0)
329 : {
330 21140 : delete poOpenInfoToDelete;
331 21140 : return FALSE;
332 : }
333 :
334 : /* -------------------------------------------------------------------- */
335 : /* Chech that it looks roughly as a XYZ dataset */
336 : /* -------------------------------------------------------------------- */
337 764 : const char* pszData = (const char*)poOpenInfo->pabyHeader;
338 :
339 : /* Skip comments line at the beginning such as in */
340 : /* http://pubs.usgs.gov/of/2003/ofr-03-230/DATA/NSLCU.XYZ */
341 764 : i=0;
342 764 : if (pszData[i] == '/')
343 : {
344 0 : nCommentLineCount ++;
345 :
346 0 : i++;
347 0 : for(;i<poOpenInfo->nHeaderBytes;i++)
348 : {
349 0 : char ch = pszData[i];
350 0 : if (ch == 13 || ch == 10)
351 : {
352 0 : if (ch == 13 && pszData[i+1] == 10)
353 0 : i++;
354 0 : if (pszData[i+1] == '/')
355 : {
356 0 : nCommentLineCount ++;
357 0 : i++;
358 : }
359 : else
360 0 : break;
361 : }
362 : }
363 : }
364 :
365 29022 : for(;i<poOpenInfo->nHeaderBytes;i++)
366 : {
367 28994 : char ch = pszData[i];
368 28994 : if (ch == 13 || ch == 10)
369 : {
370 38 : break;
371 : }
372 28956 : else if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
373 : ;
374 3180 : else if ((ch >= '0' && ch <= '9') || ch == '.' || ch == '+' ||
375 : ch == '-' || ch == 'e' || ch == 'E')
376 : ;
377 1706 : else if (ch == '"' || (ch >= 'a' && ch <= 'z') ||
378 : (ch >= 'A' && ch <= 'Z'))
379 504 : bHasHeaderLine = TRUE;
380 : else
381 : {
382 698 : delete poOpenInfoToDelete;
383 698 : return FALSE;
384 : }
385 : }
386 66 : int bHasFoundNewLine = FALSE;
387 66 : int bPrevWasSep = TRUE;
388 66 : int nCols = 0;
389 66 : int nMaxCols = 0;
390 33902 : for(;i<poOpenInfo->nHeaderBytes;i++)
391 : {
392 33838 : char ch = pszData[i];
393 35080 : if (ch == 13 || ch == 10)
394 : {
395 1242 : bHasFoundNewLine = TRUE;
396 1242 : if (!bPrevWasSep)
397 : {
398 1190 : nCols ++;
399 1190 : if (nCols > nMaxCols)
400 36 : nMaxCols = nCols;
401 : }
402 1242 : bPrevWasSep = TRUE;
403 1242 : nCols = 0;
404 : }
405 35004 : else if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
406 : {
407 2408 : if (!bPrevWasSep)
408 : {
409 2408 : nCols ++;
410 2408 : if (nCols > nMaxCols)
411 72 : nMaxCols = nCols;
412 : }
413 2408 : bPrevWasSep = TRUE;
414 : }
415 60374 : else if ((ch >= '0' && ch <= '9') || ch == '.' || ch == '+' ||
416 : ch == '-' || ch == 'e' || ch == 'E')
417 : {
418 30186 : bPrevWasSep = FALSE;
419 : }
420 : else
421 : {
422 2 : delete poOpenInfoToDelete;
423 2 : return FALSE;
424 : }
425 : }
426 :
427 64 : delete poOpenInfoToDelete;
428 64 : return bHasFoundNewLine && nMaxCols >= 3;
429 : }
430 :
431 : /************************************************************************/
432 : /* Open() */
433 : /************************************************************************/
434 :
435 3010 : GDALDataset *XYZDataset::Open( GDALOpenInfo * poOpenInfo )
436 :
437 : {
438 : int i;
439 : int bHasHeaderLine;
440 3010 : int nCommentLineCount = 0;
441 :
442 3010 : if (!IdentifyEx(poOpenInfo, bHasHeaderLine, nCommentLineCount))
443 2974 : return NULL;
444 :
445 36 : CPLString osFilename(poOpenInfo->pszFilename);
446 :
447 : /* GZipped .xyz files are common, so automagically open them */
448 : /* if the /vsigzip/ has not been explicitely passed */
449 36 : if (strlen(poOpenInfo->pszFilename) > 6 &&
450 : EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "xyz.gz") &&
451 : !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
452 : {
453 0 : osFilename = "/vsigzip/";
454 0 : osFilename += poOpenInfo->pszFilename;
455 : }
456 :
457 : /* -------------------------------------------------------------------- */
458 : /* Find dataset characteristics */
459 : /* -------------------------------------------------------------------- */
460 36 : VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "rb");
461 36 : if (fp == NULL)
462 0 : return NULL;
463 :
464 : /* For better performance of CPLReadLine2L() we create a buffered reader */
465 : /* (except for /vsigzip/ since it has one internally) */
466 36 : if (!EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
467 36 : fp = (VSILFILE*) VSICreateBufferedReaderHandle((VSIVirtualHandle*)fp);
468 :
469 : const char* pszLine;
470 36 : int nXIndex = -1, nYIndex = -1, nZIndex = -1;
471 36 : int nMinTokens = 0;
472 :
473 :
474 36 : for(i=0;i<nCommentLineCount;i++)
475 0 : CPLReadLine2L(fp, 100, NULL);
476 :
477 : /* -------------------------------------------------------------------- */
478 : /* Parse header line */
479 : /* -------------------------------------------------------------------- */
480 36 : if (bHasHeaderLine)
481 : {
482 8 : pszLine = CPLReadLine2L(fp, 100, NULL);
483 8 : if (pszLine == NULL)
484 : {
485 0 : VSIFCloseL(fp);
486 0 : return NULL;
487 : }
488 : char** papszTokens = CSLTokenizeString2( pszLine, " ,\t;",
489 8 : CSLT_HONOURSTRINGS );
490 8 : int nTokens = CSLCount(papszTokens);
491 8 : if (nTokens < 3)
492 : {
493 : CPLError(CE_Failure, CPLE_AppDefined,
494 : "At line %d, found %d tokens. Expected 3 at least",
495 0 : 1, nTokens);
496 0 : CSLDestroy(papszTokens);
497 0 : VSIFCloseL(fp);
498 0 : return NULL;
499 : }
500 : int i;
501 32 : for(i=0;i<nTokens;i++)
502 : {
503 64 : if (EQUAL(papszTokens[i], "x") ||
504 16 : EQUALN(papszTokens[i], "lon", 3) ||
505 16 : EQUALN(papszTokens[i], "east", 4))
506 8 : nXIndex = i;
507 40 : else if (EQUAL(papszTokens[i], "y") ||
508 8 : EQUALN(papszTokens[i], "lat", 3) ||
509 8 : EQUALN(papszTokens[i], "north", 5))
510 8 : nYIndex = i;
511 8 : else if (EQUAL(papszTokens[i], "z") ||
512 0 : EQUALN(papszTokens[i], "alt", 3) ||
513 0 : EQUAL(papszTokens[i], "height"))
514 8 : nZIndex = i;
515 : }
516 8 : CSLDestroy(papszTokens);
517 8 : papszTokens = NULL;
518 8 : if (nXIndex < 0 || nYIndex < 0 || nZIndex < 0)
519 : {
520 : CPLError(CE_Warning, CPLE_AppDefined,
521 0 : "Could not find one of the X, Y or Z column names in header line. Defaulting to the first 3 columns");
522 0 : nXIndex = 0;
523 0 : nYIndex = 1;
524 0 : nZIndex = 2;
525 : }
526 8 : nMinTokens = 1 + MAX(MAX(nXIndex, nYIndex), nZIndex);
527 : }
528 : else
529 : {
530 28 : nXIndex = 0;
531 28 : nYIndex = 1;
532 28 : nZIndex = 2;
533 28 : nMinTokens = 3;
534 : }
535 :
536 : /* -------------------------------------------------------------------- */
537 : /* Parse data lines */
538 : /* -------------------------------------------------------------------- */
539 :
540 36 : int nXSize = 0, nYSize = 0;
541 36 : int nLineNum = 0;
542 36 : int nDataLineNum = 0;
543 36 : double dfFirstX = 0;
544 36 : double dfX = 0, dfY = 0, dfZ = 0;
545 36 : double dfMinX = 0, dfMinY = 0, dfMaxX = 0, dfMaxY = 0;
546 36 : double dfLastX = 0, dfLastY = 0;
547 36 : double dfStepX = 0, dfStepY = 0;
548 36 : GDALDataType eDT = GDT_Byte;
549 247102 : while((pszLine = CPLReadLine2L(fp, 100, NULL)) != NULL)
550 : {
551 247030 : nLineNum ++;
552 :
553 247030 : const char* pszPtr = pszLine;
554 : char ch;
555 247030 : int nCol = 0;
556 247030 : int bLastWasSep = TRUE;
557 8483336 : while((ch = *pszPtr) != '\0')
558 : {
559 8483312 : if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
560 : {
561 494036 : if (!bLastWasSep)
562 494036 : nCol ++;
563 494036 : bLastWasSep = TRUE;
564 : }
565 : else
566 : {
567 7495240 : if (bLastWasSep)
568 : {
569 741054 : if (nCol == nXIndex)
570 247018 : dfX = CPLAtofM(pszPtr);
571 494036 : else if (nCol == nYIndex)
572 247018 : dfY = CPLAtofM(pszPtr);
573 247018 : else if (nCol == nZIndex && eDT != GDT_Float32)
574 : {
575 30262 : dfZ = CPLAtofM(pszPtr);
576 30262 : int nZ = (int)dfZ;
577 30262 : if ((double)nZ != dfZ)
578 : {
579 6 : eDT = GDT_Float32;
580 : }
581 30256 : else if ((eDT == GDT_Byte || eDT == GDT_Int16) && (nZ < 0 || nZ > 255))
582 : {
583 0 : if (nZ < -32768 || nZ > 32767)
584 0 : eDT = GDT_Int32;
585 : else
586 0 : eDT = GDT_Int16;
587 : }
588 : }
589 : }
590 7495240 : bLastWasSep = FALSE;
591 : }
592 7989276 : pszPtr ++;
593 : }
594 : /* skip empty lines */
595 247030 : if (bLastWasSep && nCol == 0)
596 : {
597 12 : continue;
598 : }
599 247018 : nDataLineNum ++;
600 247018 : nCol ++;
601 247018 : if (nCol < nMinTokens)
602 : {
603 : CPLError(CE_Failure, CPLE_AppDefined,
604 : "At line %d, found %d tokens. Expected %d at least",
605 0 : nLineNum, nCol, nMinTokens);
606 0 : VSIFCloseL(fp);
607 0 : return NULL;
608 : }
609 :
610 247018 : if (nDataLineNum == 1)
611 : {
612 36 : dfFirstX = dfMinX = dfMaxX = dfX;
613 36 : dfMinY = dfMaxY = dfY;
614 : }
615 : else
616 : {
617 246982 : if (dfX < dfMinX) dfMinX = dfX;
618 246982 : if (dfX > dfMaxX) dfMaxX = dfX;
619 246982 : if (dfY < dfMinY) dfMinY = dfY;
620 246982 : if (dfY > dfMaxY) dfMaxY = dfY;
621 : }
622 247018 : if (nDataLineNum == 2)
623 : {
624 36 : dfStepX = dfX - dfLastX;
625 36 : if (dfStepX <= 0)
626 : {
627 : CPLError(CE_Failure, CPLE_AppDefined,
628 : "Ungridded dataset: At line %d, X spacing was %f. Expected >0 value",
629 0 : nLineNum, dfStepX);
630 0 : VSIFCloseL(fp);
631 0 : return NULL;
632 : }
633 : }
634 246982 : else if (nDataLineNum > 2)
635 : {
636 246946 : double dfNewStepX = dfX - dfLastX;
637 246946 : double dfNewStepY = dfY - dfLastY;
638 246946 : if (dfNewStepY != 0)
639 : {
640 1516 : nYSize ++;
641 1516 : if (dfStepY == 0)
642 : {
643 36 : nXSize = nDataLineNum - 1;
644 36 : double dfAdjustedStepX = (dfMaxX - dfMinX) / (nXSize - 1);
645 36 : if (fabs(dfStepX - dfAdjustedStepX) > 1e-8)
646 : {
647 0 : CPLDebug("XYZ", "Adjusting stepx from %f to %f", dfStepX, dfAdjustedStepX);
648 : }
649 36 : dfStepX = dfAdjustedStepX;
650 : }
651 1516 : if (dfStepY != 0 && fabs(dfX - dfFirstX) > 1e-8)
652 : {
653 : CPLError(CE_Failure, CPLE_AppDefined,
654 : "Ungridded dataset: At line %d, X is %f, where as %f was expected",
655 0 : nLineNum, dfX, dfFirstX);
656 0 : VSIFCloseL(fp);
657 0 : return NULL;
658 : }
659 1516 : if (dfStepY != 0 && fabs(dfLastX - dfMaxX) > 1e-8)
660 : {
661 : CPLError(CE_Failure, CPLE_AppDefined,
662 : "Ungridded dataset: At line %d, X is %f, where as %f was expected",
663 0 : nLineNum - 1, dfLastX, dfMaxX);
664 0 : VSIFCloseL(fp);
665 0 : return NULL;
666 : }
667 : /*if (dfStepY != 0 && fabs(dfNewStepY - dfStepY) > 1e-8)
668 : {
669 : CPLError(CE_Failure, CPLE_AppDefined,
670 : "Ungridded dataset: At line %d, Y spacing was %f, whereas it was %f before",
671 : nLineNum, dfNewStepY, dfStepY);
672 : VSIFCloseL(fp);
673 : return NULL;
674 : }*/
675 1516 : dfStepY = dfNewStepY;
676 : }
677 : else if (dfNewStepX != 0)
678 : {
679 : /*if (dfStepX != 0 && fabs(dfNewStepX - dfStepX) > 1e-8)
680 : {
681 : CPLError(CE_Failure, CPLE_AppDefined,
682 : "At line %d, X spacing was %f, whereas it was %f before",
683 : nLineNum, dfNewStepX, dfStepX);
684 : VSIFCloseL(fp);
685 : return NULL;
686 : }*/
687 : }
688 : }
689 247018 : dfLastX = dfX;
690 247018 : dfLastY = dfY;
691 : }
692 36 : nYSize ++;
693 :
694 36 : if (dfStepX == 0)
695 : {
696 0 : CPLError(CE_Failure, CPLE_AppDefined, "Couldn't determine X spacing");
697 0 : VSIFCloseL(fp);
698 0 : return NULL;
699 : }
700 :
701 36 : if (dfStepY == 0)
702 : {
703 0 : CPLError(CE_Failure, CPLE_AppDefined, "Couldn't determine Y spacing");
704 0 : VSIFCloseL(fp);
705 0 : return NULL;
706 : }
707 :
708 36 : double dfAdjustedStepY = ((dfStepY < 0) ? -1 : 1) * (dfMaxY - dfMinY) / (nYSize - 1);
709 36 : if (fabs(dfStepY - dfAdjustedStepY) > 1e-8)
710 : {
711 0 : CPLDebug("XYZ", "Adjusting stepy from %f to %f", dfStepY, dfAdjustedStepY);
712 : }
713 36 : dfStepY = dfAdjustedStepY;
714 :
715 : //CPLDebug("XYZ", "minx=%f maxx=%f stepx=%f", dfMinX, dfMaxX, dfStepX);
716 : //CPLDebug("XYZ", "miny=%f maxy=%f stepy=%f", dfMinY, dfMaxY, dfStepY);
717 :
718 36 : if (nDataLineNum != nXSize * nYSize)
719 : {
720 : CPLError(CE_Failure, CPLE_AppDefined, "Found %d lines. Expected %d",
721 0 : nDataLineNum,nXSize * nYSize);
722 0 : VSIFCloseL(fp);
723 0 : return NULL;
724 : }
725 :
726 36 : if (poOpenInfo->eAccess == GA_Update)
727 : {
728 : CPLError( CE_Failure, CPLE_NotSupported,
729 : "The XYZ driver does not support update access to existing"
730 0 : " datasets.\n" );
731 0 : VSIFCloseL(fp);
732 0 : return NULL;
733 : }
734 :
735 : /* -------------------------------------------------------------------- */
736 : /* Create a corresponding GDALDataset. */
737 : /* -------------------------------------------------------------------- */
738 : XYZDataset *poDS;
739 :
740 36 : poDS = new XYZDataset();
741 36 : poDS->fp = fp;
742 36 : poDS->bHasHeaderLine = bHasHeaderLine;
743 36 : poDS->nCommentLineCount = nCommentLineCount;
744 36 : poDS->nXIndex = nXIndex;
745 36 : poDS->nYIndex = nYIndex;
746 36 : poDS->nZIndex = nZIndex;
747 36 : poDS->nMinTokens = nMinTokens;
748 36 : poDS->nRasterXSize = nXSize;
749 36 : poDS->nRasterYSize = nYSize;
750 36 : poDS->adfGeoTransform[0] = dfMinX - dfStepX / 2;
751 36 : poDS->adfGeoTransform[1] = dfStepX;
752 : poDS->adfGeoTransform[3] = (dfStepY < 0) ? dfMaxY - dfStepY / 2 :
753 70 : dfMinY - dfStepY / 2;
754 36 : poDS->adfGeoTransform[5] = dfStepY;
755 :
756 36 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
757 : {
758 0 : delete poDS;
759 0 : return NULL;
760 : }
761 :
762 : /* -------------------------------------------------------------------- */
763 : /* Create band information objects. */
764 : /* -------------------------------------------------------------------- */
765 36 : poDS->nBands = 1;
766 72 : for( i = 0; i < poDS->nBands; i++ )
767 36 : poDS->SetBand( i+1, new XYZRasterBand( poDS, i+1, eDT ) );
768 :
769 : /* -------------------------------------------------------------------- */
770 : /* Initialize any PAM information. */
771 : /* -------------------------------------------------------------------- */
772 36 : poDS->SetDescription( poOpenInfo->pszFilename );
773 36 : poDS->TryLoadXML();
774 :
775 : /* -------------------------------------------------------------------- */
776 : /* Support overviews. */
777 : /* -------------------------------------------------------------------- */
778 36 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
779 36 : return( poDS );
780 : }
781 :
782 : /************************************************************************/
783 : /* CreateCopy() */
784 : /************************************************************************/
785 :
786 40 : GDALDataset* XYZDataset::CreateCopy( const char * pszFilename,
787 : GDALDataset *poSrcDS,
788 : int bStrict, char ** papszOptions,
789 : GDALProgressFunc pfnProgress,
790 : void * pProgressData )
791 : {
792 : /* -------------------------------------------------------------------- */
793 : /* Some some rudimentary checks */
794 : /* -------------------------------------------------------------------- */
795 40 : int nBands = poSrcDS->GetRasterCount();
796 40 : if (nBands == 0)
797 : {
798 : CPLError( CE_Failure, CPLE_NotSupported,
799 2 : "XYZ driver does not support source dataset with zero band.\n");
800 2 : return NULL;
801 : }
802 :
803 38 : if (nBands != 1)
804 : {
805 : CPLError( (bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
806 8 : "XYZ driver only uses the first band of the dataset.\n");
807 8 : if (bStrict)
808 8 : return NULL;
809 : }
810 :
811 30 : if( pfnProgress && !pfnProgress( 0.0, NULL, pProgressData ) )
812 0 : return NULL;
813 :
814 : /* -------------------------------------------------------------------- */
815 : /* Get source dataset info */
816 : /* -------------------------------------------------------------------- */
817 :
818 30 : int nXSize = poSrcDS->GetRasterXSize();
819 30 : int nYSize = poSrcDS->GetRasterYSize();
820 : double adfGeoTransform[6];
821 30 : poSrcDS->GetGeoTransform(adfGeoTransform);
822 30 : if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
823 : {
824 : CPLError( CE_Failure, CPLE_NotSupported,
825 0 : "XYZ driver does not support CreateCopy() from skewed or rotated dataset.\n");
826 0 : return NULL;
827 : }
828 :
829 30 : GDALDataType eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
830 : GDALDataType eReqDT;
831 44 : if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 ||
832 : eSrcDT == GDT_UInt16 || eSrcDT == GDT_Int32)
833 14 : eReqDT = GDT_Int32;
834 : else
835 16 : eReqDT = GDT_Float32;
836 :
837 : /* -------------------------------------------------------------------- */
838 : /* Create target file */
839 : /* -------------------------------------------------------------------- */
840 :
841 30 : VSILFILE* fp = VSIFOpenL(pszFilename, "wb");
842 30 : if (fp == NULL)
843 : {
844 : CPLError( CE_Failure, CPLE_AppDefined,
845 4 : "Cannot create %s", pszFilename );
846 4 : return NULL;
847 : }
848 :
849 : /* -------------------------------------------------------------------- */
850 : /* Read creation options */
851 : /* -------------------------------------------------------------------- */
852 : const char* pszColSep =
853 26 : CSLFetchNameValue(papszOptions, "COLUMN_SEPARATOR");
854 26 : if (pszColSep == NULL)
855 24 : pszColSep = " ";
856 2 : else if (EQUAL(pszColSep, "COMMA"))
857 0 : pszColSep = ",";
858 2 : else if (EQUAL(pszColSep, "SPACE"))
859 0 : pszColSep = " ";
860 2 : else if (EQUAL(pszColSep, "SEMICOLON"))
861 0 : pszColSep = ";";
862 2 : else if (EQUAL(pszColSep, "\\t") || EQUAL(pszColSep, "TAB"))
863 0 : pszColSep = "\t";
864 :
865 : const char* pszAddHeaderLine =
866 26 : CSLFetchNameValue(papszOptions, "ADD_HEADER_LINE");
867 26 : if (pszAddHeaderLine != NULL && CSLTestBoolean(pszAddHeaderLine))
868 : {
869 2 : VSIFPrintfL(fp, "X%sY%sZ\n", pszColSep, pszColSep);
870 : }
871 :
872 : /* -------------------------------------------------------------------- */
873 : /* Copy imagery */
874 : /* -------------------------------------------------------------------- */
875 26 : void* pLineBuffer = (void*) CPLMalloc(nXSize * sizeof(int));
876 : int i, j;
877 26 : CPLErr eErr = CE_None;
878 26 : for(j=0;j<nYSize && eErr == CE_None;j++)
879 : {
880 : eErr = poSrcDS->GetRasterBand(1)->RasterIO(
881 : GF_Read, 0, j, nXSize, 1,
882 : pLineBuffer, nXSize, 1,
883 662 : eReqDT, 0, 0);
884 662 : if (eErr != CE_None)
885 0 : break;
886 662 : double dfY = adfGeoTransform[3] + (j + 0.5) * adfGeoTransform[5];
887 662 : CPLString osBuf;
888 84464 : for(i=0;i<nXSize;i++)
889 : {
890 : char szBuf[256];
891 83802 : double dfX = adfGeoTransform[0] + (i + 0.5) * adfGeoTransform[1];
892 83802 : if (eReqDT == GDT_Int32)
893 1600 : sprintf(szBuf, "%.18g%c%.18g%c%d\n", dfX, pszColSep[0], dfY, pszColSep[0], ((int*)pLineBuffer)[i]);
894 : else
895 82202 : sprintf(szBuf, "%.18g%c%.18g%c%.18g\n", dfX, pszColSep[0], dfY, pszColSep[0], ((float*)pLineBuffer)[i]);
896 83802 : osBuf += szBuf;
897 83802 : if( (i & 1023) == 0 || i == nXSize - 1 )
898 : {
899 1324 : if ( VSIFWriteL( osBuf, (int)osBuf.size(), 1, fp ) != 1 )
900 : {
901 0 : eErr = CE_Failure;
902 : CPLError( CE_Failure, CPLE_AppDefined,
903 0 : "Write failed, disk full?\n" );
904 0 : break;
905 : }
906 1324 : osBuf = "";
907 : }
908 : }
909 662 : if (!pfnProgress( (j+1) * 1.0 / nYSize, NULL, pProgressData))
910 : break;
911 : }
912 26 : CPLFree(pLineBuffer);
913 26 : VSIFCloseL(fp);
914 :
915 26 : if (eErr != CE_None)
916 0 : return NULL;
917 :
918 : /* If outputing to stdout, we can't reopen it, so we'll return */
919 : /* a fake dataset to make the caller happy */
920 26 : CPLPushErrorHandler(CPLQuietErrorHandler);
921 26 : GDALDataset* poDS = (GDALDataset*) GDALOpen(pszFilename, GA_ReadOnly);
922 26 : CPLPopErrorHandler();
923 26 : if (poDS)
924 26 : return poDS;
925 :
926 0 : CPLErrorReset();
927 :
928 0 : XYZDataset* poXYZ_DS = new XYZDataset();
929 0 : poXYZ_DS->nRasterXSize = nXSize;
930 0 : poXYZ_DS->nRasterYSize = nYSize;
931 0 : poXYZ_DS->nBands = 1;
932 0 : poXYZ_DS->SetBand( 1, new XYZRasterBand( poXYZ_DS, 1, eReqDT ) );
933 0 : return poXYZ_DS;
934 : }
935 :
936 : /************************************************************************/
937 : /* GetGeoTransform() */
938 : /************************************************************************/
939 :
940 2 : CPLErr XYZDataset::GetGeoTransform( double * padfTransform )
941 :
942 : {
943 2 : memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
944 :
945 2 : return( CE_None );
946 : }
947 :
948 : /************************************************************************/
949 : /* GDALRegister_XYZ() */
950 : /************************************************************************/
951 :
952 1135 : void GDALRegister_XYZ()
953 :
954 : {
955 : GDALDriver *poDriver;
956 :
957 1135 : if( GDALGetDriverByName( "XYZ" ) == NULL )
958 : {
959 1093 : poDriver = new GDALDriver();
960 :
961 1093 : poDriver->SetDescription( "XYZ" );
962 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
963 1093 : "ASCII Gridded XYZ" );
964 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
965 1093 : "frmt_xyz.html" );
966 1093 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "xyz" );
967 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
968 : "<CreationOptionList>"
969 : " <Option name='COLUMN_SEPARATOR' type='string' default=' ' description='Separator between fields.'/>"
970 : " <Option name='ADD_HEADER_LINE' type='boolean' default='false' description='Add an header line with column names.'/>"
971 1093 : "</CreationOptionList>");
972 :
973 1093 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
974 :
975 1093 : poDriver->pfnOpen = XYZDataset::Open;
976 1093 : poDriver->pfnIdentify = XYZDataset::Identify;
977 1093 : poDriver->pfnCreateCopy = XYZDataset::CreateCopy;
978 :
979 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
980 : }
981 1135 : }
982 :
|
