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