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