1 : /******************************************************************************
2 : * $Id: sagadataset.cpp 22699 2011-07-11 18:12:32Z rouault $
3 : * Project: SAGA GIS Binary Driver
4 : * Purpose: Implements the SAGA GIS Binary Grid Format.
5 : * Author: Volker Wichmann, wichmann@laserdata.at
6 : * (Based on gsbgdataset.cpp by Kevin Locke and Frank Warmerdam)
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2009, Volker Wichmann <wichmann@laserdata.at>
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_conv.h"
31 :
32 : #include <float.h>
33 : #include <limits.h>
34 : #include <assert.h>
35 :
36 : #include "gdal_pam.h"
37 :
38 : CPL_CVSID("$Id: sagadataset.cpp 22699 2011-07-11 18:12:32Z rouault $");
39 :
40 : #ifndef INT_MAX
41 : # define INT_MAX 2147483647
42 : #endif /* INT_MAX */
43 :
44 : /* NODATA Values */
45 : //#define SG_NODATA_GDT_Bit 0.0
46 : #define SG_NODATA_GDT_Byte 255
47 : #define SG_NODATA_GDT_UInt16 65535
48 : #define SG_NODATA_GDT_Int16 -32767
49 : #define SG_NODATA_GDT_UInt32 4294967295U
50 : #define SG_NODATA_GDT_Int32 -2147483647
51 : #define SG_NODATA_GDT_Float32 -99999.0
52 : #define SG_NODATA_GDT_Float64 -99999.0
53 :
54 :
55 : CPL_C_START
56 : void GDALRegister_SAGA(void);
57 : CPL_C_END
58 :
59 :
60 : /************************************************************************/
61 : /* ==================================================================== */
62 : /* SAGADataset */
63 : /* ==================================================================== */
64 : /************************************************************************/
65 :
66 : class SAGARasterBand;
67 :
68 : class SAGADataset : public GDALPamDataset
69 102 : {
70 : friend class SAGARasterBand;
71 :
72 : static CPLErr WriteHeader( CPLString osHDRFilename, GDALDataType eType,
73 : int nXSize, int nYSize,
74 : double dfMinX, double dfMinY,
75 : double dfCellsize, double dfNoData,
76 : double dfZFactor, bool bTopToBottom );
77 : VSILFILE *fp;
78 :
79 : public:
80 : ~SAGADataset();
81 :
82 : static GDALDataset *Open( GDALOpenInfo * );
83 : static GDALDataset *Create( const char * pszFilename,
84 : int nXSize, int nYSize, int nBands,
85 : GDALDataType eType,
86 : char **papszParmList );
87 : static GDALDataset *CreateCopy( const char *pszFilename,
88 : GDALDataset *poSrcDS,
89 : int bStrict, char **papszOptions,
90 : GDALProgressFunc pfnProgress,
91 : void *pProgressData );
92 :
93 : virtual char **GetFileList();
94 :
95 : CPLErr GetGeoTransform( double *padfGeoTransform );
96 : CPLErr SetGeoTransform( double *padfGeoTransform );
97 : };
98 :
99 :
100 : /************************************************************************/
101 : /* ==================================================================== */
102 : /* SAGARasterBand */
103 : /* ==================================================================== */
104 : /************************************************************************/
105 :
106 : class SAGARasterBand : public GDALPamRasterBand
107 : {
108 : friend class SAGADataset;
109 :
110 : int m_Cols;
111 : int m_Rows;
112 : double m_Xmin;
113 : double m_Ymin;
114 : double m_Cellsize;
115 : double m_NoData;
116 : int m_ByteOrder;
117 : int m_nBits;
118 :
119 : void SetDataType( GDALDataType eType );
120 : void SwapBuffer(void* pImage);
121 : public:
122 :
123 : SAGARasterBand( SAGADataset *, int );
124 : ~SAGARasterBand();
125 :
126 : CPLErr IReadBlock( int, int, void * );
127 : CPLErr IWriteBlock( int, int, void * );
128 :
129 : double GetNoDataValue( int *pbSuccess = NULL );
130 : };
131 :
132 : /************************************************************************/
133 : /* SAGARasterBand() */
134 : /************************************************************************/
135 :
136 102 : SAGARasterBand::SAGARasterBand( SAGADataset *poDS, int nBand )
137 :
138 : {
139 102 : this->poDS = poDS;
140 102 : this->nBand = nBand;
141 :
142 102 : eDataType = GDT_Float32;
143 :
144 102 : nBlockXSize = poDS->GetRasterXSize();
145 102 : nBlockYSize = 1;
146 102 : }
147 :
148 : /************************************************************************/
149 : /* ~SAGARasterBand() */
150 : /************************************************************************/
151 :
152 102 : SAGARasterBand::~SAGARasterBand( )
153 :
154 : {
155 102 : }
156 :
157 : /************************************************************************/
158 : /* SetDataType() */
159 : /************************************************************************/
160 :
161 102 : void SAGARasterBand::SetDataType( GDALDataType eType )
162 :
163 : {
164 102 : eDataType = eType;
165 : return;
166 : }
167 :
168 : /************************************************************************/
169 : /* SwapBuffer() */
170 : /************************************************************************/
171 :
172 1117 : void SAGARasterBand::SwapBuffer(void* pImage)
173 : {
174 :
175 : #ifdef CPL_LSB
176 1117 : int bSwap = ( m_ByteOrder == 1);
177 : #else
178 : int bSwap = ( m_ByteOrder == 0);
179 : #endif
180 :
181 1117 : if (bSwap)
182 : {
183 0 : if ( m_nBits == 16 )
184 : {
185 0 : short* pImage16 = (short*) pImage;
186 0 : for( int iPixel=0; iPixel<nBlockXSize; iPixel++ )
187 : {
188 0 : CPL_SWAP16PTR( pImage16 + iPixel );
189 : }
190 : }
191 0 : else if ( m_nBits == 32 )
192 : {
193 0 : int* pImage32 = (int*) pImage;
194 0 : for( int iPixel=0; iPixel<nBlockXSize; iPixel++ )
195 : {
196 0 : CPL_SWAP32PTR( pImage32 + iPixel );
197 : }
198 : }
199 0 : else if ( m_nBits == 64 )
200 : {
201 0 : double* pImage64 = (double*) pImage;
202 0 : for( int iPixel=0; iPixel<nBlockXSize; iPixel++ )
203 : {
204 0 : CPL_SWAP64PTR( pImage64 + iPixel );
205 : }
206 : }
207 : }
208 :
209 1117 : }
210 :
211 : /************************************************************************/
212 : /* IReadBlock() */
213 : /************************************************************************/
214 :
215 357 : CPLErr SAGARasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
216 : void * pImage )
217 :
218 : {
219 357 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
220 0 : return CE_Failure;
221 :
222 357 : SAGADataset *poGDS = dynamic_cast<SAGADataset *>(poDS);
223 357 : vsi_l_offset offset = (vsi_l_offset) (m_nBits / 8) * nRasterXSize * (nRasterYSize - nBlockYOff - 1);
224 :
225 357 : if( VSIFSeekL( poGDS->fp, offset, SEEK_SET ) != 0 )
226 : {
227 : CPLError( CE_Failure, CPLE_FileIO,
228 0 : "Unable to seek to beginning of grid row.\n" );
229 0 : return CE_Failure;
230 : }
231 357 : if( VSIFReadL( pImage, m_nBits / 8, nBlockXSize,
232 : poGDS->fp ) != static_cast<unsigned>(nBlockXSize) )
233 : {
234 : CPLError( CE_Failure, CPLE_FileIO,
235 0 : "Unable to read block from grid file.\n" );
236 0 : return CE_Failure;
237 : }
238 :
239 357 : SwapBuffer(pImage);
240 :
241 357 : return CE_None;
242 : }
243 :
244 : /************************************************************************/
245 : /* IWriteBlock() */
246 : /************************************************************************/
247 :
248 380 : CPLErr SAGARasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
249 : void *pImage )
250 :
251 : {
252 380 : if( eAccess == GA_ReadOnly )
253 : {
254 : CPLError( CE_Failure, CPLE_NoWriteAccess,
255 0 : "Unable to write block, dataset opened read only.\n" );
256 0 : return CE_Failure;
257 : }
258 :
259 380 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
260 0 : return CE_Failure;
261 :
262 380 : vsi_l_offset offset = (vsi_l_offset) (m_nBits / 8) * nRasterXSize * (nRasterYSize - nBlockYOff - 1);
263 380 : SAGADataset *poGDS = dynamic_cast<SAGADataset *>(poDS);
264 380 : assert( poGDS != NULL );
265 :
266 380 : if( VSIFSeekL( poGDS->fp, offset, SEEK_SET ) != 0 )
267 : {
268 : CPLError( CE_Failure, CPLE_FileIO,
269 0 : "Unable to seek to beginning of grid row.\n" );
270 0 : return CE_Failure;
271 : }
272 :
273 380 : SwapBuffer(pImage);
274 :
275 : int bSuccess = ( VSIFWriteL( pImage, m_nBits / 8, nBlockXSize,
276 380 : poGDS->fp ) == static_cast<unsigned>(nBlockXSize) );
277 :
278 380 : SwapBuffer(pImage);
279 :
280 380 : if (!bSuccess)
281 : {
282 : CPLError( CE_Failure, CPLE_FileIO,
283 0 : "Unable to write block to grid file.\n" );
284 0 : return CE_Failure;
285 : }
286 :
287 380 : return CE_None;
288 : }
289 :
290 : /************************************************************************/
291 : /* GetNoDataValue() */
292 : /************************************************************************/
293 :
294 44 : double SAGARasterBand::GetNoDataValue( int * pbSuccess )
295 : {
296 44 : if( pbSuccess )
297 44 : *pbSuccess = TRUE;
298 :
299 44 : return m_NoData;
300 : }
301 :
302 :
303 : /************************************************************************/
304 : /* ==================================================================== */
305 : /* SAGADataset */
306 : /* ==================================================================== */
307 : /************************************************************************/
308 :
309 102 : SAGADataset::~SAGADataset()
310 :
311 : {
312 102 : FlushCache();
313 102 : if( fp != NULL )
314 102 : VSIFCloseL( fp );
315 102 : }
316 :
317 :
318 : /************************************************************************/
319 : /* GetFileList() */
320 : /************************************************************************/
321 :
322 23 : char** SAGADataset::GetFileList()
323 : {
324 23 : CPLString osPath = CPLGetPath( GetDescription() );
325 23 : CPLString osName = CPLGetBasename( GetDescription() );
326 23 : CPLString osHDRFilename;
327 :
328 23 : char **papszFileList = NULL;
329 :
330 : // Main data file, etc.
331 23 : papszFileList = GDALPamDataset::GetFileList();
332 :
333 : // Header file.
334 23 : osHDRFilename = CPLFormCIFilename( osPath, osName, ".sgrd" );
335 23 : papszFileList = CSLAddString( papszFileList, osHDRFilename );
336 :
337 23 : return papszFileList;
338 : }
339 :
340 : /************************************************************************/
341 : /* Open() */
342 : /************************************************************************/
343 :
344 11735 : GDALDataset *SAGADataset::Open( GDALOpenInfo * poOpenInfo )
345 :
346 : {
347 : /* -------------------------------------------------------------------- */
348 : /* We assume the user is pointing to the binary (ie. .sdat) file. */
349 : /* -------------------------------------------------------------------- */
350 11735 : if( !EQUAL(CPLGetExtension( poOpenInfo->pszFilename ), "sdat"))
351 : {
352 11545 : return NULL;
353 : }
354 :
355 190 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
356 190 : CPLString osName = CPLGetBasename( poOpenInfo->pszFilename );
357 190 : CPLString osHDRFilename;
358 :
359 190 : osHDRFilename = CPLFormCIFilename( osPath, osName, ".sgrd" );
360 :
361 :
362 : VSILFILE *fp;
363 :
364 190 : fp = VSIFOpenL( osHDRFilename, "r" );
365 :
366 190 : if( fp == NULL )
367 : {
368 88 : return NULL;
369 : }
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Is this file a SAGA header file? Read a few lines of text */
373 : /* searching for something starting with nrows or ncols. */
374 : /* -------------------------------------------------------------------- */
375 : const char *pszLine;
376 102 : int nRows = -1, nCols = -1;
377 102 : double dXmin = 0.0, dYmin = 0.0, dCellsize = 0.0, dNoData = 0.0, dZFactor = 0.0;
378 102 : int nLineCount = 0;
379 102 : char szDataFormat[20] = "DOUBLE";
380 102 : char szByteOrderBig[10] = "FALSE";
381 102 : char szTopToBottom[10] = "FALSE";
382 102 : char **papszHDR = NULL;
383 :
384 :
385 1632 : while( (pszLine = CPLReadLineL( fp )) != NULL )
386 : {
387 : char **papszTokens;
388 :
389 1428 : nLineCount++;
390 :
391 1428 : if( nLineCount > 50 || strlen(pszLine) > 1000 )
392 0 : break;
393 :
394 1428 : papszHDR = CSLAddString( papszHDR, pszLine );
395 :
396 1428 : papszTokens = CSLTokenizeStringComplex( pszLine, " =", TRUE, FALSE );
397 1428 : if( CSLCount( papszTokens ) < 2 )
398 : {
399 204 : CSLDestroy( papszTokens );
400 204 : continue;
401 : }
402 :
403 1224 : if( EQUALN(papszTokens[0],"CELLCOUNT_X",strlen("CELLCOUNT_X")) )
404 102 : nCols = atoi(papszTokens[1]);
405 1122 : else if( EQUALN(papszTokens[0],"CELLCOUNT_Y",strlen("CELLCOUNT_Y")) )
406 102 : nRows = atoi(papszTokens[1]);
407 1020 : else if( EQUALN(papszTokens[0],"POSITION_XMIN",strlen("POSITION_XMIN")) )
408 102 : dXmin = CPLAtofM(papszTokens[1]);
409 918 : else if( EQUALN(papszTokens[0],"POSITION_YMIN",strlen("POSITION_YMIN")) )
410 102 : dYmin = CPLAtofM(papszTokens[1]);
411 816 : else if( EQUALN(papszTokens[0],"CELLSIZE",strlen("CELLSIZE")) )
412 102 : dCellsize = CPLAtofM(papszTokens[1]);
413 714 : else if( EQUALN(papszTokens[0],"NODATA_VALUE",strlen("NODATA_VALUE")) )
414 102 : dNoData = CPLAtofM(papszTokens[1]);
415 612 : else if( EQUALN(papszTokens[0],"DATAFORMAT",strlen("DATAFORMAT")) )
416 102 : strncpy( szDataFormat, papszTokens[1], sizeof(szDataFormat)-1 );
417 510 : else if( EQUALN(papszTokens[0],"BYTEORDER_BIG",strlen("BYTEORDER_BIG")) )
418 102 : strncpy( szByteOrderBig, papszTokens[1], sizeof(szByteOrderBig)-1 );
419 408 : else if( EQUALN(papszTokens[0],"TOPTOBOTTOM",strlen("TOPTOBOTTOM")) )
420 102 : strncpy( szTopToBottom, papszTokens[1], sizeof(szTopToBottom)-1 );
421 306 : else if( EQUALN(papszTokens[0],"Z_FACTOR",strlen("Z_FACTOR")) )
422 102 : dZFactor = CPLAtofM(papszTokens[1]);
423 :
424 1224 : CSLDestroy( papszTokens );
425 : }
426 :
427 102 : VSIFCloseL( fp );
428 :
429 102 : CSLDestroy( papszHDR );
430 :
431 : /* -------------------------------------------------------------------- */
432 : /* Did we get the required keywords? If not we return with */
433 : /* this never having been considered to be a match. This isn't */
434 : /* an error! */
435 : /* -------------------------------------------------------------------- */
436 102 : if( nRows == -1 || nCols == -1 )
437 : {
438 0 : return NULL;
439 : }
440 :
441 102 : if (!GDALCheckDatasetDimensions(nCols, nRows))
442 : {
443 0 : return NULL;
444 : }
445 :
446 102 : if( EQUALN(szTopToBottom,"TRUE",strlen("TRUE")) )
447 : {
448 : CPLError( CE_Failure, CPLE_AppDefined,
449 : "Currently the SAGA Binary Grid driver does not support\n"
450 0 : "SAGA grids written TOPTOBOTTOM.\n");
451 0 : return NULL;
452 : }
453 102 : if( dZFactor != 1.0)
454 : {
455 : CPLError( CE_Warning, CPLE_AppDefined,
456 : "Currently the SAGA Binary Grid driver does not support\n"
457 0 : "ZFACTORs other than 1.\n");
458 : }
459 :
460 :
461 :
462 : /* -------------------------------------------------------------------- */
463 : /* Create a corresponding GDALDataset. */
464 : /* -------------------------------------------------------------------- */
465 102 : SAGADataset *poDS = new SAGADataset();
466 :
467 102 : poDS->eAccess = poOpenInfo->eAccess;
468 102 : if( poOpenInfo->eAccess == GA_ReadOnly )
469 64 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
470 : else
471 38 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
472 :
473 102 : if( poDS->fp == NULL )
474 : {
475 0 : delete poDS;
476 : CPLError( CE_Failure, CPLE_OpenFailed,
477 : "VSIFOpenL(%s) failed unexpectedly.",
478 0 : poOpenInfo->pszFilename );
479 0 : return NULL;
480 : }
481 :
482 102 : poDS->nRasterXSize = nCols;
483 102 : poDS->nRasterYSize = nRows;
484 :
485 102 : SAGARasterBand *poBand = new SAGARasterBand( poDS, 1 );
486 :
487 :
488 : /* -------------------------------------------------------------------- */
489 : /* Figure out the byte order. */
490 : /* -------------------------------------------------------------------- */
491 102 : if( EQUALN(szByteOrderBig,"TRUE",strlen("TRUE")) )
492 0 : poBand->m_ByteOrder = 1;
493 102 : else if( EQUALN(szByteOrderBig,"FALSE",strlen("FALSE")) )
494 102 : poBand->m_ByteOrder = 0;
495 :
496 :
497 : /* -------------------------------------------------------------------- */
498 : /* Figure out the data type. */
499 : /* -------------------------------------------------------------------- */
500 102 : if( EQUAL(szDataFormat,"BIT") )
501 : {
502 0 : poBand->SetDataType(GDT_Byte);
503 0 : poBand->m_nBits = 8;
504 : }
505 102 : else if( EQUAL(szDataFormat,"BYTE_UNSIGNED") )
506 : {
507 13 : poBand->SetDataType(GDT_Byte);
508 13 : poBand->m_nBits = 8;
509 : }
510 89 : else if( EQUAL(szDataFormat,"BYTE") )
511 : {
512 0 : poBand->SetDataType(GDT_Byte);
513 0 : poBand->m_nBits = 8;
514 : }
515 89 : else if( EQUAL(szDataFormat,"SHORTINT_UNSIGNED") )
516 : {
517 13 : poBand->SetDataType(GDT_UInt16);
518 13 : poBand->m_nBits = 16;
519 : }
520 76 : else if( EQUAL(szDataFormat,"SHORTINT") )
521 : {
522 13 : poBand->SetDataType(GDT_Int16);
523 13 : poBand->m_nBits = 16;
524 : }
525 63 : else if( EQUAL(szDataFormat,"INTEGER_UNSIGNED") )
526 : {
527 13 : poBand->SetDataType(GDT_UInt32);
528 13 : poBand->m_nBits = 32;
529 : }
530 50 : else if( EQUAL(szDataFormat,"INTEGER") )
531 : {
532 13 : poBand->SetDataType(GDT_Int32);
533 13 : poBand->m_nBits = 32;
534 : }
535 37 : else if( EQUAL(szDataFormat,"FLOAT") )
536 : {
537 26 : poBand->SetDataType(GDT_Float32);
538 26 : poBand->m_nBits = 32;
539 : }
540 11 : else if( EQUAL(szDataFormat,"DOUBLE") )
541 : {
542 11 : poBand->SetDataType(GDT_Float64);
543 11 : poBand->m_nBits = 64;
544 : }
545 : else
546 : {
547 : CPLError( CE_Failure, CPLE_NotSupported,
548 : "SAGA driver does not support the dataformat %s.",
549 0 : szDataFormat );
550 0 : delete poBand;
551 0 : delete poDS;
552 0 : return NULL;
553 : }
554 :
555 : /* -------------------------------------------------------------------- */
556 : /* Save band information */
557 : /* -------------------------------------------------------------------- */
558 102 : poBand->m_Xmin = dXmin;
559 102 : poBand->m_Ymin = dYmin;
560 102 : poBand->m_NoData = dNoData;
561 102 : poBand->m_Cellsize = dCellsize;
562 102 : poBand->m_Rows = nRows;
563 102 : poBand->m_Cols = nCols;
564 :
565 102 : poDS->SetBand( 1, poBand );
566 :
567 : /* -------------------------------------------------------------------- */
568 : /* Initialize any PAM information. */
569 : /* -------------------------------------------------------------------- */
570 102 : poDS->SetDescription( poOpenInfo->pszFilename );
571 102 : poDS->TryLoadXML();
572 :
573 : /* -------------------------------------------------------------------- */
574 : /* Check for external overviews. */
575 : /* -------------------------------------------------------------------- */
576 102 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename, poOpenInfo->papszSiblingFiles );
577 :
578 102 : return poDS;
579 : }
580 :
581 : /************************************************************************/
582 : /* GetGeoTransform() */
583 : /************************************************************************/
584 :
585 11 : CPLErr SAGADataset::GetGeoTransform( double *padfGeoTransform )
586 : {
587 11 : if( padfGeoTransform == NULL )
588 0 : return CE_Failure;
589 :
590 11 : SAGARasterBand *poGRB = dynamic_cast<SAGARasterBand *>(GetRasterBand( 1 ));
591 :
592 11 : if( poGRB == NULL )
593 : {
594 0 : padfGeoTransform[0] = 0;
595 0 : padfGeoTransform[1] = 1;
596 0 : padfGeoTransform[2] = 0;
597 0 : padfGeoTransform[3] = 0;
598 0 : padfGeoTransform[4] = 0;
599 0 : padfGeoTransform[5] = 1;
600 0 : return CE_Failure;
601 : }
602 :
603 : /* check if we have a PAM GeoTransform stored */
604 11 : CPLPushErrorHandler( CPLQuietErrorHandler );
605 11 : CPLErr eErr = GDALPamDataset::GetGeoTransform( padfGeoTransform );
606 11 : CPLPopErrorHandler();
607 :
608 11 : if( eErr == CE_None )
609 0 : return CE_None;
610 :
611 11 : padfGeoTransform[1] = poGRB->m_Cellsize;
612 11 : padfGeoTransform[5] = poGRB->m_Cellsize * -1.0;
613 11 : padfGeoTransform[0] = poGRB->m_Xmin - poGRB->m_Cellsize / 2;
614 11 : padfGeoTransform[3] = poGRB->m_Ymin + (nRasterYSize - 1) * poGRB->m_Cellsize + poGRB->m_Cellsize / 2;
615 :
616 : /* tilt/rotation is not supported by SAGA grids */
617 11 : padfGeoTransform[4] = 0.0;
618 11 : padfGeoTransform[2] = 0.0;
619 :
620 11 : return CE_None;
621 : }
622 :
623 : /************************************************************************/
624 : /* SetGeoTransform() */
625 : /************************************************************************/
626 :
627 23 : CPLErr SAGADataset::SetGeoTransform( double *padfGeoTransform )
628 : {
629 :
630 23 : if( eAccess == GA_ReadOnly )
631 : {
632 : CPLError( CE_Failure, CPLE_NoWriteAccess,
633 0 : "Unable to set GeoTransform, dataset opened read only.\n" );
634 0 : return CE_Failure;
635 : }
636 :
637 23 : SAGARasterBand *poGRB = dynamic_cast<SAGARasterBand *>(GetRasterBand( 1 ));
638 :
639 23 : if( poGRB == NULL || padfGeoTransform == NULL)
640 0 : return CE_Failure;
641 :
642 23 : if( padfGeoTransform[1] != padfGeoTransform[5] * -1.0 )
643 : {
644 : CPLError( CE_Failure, CPLE_NotSupported,
645 : "Unable to set GeoTransform, SAGA binary grids only support "
646 0 : "the same cellsize in x-y.\n" );
647 0 : return CE_Failure;
648 : }
649 :
650 23 : double dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
651 : double dfMinY =
652 23 : padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
653 :
654 23 : CPLString osPath = CPLGetPath( GetDescription() );
655 23 : CPLString osName = CPLGetBasename( GetDescription() );
656 23 : CPLString osHDRFilename = CPLFormCIFilename( osPath, osName, ".sgrd" );
657 :
658 : CPLErr eErr = WriteHeader( osHDRFilename, poGRB->GetRasterDataType(),
659 : poGRB->nRasterXSize, poGRB->nRasterYSize,
660 : dfMinX, dfMinY, padfGeoTransform[1],
661 23 : poGRB->m_NoData, 1.0, false );
662 :
663 :
664 23 : if( eErr == CE_None )
665 : {
666 23 : poGRB->m_Xmin = dfMinX;
667 23 : poGRB->m_Ymin = dfMinY;
668 23 : poGRB->m_Cellsize = padfGeoTransform[1];
669 23 : poGRB->m_Cols = nRasterXSize;
670 23 : poGRB->m_Rows = nRasterYSize;
671 : }
672 :
673 23 : return eErr;
674 : }
675 :
676 : /************************************************************************/
677 : /* WriteHeader() */
678 : /************************************************************************/
679 :
680 61 : CPLErr SAGADataset::WriteHeader( CPLString osHDRFilename, GDALDataType eType,
681 : int nXSize, int nYSize,
682 : double dfMinX, double dfMinY,
683 : double dfCellsize, double dfNoData,
684 : double dfZFactor, bool bTopToBottom )
685 :
686 : {
687 : VSILFILE *fp;
688 :
689 61 : fp = VSIFOpenL( osHDRFilename, "wt" );
690 :
691 61 : if( fp == NULL )
692 : {
693 : CPLError( CE_Failure, CPLE_OpenFailed,
694 : "Failed to write .sgrd file %s.",
695 0 : osHDRFilename.c_str() );
696 0 : return CE_Failure;
697 : }
698 :
699 61 : VSIFPrintfL( fp, "NAME\t= %s\n", CPLGetBasename( osHDRFilename ) );
700 61 : VSIFPrintfL( fp, "DESCRIPTION\t=\n" );
701 61 : VSIFPrintfL( fp, "UNIT\t=\n" );
702 61 : VSIFPrintfL( fp, "DATAFILE_OFFSET\t= 0\n" );
703 :
704 61 : if( eType == GDT_Int32 )
705 8 : VSIFPrintfL( fp, "DATAFORMAT\t= INTEGER\n" );
706 53 : else if( eType == GDT_UInt32 )
707 8 : VSIFPrintfL( fp, "DATAFORMAT\t= INTEGER_UNSIGNED\n" );
708 45 : else if( eType == GDT_Int16 )
709 8 : VSIFPrintfL( fp, "DATAFORMAT\t= SHORTINT\n" );
710 37 : else if( eType == GDT_UInt16 )
711 8 : VSIFPrintfL( fp, "DATAFORMAT\t= SHORTINT_UNSIGNED\n" );
712 29 : else if( eType == GDT_Byte )
713 8 : VSIFPrintfL( fp, "DATAFORMAT\t= BYTE_UNSIGNED\n" );
714 21 : else if( eType == GDT_Float32 )
715 13 : VSIFPrintfL( fp, "DATAFORMAT\t= FLOAT\n" );
716 : else //if( eType == GDT_Float64 )
717 8 : VSIFPrintfL( fp, "DATAFORMAT\t= DOUBLE\n" );
718 : #ifdef CPL_LSB
719 61 : VSIFPrintfL( fp, "BYTEORDER_BIG\t= FALSE\n" );
720 : #else
721 : VSIFPrintfL( fp, "BYTEORDER_BIG\t= TRUE\n" );
722 : #endif
723 :
724 61 : VSIFPrintfL( fp, "POSITION_XMIN\t= %.10f\n", dfMinX );
725 61 : VSIFPrintfL( fp, "POSITION_YMIN\t= %.10f\n", dfMinY );
726 61 : VSIFPrintfL( fp, "CELLCOUNT_X\t= %d\n", nXSize );
727 61 : VSIFPrintfL( fp, "CELLCOUNT_Y\t= %d\n", nYSize );
728 61 : VSIFPrintfL( fp, "CELLSIZE\t= %.10f\n", dfCellsize );
729 61 : VSIFPrintfL( fp, "Z_FACTOR\t= %f\n", dfZFactor );
730 61 : VSIFPrintfL( fp, "NODATA_VALUE\t= %f\n", dfNoData );
731 61 : if (bTopToBottom)
732 0 : VSIFPrintfL( fp, "TOPTOBOTTOM\t= TRUE\n" );
733 : else
734 61 : VSIFPrintfL( fp, "TOPTOBOTTOM\t= FALSE\n" );
735 :
736 :
737 61 : VSIFCloseL( fp );
738 :
739 61 : return CE_None;
740 : }
741 :
742 :
743 : /************************************************************************/
744 : /* Create() */
745 : /************************************************************************/
746 :
747 70 : GDALDataset *SAGADataset::Create( const char * pszFilename,
748 : int nXSize, int nYSize, int nBands,
749 : GDALDataType eType,
750 : char **papszParmList )
751 :
752 : {
753 70 : if( nXSize <= 0 || nYSize <= 0 )
754 : {
755 : CPLError( CE_Failure, CPLE_IllegalArg,
756 : "Unable to create grid, both X and Y size must be "
757 0 : "non-negative.\n" );
758 :
759 0 : return NULL;
760 : }
761 :
762 70 : if( nBands != 1 )
763 : {
764 : CPLError( CE_Failure, CPLE_IllegalArg,
765 15 : "SAGA Binary Grid only supports 1 band" );
766 15 : return NULL;
767 : }
768 :
769 55 : if( eType != GDT_Byte && eType != GDT_UInt16 && eType != GDT_Int16
770 : && eType != GDT_UInt32 && eType != GDT_Int32 && eType != GDT_Float32
771 : && eType != GDT_Float64 )
772 : {
773 : CPLError( CE_Failure, CPLE_AppDefined,
774 : "SAGA Binary Grid only supports Byte, UInt16, Int16, "
775 : "UInt32, Int32, Float32 and Float64 datatypes. Unable to "
776 8 : "create with type %s.\n", GDALGetDataTypeName( eType ) );
777 :
778 8 : return NULL;
779 : }
780 :
781 47 : VSILFILE *fp = VSIFOpenL( pszFilename, "w+b" );
782 :
783 47 : if( fp == NULL )
784 : {
785 : CPLError( CE_Failure, CPLE_OpenFailed,
786 : "Attempt to create file '%s' failed.\n",
787 9 : pszFilename );
788 9 : return NULL;
789 : }
790 :
791 : char abyNoData[8];
792 38 : double dfNoDataVal = 0.0;
793 :
794 38 : const char* pszNoDataValue = CSLFetchNameValue(papszParmList, "NODATA_VALUE");
795 38 : if (pszNoDataValue)
796 : {
797 2 : dfNoDataVal = CPLAtofM(pszNoDataValue);
798 : }
799 : else
800 : {
801 36 : switch (eType) /* GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32 */
802 : { /* GDT_Int32, GDT_Float32, GDT_Float64 */
803 : case (GDT_Byte):
804 : {
805 5 : dfNoDataVal = SG_NODATA_GDT_Byte;
806 5 : break;
807 : }
808 : case (GDT_UInt16):
809 : {
810 5 : dfNoDataVal = SG_NODATA_GDT_UInt16;
811 5 : break;
812 : }
813 : case (GDT_Int16):
814 : {
815 5 : dfNoDataVal = SG_NODATA_GDT_Int16;
816 5 : break;
817 : }
818 : case (GDT_UInt32):
819 : {
820 5 : dfNoDataVal = SG_NODATA_GDT_UInt32;
821 5 : break;
822 : }
823 : case (GDT_Int32):
824 : {
825 5 : dfNoDataVal = SG_NODATA_GDT_Int32;
826 5 : break;
827 : }
828 : default:
829 : case (GDT_Float32):
830 : {
831 6 : dfNoDataVal = SG_NODATA_GDT_Float32;
832 6 : break;
833 : }
834 : case (GDT_Float64):
835 : {
836 5 : dfNoDataVal = SG_NODATA_GDT_Float64;
837 : break;
838 : }
839 : }
840 : }
841 :
842 : GDALCopyWords(&dfNoDataVal, GDT_Float64, 0,
843 38 : abyNoData, eType, 0, 1);
844 :
845 38 : CPLString osHdrFilename = CPLResetExtension( pszFilename, "sgrd" );
846 : CPLErr eErr = WriteHeader( osHdrFilename, eType,
847 : nXSize, nYSize,
848 : 0.0, 0.0, 1.0,
849 38 : dfNoDataVal, 1.0, false );
850 :
851 38 : if( eErr != CE_None )
852 : {
853 0 : VSIFCloseL( fp );
854 0 : return NULL;
855 : }
856 :
857 38 : if (CSLFetchBoolean( papszParmList , "FILL_NODATA", TRUE ))
858 : {
859 22 : int nDataTypeSize = GDALGetDataTypeSize(eType) / 8;
860 22 : GByte* pabyNoDataBuf = (GByte*) VSIMalloc2(nDataTypeSize, nXSize);
861 22 : if (pabyNoDataBuf == NULL)
862 : {
863 0 : VSIFCloseL( fp );
864 0 : return NULL;
865 : }
866 :
867 886 : for( int iCol = 0; iCol < nXSize; iCol++)
868 : {
869 864 : memcpy(pabyNoDataBuf + iCol * nDataTypeSize, abyNoData, nDataTypeSize);
870 : }
871 :
872 886 : for( int iRow = 0; iRow < nYSize; iRow++ )
873 : {
874 864 : if( VSIFWriteL( pabyNoDataBuf, nDataTypeSize, nXSize, fp ) != (unsigned)nXSize )
875 : {
876 0 : VSIFCloseL( fp );
877 0 : VSIFree(pabyNoDataBuf);
878 : CPLError( CE_Failure, CPLE_FileIO,
879 0 : "Unable to write grid cell. Disk full?\n" );
880 0 : return NULL;
881 : }
882 : }
883 :
884 22 : VSIFree(pabyNoDataBuf);
885 : }
886 :
887 38 : VSIFCloseL( fp );
888 :
889 38 : return (GDALDataset *)GDALOpen( pszFilename, GA_Update );
890 : }
891 :
892 : /************************************************************************/
893 : /* CreateCopy() */
894 : /************************************************************************/
895 :
896 34 : GDALDataset *SAGADataset::CreateCopy( const char *pszFilename,
897 : GDALDataset *poSrcDS,
898 : int bStrict, char **papszOptions,
899 : GDALProgressFunc pfnProgress,
900 : void *pProgressData )
901 : {
902 34 : if( pfnProgress == NULL )
903 0 : pfnProgress = GDALDummyProgress;
904 :
905 34 : int nBands = poSrcDS->GetRasterCount();
906 34 : if (nBands == 0)
907 : {
908 : CPLError( CE_Failure, CPLE_NotSupported,
909 1 : "SAGA driver does not support source dataset with zero band.\n");
910 1 : return NULL;
911 : }
912 33 : else if (nBands > 1)
913 : {
914 4 : if( bStrict )
915 : {
916 : CPLError( CE_Failure, CPLE_NotSupported,
917 : "Unable to create copy, SAGA Binary Grid "
918 4 : "format only supports one raster band.\n" );
919 4 : return NULL;
920 : }
921 : else
922 : CPLError( CE_Warning, CPLE_NotSupported,
923 : "SAGA Binary Grid format only supports one "
924 0 : "raster band, first band will be copied.\n" );
925 : }
926 :
927 29 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );
928 :
929 29 : char** papszCreateOptions = NULL;
930 29 : papszCreateOptions = CSLSetNameValue(papszCreateOptions, "FILL_NODATA", "NO");
931 :
932 29 : int bHasNoDataValue = FALSE;
933 29 : double dfNoDataValue = poSrcBand->GetNoDataValue(&bHasNoDataValue);
934 29 : if (bHasNoDataValue)
935 : papszCreateOptions = CSLSetNameValue(papszCreateOptions, "NODATA_VALUE",
936 2 : CPLSPrintf("%.16g", dfNoDataValue));
937 :
938 : GDALDataset* poDstDS =
939 : Create(pszFilename, poSrcBand->GetXSize(), poSrcBand->GetYSize(),
940 29 : 1, poSrcBand->GetRasterDataType(), papszCreateOptions);
941 29 : CSLDestroy(papszCreateOptions);
942 :
943 29 : if (poDstDS == NULL)
944 13 : return NULL;
945 :
946 : /* -------------------------------------------------------------------- */
947 : /* Copy band data. */
948 : /* -------------------------------------------------------------------- */
949 :
950 : CPLErr eErr;
951 :
952 : eErr = GDALDatasetCopyWholeRaster( (GDALDatasetH) poSrcDS,
953 : (GDALDatasetH) poDstDS,
954 : NULL,
955 16 : pfnProgress, pProgressData );
956 :
957 16 : if (eErr == CE_Failure)
958 : {
959 0 : delete poDstDS;
960 0 : return NULL;
961 : }
962 :
963 : double adfGeoTransform[6];
964 :
965 16 : poSrcDS->GetGeoTransform( adfGeoTransform );
966 16 : poDstDS->SetGeoTransform( adfGeoTransform );
967 :
968 16 : return poDstDS;
969 : }
970 :
971 : /************************************************************************/
972 : /* GDALRegister_SAGA() */
973 : /************************************************************************/
974 :
975 582 : void GDALRegister_SAGA()
976 :
977 : {
978 : GDALDriver *poDriver;
979 :
980 582 : if( GDALGetDriverByName( "SAGA" ) == NULL )
981 : {
982 561 : poDriver = new GDALDriver();
983 :
984 561 : poDriver->SetDescription( "SAGA" );
985 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
986 561 : "SAGA GIS Binary Grid (.sdat)" );
987 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
988 561 : "frmt_various.html#SAGA" );
989 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "sdat" );
990 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
991 561 : "Byte Int16 UInt16 Int32 UInt32 Float32 Float64" );
992 :
993 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
994 :
995 561 : poDriver->pfnOpen = SAGADataset::Open;
996 561 : poDriver->pfnCreate = SAGADataset::Create;
997 561 : poDriver->pfnCreateCopy = SAGADataset::CreateCopy;
998 :
999 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1000 : }
1001 582 : }
|
