1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // //
10 // io_gdal //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // gdal_driver.cpp //
15 // //
16 // Copyright (C) 2007 by O. Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.de //
40 // //
41 // contact: Olaf Conrad //
42 // Bundesstr. 55 //
43 // D-20146 Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "gdal_driver.h"
50
51 #include <gdal_vrt.h>
52
53 #include <ogr_srs_api.h>
54
55 #include <cpl_string.h>
56 #include <cpl_error.h>
57
58
59 ///////////////////////////////////////////////////////////
60 // //
61 // //
62 // //
63 ///////////////////////////////////////////////////////////
64
65 //---------------------------------------------------------
66 CSG_GDAL_Drivers gSG_GDAL_Drivers;
67
SG_Get_GDAL_Drivers(void)68 const CSG_GDAL_Drivers & SG_Get_GDAL_Drivers (void)
69 {
70 return( gSG_GDAL_Drivers );
71 }
72
73
74 ///////////////////////////////////////////////////////////
75 // //
76 // //
77 // //
78 ///////////////////////////////////////////////////////////
79
80 //---------------------------------------------------------
CSG_GDAL_Drivers(void)81 CSG_GDAL_Drivers::CSG_GDAL_Drivers(void)
82 {
83 GDALAllRegister();
84
85 // affects Windows only, might be appropriate for applications
86 // that treat filenames as being in the local encoding.
87 // for more info see: http://trac.osgeo.org/gdal/wiki/ConfigOptions
88 CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
89 }
90
91 //---------------------------------------------------------
~CSG_GDAL_Drivers(void)92 CSG_GDAL_Drivers::~CSG_GDAL_Drivers(void)
93 {
94 GDALDestroyDriverManager();
95 }
96
97
98 ///////////////////////////////////////////////////////////
99 // //
100 ///////////////////////////////////////////////////////////
101
102 //---------------------------------------------------------
Get_Version(void) const103 CSG_String CSG_GDAL_Drivers::Get_Version(void) const
104 {
105 return( GDALVersionInfo("RELEASE_NAME") );
106 }
107
108 //---------------------------------------------------------
Get_Count(void) const109 int CSG_GDAL_Drivers::Get_Count(void) const
110 {
111 return( GDALGetDriverCount() );
112 }
113
114 //---------------------------------------------------------
Get_Driver(const CSG_String & Name) const115 GDALDriverH CSG_GDAL_Drivers::Get_Driver(const CSG_String &Name) const
116 {
117 return( GDALGetDriverByName(Name) );
118 }
119
120 //---------------------------------------------------------
Get_Driver(int Index) const121 GDALDriverH CSG_GDAL_Drivers::Get_Driver(int Index) const
122 {
123 return( GDALGetDriver(Index) );
124 }
125
126 //---------------------------------------------------------
Get_Name(int Index) const127 CSG_String CSG_GDAL_Drivers::Get_Name(int Index) const
128 {
129 const char *s = GDALGetMetadataItem(Get_Driver(Index), GDAL_DMD_LONGNAME, "");
130
131 return( s ? s : "" );
132 }
133
134 //---------------------------------------------------------
Get_Description(int Index) const135 CSG_String CSG_GDAL_Drivers::Get_Description(int Index) const
136 {
137 const char *s = GDALGetDescription(Get_Driver(Index));
138
139 return( s ? s : "" );
140 }
141
142 //---------------------------------------------------------
Get_Extension(int Index) const143 CSG_String CSG_GDAL_Drivers::Get_Extension(int Index) const
144 {
145 const char *s = GDALGetMetadataItem(Get_Driver(Index), GDAL_DMD_EXTENSION, "");
146
147 return( s ? s : "" );
148 }
149
150 //---------------------------------------------------------
has_Capability(GDALDriverH pDriver,const char * Capapility)151 bool CSG_GDAL_Drivers::has_Capability(GDALDriverH pDriver, const char *Capapility)
152 {
153 const char *s = GDALGetMetadataItem(pDriver, Capapility, "");
154
155 return( s && SG_STR_CMP(s, "YES") == 0 );
156 }
157
158 //---------------------------------------------------------
is_Raster(int Index) const159 bool CSG_GDAL_Drivers::is_Raster(int Index) const
160 {
161 #ifdef GDAL_V2_0_OR_NEWER
162 return( has_Capability(Get_Driver(Index), GDAL_DCAP_RASTER) );
163 #else
164 return( true );
165 #endif
166 }
167
168 //---------------------------------------------------------
Can_Read(int Index) const169 bool CSG_GDAL_Drivers::Can_Read(int Index) const
170 {
171 #ifdef GDAL_DCAP_OPEN
172 return( has_Capability(Get_Driver(Index), GDAL_DCAP_OPEN) );
173 #else
174 return( true );
175 #endif
176 }
177
178 //---------------------------------------------------------
Can_Write(int Index) const179 bool CSG_GDAL_Drivers::Can_Write(int Index) const
180 {
181 return( has_Capability(Get_Driver(Index), GDAL_DCAP_CREATE) );
182 }
183
184 //---------------------------------------------------------
Can_Copy(int Index) const185 bool CSG_GDAL_Drivers::Can_Copy(int Index) const
186 {
187 return( has_Capability(Get_Driver(Index), GDAL_DCAP_CREATECOPY) );
188 }
189
190
191 ///////////////////////////////////////////////////////////
192 // //
193 ///////////////////////////////////////////////////////////
194
195 //---------------------------------------------------------
Get_GDAL_Type(TSG_Data_Type Type)196 int CSG_GDAL_Drivers::Get_GDAL_Type(TSG_Data_Type Type)
197 {
198 switch( Type )
199 {
200 case SG_DATATYPE_Bit : return( GDT_Byte ); // Eight bit unsigned integer
201 case SG_DATATYPE_Byte : return( GDT_Byte ); // Eight bit unsigned integer
202 case SG_DATATYPE_Char : return( GDT_Byte ); // Eight bit unsigned integer
203 case SG_DATATYPE_Word : return( GDT_UInt16 ); // Sixteen bit unsigned integer
204 case SG_DATATYPE_Short : return( GDT_Int16 ); // Sixteen bit signed integer
205 case SG_DATATYPE_DWord : return( GDT_UInt32 ); // Thirty two bit unsigned integer
206 case SG_DATATYPE_Int : return( GDT_Int32 ); // Thirty two bit signed integer
207 case SG_DATATYPE_Float : return( GDT_Float32 ); // Thirty two bit floating point
208 case SG_DATATYPE_Double: return( GDT_Float64 ); // Sixty four bit floating point
209
210 default : return( GDT_Float64 );
211 }
212 }
213
214 //---------------------------------------------------------
Get_SAGA_Type(int Type)215 TSG_Data_Type CSG_GDAL_Drivers::Get_SAGA_Type(int Type)
216 {
217 switch( Type )
218 {
219 case GDT_Byte : return( SG_DATATYPE_Byte ); // Eight bit unsigned integer
220 case GDT_UInt16 : return( SG_DATATYPE_Word ); // Sixteen bit unsigned integer
221 case GDT_Int16 : return( SG_DATATYPE_Short ); // Sixteen bit signed integer
222 case GDT_UInt32 : return( SG_DATATYPE_DWord ); // Thirty two bit unsigned integer
223 case GDT_Int32 : return( SG_DATATYPE_Int ); // Thirty two bit signed integer
224 case GDT_Float32 : return( SG_DATATYPE_Float ); // Thirty two bit floating point
225 case GDT_Float64 : return( SG_DATATYPE_Double ); // Sixty four bit floating point
226
227 case GDT_CInt16 : return( SG_DATATYPE_Undefined ); // Complex Int16
228 case GDT_CInt32 : return( SG_DATATYPE_Undefined ); // Complex Int32
229 case GDT_CFloat32: return( SG_DATATYPE_Undefined ); // Complex Float32
230 case GDT_CFloat64: return( SG_DATATYPE_Undefined ); // Complex Float64
231
232 default : return( SG_DATATYPE_Undefined );
233 }
234 }
235
236
237 ///////////////////////////////////////////////////////////
238 // //
239 // //
240 // //
241 ///////////////////////////////////////////////////////////
242
243 //---------------------------------------------------------
CSG_GDAL_DataSet(void)244 CSG_GDAL_DataSet::CSG_GDAL_DataSet(void)
245 {
246 m_pDataSet = m_pVrtSource = NULL;
247
248 m_TF_A.Create(2);
249 m_TF_B.Create(2, 2);
250 }
251
252 //---------------------------------------------------------
CSG_GDAL_DataSet(const CSG_String & File_Name)253 CSG_GDAL_DataSet::CSG_GDAL_DataSet(const CSG_String &File_Name)
254 {
255 m_pDataSet = m_pVrtSource = NULL;
256
257 m_TF_A.Create(2);
258 m_TF_B.Create(2, 2);
259
260 Open_Read(File_Name);
261 }
262
263 //---------------------------------------------------------
~CSG_GDAL_DataSet(void)264 CSG_GDAL_DataSet::~CSG_GDAL_DataSet(void)
265 {
266 Close();
267 }
268
269
270 ///////////////////////////////////////////////////////////
271 // //
272 ///////////////////////////////////////////////////////////
273
274 //---------------------------------------------------------
Open_Read(const CSG_String & File_Name,const char * Drivers[])275 bool CSG_GDAL_DataSet::Open_Read(const CSG_String &File_Name, const char *Drivers[])
276 {
277 Close();
278
279 #ifdef GDAL_V2_0_OR_NEWER
280 if( Drivers )
281 {
282 m_pDataSet = GDALOpenEx(File_Name, GA_ReadOnly, Drivers, NULL, NULL);
283 }
284 #else
285 m_pDataSet = NULL;
286 #endif
287
288 if( !m_pDataSet && (m_pDataSet = GDALOpen(File_Name, GA_ReadOnly)) == NULL )
289 {
290 return( false );
291 }
292
293 m_File_Name = File_Name;
294
295 m_Access = SG_GDAL_IO_READ;
296
297 return( _Set_Transformation() );
298 }
299
300 //---------------------------------------------------------
Open_Read(const CSG_String & File_Name,const CSG_Grid_System & System)301 bool CSG_GDAL_DataSet::Open_Read(const CSG_String &File_Name, const CSG_Grid_System &System)
302 {
303 Close();
304
305 if( (m_pVrtSource = GDALOpen(File_Name, GA_ReadOnly)) == NULL )
306 {
307 return( false );
308 }
309
310 //-----------------------------------------------------
311 if( (m_pDataSet = VRTCreate(System.Get_NX(), System.Get_NY())) == NULL )
312 {
313 Close();
314
315 return( false );
316 }
317
318 GDALSetProjection(m_pDataSet, GDALGetProjectionRef(m_pVrtSource));
319
320 double Transform[6] =
321 {
322 System.Get_XMin(true), System.Get_Cellsize(), 0.,
323 System.Get_YMax(true), 0., -System.Get_Cellsize()
324 };
325
326 GDALSetGeoTransform(m_pDataSet, Transform);
327
328 //-----------------------------------------------------
329 GDALGetGeoTransform(m_pVrtSource, Transform);
330
331 if( Transform[2] != 0. || Transform[4] != 0. )
332 {
333 return( false ); // geotransform is rotated, this configuration is not supported...
334 }
335
336 int xOff = (int)floor((System.Get_XMin (true) - Transform[0]) / Transform[1] + 0.001);
337 int yOff = (int)floor((System.Get_YMax (true) - Transform[3]) / Transform[5] + 0.001);
338 int xSize = (int) ( System.Get_XRange(true) / Transform[1] + 0.5 );
339 int ySize = (int) ( System.Get_YRange(true) / fabs(Transform[5]) + 0.5 );
340
341 //-----------------------------------------------------
342 for(int i=0; i<GDALGetRasterCount(m_pVrtSource); i++)
343 {
344 GDALRasterBandH pSrcBand = GDALGetRasterBand(m_pVrtSource, i + 1);
345
346 GDALAddBand(m_pDataSet, GDALGetRasterDataType(pSrcBand), NULL);
347
348 VRTSourcedRasterBandH pVrtBand = (VRTSourcedRasterBandH)GDALGetRasterBand(m_pDataSet, i + 1);
349
350 VRTAddSimpleSource(pVrtBand, pSrcBand,
351 xOff, yOff, xSize, ySize, 0, 0, System.Get_NX(), System.Get_NY(), "near", VRT_NODATA_UNSET
352 );
353
354 int bSuccess; double zNoData = GDALGetRasterNoDataValue(pSrcBand, &bSuccess);
355
356 if( bSuccess )
357 {
358 GDALSetRasterNoDataValue(pVrtBand, zNoData);
359 }
360
361 //#if GDAL_VERSION_MAJOR >= 2 // instead of pVrtBand->AddSimpleSource(...)
362 // VRTSimpleSource *pSrcSimple = new VRTSimpleSource();
363 //
364 // // pSrcSimple->SetResampling(pszResampling);
365 //
366 // pVrtBand->ConfigureSource(pSrcSimple, pSrcBand, 0,
367 // xOff, yOff, xSize, ySize, 0, 0, System.Get_NX(), System.Get_NY()
368 // );
369 //
370 // pVrtBand->AddSource(pSrcSimple);
371 //#endif
372 }
373
374 //-----------------------------------------------------
375 m_File_Name = File_Name;
376
377 m_Access = SG_GDAL_IO_READ;
378
379 return( _Set_Transformation() );
380 }
381
382 //---------------------------------------------------------
Open_Read(const CSG_String & File_Name,const TSG_Rect & Extent)383 bool CSG_GDAL_DataSet::Open_Read(const CSG_String &File_Name, const TSG_Rect &Extent)
384 {
385 CSG_GDAL_DataSet DataSet;
386
387 if( DataSet.Open_Read(File_Name) == false )
388 {
389 return( false );
390 }
391
392 double c = DataSet.Get_System().Get_Cellsize();
393 TSG_Rect r = DataSet.Get_System().Get_Extent(true);
394
395 r.xMin = r.xMin + (floor((Extent.xMin - r.xMin) / c) + 0.5) * c;
396 r.xMax = r.xMax + (ceil ((Extent.xMax - r.xMax) / c) - 0.5) * c;
397 r.yMin = r.yMin + (floor((Extent.yMin - r.yMin) / c) + 0.5) * c;
398 r.yMax = r.yMax + (ceil ((Extent.yMax - r.yMax) / c) - 0.5) * c;
399
400 CSG_Grid_System System(c, r);
401
402 return( System.is_Valid() && System.Get_Extent(true).Intersects(DataSet.Get_System().Get_Extent(true)) && Open_Read(File_Name, System) );
403 }
404
405 //---------------------------------------------------------
Get_SubDataSets(bool bDescription) const406 CSG_Strings CSG_GDAL_DataSet::Get_SubDataSets(bool bDescription) const
407 {
408 CSG_MetaData MetaData; Get_MetaData(MetaData, "SUBDATASETS");
409
410 CSG_Strings SubDataSets;
411
412 const SG_Char *Type = bDescription ? SG_T("DESC") : SG_T("NAME");
413
414 for(int i=0; i==SubDataSets.Get_Count(); i++)
415 {
416 CSG_MetaData *pSubDataSet = MetaData(CSG_String::Format("SUBDATASET_%d_%s", i + 1, Type));
417
418 if( pSubDataSet )
419 {
420 SubDataSets += pSubDataSet->Get_Content();
421 }
422 }
423
424 return( SubDataSets );
425 }
426
427 //---------------------------------------------------------
_Get_Transformation(double Transform[6])428 bool CSG_GDAL_DataSet::_Get_Transformation(double Transform[6])
429 {
430 if( GDALGetGeoTransform(m_pDataSet, Transform) == CE_None )
431 {
432 return( true );
433 }
434
435 //-----------------------------------------------------
436 Transform[0] = 0.;
437 Transform[1] = 1.;
438 Transform[2] = 0.;
439 Transform[3] = 0.;
440 Transform[4] = 0.;
441 Transform[5] = 1.;
442
443 bool bResult = false; CSG_String Data;
444
445 if( Get_MetaData_Item(Data, "XORIG") && Data.asDouble(Transform[0]) ) { bResult = true; }
446 if( Get_MetaData_Item(Data, "XCELL") && Data.asDouble(Transform[1]) ) { bResult = true; }
447 if( Get_MetaData_Item(Data, "YORIG") && Data.asDouble(Transform[3]) ) { bResult = true; }
448 if( Get_MetaData_Item(Data, "YCELL") && Data.asDouble(Transform[5]) ) { bResult = true; }
449
450 return( bResult );
451 }
452
453 //---------------------------------------------------------
_Set_Transformation(void)454 bool CSG_GDAL_DataSet::_Set_Transformation(void)
455 {
456 if( !m_pDataSet )
457 {
458 return( false );
459 }
460
461 double Transform[6];
462
463 m_NX = GDALGetRasterXSize(m_pDataSet);
464 m_NY = GDALGetRasterYSize(m_pDataSet);
465
466 if( _Get_Transformation(Transform) == false )
467 {
468 m_bTransform = false;
469 m_Cellsize = 1.;
470 m_xMin = 0.;
471 m_yMin = 0.;
472 }
473 else if( Transform[1] == -Transform[5] && Transform[2] == 0. && Transform[4] == 0. ) // nothing to transform
474 {
475 m_bTransform = false;
476 m_Cellsize = Transform[1]; // pixel width (== pixel height)
477 m_xMin = Transform[0] + m_Cellsize * 0.5; // center (x) of left edge pixels
478 m_yMin = Transform[3] + m_Cellsize * (0.5 - m_NY); // center (y) of lower edge pixels
479 }
480 else
481 {
482 m_bTransform = true;
483 m_Cellsize = 1.;
484 m_xMin = 0.5;
485 m_yMin = 0.5;
486 }
487
488 m_TF_A[0] = Transform[0];
489 m_TF_A[1] = Transform[3];
490 m_TF_B[0][0] = Transform[1];
491 m_TF_B[0][1] = Transform[2];
492 m_TF_B[1][0] = Transform[4];
493 m_TF_B[1][1] = Transform[5];
494 m_TF_BInv = m_TF_B.Get_Inverse();
495
496 return( true );
497 }
498
499
500 ///////////////////////////////////////////////////////////
501 // //
502 ///////////////////////////////////////////////////////////
503
504 //---------------------------------------------------------
Open_Write(const CSG_String & File_Name,const CSG_String & Driver,const CSG_String & Options,TSG_Data_Type Type,int NBands,const CSG_Grid_System & System,const CSG_Projection & Projection)505 bool CSG_GDAL_DataSet::Open_Write(const CSG_String &File_Name, const CSG_String &Driver, const CSG_String &Options, TSG_Data_Type Type, int NBands, const CSG_Grid_System &System, const CSG_Projection &Projection)
506 {
507 Close();
508
509 //--------------------------------------------------------
510 GDALDriverH pDriver;
511
512 if( (pDriver = gSG_GDAL_Drivers.Get_Driver(Driver)) == NULL )
513 {
514 SG_UI_Msg_Add_Error(CSG_String::Format("%s: %s", _TL("driver not found."), Driver.c_str()));
515
516 return( false );
517 }
518
519 if( !CSG_GDAL_Drivers::has_Capability(pDriver, GDAL_DCAP_CREATE) )
520 {
521 SG_UI_Msg_Add_Error(_TL("Driver does not support file creation."));
522
523 return( false );
524 }
525
526 //--------------------------------------------------------
527 char **pOptions = Options.is_Empty() ? NULL : CSLTokenizeString2(Options, " ", CSLT_STRIPLEADSPACES);
528
529 if( !GDALValidateCreationOptions(pDriver, pOptions) )
530 {
531 SG_UI_Msg_Add_Error(CSG_String::Format("%s: %s", _TL("Creation option(s) not supported by the driver"), Options.c_str()));
532
533 CSLDestroy(pOptions);
534
535 return( false );
536 }
537
538 if( (m_pDataSet = GDALCreate(pDriver, File_Name, System.Get_NX(), System.Get_NY(), NBands, (GDALDataType)gSG_GDAL_Drivers.Get_GDAL_Type(Type), pOptions)) == NULL )
539 {
540 SG_UI_Msg_Add_Error(_TL("Could not create dataset."));
541
542 CSLDestroy(pOptions);
543
544 return( false );
545 }
546
547 CSLDestroy(pOptions);
548
549 //--------------------------------------------------------
550 m_File_Name = File_Name;
551
552 m_Access = SG_GDAL_IO_WRITE;
553
554 if( Projection.is_Okay() )
555 {
556 GDALSetProjection(m_pDataSet, Projection.Get_WKT());
557 }
558
559 double Transform[6] =
560 {
561 System.Get_XMin() - 0.5 * System.Get_Cellsize(), System.Get_Cellsize(), 0.,
562 System.Get_YMax() + 0.5 * System.Get_Cellsize(), 0., -System.Get_Cellsize()
563 };
564
565 GDALSetGeoTransform(m_pDataSet, Transform);
566
567 m_NX = GDALGetRasterXSize(m_pDataSet);
568 m_NY = GDALGetRasterYSize(m_pDataSet);
569
570 m_bTransform = false;
571 m_Cellsize = 1.;
572 m_xMin = 0.5;
573 m_yMin = 0.5;
574
575 return( true );
576 }
577
578 //---------------------------------------------------------
Close(void)579 bool CSG_GDAL_DataSet::Close(void)
580 {
581 if( m_pVrtSource )
582 {
583 GDALClose(m_pVrtSource); m_pVrtSource = NULL;
584
585 if( m_pDataSet )
586 {
587 // GDALClose(m_pDataSet); // this crashes in debug mode, gdal2.0dev!!!(???)
588 m_pDataSet = NULL;
589 }
590 }
591
592 if( m_pDataSet )
593 {
594 GDALClose(m_pDataSet); m_pDataSet = NULL;
595 }
596
597 m_File_Name.Clear();
598
599 m_Access = SG_GDAL_IO_CLOSED;
600
601 if( strlen(CPLGetLastErrorMsg()) > 3 )
602 {
603 SG_UI_Msg_Add_Error(CSG_String::Format("%s: %s", _TL("Dataset creation failed"), SG_STR_MBTOSG(CPLGetLastErrorMsg())));
604
605 CPLErrorReset();
606
607 return( false );
608 }
609
610 return( true );
611 }
612
613
614 ///////////////////////////////////////////////////////////
615 // //
616 ///////////////////////////////////////////////////////////
617
618 //---------------------------------------------------------
Get_Driver(void) const619 GDALDriverH CSG_GDAL_DataSet::Get_Driver(void) const
620 {
621 return( m_pDataSet ? GDALGetDatasetDriver(m_pDataSet) : NULL );
622 }
623
624 //---------------------------------------------------------
Get_DriverID(void) const625 CSG_String CSG_GDAL_DataSet::Get_DriverID(void) const
626 {
627 const char *s = GDALGetDescription(Get_Driver());
628
629 return( s ? s : "" );
630 }
631
632 //---------------------------------------------------------
Get_Projection(void) const633 const char * CSG_GDAL_DataSet::Get_Projection(void) const
634 {
635 const char *s = GDALGetProjectionRef(m_pDataSet);
636
637 return( s ? s : "" );
638 }
639
640 //---------------------------------------------------------
Get_Name(void) const641 CSG_String CSG_GDAL_DataSet::Get_Name(void) const
642 {
643 const char *s = GDALGetMetadataItem(m_pDataSet, GDAL_DMD_LONGNAME, 0);
644
645 return( s ? s : "" );
646 }
647
648 //---------------------------------------------------------
Get_Description(void) const649 CSG_String CSG_GDAL_DataSet::Get_Description(void) const
650 {
651 const char *s = GDALGetDescription(m_pDataSet);
652
653 return( s ? s : "" );
654 }
655
656 //---------------------------------------------------------
Get_File_Name(void) const657 CSG_String CSG_GDAL_DataSet::Get_File_Name(void) const
658 {
659 return( m_File_Name );
660 }
661
662 //---------------------------------------------------------
Get_MetaData_Item(const char * pszName,const char * pszDomain) const663 const char * CSG_GDAL_DataSet::Get_MetaData_Item(const char *pszName, const char *pszDomain) const
664 {
665 const char *s = GDALGetMetadataItem(m_pDataSet, pszName, pszDomain);
666
667 return( s ? s : "" );
668 }
669
670 //---------------------------------------------------------
Get_MetaData(const char * pszDomain) const671 const char ** CSG_GDAL_DataSet::Get_MetaData(const char *pszDomain) const
672 {
673 return( m_pDataSet ? (const char **)GDALGetMetadata(m_pDataSet, pszDomain) : NULL );
674 }
675
676 //---------------------------------------------------------
Get_MetaData_Item(CSG_String & MetaData,const char * pszName,const char * pszDomain) const677 bool CSG_GDAL_DataSet::Get_MetaData_Item(CSG_String &MetaData, const char *pszName, const char *pszDomain) const
678 {
679 const char *s = Get_MetaData_Item(pszName, pszDomain);
680
681 if( s && *s )
682 {
683 MetaData = s;
684
685 return( true );
686 }
687
688 return( false );
689 }
690
691 //---------------------------------------------------------
Get_MetaData_Domains(void) const692 CSG_Strings CSG_GDAL_DataSet::Get_MetaData_Domains(void) const
693 {
694 CSG_Strings Domains;
695
696 if( m_pDataSet && is_Reading() )
697 {
698 char **pDomains = GDALGetMetadataDomainList(m_pDataSet);
699
700 if( pDomains )
701 {
702 while( *pDomains )
703 {
704 Domains += *pDomains;
705
706 pDomains++;
707 }
708 }
709 }
710
711 return( Domains );
712 }
713
714 //---------------------------------------------------------
Get_MetaData(CSG_MetaData & MetaData) const715 bool CSG_GDAL_DataSet::Get_MetaData(CSG_MetaData &MetaData) const
716 {
717 if( m_pDataSet && is_Reading() )
718 {
719 char **pMetaData = GDALGetMetadata(m_pDataSet, 0);
720
721 if( pMetaData )
722 {
723 while( *pMetaData )
724 {
725 CSG_String s(*pMetaData);
726
727 CSG_String Key = s.BeforeFirst('='); Key.Replace("#", "."); // '#' in tag is not XML conform
728 CSG_String Value = s. AfterFirst('=');
729
730 MetaData.Add_Child(Key, Value);
731
732 pMetaData++;
733 }
734
735 return( true );
736 }
737 }
738
739 return( false );
740 }
741
742 //---------------------------------------------------------
Get_MetaData(CSG_MetaData & MetaData,const char * pszDomain) const743 bool CSG_GDAL_DataSet::Get_MetaData(CSG_MetaData &MetaData, const char *pszDomain) const
744 {
745 if( m_pDataSet && is_Reading() )
746 {
747 char **pMetaData = GDALGetMetadata(m_pDataSet, pszDomain);
748
749 if( pMetaData )
750 {
751 while( *pMetaData )
752 {
753 CSG_String s(*pMetaData);
754
755 CSG_String Key = s.BeforeFirst('='); Key.Replace("#", "."); // '#' in tag is not XML conform
756 CSG_String Value = s. AfterFirst('=');
757
758 MetaData.Add_Child(Key, Value);
759
760 pMetaData++;
761 }
762
763 return( true );
764 }
765 }
766
767 return( false );
768 }
769
770
771 ///////////////////////////////////////////////////////////
772 // //
773 ///////////////////////////////////////////////////////////
774
775 //---------------------------------------------------------
Get_Count(void) const776 int CSG_GDAL_DataSet::Get_Count(void) const
777 {
778 return( m_pDataSet ? GDALGetRasterCount(m_pDataSet) : 0 );
779 }
780
781 //---------------------------------------------------------
Get_Name(int i) const782 CSG_String CSG_GDAL_DataSet::Get_Name(int i) const
783 {
784 CSG_String Name;
785
786 GDALRasterBandH pBand;
787
788 if( is_Reading() && (pBand = GDALGetRasterBand(m_pDataSet, i + 1)) != NULL )
789 {
790 const char *s;
791
792 //-------------------------------------------------
793 if( !Get_DriverID().Cmp("GRIB") )
794 {
795 if( (s = GDALGetMetadataItem(pBand, "GRIB_COMMENT", 0)) != NULL && *s )
796 {
797 Name = s; CSG_DateTime d;
798
799 if( (s = GDALGetMetadataItem(pBand, "GRIB_ELEMENT" , 0)) != NULL && *s ) { Name += "["; Name += s; Name += "]"; }
800 if( (s = GDALGetMetadataItem(pBand, "GRIB_SHORT_NAME", 0)) != NULL && *s ) { Name += "["; Name += s; Name += "]"; }
801 if( (s = GDALGetMetadataItem(pBand, "GRIB_VALID_TIME", 0)) != NULL && *s ) { d.Set_Unix_Time(atoi(s)); Name += "[" + d.Format_ISOCombined() + "]"; }
802 }
803 }
804
805 //-------------------------------------------------
806 if( !Get_DriverID().Cmp("netCDF") )
807 {
808 if( (s = GDALGetMetadataItem(pBand, "NETCDF_VARNAME" , 0)) != NULL && *s ) { Name += "["; Name += s; Name += "]"; }
809 if( (s = GDALGetMetadataItem(pBand, "NETCDF_DIMENSION_time" , 0)) != NULL && *s ) { Name += "["; Name += s; Name += "]"; }
810 if( (s = GDALGetMetadataItem(pBand, "NETCDF_DIMENSION_level", 0)) != NULL && *s ) { Name += "["; Name += s; Name += "]"; }
811 }
812
813 //-------------------------------------------------
814 if( Name.is_Empty() )
815 {
816 Name = Get_Name();
817
818 if( Name.is_Empty() )
819 {
820 Name = _TL("Band");
821 }
822
823 Name += CSG_String::Format(" %0*d", SG_Get_Digit_Count(Get_Count() + 1), i + 1);
824 }
825 }
826
827 return( Name );
828 }
829
830 //---------------------------------------------------------
Get_Description(int i) const831 CSG_String CSG_GDAL_DataSet::Get_Description(int i) const
832 {
833 CSG_String Description;
834
835 GDALRasterBandH pBand;
836
837 if( is_Reading() && (pBand = GDALGetRasterBand(m_pDataSet, i + 1)) != NULL )
838 {
839 char **pMetaData = GDALGetMetadata(pBand, 0);
840
841 if( pMetaData )
842 {
843 while( *pMetaData )
844 {
845 CSG_String s(*pMetaData);
846
847 Description += s + "\n";
848
849 pMetaData++;
850 }
851 }
852 }
853
854 return( Description );
855 }
856
857 //---------------------------------------------------------
Get_MetaData(int i,CSG_MetaData & MetaData) const858 bool CSG_GDAL_DataSet::Get_MetaData(int i, CSG_MetaData &MetaData) const
859 {
860 GDALRasterBandH pBand;
861
862 if( is_Reading() && (pBand = GDALGetRasterBand(m_pDataSet, i + 1)) != NULL )
863 {
864 char **pMetaData = GDALGetMetadata(pBand, 0);
865
866 if( pMetaData )
867 {
868 while( *pMetaData )
869 {
870 CSG_String s(*pMetaData);
871
872 MetaData.Add_Child(s.BeforeFirst('='), s.AfterFirst('='));
873
874 pMetaData++;
875 }
876
877 return( true );
878 }
879 }
880
881 return( false );
882 }
883
884 //---------------------------------------------------------
Get_MetaData_Item(int i,const char * pszName) const885 const char * CSG_GDAL_DataSet::Get_MetaData_Item(int i, const char *pszName) const
886 {
887 GDALRasterBandH pBand = GDALGetRasterBand(m_pDataSet, i + 1);
888
889 return( pBand ? GDALGetMetadataItem(pBand, pszName, 0) : "" );
890 }
891
Get_MetaData_Item(int i,const char * pszName,CSG_String & MetaData) const892 bool CSG_GDAL_DataSet::Get_MetaData_Item(int i, const char *pszName, CSG_String &MetaData) const
893 {
894 GDALRasterBandH pBand = GDALGetRasterBand(m_pDataSet, i + 1);
895
896 if( pBand != NULL )
897 {
898 const char *pMetaData = GDALGetMetadataItem(pBand, pszName, 0);
899
900 if( pMetaData && *pMetaData )
901 {
902 MetaData = pMetaData;
903
904 return( true );
905 }
906 }
907
908 return( false );
909 }
910
911 //---------------------------------------------------------
Read(int i)912 CSG_Grid * CSG_GDAL_DataSet::Read(int i)
913 {
914 if( !is_Reading() )
915 {
916 return( NULL );
917 }
918
919 //-------------------------------------------------
920 GDALRasterBandH pBand = GDALGetRasterBand(m_pDataSet, i + 1);
921
922 if( !pBand )
923 {
924 return( NULL );
925 }
926
927 //-------------------------------------------------
928 TSG_Data_Type Type = gSG_GDAL_Drivers.Get_SAGA_Type(GDALGetRasterDataType(pBand));
929
930 CSG_Grid *pGrid = SG_Create_Grid(Type, Get_NX(), Get_NY(), Get_Cellsize(), Get_xMin(), Get_yMin());
931
932 if( !pGrid )
933 {
934 return( NULL );
935 }
936
937 //-------------------------------------------------
938 int bSuccess;
939
940 double zScale = GDALGetRasterScale (pBand, &bSuccess); if( !bSuccess || !zScale ) zScale = 1.;
941 double zOffset = GDALGetRasterOffset(pBand, &bSuccess); if( !bSuccess ) zOffset = 0.;
942
943 pGrid->Set_Name (Get_Name (i));
944 pGrid->Set_Description (Get_Description(i));
945 pGrid->Set_Unit (CSG_String(GDALGetRasterUnitType(pBand)));
946 pGrid->Set_Scaling (zScale, zOffset);
947
948 //-------------------------------------------------
949 OGRSpatialReferenceH SRef = OSRNewSpatialReference(Get_Projection());
950
951 char *Proj4 = NULL;
952
953 if( OSRExportToProj4(SRef, &Proj4) == OGRERR_NONE )
954 {
955 pGrid->Get_Projection().Create(Get_Projection(), Proj4);
956
957 CPLFree(Proj4);
958 }
959 else
960 {
961 pGrid->Get_Projection().Create(Get_Projection());
962 }
963
964 CPLFree(SRef);
965
966 if( pGrid->Get_Projection().is_Okay() == false )
967 {
968 CSG_String Data; int EPSG;
969
970 if( !Get_MetaData_Item(Data, "EPSG") || !Data.asInt(EPSG) || !pGrid->Get_Projection().Create(EPSG) )
971 {
972 if( Get_MetaData_Item(Data, "proj4_string") )
973 {
974 pGrid->Get_Projection().Create(Data, SG_PROJ_FMT_Proj4);
975 }
976 }
977 }
978
979 //-------------------------------------------------
980 CSG_MetaData &MetaData = pGrid->Get_MetaData();
981
982 MetaData.Add_Child("GDAL_DRIVER", Get_DriverID());
983
984 Get_MetaData(MetaData);
985
986 Get_MetaData(i, *MetaData.Add_Child("Band"));
987
988 //-------------------------------------------------
989 double zNoData = GDALGetRasterNoDataValue(pBand, &bSuccess);
990
991 if( bSuccess )
992 {
993 switch( Type )
994 {
995 default : pGrid->Set_NoData_Value( (int)zNoData); break;
996 case SG_DATATYPE_Float : pGrid->Set_NoData_Value( (float)zNoData); break;
997 case SG_DATATYPE_Double: pGrid->Set_NoData_Value((double)zNoData); break;
998 }
999 }
1000
1001 //-------------------------------------------------
1002 void *zLine; GDALDataType zType;
1003
1004 switch( Type )
1005 {
1006 default : zLine = SG_Malloc(Get_NX() * sizeof( int)); zType = GDT_Int32 ; break;
1007 case SG_DATATYPE_Float : zLine = SG_Malloc(Get_NX() * sizeof( float)); zType = GDT_Float32; break;
1008 case SG_DATATYPE_Double: zLine = SG_Malloc(Get_NX() * sizeof(double)); zType = GDT_Float64; break;
1009 }
1010
1011 for(int y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
1012 {
1013 int yy = m_bTransform ? y : Get_NY() - 1 - y;
1014
1015 if( GDALRasterIO(pBand, GF_Read, 0, y, Get_NX(), 1, zLine, Get_NX(), 1, zType, 0, 0) == CE_None )
1016 {
1017 for(int x=0; x<Get_NX(); x++)
1018 {
1019 switch( Type )
1020 {
1021 default : pGrid->Set_Value(x, yy, ((int *)zLine)[x], false); break;
1022 case SG_DATATYPE_Float : pGrid->Set_Value(x, yy, ((float *)zLine)[x], false); break;
1023 case SG_DATATYPE_Double: pGrid->Set_Value(x, yy, ((double *)zLine)[x], false); break;
1024 }
1025 }
1026 }
1027 }
1028
1029 SG_Free(zLine);
1030
1031 return( pGrid );
1032 }
1033
1034 //---------------------------------------------------------
Write(int i,CSG_Grid * pGrid,double noDataValue)1035 bool CSG_GDAL_DataSet::Write(int i, CSG_Grid *pGrid, double noDataValue)
1036 {
1037 if( !m_pDataSet || !pGrid || pGrid->Get_NX() != Get_NX() || pGrid->Get_NY() != Get_NY() || i < 0 || i >= Get_Count() )
1038 {
1039 return( false );
1040 }
1041
1042 GDALRasterBandH pBand = GDALGetRasterBand(m_pDataSet, i + 1);
1043
1044 //-----------------------------------------------------
1045 CPLErr Error = CE_None;
1046
1047 double *zLine = (double *)SG_Malloc(Get_NX() * sizeof(double));
1048
1049 for(int y=0, yy=Get_NY()-1; Error==CE_None && y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++, yy--)
1050 {
1051 for(int x=0; x<Get_NX(); x++)
1052 {
1053 zLine[x] = pGrid->is_NoData(x, yy) ? noDataValue : pGrid->asDouble(x, yy);
1054 }
1055
1056 Error = GDALRasterIO(pBand, GF_Write, 0, y, Get_NX(), 1, zLine, Get_NX(), 1, GDT_Float64, 0, 0);
1057 }
1058
1059 SG_Free(zLine);
1060
1061 //-----------------------------------------------------
1062 if( Error != CE_None )
1063 {
1064 SG_UI_Msg_Add_Error(CSG_String::Format("%s", _TL("Writing dataset failed.")));
1065
1066 return( false );
1067 }
1068
1069 //-----------------------------------------------------
1070 GDALSetRasterNoDataValue(pBand, noDataValue);
1071 GDALSetRasterStatistics (pBand, pGrid->Get_Min(), pGrid->Get_Max(), pGrid->Get_Mean(), pGrid->Get_StdDev());
1072
1073 return( true );
1074 }
1075
1076 //---------------------------------------------------------
Write(int i,CSG_Grid * pGrid)1077 bool CSG_GDAL_DataSet::Write(int i, CSG_Grid *pGrid)
1078 {
1079 return (CSG_GDAL_DataSet::Write (i, pGrid, pGrid->Get_NoData_Value()));
1080 }
1081
1082
1083 ///////////////////////////////////////////////////////////
1084 // //
1085 ///////////////////////////////////////////////////////////
1086
1087 //---------------------------------------------------------
Get_Extent(bool bTransform) const1088 CSG_Rect CSG_GDAL_DataSet::Get_Extent(bool bTransform) const
1089 {
1090 if( bTransform && Needs_Transformation() )
1091 {
1092 CSG_Grid_System System;
1093
1094 if( Get_Transformation(System, false) )
1095 {
1096 return( System.Get_Extent() );
1097 }
1098 }
1099
1100 return( CSG_Rect(Get_xMin(), Get_yMin(), Get_xMax(), Get_yMax()) );
1101 }
1102
1103 //---------------------------------------------------------
Get_System(void) const1104 CSG_Grid_System CSG_GDAL_DataSet::Get_System(void) const
1105 {
1106 CSG_Grid_System System;
1107
1108 if( !Needs_Transformation() || !Get_Transformation(System, false) )
1109 {
1110 System.Assign(Get_Cellsize(), Get_xMin(), Get_yMin(), Get_NX(), Get_NY());
1111 }
1112
1113 return( System );
1114 }
1115
1116
1117 ///////////////////////////////////////////////////////////
1118 // //
1119 ///////////////////////////////////////////////////////////
1120
1121 //---------------------------------------------------------
Get_Transformation(CSG_Grid_System & System,bool bVerbose) const1122 bool CSG_GDAL_DataSet::Get_Transformation(CSG_Grid_System &System, bool bVerbose) const
1123 {
1124 CSG_Vector A;
1125 CSG_Matrix B;
1126
1127 Get_Transformation(A, B);
1128
1129 //-----------------------------------------------------
1130 if( Needs_Transformation() )
1131 {
1132 CSG_Vector v(2);
1133 CSG_Rect r;
1134
1135 v[0] = Get_xMin(); v[1] = Get_yMin(); v = B * v + A; r.Assign(v[0], v[1], v[0], v[1]);
1136 v[0] = Get_xMin(); v[1] = Get_yMax(); v = B * v + A; r.Union(CSG_Point(v[0], v[1]));
1137 v[0] = Get_xMax(); v[1] = Get_yMax(); v = B * v + A; r.Union(CSG_Point(v[0], v[1]));
1138 v[0] = Get_xMax(); v[1] = Get_yMin(); v = B * v + A; r.Union(CSG_Point(v[0], v[1]));
1139
1140 v[0] = 1; v[1] = 0; v = B * v; double dx = v.Get_Length();
1141 v[0] = 0; v[1] = 1; v = B * v; double dy = v.Get_Length();
1142
1143 if( dx != dy )
1144 {
1145 if( bVerbose )
1146 {
1147 SG_UI_Msg_Add_Execution(CSG_String::Format("\n%s: %s\n\t%s: %f",
1148 _TL("warning"), _TL("top-to-bottom and left-to-right cell sizes differ."), _TL("Difference"), fabs(dy - dx)), false
1149 );
1150 }
1151
1152 if( dx > dy )
1153 {
1154 dx = dy;
1155 }
1156
1157 if( bVerbose )
1158 {
1159 SG_UI_Msg_Add_Execution(CSG_String::Format("\n\t%s: %f\n", _TL("using cellsize"), dx), false);
1160 }
1161 }
1162
1163 return( System.Assign(dx, r) );
1164 }
1165
1166 //-----------------------------------------------------
1167 return( false );
1168 }
1169
1170 //---------------------------------------------------------
Get_Transformation(CSG_Grid ** ppGrid,TSG_Grid_Resampling Interpolation,bool bVerbose) const1171 bool CSG_GDAL_DataSet::Get_Transformation(CSG_Grid **ppGrid, TSG_Grid_Resampling Interpolation, bool bVerbose) const
1172 {
1173 CSG_Grid_System System;
1174
1175 if( Get_Transformation(System, bVerbose) )
1176 {
1177 return( Get_Transformation(ppGrid, Interpolation, System, bVerbose) );
1178 }
1179
1180 return( false );
1181 }
1182
1183 //---------------------------------------------------------
Get_Transformation(CSG_Grid ** ppGrid,TSG_Grid_Resampling Interpolation,const CSG_Grid_System & System,bool bVerbose) const1184 bool CSG_GDAL_DataSet::Get_Transformation(CSG_Grid **ppGrid, TSG_Grid_Resampling Interpolation, const CSG_Grid_System &System, bool bVerbose) const
1185 {
1186 if( !System.is_Valid() )
1187 {
1188 return( false );
1189 }
1190
1191 //-----------------------------------------------------
1192 CSG_Vector A;
1193 CSG_Matrix B, BInv;
1194
1195 Get_Transformation(A, B);
1196
1197 BInv = B.Get_Inverse();
1198
1199 //-----------------------------------------------------
1200 CSG_Grid *pImage = *ppGrid;
1201 CSG_Grid *pWorld = SG_Create_Grid(System, pImage->Get_Type());
1202
1203 if( !pWorld )
1204 {
1205 return( false );
1206 }
1207
1208 *ppGrid = pWorld;
1209
1210 pWorld->Set_Name (pImage->Get_Name ());
1211 pWorld->Set_Description (pImage->Get_Description ());
1212 pWorld->Set_Unit (pImage->Get_Unit ());
1213 pWorld->Set_Scaling (pImage->Get_Scaling (), pImage->Get_Offset());
1214 pWorld->Set_NoData_Value_Range(pImage->Get_NoData_Value(), pImage->Get_NoData_Value(true));
1215 pWorld->Get_MetaData() = pImage->Get_MetaData ();
1216 pWorld->Get_Projection() = pImage->Get_Projection ();
1217
1218 //-----------------------------------------------------
1219 // #pragma omp parallel for
1220 // for(int y=0; y<pWorld->Get_NY(); y++)
1221 // {
1222 // Process_Get_Okay();
1223
1224 for(int y=0; y<pWorld->Get_NY() && SG_UI_Process_Set_Progress(y, pWorld->Get_NY()); y++)
1225 {
1226 #pragma omp parallel for
1227 for(int x=0; x<pWorld->Get_NX(); x++)
1228 {
1229 double z;
1230 CSG_Vector vWorld(2), vImage;
1231
1232 vWorld[0] = pWorld->Get_XMin() + x * pWorld->Get_Cellsize();
1233 vWorld[1] = pWorld->Get_YMin() + y * pWorld->Get_Cellsize();
1234
1235 vImage = BInv * (vWorld - A);
1236
1237 if( pImage->Get_Value(vImage[0], vImage[1], z, Interpolation, true) )
1238 {
1239 pWorld->Set_Value(x, y, z);
1240 }
1241 else
1242 {
1243 pWorld->Set_NoData(x, y);
1244 }
1245 }
1246 }
1247
1248 delete(pImage);
1249
1250 return( true );
1251 }
1252
1253
1254 ///////////////////////////////////////////////////////////
1255 // //
1256 // //
1257 // //
1258 ///////////////////////////////////////////////////////////
1259
1260 //---------------------------------------------------------
SG_Get_Grid_Type(CSG_Parameter_Grid_List * pGrids)1261 TSG_Data_Type SG_Get_Grid_Type (CSG_Parameter_Grid_List *pGrids)
1262 {
1263 TSG_Data_Type Type = SG_DATATYPE_Byte;
1264
1265 if( pGrids )
1266 {
1267 for(int i=0; i<pGrids->Get_Grid_Count(); i++)
1268 {
1269 if( SG_Data_Type_Get_Size(Type) <= SG_Data_Type_Get_Size(pGrids->Get_Grid(i)->Get_Type()) )
1270 {
1271 Type = pGrids->Get_Grid(i)->Get_Type();
1272 }
1273 }
1274 }
1275
1276 return( Type );
1277 }
1278
1279
1280 ///////////////////////////////////////////////////////////
1281 // //
1282 // //
1283 // //
1284 ///////////////////////////////////////////////////////////
1285
1286 //---------------------------------------------------------
1287