1 /******************************************************************************
2 * $Id: gdal_misc.cpp 29326 2015-06-10 20:36:31Z rouault $
3 *
4 * Project: GDAL Core
5 * Purpose: Free standing functions for GDAL.
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 *
8 ******************************************************************************
9 * Copyright (c) 1999, Frank Warmerdam
10 * Copyright (c) 2007-2013, Even Rouault <even dot rouault at mines-paris dot org>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ****************************************************************************/
30
31 #include "gdal_priv.h"
32 #include "cpl_string.h"
33 #include "cpl_minixml.h"
34 #include "cpl_multiproc.h"
35 #include <ctype.h>
36 #include <string>
37
38 CPL_CVSID("$Id: gdal_misc.cpp 29326 2015-06-10 20:36:31Z rouault $");
39
40 #include "ogr_spatialref.h"
41 #include "gdal_mdreader.h"
42
43 /************************************************************************/
44 /* __pure_virtual() */
45 /* */
46 /* The following is a gross hack to remove the last remaining */
47 /* dependency on the GNU C++ standard library. */
48 /************************************************************************/
49
50 #ifdef __GNUC__
51
52 extern "C"
__pure_virtual()53 void __pure_virtual()
54
55 {
56 }
57
58 #endif
59
60 /************************************************************************/
61 /* GDALDataTypeUnion() */
62 /************************************************************************/
63
64 /**
65 * \brief Return the smallest data type that can fully express both input data
66 * types.
67 *
68 * @param eType1 first data type.
69 * @param eType2 second data type.
70 *
71 * @return a data type able to express eType1 and eType2.
72 */
73
74 GDALDataType CPL_STDCALL
GDALDataTypeUnion(GDALDataType eType1,GDALDataType eType2)75 GDALDataTypeUnion( GDALDataType eType1, GDALDataType eType2 )
76
77 {
78 int bFloating, bComplex, nBits, bSigned;
79
80 bComplex = GDALDataTypeIsComplex(eType1) | GDALDataTypeIsComplex(eType2);
81
82 switch( eType1 )
83 {
84 case GDT_Byte:
85 nBits = 8;
86 bSigned = FALSE;
87 bFloating = FALSE;
88 break;
89
90 case GDT_Int16:
91 case GDT_CInt16:
92 nBits = 16;
93 bSigned = TRUE;
94 bFloating = FALSE;
95 break;
96
97 case GDT_UInt16:
98 nBits = 16;
99 bSigned = FALSE;
100 bFloating = FALSE;
101 break;
102
103 case GDT_Int32:
104 case GDT_CInt32:
105 nBits = 32;
106 bSigned = TRUE;
107 bFloating = FALSE;
108 break;
109
110 case GDT_UInt32:
111 nBits = 32;
112 bSigned = FALSE;
113 bFloating = FALSE;
114 break;
115
116 case GDT_Float32:
117 case GDT_CFloat32:
118 nBits = 32;
119 bSigned = TRUE;
120 bFloating = TRUE;
121 break;
122
123 case GDT_Float64:
124 case GDT_CFloat64:
125 nBits = 64;
126 bSigned = TRUE;
127 bFloating = TRUE;
128 break;
129
130 default:
131 CPLAssert( FALSE );
132 return GDT_Unknown;
133 }
134
135 switch( eType2 )
136 {
137 case GDT_Byte:
138 break;
139
140 case GDT_Int16:
141 case GDT_CInt16:
142 nBits = MAX(nBits,16);
143 bSigned = TRUE;
144 break;
145
146 case GDT_UInt16:
147 nBits = MAX(nBits,16);
148 break;
149
150 case GDT_Int32:
151 case GDT_CInt32:
152 nBits = MAX(nBits,32);
153 bSigned = TRUE;
154 break;
155
156 case GDT_UInt32:
157 nBits = MAX(nBits,32);
158 break;
159
160 case GDT_Float32:
161 case GDT_CFloat32:
162 nBits = MAX(nBits,32);
163 bSigned = TRUE;
164 bFloating = TRUE;
165 break;
166
167 case GDT_Float64:
168 case GDT_CFloat64:
169 nBits = MAX(nBits,64);
170 bSigned = TRUE;
171 bFloating = TRUE;
172 break;
173
174 default:
175 CPLAssert( FALSE );
176 return GDT_Unknown;
177 }
178
179 if( nBits == 8 )
180 return GDT_Byte;
181 else if( nBits == 16 && bComplex )
182 return GDT_CInt16;
183 else if( nBits == 16 && bSigned )
184 return GDT_Int16;
185 else if( nBits == 16 && !bSigned )
186 return GDT_UInt16;
187 else if( nBits == 32 && bFloating && bComplex )
188 return GDT_CFloat32;
189 else if( nBits == 32 && bFloating )
190 return GDT_Float32;
191 else if( nBits == 32 && bComplex )
192 return GDT_CInt32;
193 else if( nBits == 32 && bSigned )
194 return GDT_Int32;
195 else if( nBits == 32 && !bSigned )
196 return GDT_UInt32;
197 else if( nBits == 64 && bComplex )
198 return GDT_CFloat64;
199 else
200 return GDT_Float64;
201 }
202
203
204 /************************************************************************/
205 /* GDALGetDataTypeSize() */
206 /************************************************************************/
207
208 /**
209 * \brief Get data type size in bits.
210 *
211 * Returns the size of a a GDT_* type in bits, <b>not bytes</b>!
212 *
213 * @param eDataType type, such as GDT_Byte.
214 * @return the number of bits or zero if it is not recognised.
215 */
216
GDALGetDataTypeSize(GDALDataType eDataType)217 int CPL_STDCALL GDALGetDataTypeSize( GDALDataType eDataType )
218
219 {
220 switch( eDataType )
221 {
222 case GDT_Byte:
223 return 8;
224
225 case GDT_UInt16:
226 case GDT_Int16:
227 return 16;
228
229 case GDT_UInt32:
230 case GDT_Int32:
231 case GDT_Float32:
232 case GDT_CInt16:
233 return 32;
234
235 case GDT_Float64:
236 case GDT_CInt32:
237 case GDT_CFloat32:
238 return 64;
239
240 case GDT_CFloat64:
241 return 128;
242
243 default:
244 return 0;
245 }
246 }
247
248 /************************************************************************/
249 /* GDALDataTypeIsComplex() */
250 /************************************************************************/
251
252 /**
253 * \brief Is data type complex?
254 *
255 * @return TRUE if the passed type is complex (one of GDT_CInt16, GDT_CInt32,
256 * GDT_CFloat32 or GDT_CFloat64), that is it consists of a real and imaginary
257 * component.
258 */
259
GDALDataTypeIsComplex(GDALDataType eDataType)260 int CPL_STDCALL GDALDataTypeIsComplex( GDALDataType eDataType )
261
262 {
263 switch( eDataType )
264 {
265 case GDT_CInt16:
266 case GDT_CInt32:
267 case GDT_CFloat32:
268 case GDT_CFloat64:
269 return TRUE;
270
271 default:
272 return FALSE;
273 }
274 }
275
276 /************************************************************************/
277 /* GDALGetDataTypeName() */
278 /************************************************************************/
279
280 /**
281 * \brief Get name of data type.
282 *
283 * Returns a symbolic name for the data type. This is essentially the
284 * the enumerated item name with the GDT_ prefix removed. So GDT_Byte returns
285 * "Byte". The returned strings are static strings and should not be modified
286 * or freed by the application. These strings are useful for reporting
287 * datatypes in debug statements, errors and other user output.
288 *
289 * @param eDataType type to get name of.
290 * @return string corresponding to existing data type
291 * or NULL pointer if invalid type given.
292 */
293
GDALGetDataTypeName(GDALDataType eDataType)294 const char * CPL_STDCALL GDALGetDataTypeName( GDALDataType eDataType )
295
296 {
297 switch( eDataType )
298 {
299 case GDT_Unknown:
300 return "Unknown";
301
302 case GDT_Byte:
303 return "Byte";
304
305 case GDT_UInt16:
306 return "UInt16";
307
308 case GDT_Int16:
309 return "Int16";
310
311 case GDT_UInt32:
312 return "UInt32";
313
314 case GDT_Int32:
315 return "Int32";
316
317 case GDT_Float32:
318 return "Float32";
319
320 case GDT_Float64:
321 return "Float64";
322
323 case GDT_CInt16:
324 return "CInt16";
325
326 case GDT_CInt32:
327 return "CInt32";
328
329 case GDT_CFloat32:
330 return "CFloat32";
331
332 case GDT_CFloat64:
333 return "CFloat64";
334
335 default:
336 return NULL;
337 }
338 }
339
340 /************************************************************************/
341 /* GDALGetDataTypeByName() */
342 /************************************************************************/
343
344 /**
345 * \brief Get data type by symbolic name.
346 *
347 * Returns a data type corresponding to the given symbolic name. This
348 * function is opposite to the GDALGetDataTypeName().
349 *
350 * @param pszName string containing the symbolic name of the type.
351 *
352 * @return GDAL data type.
353 */
354
GDALGetDataTypeByName(const char * pszName)355 GDALDataType CPL_STDCALL GDALGetDataTypeByName( const char *pszName )
356
357 {
358 VALIDATE_POINTER1( pszName, "GDALGetDataTypeByName", GDT_Unknown );
359
360 int iType;
361
362 for( iType = 1; iType < GDT_TypeCount; iType++ )
363 {
364 if( GDALGetDataTypeName((GDALDataType)iType) != NULL
365 && EQUAL(GDALGetDataTypeName((GDALDataType)iType), pszName) )
366 {
367 return (GDALDataType)iType;
368 }
369 }
370
371 return GDT_Unknown;
372 }
373
374 /************************************************************************/
375 /* GDALGetAsyncStatusTypeByName() */
376 /************************************************************************/
377 /**
378 * Get AsyncStatusType by symbolic name.
379 *
380 * Returns a data type corresponding to the given symbolic name. This
381 * function is opposite to the GDALGetAsyncStatusTypeName().
382 *
383 * @param pszName string containing the symbolic name of the type.
384 *
385 * @return GDAL AsyncStatus type.
386 */
GDALGetAsyncStatusTypeByName(const char * pszName)387 GDALAsyncStatusType CPL_DLL CPL_STDCALL GDALGetAsyncStatusTypeByName( const char *pszName )
388 {
389 VALIDATE_POINTER1( pszName, "GDALGetAsyncStatusTypeByName", GARIO_ERROR);
390
391 int iType;
392
393 for( iType = 1; iType < GARIO_TypeCount; iType++ )
394 {
395 if( GDALGetAsyncStatusTypeName((GDALAsyncStatusType)iType) != NULL
396 && EQUAL(GDALGetAsyncStatusTypeName((GDALAsyncStatusType)iType), pszName) )
397 {
398 return (GDALAsyncStatusType)iType;
399 }
400 }
401
402 return GARIO_ERROR;
403 }
404
405
406 /************************************************************************/
407 /* GDALGetAsyncStatusTypeName() */
408 /************************************************************************/
409
410 /**
411 * Get name of AsyncStatus data type.
412 *
413 * Returns a symbolic name for the AsyncStatus data type. This is essentially the
414 * the enumerated item name with the GARIO_ prefix removed. So GARIO_COMPLETE returns
415 * "COMPLETE". The returned strings are static strings and should not be modified
416 * or freed by the application. These strings are useful for reporting
417 * datatypes in debug statements, errors and other user output.
418 *
419 * @param eAsyncStatusType type to get name of.
420 * @return string corresponding to type.
421 */
422
GDALGetAsyncStatusTypeName(GDALAsyncStatusType eAsyncStatusType)423 const char * CPL_STDCALL GDALGetAsyncStatusTypeName( GDALAsyncStatusType eAsyncStatusType )
424
425 {
426 switch( eAsyncStatusType )
427 {
428 case GARIO_PENDING:
429 return "PENDING";
430
431 case GARIO_UPDATE:
432 return "UPDATE";
433
434 case GARIO_ERROR:
435 return "ERROR";
436
437 case GARIO_COMPLETE:
438 return "COMPLETE";
439 default:
440 return NULL;
441 }
442 }
443
444 /************************************************************************/
445 /* GDALGetPaletteInterpretationName() */
446 /************************************************************************/
447
448 /**
449 * \brief Get name of palette interpretation
450 *
451 * Returns a symbolic name for the palette interpretation. This is the
452 * the enumerated item name with the GPI_ prefix removed. So GPI_Gray returns
453 * "Gray". The returned strings are static strings and should not be modified
454 * or freed by the application.
455 *
456 * @param eInterp palette interpretation to get name of.
457 * @return string corresponding to palette interpretation.
458 */
459
GDALGetPaletteInterpretationName(GDALPaletteInterp eInterp)460 const char *GDALGetPaletteInterpretationName( GDALPaletteInterp eInterp )
461
462 {
463 switch( eInterp )
464 {
465 case GPI_Gray:
466 return "Gray";
467
468 case GPI_RGB:
469 return "RGB";
470
471 case GPI_CMYK:
472 return "CMYK";
473
474 case GPI_HLS:
475 return "HLS";
476
477 default:
478 return "Unknown";
479 }
480 }
481
482 /************************************************************************/
483 /* GDALGetColorInterpretationName() */
484 /************************************************************************/
485
486 /**
487 * \brief Get name of color interpretation
488 *
489 * Returns a symbolic name for the color interpretation. This is derived from
490 * the enumerated item name with the GCI_ prefix removed, but there are some
491 * variations. So GCI_GrayIndex returns "Gray" and GCI_RedBand returns "Red".
492 * The returned strings are static strings and should not be modified
493 * or freed by the application.
494 *
495 * @param eInterp color interpretation to get name of.
496 * @return string corresponding to color interpretation
497 * or NULL pointer if invalid enumerator given.
498 */
499
GDALGetColorInterpretationName(GDALColorInterp eInterp)500 const char *GDALGetColorInterpretationName( GDALColorInterp eInterp )
501
502 {
503 switch( eInterp )
504 {
505 case GCI_Undefined:
506 return "Undefined";
507
508 case GCI_GrayIndex:
509 return "Gray";
510
511 case GCI_PaletteIndex:
512 return "Palette";
513
514 case GCI_RedBand:
515 return "Red";
516
517 case GCI_GreenBand:
518 return "Green";
519
520 case GCI_BlueBand:
521 return "Blue";
522
523 case GCI_AlphaBand:
524 return "Alpha";
525
526 case GCI_HueBand:
527 return "Hue";
528
529 case GCI_SaturationBand:
530 return "Saturation";
531
532 case GCI_LightnessBand:
533 return "Lightness";
534
535 case GCI_CyanBand:
536 return "Cyan";
537
538 case GCI_MagentaBand:
539 return "Magenta";
540
541 case GCI_YellowBand:
542 return "Yellow";
543
544 case GCI_BlackBand:
545 return "Black";
546
547 case GCI_YCbCr_YBand:
548 return "YCbCr_Y";
549
550 case GCI_YCbCr_CbBand:
551 return "YCbCr_Cb";
552
553 case GCI_YCbCr_CrBand:
554 return "YCbCr_Cr";
555
556 default:
557 return "Unknown";
558 }
559 }
560
561 /************************************************************************/
562 /* GDALGetColorInterpretationByName() */
563 /************************************************************************/
564
565 /**
566 * \brief Get color interpreation by symbolic name.
567 *
568 * Returns a color interpreation corresponding to the given symbolic name. This
569 * function is opposite to the GDALGetColorInterpretationName().
570 *
571 * @param pszName string containing the symbolic name of the color interpretation.
572 *
573 * @return GDAL color interpretation.
574 *
575 * @since GDAL 1.7.0
576 */
577
GDALGetColorInterpretationByName(const char * pszName)578 GDALColorInterp GDALGetColorInterpretationByName( const char *pszName )
579
580 {
581 VALIDATE_POINTER1( pszName, "GDALGetColorInterpretationByName", GCI_Undefined );
582
583 int iType;
584
585 for( iType = 0; iType <= GCI_Max; iType++ )
586 {
587 if( EQUAL(GDALGetColorInterpretationName((GDALColorInterp)iType), pszName) )
588 {
589 return (GDALColorInterp)iType;
590 }
591 }
592
593 return GCI_Undefined;
594 }
595
596 /************************************************************************/
597 /* GDALGetRandomRasterSample() */
598 /************************************************************************/
599
600 int CPL_STDCALL
GDALGetRandomRasterSample(GDALRasterBandH hBand,int nSamples,float * pafSampleBuf)601 GDALGetRandomRasterSample( GDALRasterBandH hBand, int nSamples,
602 float *pafSampleBuf )
603
604 {
605 VALIDATE_POINTER1( hBand, "GDALGetRandomRasterSample", 0 );
606
607 GDALRasterBand *poBand;
608
609 poBand = (GDALRasterBand *) GDALGetRasterSampleOverview( hBand, nSamples );
610 CPLAssert( NULL != poBand );
611
612 /* -------------------------------------------------------------------- */
613 /* Figure out the ratio of blocks we will read to get an */
614 /* approximate value. */
615 /* -------------------------------------------------------------------- */
616 int nBlockXSize, nBlockYSize;
617 int nBlocksPerRow, nBlocksPerColumn;
618 int nSampleRate;
619 int bGotNoDataValue;
620 double dfNoDataValue;
621 int nActualSamples = 0;
622 int nBlockSampleRate;
623 int nBlockPixels, nBlockCount;
624
625 dfNoDataValue = poBand->GetNoDataValue( &bGotNoDataValue );
626
627 poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
628
629 nBlocksPerRow = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
630 nBlocksPerColumn = (poBand->GetYSize() + nBlockYSize - 1) / nBlockYSize;
631
632 nBlockPixels = nBlockXSize * nBlockYSize;
633 nBlockCount = nBlocksPerRow * nBlocksPerColumn;
634
635 if( nBlocksPerRow == 0 || nBlocksPerColumn == 0 || nBlockPixels == 0
636 || nBlockCount == 0 )
637 {
638 CPLError( CE_Failure, CPLE_AppDefined,
639 "GDALGetRandomRasterSample(): returning because band"
640 " appears degenerate." );
641
642 return FALSE;
643 }
644
645 nSampleRate = (int) MAX(1,sqrt((double) nBlockCount)-2.0);
646
647 if( nSampleRate == nBlocksPerRow && nSampleRate > 1 )
648 nSampleRate--;
649
650 while( nSampleRate > 1
651 && ((nBlockCount-1) / nSampleRate + 1) * nBlockPixels < nSamples )
652 nSampleRate--;
653
654 if ((nSamples / ((nBlockCount-1) / nSampleRate + 1)) == 0)
655 nBlockSampleRate = 1;
656 else
657 nBlockSampleRate =
658 MAX(1,nBlockPixels / (nSamples / ((nBlockCount-1) / nSampleRate + 1)));
659
660 for( int iSampleBlock = 0;
661 iSampleBlock < nBlockCount;
662 iSampleBlock += nSampleRate )
663 {
664 double dfValue = 0.0, dfReal, dfImag;
665 int iXBlock, iYBlock, iX, iY, iXValid, iYValid, iRemainder = 0;
666 GDALRasterBlock *poBlock;
667
668 iYBlock = iSampleBlock / nBlocksPerRow;
669 iXBlock = iSampleBlock - nBlocksPerRow * iYBlock;
670
671 poBlock = poBand->GetLockedBlockRef( iXBlock, iYBlock );
672 if( poBlock == NULL )
673 continue;
674 if( poBlock->GetDataRef() == NULL )
675 {
676 poBlock->DropLock();
677 continue;
678 }
679
680 if( (iXBlock + 1) * nBlockXSize > poBand->GetXSize() )
681 iXValid = poBand->GetXSize() - iXBlock * nBlockXSize;
682 else
683 iXValid = nBlockXSize;
684
685 if( (iYBlock + 1) * nBlockYSize > poBand->GetYSize() )
686 iYValid = poBand->GetYSize() - iYBlock * nBlockYSize;
687 else
688 iYValid = nBlockYSize;
689
690 for( iY = 0; iY < iYValid; iY++ )
691 {
692 for( iX = iRemainder; iX < iXValid; iX += nBlockSampleRate )
693 {
694 int iOffset;
695
696 iOffset = iX + iY * nBlockXSize;
697 switch( poBlock->GetDataType() )
698 {
699 case GDT_Byte:
700 dfValue = ((GByte *) poBlock->GetDataRef())[iOffset];
701 break;
702 case GDT_UInt16:
703 dfValue = ((GUInt16 *) poBlock->GetDataRef())[iOffset];
704 break;
705 case GDT_Int16:
706 dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset];
707 break;
708 case GDT_UInt32:
709 dfValue = ((GUInt32 *) poBlock->GetDataRef())[iOffset];
710 break;
711 case GDT_Int32:
712 dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset];
713 break;
714 case GDT_Float32:
715 dfValue = ((float *) poBlock->GetDataRef())[iOffset];
716 break;
717 case GDT_Float64:
718 dfValue = ((double *) poBlock->GetDataRef())[iOffset];
719 break;
720 case GDT_CInt16:
721 dfReal = ((GInt16 *) poBlock->GetDataRef())[iOffset*2];
722 dfImag = ((GInt16 *) poBlock->GetDataRef())[iOffset*2+1];
723 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
724 break;
725 case GDT_CInt32:
726 dfReal = ((GInt32 *) poBlock->GetDataRef())[iOffset*2];
727 dfImag = ((GInt32 *) poBlock->GetDataRef())[iOffset*2+1];
728 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
729 break;
730 case GDT_CFloat32:
731 dfReal = ((float *) poBlock->GetDataRef())[iOffset*2];
732 dfImag = ((float *) poBlock->GetDataRef())[iOffset*2+1];
733 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
734 break;
735 case GDT_CFloat64:
736 dfReal = ((double *) poBlock->GetDataRef())[iOffset*2];
737 dfImag = ((double *) poBlock->GetDataRef())[iOffset*2+1];
738 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
739 break;
740 default:
741 CPLAssert( FALSE );
742 }
743
744 if( bGotNoDataValue && dfValue == dfNoDataValue )
745 continue;
746
747 if( nActualSamples < nSamples )
748 pafSampleBuf[nActualSamples++] = (float) dfValue;
749 }
750
751 iRemainder = iX - iXValid;
752 }
753
754 poBlock->DropLock();
755 }
756
757 return nActualSamples;
758 }
759
760 /************************************************************************/
761 /* GDALInitGCPs() */
762 /************************************************************************/
763
GDALInitGCPs(int nCount,GDAL_GCP * psGCP)764 void CPL_STDCALL GDALInitGCPs( int nCount, GDAL_GCP * psGCP )
765
766 {
767 if( nCount > 0 )
768 {
769 VALIDATE_POINTER0( psGCP, "GDALInitGCPs" );
770 }
771
772 for( int iGCP = 0; iGCP < nCount; iGCP++ )
773 {
774 memset( psGCP, 0, sizeof(GDAL_GCP) );
775 psGCP->pszId = CPLStrdup("");
776 psGCP->pszInfo = CPLStrdup("");
777 psGCP++;
778 }
779 }
780
781 /************************************************************************/
782 /* GDALDeinitGCPs() */
783 /************************************************************************/
784
GDALDeinitGCPs(int nCount,GDAL_GCP * psGCP)785 void CPL_STDCALL GDALDeinitGCPs( int nCount, GDAL_GCP * psGCP )
786
787 {
788 if ( nCount > 0 )
789 {
790 VALIDATE_POINTER0( psGCP, "GDALDeinitGCPs" );
791 }
792
793 for( int iGCP = 0; iGCP < nCount; iGCP++ )
794 {
795 CPLFree( psGCP->pszId );
796 CPLFree( psGCP->pszInfo );
797 psGCP++;
798 }
799 }
800
801 /************************************************************************/
802 /* GDALDuplicateGCPs() */
803 /************************************************************************/
804
805 GDAL_GCP * CPL_STDCALL
GDALDuplicateGCPs(int nCount,const GDAL_GCP * pasGCPList)806 GDALDuplicateGCPs( int nCount, const GDAL_GCP *pasGCPList )
807
808 {
809 GDAL_GCP *pasReturn;
810
811 pasReturn = (GDAL_GCP *) CPLMalloc(sizeof(GDAL_GCP) * nCount);
812 GDALInitGCPs( nCount, pasReturn );
813
814 for( int iGCP = 0; iGCP < nCount; iGCP++ )
815 {
816 CPLFree( pasReturn[iGCP].pszId );
817 pasReturn[iGCP].pszId = CPLStrdup( pasGCPList[iGCP].pszId );
818
819 CPLFree( pasReturn[iGCP].pszInfo );
820 pasReturn[iGCP].pszInfo = CPLStrdup( pasGCPList[iGCP].pszInfo );
821
822 pasReturn[iGCP].dfGCPPixel = pasGCPList[iGCP].dfGCPPixel;
823 pasReturn[iGCP].dfGCPLine = pasGCPList[iGCP].dfGCPLine;
824 pasReturn[iGCP].dfGCPX = pasGCPList[iGCP].dfGCPX;
825 pasReturn[iGCP].dfGCPY = pasGCPList[iGCP].dfGCPY;
826 pasReturn[iGCP].dfGCPZ = pasGCPList[iGCP].dfGCPZ;
827 }
828
829 return pasReturn;
830 }
831
832 /************************************************************************/
833 /* GDALFindAssociatedFile() */
834 /************************************************************************/
835
836 /**
837 * Find file with alternate extension.
838 *
839 * Finds the file with the indicated extension, substituting it in place
840 * of the extension of the base filename. Generally used to search for
841 * associated files like world files .RPB files, etc. If necessary, the
842 * extension will be tried in both upper and lower case. If a sibling file
843 * list is available it will be used instead of doing VSIStatExL() calls to
844 * probe the file system.
845 *
846 * Note that the result is a dynamic CPLString so this method should not
847 * be used in a situation where there could be cross heap issues. It is
848 * generally imprudent for application built on GDAL to use this function
849 * unless they are sure they will always use the same runtime heap as GDAL.
850 *
851 * @param pszBaseFilename the filename relative to which to search.
852 * @param pszExt the target extension in either upper or lower case.
853 * @param papszSiblingFiles the list of files in the same directory as
854 * pszBaseFilename or NULL if they are not known.
855 * @param nFlags special options controlling search. None defined yet, just
856 * pass 0.
857 *
858 * @return an empty string if the target is not found, otherwise the target
859 * file with similar path style as the pszBaseFilename.
860 */
861
GDALFindAssociatedFile(const char * pszBaseFilename,const char * pszExt,char ** papszSiblingFiles,int nFlags)862 CPLString GDALFindAssociatedFile( const char *pszBaseFilename,
863 const char *pszExt,
864 char **papszSiblingFiles,
865 int nFlags )
866
867 {
868 (void) nFlags;
869
870 CPLString osTarget = CPLResetExtension( pszBaseFilename, pszExt );
871
872 if( papszSiblingFiles == NULL )
873 {
874 VSIStatBufL sStatBuf;
875
876 if( VSIStatExL( osTarget, &sStatBuf, VSI_STAT_EXISTS_FLAG ) != 0 )
877 {
878 CPLString osAltExt = pszExt;
879
880 if( islower( pszExt[0] ) )
881 osAltExt.toupper();
882 else
883 osAltExt.tolower();
884
885 osTarget = CPLResetExtension( pszBaseFilename, osAltExt );
886
887 if( VSIStatExL( osTarget, &sStatBuf, VSI_STAT_EXISTS_FLAG ) != 0 )
888 return "";
889 }
890 }
891 else
892 {
893 int iSibling = CSLFindString( papszSiblingFiles,
894 CPLGetFilename(osTarget) );
895 if( iSibling < 0 )
896 return "";
897
898 osTarget.resize(osTarget.size() - strlen(papszSiblingFiles[iSibling]));
899 osTarget += papszSiblingFiles[iSibling];
900 }
901
902 return osTarget;
903 }
904
905 /************************************************************************/
906 /* GDALLoadOziMapFile() */
907 /************************************************************************/
908
909 #define MAX_GCP 30
910
GDALLoadOziMapFile(const char * pszFilename,double * padfGeoTransform,char ** ppszWKT,int * pnGCPCount,GDAL_GCP ** ppasGCPs)911 int CPL_STDCALL GDALLoadOziMapFile( const char *pszFilename,
912 double *padfGeoTransform, char **ppszWKT,
913 int *pnGCPCount, GDAL_GCP **ppasGCPs )
914
915
916 {
917 char **papszLines;
918 int iLine, nLines=0;
919 int nCoordinateCount = 0;
920 GDAL_GCP asGCPs[MAX_GCP];
921
922 VALIDATE_POINTER1( pszFilename, "GDALLoadOziMapFile", FALSE );
923 VALIDATE_POINTER1( padfGeoTransform, "GDALLoadOziMapFile", FALSE );
924 VALIDATE_POINTER1( pnGCPCount, "GDALLoadOziMapFile", FALSE );
925 VALIDATE_POINTER1( ppasGCPs, "GDALLoadOziMapFile", FALSE );
926
927 papszLines = CSLLoad2( pszFilename, 1000, 200, NULL );
928
929 if ( !papszLines )
930 return FALSE;
931
932 nLines = CSLCount( papszLines );
933
934 // Check the OziExplorer Map file signature
935 if ( nLines < 5
936 || !EQUALN(papszLines[0], "OziExplorer Map Data File Version ", 34) )
937 {
938 CPLError( CE_Failure, CPLE_AppDefined,
939 "GDALLoadOziMapFile(): file \"%s\" is not in OziExplorer Map format.",
940 pszFilename );
941 CSLDestroy( papszLines );
942 return FALSE;
943 }
944
945 OGRSpatialReference oSRS;
946 OGRErr eErr = OGRERR_NONE;
947
948 /* The Map Scale Factor has been introduced recently on the 6th line */
949 /* and is a trick that is used to just change that line without changing */
950 /* the rest of the MAP file but providing an imagery that is smaller or larger */
951 /* so we have to correct the pixel/line values read in the .MAP file so they */
952 /* match the actual imagery dimension. Well, this is a bad summary of what */
953 /* is explained at http://tech.groups.yahoo.com/group/OziUsers-L/message/12484 */
954 double dfMSF = 1;
955
956 for ( iLine = 5; iLine < nLines; iLine++ )
957 {
958 if ( EQUALN(papszLines[iLine], "MSF,", 4) )
959 {
960 dfMSF = CPLAtof(papszLines[iLine] + 4);
961 if (dfMSF <= 0.01) /* Suspicious values */
962 {
963 CPLDebug("OZI", "Suspicious MSF value : %s", papszLines[iLine]);
964 dfMSF = 1;
965 }
966 }
967 }
968
969 eErr = oSRS.importFromOzi( papszLines );
970 if ( eErr == OGRERR_NONE )
971 {
972 if ( ppszWKT != NULL )
973 oSRS.exportToWkt( ppszWKT );
974 }
975
976 // Iterate all lines in the MAP-file
977 for ( iLine = 5; iLine < nLines; iLine++ )
978 {
979 char **papszTok = NULL;
980
981 papszTok = CSLTokenizeString2( papszLines[iLine], ",",
982 CSLT_ALLOWEMPTYTOKENS
983 | CSLT_STRIPLEADSPACES
984 | CSLT_STRIPENDSPACES );
985
986 if ( CSLCount(papszTok) < 12 )
987 {
988 CSLDestroy(papszTok);
989 continue;
990 }
991
992 if ( CSLCount(papszTok) >= 17
993 && EQUALN(papszTok[0], "Point", 5)
994 && !EQUAL(papszTok[2], "")
995 && !EQUAL(papszTok[3], "")
996 && nCoordinateCount < MAX_GCP )
997 {
998 int bReadOk = FALSE;
999 double dfLon = 0., dfLat = 0.;
1000
1001 if ( !EQUAL(papszTok[6], "")
1002 && !EQUAL(papszTok[7], "")
1003 && !EQUAL(papszTok[9], "")
1004 && !EQUAL(papszTok[10], "") )
1005 {
1006 // Set geographical coordinates of the pixels
1007 dfLon = CPLAtofM(papszTok[9]) + CPLAtofM(papszTok[10]) / 60.0;
1008 dfLat = CPLAtofM(papszTok[6]) + CPLAtofM(papszTok[7]) / 60.0;
1009 if ( EQUAL(papszTok[11], "W") )
1010 dfLon = -dfLon;
1011 if ( EQUAL(papszTok[8], "S") )
1012 dfLat = -dfLat;
1013
1014 // Transform from the geographical coordinates into projected
1015 // coordinates.
1016 if ( eErr == OGRERR_NONE )
1017 {
1018 OGRSpatialReference *poLatLong = NULL;
1019 OGRCoordinateTransformation *poTransform = NULL;
1020
1021 poLatLong = oSRS.CloneGeogCS();
1022 if ( poLatLong )
1023 {
1024 poTransform = OGRCreateCoordinateTransformation( poLatLong, &oSRS );
1025 if ( poTransform )
1026 {
1027 bReadOk = poTransform->Transform( 1, &dfLon, &dfLat );
1028 delete poTransform;
1029 }
1030 delete poLatLong;
1031 }
1032 }
1033 }
1034 else if ( !EQUAL(papszTok[14], "")
1035 && !EQUAL(papszTok[15], "") )
1036 {
1037 // Set cartesian coordinates of the pixels.
1038 dfLon = CPLAtofM(papszTok[14]);
1039 dfLat = CPLAtofM(papszTok[15]);
1040 bReadOk = TRUE;
1041
1042 //if ( EQUAL(papszTok[16], "S") )
1043 // dfLat = -dfLat;
1044 }
1045
1046 if ( bReadOk )
1047 {
1048 GDALInitGCPs( 1, asGCPs + nCoordinateCount );
1049
1050 // Set pixel/line part
1051 asGCPs[nCoordinateCount].dfGCPPixel = CPLAtofM(papszTok[2]) / dfMSF;
1052 asGCPs[nCoordinateCount].dfGCPLine = CPLAtofM(papszTok[3]) / dfMSF;
1053
1054 asGCPs[nCoordinateCount].dfGCPX = dfLon;
1055 asGCPs[nCoordinateCount].dfGCPY = dfLat;
1056
1057 nCoordinateCount++;
1058 }
1059 }
1060
1061 CSLDestroy( papszTok );
1062 }
1063
1064 CSLDestroy( papszLines );
1065
1066 if ( nCoordinateCount == 0 )
1067 {
1068 CPLDebug( "GDAL", "GDALLoadOziMapFile(\"%s\") did read no GCPs.",
1069 pszFilename );
1070 return FALSE;
1071 }
1072
1073 /* -------------------------------------------------------------------- */
1074 /* Try to convert the GCPs into a geotransform definition, if */
1075 /* possible. Otherwise we will need to use them as GCPs. */
1076 /* -------------------------------------------------------------------- */
1077 if( !GDALGCPsToGeoTransform( nCoordinateCount, asGCPs, padfGeoTransform,
1078 CSLTestBoolean(CPLGetConfigOption("OZI_APPROX_GEOTRANSFORM", "NO")) ) )
1079 {
1080 if ( pnGCPCount && ppasGCPs )
1081 {
1082 CPLDebug( "GDAL",
1083 "GDALLoadOziMapFile(%s) found file, wasn't able to derive a\n"
1084 "first order geotransform. Using points as GCPs.",
1085 pszFilename );
1086
1087 *ppasGCPs = (GDAL_GCP *)
1088 CPLCalloc( sizeof(GDAL_GCP),nCoordinateCount );
1089 memcpy( *ppasGCPs, asGCPs, sizeof(GDAL_GCP) * nCoordinateCount );
1090 *pnGCPCount = nCoordinateCount;
1091 }
1092 }
1093 else
1094 {
1095 GDALDeinitGCPs( nCoordinateCount, asGCPs );
1096 }
1097
1098 return TRUE;
1099 }
1100
1101 #undef MAX_GCP
1102
1103 /************************************************************************/
1104 /* GDALReadOziMapFile() */
1105 /************************************************************************/
1106
GDALReadOziMapFile(const char * pszBaseFilename,double * padfGeoTransform,char ** ppszWKT,int * pnGCPCount,GDAL_GCP ** ppasGCPs)1107 int CPL_STDCALL GDALReadOziMapFile( const char * pszBaseFilename,
1108 double *padfGeoTransform, char **ppszWKT,
1109 int *pnGCPCount, GDAL_GCP **ppasGCPs )
1110
1111
1112 {
1113 const char *pszOzi;
1114 FILE *fpOzi;
1115
1116 /* -------------------------------------------------------------------- */
1117 /* Try lower case, then upper case. */
1118 /* -------------------------------------------------------------------- */
1119 pszOzi = CPLResetExtension( pszBaseFilename, "map" );
1120
1121 fpOzi = VSIFOpen( pszOzi, "rt" );
1122
1123 if ( fpOzi == NULL && VSIIsCaseSensitiveFS(pszOzi) )
1124 {
1125 pszOzi = CPLResetExtension( pszBaseFilename, "MAP" );
1126 fpOzi = VSIFOpen( pszOzi, "rt" );
1127 }
1128
1129 if ( fpOzi == NULL )
1130 return FALSE;
1131
1132 VSIFClose( fpOzi );
1133
1134 /* -------------------------------------------------------------------- */
1135 /* We found the file, now load and parse it. */
1136 /* -------------------------------------------------------------------- */
1137 return GDALLoadOziMapFile( pszOzi, padfGeoTransform, ppszWKT,
1138 pnGCPCount, ppasGCPs );
1139 }
1140
1141 /************************************************************************/
1142 /* GDALLoadTabFile() */
1143 /* */
1144 /* Helper function for translator implementators wanting */
1145 /* support for MapInfo .tab-files. */
1146 /************************************************************************/
1147
1148 #define MAX_GCP 256
1149
GDALLoadTabFile(const char * pszFilename,double * padfGeoTransform,char ** ppszWKT,int * pnGCPCount,GDAL_GCP ** ppasGCPs)1150 int CPL_STDCALL GDALLoadTabFile( const char *pszFilename,
1151 double *padfGeoTransform, char **ppszWKT,
1152 int *pnGCPCount, GDAL_GCP **ppasGCPs )
1153
1154
1155 {
1156 char **papszLines;
1157 char **papszTok=NULL;
1158 int bTypeRasterFound = FALSE;
1159 int bInsideTableDef = FALSE;
1160 int iLine, numLines=0;
1161 int nCoordinateCount = 0;
1162 GDAL_GCP asGCPs[MAX_GCP];
1163
1164 papszLines = CSLLoad2( pszFilename, 1000, 200, NULL );
1165
1166 if ( !papszLines )
1167 return FALSE;
1168
1169 numLines = CSLCount(papszLines);
1170
1171 // Iterate all lines in the TAB-file
1172 for(iLine=0; iLine<numLines; iLine++)
1173 {
1174 CSLDestroy(papszTok);
1175 papszTok = CSLTokenizeStringComplex(papszLines[iLine], " \t(),;",
1176 TRUE, FALSE);
1177
1178 if (CSLCount(papszTok) < 2)
1179 continue;
1180
1181 // Did we find table definition
1182 if (EQUAL(papszTok[0], "Definition") && EQUAL(papszTok[1], "Table") )
1183 {
1184 bInsideTableDef = TRUE;
1185 }
1186 else if (bInsideTableDef && (EQUAL(papszTok[0], "Type")) )
1187 {
1188 // Only RASTER-type will be handled
1189 if (EQUAL(papszTok[1], "RASTER"))
1190 {
1191 bTypeRasterFound = TRUE;
1192 }
1193 else
1194 {
1195 CSLDestroy(papszTok);
1196 CSLDestroy(papszLines);
1197 return FALSE;
1198 }
1199 }
1200 else if (bTypeRasterFound && bInsideTableDef
1201 && CSLCount(papszTok) > 4
1202 && EQUAL(papszTok[4], "Label")
1203 && nCoordinateCount < MAX_GCP )
1204 {
1205 GDALInitGCPs( 1, asGCPs + nCoordinateCount );
1206
1207 asGCPs[nCoordinateCount].dfGCPPixel = CPLAtofM(papszTok[2]);
1208 asGCPs[nCoordinateCount].dfGCPLine = CPLAtofM(papszTok[3]);
1209 asGCPs[nCoordinateCount].dfGCPX = CPLAtofM(papszTok[0]);
1210 asGCPs[nCoordinateCount].dfGCPY = CPLAtofM(papszTok[1]);
1211 if( papszTok[5] != NULL )
1212 {
1213 CPLFree( asGCPs[nCoordinateCount].pszId );
1214 asGCPs[nCoordinateCount].pszId = CPLStrdup(papszTok[5]);
1215 }
1216
1217 nCoordinateCount++;
1218 }
1219 else if( bTypeRasterFound && bInsideTableDef
1220 && EQUAL(papszTok[0],"CoordSys")
1221 && ppszWKT != NULL )
1222 {
1223 OGRSpatialReference oSRS;
1224
1225 if( oSRS.importFromMICoordSys( papszLines[iLine] ) == OGRERR_NONE )
1226 oSRS.exportToWkt( ppszWKT );
1227 }
1228 else if( EQUAL(papszTok[0],"Units")
1229 && CSLCount(papszTok) > 1
1230 && EQUAL(papszTok[1],"degree") )
1231 {
1232 /*
1233 ** If we have units of "degree", but a projected coordinate
1234 ** system we need to convert it to geographic. See to01_02.TAB.
1235 */
1236 if( ppszWKT != NULL && *ppszWKT != NULL
1237 && EQUALN(*ppszWKT,"PROJCS",6) )
1238 {
1239 OGRSpatialReference oSRS, oSRSGeogCS;
1240 char *pszSrcWKT = *ppszWKT;
1241
1242 oSRS.importFromWkt( &pszSrcWKT );
1243 oSRSGeogCS.CopyGeogCSFrom( &oSRS );
1244 CPLFree( *ppszWKT );
1245 oSRSGeogCS.exportToWkt( ppszWKT );
1246 }
1247 }
1248
1249 }
1250
1251 CSLDestroy(papszTok);
1252 CSLDestroy(papszLines);
1253
1254 if( nCoordinateCount == 0 )
1255 {
1256 CPLDebug( "GDAL", "GDALLoadTabFile(%s) did not get any GCPs.",
1257 pszFilename );
1258 return FALSE;
1259 }
1260
1261 /* -------------------------------------------------------------------- */
1262 /* Try to convert the GCPs into a geotransform definition, if */
1263 /* possible. Otherwise we will need to use them as GCPs. */
1264 /* -------------------------------------------------------------------- */
1265 if( !GDALGCPsToGeoTransform( nCoordinateCount, asGCPs, padfGeoTransform,
1266 CSLTestBoolean(CPLGetConfigOption("TAB_APPROX_GEOTRANSFORM", "NO")) ) )
1267 {
1268 if (pnGCPCount && ppasGCPs)
1269 {
1270 CPLDebug( "GDAL",
1271 "GDALLoadTabFile(%s) found file, wasn't able to derive a\n"
1272 "first order geotransform. Using points as GCPs.",
1273 pszFilename );
1274
1275 *ppasGCPs = (GDAL_GCP *)
1276 CPLCalloc( sizeof(GDAL_GCP),nCoordinateCount );
1277 memcpy( *ppasGCPs, asGCPs, sizeof(GDAL_GCP) * nCoordinateCount );
1278 *pnGCPCount = nCoordinateCount;
1279 }
1280 }
1281 else
1282 {
1283 GDALDeinitGCPs( nCoordinateCount, asGCPs );
1284 }
1285
1286 return TRUE;
1287 }
1288
1289 #undef MAX_GCP
1290
1291 /************************************************************************/
1292 /* GDALReadTabFile() */
1293 /* */
1294 /* Helper function for translator implementators wanting */
1295 /* support for MapInfo .tab-files. */
1296 /************************************************************************/
1297
GDALReadTabFile(const char * pszBaseFilename,double * padfGeoTransform,char ** ppszWKT,int * pnGCPCount,GDAL_GCP ** ppasGCPs)1298 int CPL_STDCALL GDALReadTabFile( const char * pszBaseFilename,
1299 double *padfGeoTransform, char **ppszWKT,
1300 int *pnGCPCount, GDAL_GCP **ppasGCPs )
1301
1302
1303 {
1304 return GDALReadTabFile2(pszBaseFilename, padfGeoTransform,
1305 ppszWKT, pnGCPCount, ppasGCPs,
1306 NULL, NULL);
1307 }
1308
1309
GDALReadTabFile2(const char * pszBaseFilename,double * padfGeoTransform,char ** ppszWKT,int * pnGCPCount,GDAL_GCP ** ppasGCPs,char ** papszSiblingFiles,char ** ppszTabFileNameOut)1310 int GDALReadTabFile2( const char * pszBaseFilename,
1311 double *padfGeoTransform, char **ppszWKT,
1312 int *pnGCPCount, GDAL_GCP **ppasGCPs,
1313 char** papszSiblingFiles, char** ppszTabFileNameOut )
1314 {
1315 const char *pszTAB;
1316 VSILFILE *fpTAB;
1317
1318 if (ppszTabFileNameOut)
1319 *ppszTabFileNameOut = NULL;
1320
1321 if( !GDALCanFileAcceptSidecarFile(pszBaseFilename) )
1322 return FALSE;
1323
1324 pszTAB = CPLResetExtension( pszBaseFilename, "tab" );
1325
1326 if (papszSiblingFiles)
1327 {
1328 int iSibling = CSLFindString(papszSiblingFiles, CPLGetFilename(pszTAB));
1329 if (iSibling >= 0)
1330 {
1331 CPLString osTabFilename = pszBaseFilename;
1332 osTabFilename.resize(strlen(pszBaseFilename) -
1333 strlen(CPLGetFilename(pszBaseFilename)));
1334 osTabFilename += papszSiblingFiles[iSibling];
1335 if ( GDALLoadTabFile(osTabFilename, padfGeoTransform, ppszWKT,
1336 pnGCPCount, ppasGCPs ) )
1337 {
1338 if (ppszTabFileNameOut)
1339 *ppszTabFileNameOut = CPLStrdup(osTabFilename);
1340 return TRUE;
1341 }
1342 }
1343 return FALSE;
1344 }
1345
1346 /* -------------------------------------------------------------------- */
1347 /* Try lower case, then upper case. */
1348 /* -------------------------------------------------------------------- */
1349
1350 fpTAB = VSIFOpenL( pszTAB, "rt" );
1351
1352 if( fpTAB == NULL && VSIIsCaseSensitiveFS(pszTAB) )
1353 {
1354 pszTAB = CPLResetExtension( pszBaseFilename, "TAB" );
1355 fpTAB = VSIFOpenL( pszTAB, "rt" );
1356 }
1357
1358 if( fpTAB == NULL )
1359 return FALSE;
1360
1361 VSIFCloseL( fpTAB );
1362
1363 /* -------------------------------------------------------------------- */
1364 /* We found the file, now load and parse it. */
1365 /* -------------------------------------------------------------------- */
1366 if (GDALLoadTabFile( pszTAB, padfGeoTransform, ppszWKT,
1367 pnGCPCount, ppasGCPs ) )
1368 {
1369 if (ppszTabFileNameOut)
1370 *ppszTabFileNameOut = CPLStrdup(pszTAB);
1371 return TRUE;
1372 }
1373 return FALSE;
1374 }
1375
1376 /************************************************************************/
1377 /* GDALLoadWorldFile() */
1378 /************************************************************************/
1379
1380 /**
1381 * \brief Read ESRI world file.
1382 *
1383 * This function reads an ESRI style world file, and formats a geotransform
1384 * from its contents.
1385 *
1386 * The world file contains an affine transformation with the parameters
1387 * in a different order than in a geotransform array.
1388 *
1389 * <ul>
1390 * <li> geotransform[1] : width of pixel
1391 * <li> geotransform[4] : rotational coefficient, zero for north up images.
1392 * <li> geotransform[2] : rotational coefficient, zero for north up images.
1393 * <li> geotransform[5] : height of pixel (but negative)
1394 * <li> geotransform[0] + 0.5 * geotransform[1] + 0.5 * geotransform[2] : x offset to center of top left pixel.
1395 * <li> geotransform[3] + 0.5 * geotransform[4] + 0.5 * geotransform[5] : y offset to center of top left pixel.
1396 * </ul>
1397 *
1398 * @param pszFilename the world file name.
1399 * @param padfGeoTransform the six double array into which the
1400 * geotransformation should be placed.
1401 *
1402 * @return TRUE on success or FALSE on failure.
1403 */
1404
1405 int CPL_STDCALL
GDALLoadWorldFile(const char * pszFilename,double * padfGeoTransform)1406 GDALLoadWorldFile( const char *pszFilename, double *padfGeoTransform )
1407
1408 {
1409 char **papszLines;
1410
1411 VALIDATE_POINTER1( pszFilename, "GDALLoadWorldFile", FALSE );
1412 VALIDATE_POINTER1( padfGeoTransform, "GDALLoadWorldFile", FALSE );
1413
1414 papszLines = CSLLoad2( pszFilename, 100, 100, NULL );
1415
1416 if ( !papszLines )
1417 return FALSE;
1418
1419 double world[6];
1420 // reads the first 6 non-empty lines
1421 int nLines = 0;
1422 int nLinesCount = CSLCount(papszLines);
1423 for( int i = 0; i < nLinesCount && nLines < 6; ++i )
1424 {
1425 CPLString line(papszLines[i]);
1426 if( line.Trim().empty() )
1427 continue;
1428
1429 world[nLines] = CPLAtofM(line);
1430 ++nLines;
1431 }
1432
1433 if( nLines == 6
1434 && (world[0] != 0.0 || world[2] != 0.0)
1435 && (world[3] != 0.0 || world[1] != 0.0) )
1436 {
1437 padfGeoTransform[0] = world[4];
1438 padfGeoTransform[1] = world[0];
1439 padfGeoTransform[2] = world[2];
1440 padfGeoTransform[3] = world[5];
1441 padfGeoTransform[4] = world[1];
1442 padfGeoTransform[5] = world[3];
1443
1444 // correct for center of pixel vs. top left of pixel
1445 padfGeoTransform[0] -= 0.5 * padfGeoTransform[1];
1446 padfGeoTransform[0] -= 0.5 * padfGeoTransform[2];
1447 padfGeoTransform[3] -= 0.5 * padfGeoTransform[4];
1448 padfGeoTransform[3] -= 0.5 * padfGeoTransform[5];
1449
1450 CSLDestroy(papszLines);
1451
1452 return TRUE;
1453 }
1454 else
1455 {
1456 CPLDebug( "GDAL",
1457 "GDALLoadWorldFile(%s) found file, but it was corrupt.",
1458 pszFilename );
1459 CSLDestroy(papszLines);
1460 return FALSE;
1461 }
1462 }
1463
1464 /************************************************************************/
1465 /* GDALReadWorldFile() */
1466 /************************************************************************/
1467
1468 /**
1469 * \brief Read ESRI world file.
1470 *
1471 * This function reads an ESRI style world file, and formats a geotransform
1472 * from its contents. It does the same as GDALLoadWorldFile() function, but
1473 * it will form the filename for the worldfile from the filename of the raster
1474 * file referred and the suggested extension. If no extension is provided,
1475 * the code will internally try the unix style and windows style world file
1476 * extensions (eg. for .tif these would be .tfw and .tifw).
1477 *
1478 * The world file contains an affine transformation with the parameters
1479 * in a different order than in a geotransform array.
1480 *
1481 * <ul>
1482 * <li> geotransform[1] : width of pixel
1483 * <li> geotransform[4] : rotational coefficient, zero for north up images.
1484 * <li> geotransform[2] : rotational coefficient, zero for north up images.
1485 * <li> geotransform[5] : height of pixel (but negative)
1486 * <li> geotransform[0] + 0.5 * geotransform[1] + 0.5 * geotransform[2] : x offset to center of top left pixel.
1487 * <li> geotransform[3] + 0.5 * geotransform[4] + 0.5 * geotransform[5] : y offset to center of top left pixel.
1488 * </ul>
1489 *
1490 * @param pszBaseFilename the target raster file.
1491 * @param pszExtension the extension to use (ie. ".wld") or NULL to derive it
1492 * from the pszBaseFilename
1493 * @param padfGeoTransform the six double array into which the
1494 * geotransformation should be placed.
1495 *
1496 * @return TRUE on success or FALSE on failure.
1497 */
1498
1499 int CPL_STDCALL
GDALReadWorldFile(const char * pszBaseFilename,const char * pszExtension,double * padfGeoTransform)1500 GDALReadWorldFile( const char *pszBaseFilename, const char *pszExtension,
1501 double *padfGeoTransform )
1502
1503 {
1504 return GDALReadWorldFile2(pszBaseFilename, pszExtension,
1505 padfGeoTransform, NULL, NULL);
1506 }
1507
GDALReadWorldFile2(const char * pszBaseFilename,const char * pszExtension,double * padfGeoTransform,char ** papszSiblingFiles,char ** ppszWorldFileNameOut)1508 int GDALReadWorldFile2( const char *pszBaseFilename, const char *pszExtension,
1509 double *padfGeoTransform, char** papszSiblingFiles,
1510 char** ppszWorldFileNameOut )
1511 {
1512 const char *pszTFW;
1513 char szExtUpper[32], szExtLower[32];
1514 int i;
1515
1516 VALIDATE_POINTER1( pszBaseFilename, "GDALReadWorldFile", FALSE );
1517 VALIDATE_POINTER1( padfGeoTransform, "GDALReadWorldFile", FALSE );
1518
1519 if (ppszWorldFileNameOut)
1520 *ppszWorldFileNameOut = NULL;
1521
1522 if( !GDALCanFileAcceptSidecarFile(pszBaseFilename) )
1523 return FALSE;
1524
1525 /* -------------------------------------------------------------------- */
1526 /* If we aren't given an extension, try both the unix and */
1527 /* windows style extensions. */
1528 /* -------------------------------------------------------------------- */
1529 if( pszExtension == NULL )
1530 {
1531 char szDerivedExtension[100];
1532 std::string oBaseExt = CPLGetExtension( pszBaseFilename );
1533
1534 if( oBaseExt.length() < 2 )
1535 return FALSE;
1536
1537 // windows version - first + last + 'w'
1538 szDerivedExtension[0] = oBaseExt[0];
1539 szDerivedExtension[1] = oBaseExt[oBaseExt.length()-1];
1540 szDerivedExtension[2] = 'w';
1541 szDerivedExtension[3] = '\0';
1542
1543 if( GDALReadWorldFile2( pszBaseFilename, szDerivedExtension,
1544 padfGeoTransform, papszSiblingFiles,
1545 ppszWorldFileNameOut ) )
1546 return TRUE;
1547
1548 // unix version - extension + 'w'
1549 if( oBaseExt.length() > sizeof(szDerivedExtension)-2 )
1550 return FALSE;
1551
1552 strcpy( szDerivedExtension, oBaseExt.c_str() );
1553 strcat( szDerivedExtension, "w" );
1554 return GDALReadWorldFile2( pszBaseFilename, szDerivedExtension,
1555 padfGeoTransform, papszSiblingFiles,
1556 ppszWorldFileNameOut );
1557 }
1558
1559 /* -------------------------------------------------------------------- */
1560 /* Skip the leading period in the extension if there is one. */
1561 /* -------------------------------------------------------------------- */
1562 if( *pszExtension == '.' )
1563 pszExtension++;
1564
1565 /* -------------------------------------------------------------------- */
1566 /* Generate upper and lower case versions of the extension. */
1567 /* -------------------------------------------------------------------- */
1568 CPLStrlcpy( szExtUpper, pszExtension, sizeof(szExtUpper) );
1569 CPLStrlcpy( szExtLower, pszExtension, sizeof(szExtLower) );
1570
1571 for( i = 0; szExtUpper[i] != '\0'; i++ )
1572 {
1573 szExtUpper[i] = (char) toupper(szExtUpper[i]);
1574 szExtLower[i] = (char) tolower(szExtLower[i]);
1575 }
1576
1577 VSIStatBufL sStatBuf;
1578 int bGotTFW;
1579
1580 pszTFW = CPLResetExtension( pszBaseFilename, szExtLower );
1581
1582 if (papszSiblingFiles)
1583 {
1584 int iSibling = CSLFindString(papszSiblingFiles, CPLGetFilename(pszTFW));
1585 if (iSibling >= 0)
1586 {
1587 CPLString osTFWFilename = pszBaseFilename;
1588 osTFWFilename.resize(strlen(pszBaseFilename) -
1589 strlen(CPLGetFilename(pszBaseFilename)));
1590 osTFWFilename += papszSiblingFiles[iSibling];
1591 if (GDALLoadWorldFile( osTFWFilename, padfGeoTransform ))
1592 {
1593 if (ppszWorldFileNameOut)
1594 *ppszWorldFileNameOut = CPLStrdup(osTFWFilename);
1595 return TRUE;
1596 }
1597 }
1598 return FALSE;
1599 }
1600
1601 /* -------------------------------------------------------------------- */
1602 /* Try lower case, then upper case. */
1603 /* -------------------------------------------------------------------- */
1604
1605 bGotTFW = VSIStatExL( pszTFW, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0;
1606
1607 if( !bGotTFW && VSIIsCaseSensitiveFS(pszTFW) )
1608 {
1609 pszTFW = CPLResetExtension( pszBaseFilename, szExtUpper );
1610 bGotTFW = VSIStatExL( pszTFW, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0;
1611 }
1612
1613 if( !bGotTFW )
1614 return FALSE;
1615
1616 /* -------------------------------------------------------------------- */
1617 /* We found the file, now load and parse it. */
1618 /* -------------------------------------------------------------------- */
1619 if (GDALLoadWorldFile( pszTFW, padfGeoTransform ))
1620 {
1621 if (ppszWorldFileNameOut)
1622 *ppszWorldFileNameOut = CPLStrdup(pszTFW);
1623 return TRUE;
1624 }
1625 return FALSE;
1626 }
1627
1628 /************************************************************************/
1629 /* GDALWriteWorldFile() */
1630 /* */
1631 /* Helper function for translator implementators wanting */
1632 /* support for ESRI world files. */
1633 /************************************************************************/
1634
1635 /**
1636 * \brief Write ESRI world file.
1637 *
1638 * This function writes an ESRI style world file from the passed geotransform.
1639 *
1640 * The world file contains an affine transformation with the parameters
1641 * in a different order than in a geotransform array.
1642 *
1643 * <ul>
1644 * <li> geotransform[1] : width of pixel
1645 * <li> geotransform[4] : rotational coefficient, zero for north up images.
1646 * <li> geotransform[2] : rotational coefficient, zero for north up images.
1647 * <li> geotransform[5] : height of pixel (but negative)
1648 * <li> geotransform[0] + 0.5 * geotransform[1] + 0.5 * geotransform[2] : x offset to center of top left pixel.
1649 * <li> geotransform[3] + 0.5 * geotransform[4] + 0.5 * geotransform[5] : y offset to center of top left pixel.
1650 * </ul>
1651 *
1652 * @param pszBaseFilename the target raster file.
1653 * @param pszExtension the extension to use (ie. ".wld"). Must not be NULL
1654 * @param padfGeoTransform the six double array from which the
1655 * geotransformation should be read.
1656 *
1657 * @return TRUE on success or FALSE on failure.
1658 */
1659
1660 int CPL_STDCALL
GDALWriteWorldFile(const char * pszBaseFilename,const char * pszExtension,double * padfGeoTransform)1661 GDALWriteWorldFile( const char * pszBaseFilename, const char *pszExtension,
1662 double *padfGeoTransform )
1663
1664 {
1665 VALIDATE_POINTER1( pszBaseFilename, "GDALWriteWorldFile", FALSE );
1666 VALIDATE_POINTER1( pszExtension, "GDALWriteWorldFile", FALSE );
1667 VALIDATE_POINTER1( padfGeoTransform, "GDALWriteWorldFile", FALSE );
1668
1669 /* -------------------------------------------------------------------- */
1670 /* Prepare the text to write to the file. */
1671 /* -------------------------------------------------------------------- */
1672 CPLString osTFWText;
1673
1674 osTFWText.Printf( "%.10f\n%.10f\n%.10f\n%.10f\n%.10f\n%.10f\n",
1675 padfGeoTransform[1],
1676 padfGeoTransform[4],
1677 padfGeoTransform[2],
1678 padfGeoTransform[5],
1679 padfGeoTransform[0]
1680 + 0.5 * padfGeoTransform[1]
1681 + 0.5 * padfGeoTransform[2],
1682 padfGeoTransform[3]
1683 + 0.5 * padfGeoTransform[4]
1684 + 0.5 * padfGeoTransform[5] );
1685
1686 /* -------------------------------------------------------------------- */
1687 /* Update extention, and write to disk. */
1688 /* -------------------------------------------------------------------- */
1689 const char *pszTFW;
1690 VSILFILE *fpTFW;
1691
1692 pszTFW = CPLResetExtension( pszBaseFilename, pszExtension );
1693 fpTFW = VSIFOpenL( pszTFW, "wt" );
1694 if( fpTFW == NULL )
1695 return FALSE;
1696
1697 VSIFWriteL( (void *) osTFWText.c_str(), 1, osTFWText.size(), fpTFW );
1698 VSIFCloseL( fpTFW );
1699
1700 return TRUE;
1701 }
1702
1703 /************************************************************************/
1704 /* GDALVersionInfo() */
1705 /************************************************************************/
1706
1707 /**
1708 * \brief Get runtime version information.
1709 *
1710 * Available pszRequest values:
1711 * <ul>
1712 * <li> "VERSION_NUM": Returns GDAL_VERSION_NUM formatted as a string. ie. "1170"
1713 * Note: starting with GDAL 1.10, this string will be longer than 4 characters.
1714 * <li> "RELEASE_DATE": Returns GDAL_RELEASE_DATE formatted as a string.
1715 * ie. "20020416".
1716 * <li> "RELEASE_NAME": Returns the GDAL_RELEASE_NAME. ie. "1.1.7"
1717 * <li> "--version": Returns one line version message suitable for use in
1718 * response to --version requests. ie. "GDAL 1.1.7, released 2002/04/16"
1719 * <li> "LICENSE": Returns the content of the LICENSE.TXT file from the GDAL_DATA directory.
1720 * Before GDAL 1.7.0, the returned string was leaking memory but this is now resolved.
1721 * So the result should not been freed by the caller.
1722 * <li> "BUILD_INFO": List of NAME=VALUE pairs separated by newlines with
1723 * information on build time options.
1724 * </ul>
1725 *
1726 * @param pszRequest the type of version info desired, as listed above.
1727 *
1728 * @return an internal string containing the requested information.
1729 */
1730
GDALVersionInfo(const char * pszRequest)1731 const char * CPL_STDCALL GDALVersionInfo( const char *pszRequest )
1732
1733 {
1734 /* -------------------------------------------------------------------- */
1735 /* Try to capture as much build information as practical. */
1736 /* -------------------------------------------------------------------- */
1737 if( pszRequest != NULL && EQUAL(pszRequest,"BUILD_INFO") )
1738 {
1739 CPLString osBuildInfo;
1740
1741 #ifdef ESRI_BUILD
1742 osBuildInfo += "ESRI_BUILD=YES\n";
1743 #endif
1744 #ifdef PAM_ENABLED
1745 osBuildInfo += "PAM_ENABLED=YES\n";
1746 #endif
1747 #ifdef OGR_ENABLED
1748 osBuildInfo += "OGR_ENABLED=YES\n";
1749 #endif
1750
1751 CPLFree(CPLGetTLS(CTLS_VERSIONINFO));
1752 CPLSetTLS(CTLS_VERSIONINFO, CPLStrdup(osBuildInfo), TRUE );
1753 return (char *) CPLGetTLS(CTLS_VERSIONINFO);
1754 }
1755
1756 /* -------------------------------------------------------------------- */
1757 /* LICENSE is a special case. We try to find and read the */
1758 /* LICENSE.TXT file from the GDAL_DATA directory and return it */
1759 /* -------------------------------------------------------------------- */
1760 if( pszRequest != NULL && EQUAL(pszRequest,"LICENSE") )
1761 {
1762 char* pszResultLicence = (char*) CPLGetTLS( CTLS_VERSIONINFO_LICENCE );
1763 if( pszResultLicence != NULL )
1764 {
1765 return pszResultLicence;
1766 }
1767
1768 const char *pszFilename = CPLFindFile( "etc", "LICENSE.TXT" );
1769 VSILFILE *fp = NULL;
1770 int nLength;
1771
1772 if( pszFilename != NULL )
1773 fp = VSIFOpenL( pszFilename, "r" );
1774
1775 if( fp != NULL )
1776 {
1777 VSIFSeekL( fp, 0, SEEK_END );
1778 nLength = (int) VSIFTellL( fp ) + 1;
1779 VSIFSeekL( fp, SEEK_SET, 0 );
1780
1781 pszResultLicence = (char *) VSICalloc(1,nLength);
1782 if (pszResultLicence)
1783 VSIFReadL( pszResultLicence, 1, nLength-1, fp );
1784
1785 VSIFCloseL( fp );
1786 }
1787
1788 if (!pszResultLicence)
1789 {
1790 pszResultLicence = CPLStrdup(
1791 "GDAL/OGR is released under the MIT/X license.\n"
1792 "The LICENSE.TXT distributed with GDAL/OGR should\n"
1793 "contain additional details.\n" );
1794 }
1795
1796 CPLSetTLS( CTLS_VERSIONINFO_LICENCE, pszResultLicence, TRUE );
1797 return pszResultLicence;
1798 }
1799
1800 /* -------------------------------------------------------------------- */
1801 /* All other strings are fairly small. */
1802 /* -------------------------------------------------------------------- */
1803 CPLString osVersionInfo;
1804
1805 if( pszRequest == NULL || EQUAL(pszRequest,"VERSION_NUM") )
1806 osVersionInfo.Printf( "%d", GDAL_VERSION_NUM );
1807 else if( EQUAL(pszRequest,"RELEASE_DATE") )
1808 osVersionInfo.Printf( "%d", GDAL_RELEASE_DATE );
1809 else if( EQUAL(pszRequest,"RELEASE_NAME") )
1810 osVersionInfo.Printf( GDAL_RELEASE_NAME );
1811 else // --version
1812 osVersionInfo.Printf( "GDAL %s, released %d/%02d/%02d",
1813 GDAL_RELEASE_NAME,
1814 GDAL_RELEASE_DATE / 10000,
1815 (GDAL_RELEASE_DATE % 10000) / 100,
1816 GDAL_RELEASE_DATE % 100 );
1817
1818 CPLFree(CPLGetTLS(CTLS_VERSIONINFO)); // clear old value.
1819 CPLSetTLS(CTLS_VERSIONINFO, CPLStrdup(osVersionInfo), TRUE );
1820 return (char *) CPLGetTLS(CTLS_VERSIONINFO);
1821 }
1822
1823 /************************************************************************/
1824 /* GDALCheckVersion() */
1825 /************************************************************************/
1826
1827 /** Return TRUE if GDAL library version at runtime matches nVersionMajor.nVersionMinor.
1828
1829 The purpose of this method is to ensure that calling code will run with the GDAL
1830 version it is compiled for. It is primarly intented for external plugins.
1831
1832 @param nVersionMajor Major version to be tested against
1833 @param nVersionMinor Minor version to be tested against
1834 @param pszCallingComponentName If not NULL, in case of version mismatch, the method
1835 will issue a failure mentionning the name of
1836 the calling component.
1837
1838 @return TRUE if GDAL library version at runtime matches nVersionMajor.nVersionMinor, FALSE otherwise.
1839 */
GDALCheckVersion(int nVersionMajor,int nVersionMinor,const char * pszCallingComponentName)1840 int CPL_STDCALL GDALCheckVersion( int nVersionMajor, int nVersionMinor,
1841 const char* pszCallingComponentName)
1842 {
1843 if (nVersionMajor == GDAL_VERSION_MAJOR &&
1844 nVersionMinor == GDAL_VERSION_MINOR)
1845 return TRUE;
1846
1847 if (pszCallingComponentName)
1848 {
1849 CPLError( CE_Failure, CPLE_AppDefined,
1850 "%s was compiled against GDAL %d.%d but current library version is %d.%d\n",
1851 pszCallingComponentName, nVersionMajor, nVersionMinor,
1852 GDAL_VERSION_MAJOR, GDAL_VERSION_MINOR);
1853 }
1854 return FALSE;
1855 }
1856
1857 /************************************************************************/
1858 /* GDALDecToDMS() */
1859 /* */
1860 /* Translate a decimal degrees value to a DMS string with */
1861 /* hemisphere. */
1862 /************************************************************************/
1863
GDALDecToDMS(double dfAngle,const char * pszAxis,int nPrecision)1864 const char * CPL_STDCALL GDALDecToDMS( double dfAngle, const char * pszAxis,
1865 int nPrecision )
1866
1867 {
1868 return CPLDecToDMS( dfAngle, pszAxis, nPrecision );
1869 }
1870
1871 /************************************************************************/
1872 /* GDALPackedDMSToDec() */
1873 /************************************************************************/
1874
1875 /**
1876 * \brief Convert a packed DMS value (DDDMMMSSS.SS) into decimal degrees.
1877 *
1878 * See CPLPackedDMSToDec().
1879 */
1880
GDALPackedDMSToDec(double dfPacked)1881 double CPL_STDCALL GDALPackedDMSToDec( double dfPacked )
1882
1883 {
1884 return CPLPackedDMSToDec( dfPacked );
1885 }
1886
1887 /************************************************************************/
1888 /* GDALDecToPackedDMS() */
1889 /************************************************************************/
1890
1891 /**
1892 * \brief Convert decimal degrees into packed DMS value (DDDMMMSSS.SS).
1893 *
1894 * See CPLDecToPackedDMS().
1895 */
1896
GDALDecToPackedDMS(double dfDec)1897 double CPL_STDCALL GDALDecToPackedDMS( double dfDec )
1898
1899 {
1900 return CPLDecToPackedDMS( dfDec );
1901 }
1902
1903 /************************************************************************/
1904 /* GDALGCPsToGeoTransform() */
1905 /************************************************************************/
1906
1907 /**
1908 * \brief Generate Geotransform from GCPs.
1909 *
1910 * Given a set of GCPs perform first order fit as a geotransform.
1911 *
1912 * Due to imprecision in the calculations the fit algorithm will often
1913 * return non-zero rotational coefficients even if given perfectly non-rotated
1914 * inputs. A special case has been implemented for corner corner coordinates
1915 * given in TL, TR, BR, BL order. So when using this to get a geotransform
1916 * from 4 corner coordinates, pass them in this order.
1917 *
1918 * @param nGCPCount the number of GCPs being passed in.
1919 * @param pasGCPs the list of GCP structures.
1920 * @param padfGeoTransform the six double array in which the affine
1921 * geotransformation will be returned.
1922 * @param bApproxOK If FALSE the function will fail if the geotransform is not
1923 * essentially an exact fit (within 0.25 pixel) for all GCPs.
1924 *
1925 * @return TRUE on success or FALSE if there aren't enough points to prepare a
1926 * geotransform, the pointers are ill-determined or if bApproxOK is FALSE
1927 * and the fit is poor.
1928 */
1929
1930 int CPL_STDCALL
GDALGCPsToGeoTransform(int nGCPCount,const GDAL_GCP * pasGCPs,double * padfGeoTransform,int bApproxOK)1931 GDALGCPsToGeoTransform( int nGCPCount, const GDAL_GCP *pasGCPs,
1932 double *padfGeoTransform, int bApproxOK )
1933
1934 {
1935 int i;
1936
1937 /* -------------------------------------------------------------------- */
1938 /* Recognise a few special cases. */
1939 /* -------------------------------------------------------------------- */
1940 if( nGCPCount < 2 )
1941 return FALSE;
1942
1943 if( nGCPCount == 2 )
1944 {
1945 if( pasGCPs[1].dfGCPPixel == pasGCPs[0].dfGCPPixel
1946 || pasGCPs[1].dfGCPLine == pasGCPs[0].dfGCPLine )
1947 return FALSE;
1948
1949 padfGeoTransform[1] = (pasGCPs[1].dfGCPX - pasGCPs[0].dfGCPX)
1950 / (pasGCPs[1].dfGCPPixel - pasGCPs[0].dfGCPPixel);
1951 padfGeoTransform[2] = 0.0;
1952
1953 padfGeoTransform[4] = 0.0;
1954 padfGeoTransform[5] = (pasGCPs[1].dfGCPY - pasGCPs[0].dfGCPY)
1955 / (pasGCPs[1].dfGCPLine - pasGCPs[0].dfGCPLine);
1956
1957 padfGeoTransform[0] = pasGCPs[0].dfGCPX
1958 - pasGCPs[0].dfGCPPixel * padfGeoTransform[1]
1959 - pasGCPs[0].dfGCPLine * padfGeoTransform[2];
1960
1961 padfGeoTransform[3] = pasGCPs[0].dfGCPY
1962 - pasGCPs[0].dfGCPPixel * padfGeoTransform[4]
1963 - pasGCPs[0].dfGCPLine * padfGeoTransform[5];
1964
1965 return TRUE;
1966 }
1967
1968 /* -------------------------------------------------------------------- */
1969 /* Special case of 4 corner coordinates of a non-rotated */
1970 /* image. The points must be in TL-TR-BR-BL order for now. */
1971 /* This case helps avoid some imprecision in the general */
1972 /* calcuations. */
1973 /* -------------------------------------------------------------------- */
1974 if( nGCPCount == 4
1975 && pasGCPs[0].dfGCPLine == pasGCPs[1].dfGCPLine
1976 && pasGCPs[2].dfGCPLine == pasGCPs[3].dfGCPLine
1977 && pasGCPs[0].dfGCPPixel == pasGCPs[3].dfGCPPixel
1978 && pasGCPs[1].dfGCPPixel == pasGCPs[2].dfGCPPixel
1979 && pasGCPs[0].dfGCPLine != pasGCPs[2].dfGCPLine
1980 && pasGCPs[0].dfGCPPixel != pasGCPs[1].dfGCPPixel
1981 && pasGCPs[0].dfGCPY == pasGCPs[1].dfGCPY
1982 && pasGCPs[2].dfGCPY == pasGCPs[3].dfGCPY
1983 && pasGCPs[0].dfGCPX == pasGCPs[3].dfGCPX
1984 && pasGCPs[1].dfGCPX == pasGCPs[2].dfGCPX
1985 && pasGCPs[0].dfGCPY != pasGCPs[2].dfGCPY
1986 && pasGCPs[0].dfGCPX != pasGCPs[1].dfGCPX )
1987 {
1988 padfGeoTransform[1] = (pasGCPs[1].dfGCPX - pasGCPs[0].dfGCPX)
1989 / (pasGCPs[1].dfGCPPixel - pasGCPs[0].dfGCPPixel);
1990 padfGeoTransform[2] = 0.0;
1991 padfGeoTransform[4] = 0.0;
1992 padfGeoTransform[5] = (pasGCPs[2].dfGCPY - pasGCPs[1].dfGCPY)
1993 / (pasGCPs[2].dfGCPLine - pasGCPs[1].dfGCPLine);
1994
1995 padfGeoTransform[0] =
1996 pasGCPs[0].dfGCPX - pasGCPs[0].dfGCPPixel * padfGeoTransform[1];
1997 padfGeoTransform[3] =
1998 pasGCPs[0].dfGCPY - pasGCPs[0].dfGCPLine * padfGeoTransform[5];
1999 return TRUE;
2000 }
2001
2002 /* -------------------------------------------------------------------- */
2003 /* Compute source and destination ranges so we can normalize */
2004 /* the values to make the least squares computation more stable. */
2005 /* -------------------------------------------------------------------- */
2006 double min_pixel = pasGCPs[0].dfGCPPixel;
2007 double max_pixel = pasGCPs[0].dfGCPPixel;
2008 double min_line = pasGCPs[0].dfGCPLine;
2009 double max_line = pasGCPs[0].dfGCPLine;
2010 double min_geox = pasGCPs[0].dfGCPX;
2011 double max_geox = pasGCPs[0].dfGCPX;
2012 double min_geoy = pasGCPs[0].dfGCPY;
2013 double max_geoy = pasGCPs[0].dfGCPY;
2014
2015 for (i = 1; i < nGCPCount; ++i) {
2016 min_pixel = MIN(min_pixel, pasGCPs[i].dfGCPPixel);
2017 max_pixel = MAX(max_pixel, pasGCPs[i].dfGCPPixel);
2018 min_line = MIN(min_line, pasGCPs[i].dfGCPLine);
2019 max_line = MAX(max_line, pasGCPs[i].dfGCPLine);
2020 min_geox = MIN(min_geox, pasGCPs[i].dfGCPX);
2021 max_geox = MAX(max_geox, pasGCPs[i].dfGCPX);
2022 min_geoy = MIN(min_geoy, pasGCPs[i].dfGCPY);
2023 max_geoy = MAX(max_geoy, pasGCPs[i].dfGCPY);
2024 }
2025
2026 double EPS = 1.0e-12;
2027
2028 if( ABS(max_pixel - min_pixel) < EPS
2029 || ABS(max_line - min_line) < EPS
2030 || ABS(max_geox - min_geox) < EPS
2031 || ABS(max_geoy - min_geoy) < EPS)
2032 {
2033 return FALSE; // degenerate in at least one dimension.
2034 }
2035
2036 double pl_normalize[6], geo_normalize[6];
2037
2038 pl_normalize[0] = -min_pixel / (max_pixel - min_pixel);
2039 pl_normalize[1] = 1.0 / (max_pixel - min_pixel);
2040 pl_normalize[2] = 0.0;
2041 pl_normalize[3] = -min_line / (max_line - min_line);
2042 pl_normalize[4] = 0.0;
2043 pl_normalize[5] = 1.0 / (max_line - min_line);
2044
2045 geo_normalize[0] = -min_geox / (max_geox - min_geox);
2046 geo_normalize[1] = 1.0 / (max_geox - min_geox);
2047 geo_normalize[2] = 0.0;
2048 geo_normalize[3] = -min_geoy / (max_geoy - min_geoy);
2049 geo_normalize[4] = 0.0;
2050 geo_normalize[5] = 1.0 / (max_geoy - min_geoy);
2051
2052 /* -------------------------------------------------------------------- */
2053 /* In the general case, do a least squares error approximation by */
2054 /* solving the equation Sum[(A - B*x + C*y - Lon)^2] = minimum */
2055 /* -------------------------------------------------------------------- */
2056
2057 double sum_x = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_xx = 0.0, sum_yy = 0.0;
2058 double sum_Lon = 0.0, sum_Lonx = 0.0, sum_Lony = 0.0;
2059 double sum_Lat = 0.0, sum_Latx = 0.0, sum_Laty = 0.0;
2060 double divisor;
2061
2062 for (i = 0; i < nGCPCount; ++i) {
2063 double pixel, line, geox, geoy;
2064
2065 GDALApplyGeoTransform(pl_normalize,
2066 pasGCPs[i].dfGCPPixel,
2067 pasGCPs[i].dfGCPLine,
2068 &pixel, &line);
2069 GDALApplyGeoTransform(geo_normalize,
2070 pasGCPs[i].dfGCPX,
2071 pasGCPs[i].dfGCPY,
2072 &geox, &geoy);
2073
2074 sum_x += pixel;
2075 sum_y += line;
2076 sum_xy += pixel * line;
2077 sum_xx += pixel * pixel;
2078 sum_yy += line * line;
2079 sum_Lon += geox;
2080 sum_Lonx += geox * pixel;
2081 sum_Lony += geox * line;
2082 sum_Lat += geoy;
2083 sum_Latx += geoy * pixel;
2084 sum_Laty += geoy * line;
2085 }
2086
2087 divisor = nGCPCount * (sum_xx * sum_yy - sum_xy * sum_xy)
2088 + 2 * sum_x * sum_y * sum_xy - sum_y * sum_y * sum_xx
2089 - sum_x * sum_x * sum_yy;
2090
2091 /* -------------------------------------------------------------------- */
2092 /* If the divisor is zero, there is no valid solution. */
2093 /* -------------------------------------------------------------------- */
2094 if (divisor == 0.0)
2095 return FALSE;
2096
2097 /* -------------------------------------------------------------------- */
2098 /* Compute top/left origin. */
2099 /* -------------------------------------------------------------------- */
2100 double gt_normalized[6];
2101 gt_normalized[0] = (sum_Lon * (sum_xx * sum_yy - sum_xy * sum_xy)
2102 + sum_Lonx * (sum_y * sum_xy - sum_x * sum_yy)
2103 + sum_Lony * (sum_x * sum_xy - sum_y * sum_xx))
2104 / divisor;
2105
2106 gt_normalized[3] = (sum_Lat * (sum_xx * sum_yy - sum_xy * sum_xy)
2107 + sum_Latx * (sum_y * sum_xy - sum_x * sum_yy)
2108 + sum_Laty * (sum_x * sum_xy - sum_y * sum_xx))
2109 / divisor;
2110
2111 /* -------------------------------------------------------------------- */
2112 /* Compute X related coefficients. */
2113 /* -------------------------------------------------------------------- */
2114 gt_normalized[1] = (sum_Lon * (sum_y * sum_xy - sum_x * sum_yy)
2115 + sum_Lonx * (nGCPCount * sum_yy - sum_y * sum_y)
2116 + sum_Lony * (sum_x * sum_y - sum_xy * nGCPCount))
2117 / divisor;
2118
2119 gt_normalized[2] = (sum_Lon * (sum_x * sum_xy - sum_y * sum_xx)
2120 + sum_Lonx * (sum_x * sum_y - nGCPCount * sum_xy)
2121 + sum_Lony * (nGCPCount * sum_xx - sum_x * sum_x))
2122 / divisor;
2123
2124 /* -------------------------------------------------------------------- */
2125 /* Compute Y related coefficients. */
2126 /* -------------------------------------------------------------------- */
2127 gt_normalized[4] = (sum_Lat * (sum_y * sum_xy - sum_x * sum_yy)
2128 + sum_Latx * (nGCPCount * sum_yy - sum_y * sum_y)
2129 + sum_Laty * (sum_x * sum_y - sum_xy * nGCPCount))
2130 / divisor;
2131
2132 gt_normalized[5] = (sum_Lat * (sum_x * sum_xy - sum_y * sum_xx)
2133 + sum_Latx * (sum_x * sum_y - nGCPCount * sum_xy)
2134 + sum_Laty * (nGCPCount * sum_xx - sum_x * sum_x))
2135 / divisor;
2136
2137 /* -------------------------------------------------------------------- */
2138 /* Compose the resulting transformation with the normalization */
2139 /* geotransformations. */
2140 /* -------------------------------------------------------------------- */
2141 double gt1p2[6], inv_geo_normalize[6];
2142 if( !GDALInvGeoTransform(geo_normalize, inv_geo_normalize))
2143 return FALSE;
2144
2145 GDALComposeGeoTransforms(pl_normalize, gt_normalized, gt1p2);
2146 GDALComposeGeoTransforms(gt1p2, inv_geo_normalize, padfGeoTransform);
2147
2148 /* -------------------------------------------------------------------- */
2149 /* Now check if any of the input points fit this poorly. */
2150 /* -------------------------------------------------------------------- */
2151 if( !bApproxOK )
2152 {
2153 // FIXME? Not sure if it is the more accurate way of computing
2154 // pixel size
2155 double dfPixelSize = 0.5 * (ABS(padfGeoTransform[1])
2156 + ABS(padfGeoTransform[2])
2157 + ABS(padfGeoTransform[4])
2158 + ABS(padfGeoTransform[5]));
2159
2160 for( i = 0; i < nGCPCount; i++ )
2161 {
2162 double dfErrorX, dfErrorY;
2163
2164 dfErrorX =
2165 (pasGCPs[i].dfGCPPixel * padfGeoTransform[1]
2166 + pasGCPs[i].dfGCPLine * padfGeoTransform[2]
2167 + padfGeoTransform[0])
2168 - pasGCPs[i].dfGCPX;
2169 dfErrorY =
2170 (pasGCPs[i].dfGCPPixel * padfGeoTransform[4]
2171 + pasGCPs[i].dfGCPLine * padfGeoTransform[5]
2172 + padfGeoTransform[3])
2173 - pasGCPs[i].dfGCPY;
2174
2175 if( ABS(dfErrorX) > 0.25 * dfPixelSize
2176 || ABS(dfErrorY) > 0.25 * dfPixelSize )
2177 {
2178 CPLDebug("GDAL", "dfErrorX/dfPixelSize = %.2f, dfErrorY/dfPixelSize = %.2f",
2179 ABS(dfErrorX)/dfPixelSize, ABS(dfErrorY)/dfPixelSize);
2180 return FALSE;
2181 }
2182 }
2183 }
2184
2185 return TRUE;
2186 }
2187
2188 /************************************************************************/
2189 /* GDALComposeGeoTransforms() */
2190 /************************************************************************/
2191
2192 /**
2193 * \brief Compose two geotransforms.
2194 *
2195 * The resulting geotransform is the equivelent to padfGT1 and then padfGT2
2196 * being applied to a point.
2197 *
2198 * @param padfGT1 the first geotransform, six values.
2199 * @param padfGT2 the second geotransform, six values.
2200 * @param padfGTOut the output geotransform, six values, may safely be the same
2201 * array as padfGT1 or padfGT2.
2202 */
2203
GDALComposeGeoTransforms(const double * padfGT1,const double * padfGT2,double * padfGTOut)2204 void GDALComposeGeoTransforms(const double *padfGT1, const double *padfGT2,
2205 double *padfGTOut)
2206
2207 {
2208 double gtwrk[6];
2209 // We need to think of the geotransform in a more normal form to do
2210 // the matrix multiple:
2211 //
2212 // __ __
2213 // | gt[1] gt[2] gt[0] |
2214 // | gt[4] gt[5] gt[3] |
2215 // | 0.0 0.0 1.0 |
2216 // -- --
2217 //
2218 // Then we can use normal matrix multiplication to produce the
2219 // composed transformation. I don't actually reform the matrix
2220 // explicitly which is why the following may seem kind of spagettish.
2221
2222 gtwrk[1] =
2223 padfGT2[1] * padfGT1[1]
2224 + padfGT2[2] * padfGT1[4];
2225 gtwrk[2] =
2226 padfGT2[1] * padfGT1[2]
2227 + padfGT2[2] * padfGT1[5];
2228 gtwrk[0] =
2229 padfGT2[1] * padfGT1[0]
2230 + padfGT2[2] * padfGT1[3]
2231 + padfGT2[0] * 1.0;
2232
2233 gtwrk[4] =
2234 padfGT2[4] * padfGT1[1]
2235 + padfGT2[5] * padfGT1[4];
2236 gtwrk[5] =
2237 padfGT2[4] * padfGT1[2]
2238 + padfGT2[5] * padfGT1[5];
2239 gtwrk[3] =
2240 padfGT2[4] * padfGT1[0]
2241 + padfGT2[5] * padfGT1[3]
2242 + padfGT2[3] * 1.0;
2243 memcpy(padfGTOut, gtwrk, sizeof(double) * 6);
2244 }
2245
2246 /************************************************************************/
2247 /* GDALGeneralCmdLineProcessor() */
2248 /************************************************************************/
2249
2250 /**
2251 * \brief General utility option processing.
2252 *
2253 * This function is intended to provide a variety of generic commandline
2254 * options for all GDAL commandline utilities. It takes care of the following
2255 * commandline options:
2256 *
2257 * --version: report version of GDAL in use.
2258 * --build: report build info about GDAL in use.
2259 * --license: report GDAL license info.
2260 * --formats: report all format drivers configured.
2261 * --format [format]: report details of one format driver.
2262 * --optfile filename: expand an option file into the argument list.
2263 * --config key value: set system configuration option.
2264 * --debug [on/off/value]: set debug level.
2265 * --mempreload dir: preload directory contents into /vsimem
2266 * --pause: Pause for user input (allows time to attach debugger)
2267 * --locale [locale]: Install a locale using setlocale() (debugging)
2268 * --help-general: report detailed help on general options.
2269 *
2270 * The argument array is replaced "in place" and should be freed with
2271 * CSLDestroy() when no longer needed. The typical usage looks something
2272 * like the following. Note that the formats should be registered so that
2273 * the --formats and --format options will work properly.
2274 *
2275 * int main( int argc, char ** argv )
2276 * {
2277 * GDALAllRegister();
2278 *
2279 * argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
2280 * if( argc < 1 )
2281 * exit( -argc );
2282 *
2283 * @param nArgc number of values in the argument list.
2284 * @param ppapszArgv pointer to the argument list array (will be updated in place).
2285 * @param nOptions a or-able combination of GDAL_OF_RASTER and GDAL_OF_VECTOR
2286 * to determine which drivers should be displayed by --formats.
2287 * If set to 0, GDAL_OF_RASTER is assumed.
2288 *
2289 * @return updated nArgc argument count. Return of 0 requests terminate
2290 * without error, return of -1 requests exit with error code.
2291 */
2292
2293 int CPL_STDCALL
GDALGeneralCmdLineProcessor(int nArgc,char *** ppapszArgv,int nOptions)2294 GDALGeneralCmdLineProcessor( int nArgc, char ***ppapszArgv, int nOptions )
2295
2296 {
2297 char **papszReturn = NULL;
2298 int iArg;
2299 char **papszArgv = *ppapszArgv;
2300
2301 /* -------------------------------------------------------------------- */
2302 /* Preserve the program name. */
2303 /* -------------------------------------------------------------------- */
2304 papszReturn = CSLAddString( papszReturn, papszArgv[0] );
2305
2306 /* ==================================================================== */
2307 /* Loop over all arguments. */
2308 /* ==================================================================== */
2309 for( iArg = 1; iArg < nArgc; iArg++ )
2310 {
2311 /* -------------------------------------------------------------------- */
2312 /* --version */
2313 /* -------------------------------------------------------------------- */
2314 if( EQUAL(papszArgv[iArg],"--version") )
2315 {
2316 printf( "%s\n", GDALVersionInfo( "--version" ) );
2317 CSLDestroy( papszReturn );
2318 return 0;
2319 }
2320
2321 /* -------------------------------------------------------------------- */
2322 /* --build */
2323 /* -------------------------------------------------------------------- */
2324 else if( EQUAL(papszArgv[iArg],"--build") )
2325 {
2326 printf( "%s", GDALVersionInfo( "BUILD_INFO" ) );
2327 CSLDestroy( papszReturn );
2328 return 0;
2329 }
2330
2331 /* -------------------------------------------------------------------- */
2332 /* --license */
2333 /* -------------------------------------------------------------------- */
2334 else if( EQUAL(papszArgv[iArg],"--license") )
2335 {
2336 printf( "%s\n", GDALVersionInfo( "LICENSE" ) );
2337 CSLDestroy( papszReturn );
2338 return 0;
2339 }
2340
2341 /* -------------------------------------------------------------------- */
2342 /* --config */
2343 /* -------------------------------------------------------------------- */
2344 else if( EQUAL(papszArgv[iArg],"--config") )
2345 {
2346 if( iArg + 2 >= nArgc )
2347 {
2348 CPLError( CE_Failure, CPLE_AppDefined,
2349 "--config option given without a key and value argument." );
2350 CSLDestroy( papszReturn );
2351 return -1;
2352 }
2353
2354 CPLSetConfigOption( papszArgv[iArg+1], papszArgv[iArg+2] );
2355
2356 iArg += 2;
2357 }
2358
2359 /* -------------------------------------------------------------------- */
2360 /* --mempreload */
2361 /* -------------------------------------------------------------------- */
2362 else if( EQUAL(papszArgv[iArg],"--mempreload") )
2363 {
2364 int i;
2365
2366 if( iArg + 1 >= nArgc )
2367 {
2368 CPLError( CE_Failure, CPLE_AppDefined,
2369 "--mempreload option given without directory path.");
2370 CSLDestroy( papszReturn );
2371 return -1;
2372 }
2373
2374 char **papszFiles = CPLReadDir( papszArgv[iArg+1] );
2375 if( CSLCount(papszFiles) == 0 )
2376 {
2377 CPLError( CE_Failure, CPLE_AppDefined,
2378 "--mempreload given invalid or empty directory.");
2379 CSLDestroy( papszReturn );
2380 return -1;
2381 }
2382
2383 for( i = 0; papszFiles[i] != NULL; i++ )
2384 {
2385 CPLString osOldPath, osNewPath;
2386 VSIStatBufL sStatBuf;
2387
2388 if( EQUAL(papszFiles[i],".") || EQUAL(papszFiles[i],"..") )
2389 continue;
2390
2391 osOldPath = CPLFormFilename( papszArgv[iArg+1],
2392 papszFiles[i], NULL );
2393 osNewPath.Printf( "/vsimem/%s", papszFiles[i] );
2394
2395 if( VSIStatL( osOldPath, &sStatBuf ) != 0
2396 || VSI_ISDIR( sStatBuf.st_mode ) )
2397 {
2398 CPLDebug( "VSI", "Skipping preload of %s.",
2399 osOldPath.c_str() );
2400 continue;
2401 }
2402
2403 CPLDebug( "VSI", "Preloading %s to %s.",
2404 osOldPath.c_str(), osNewPath.c_str() );
2405
2406 if( CPLCopyFile( osNewPath, osOldPath ) != 0 )
2407 {
2408 CPLError( CE_Failure, CPLE_AppDefined,
2409 "Failed to copy %s to /vsimem",
2410 osOldPath.c_str() );
2411 return -1;
2412 }
2413 }
2414
2415 CSLDestroy( papszFiles );
2416 iArg += 1;
2417 }
2418
2419 /* -------------------------------------------------------------------- */
2420 /* --debug */
2421 /* -------------------------------------------------------------------- */
2422 else if( EQUAL(papszArgv[iArg],"--debug") )
2423 {
2424 if( iArg + 1 >= nArgc )
2425 {
2426 CPLError( CE_Failure, CPLE_AppDefined,
2427 "--debug option given without debug level." );
2428 CSLDestroy( papszReturn );
2429 return -1;
2430 }
2431
2432 CPLSetConfigOption( "CPL_DEBUG", papszArgv[iArg+1] );
2433 iArg += 1;
2434 }
2435
2436 /* -------------------------------------------------------------------- */
2437 /* --optfile */
2438 /* */
2439 /* Annoyingly the options inserted by --optfile will *not* be */
2440 /* processed properly if they are general options. */
2441 /* -------------------------------------------------------------------- */
2442 else if( EQUAL(papszArgv[iArg],"--optfile") )
2443 {
2444 const char *pszLine;
2445 FILE *fpOptFile;
2446
2447 if( iArg + 1 >= nArgc )
2448 {
2449 CPLError( CE_Failure, CPLE_AppDefined,
2450 "--optfile option given without filename." );
2451 CSLDestroy( papszReturn );
2452 return -1;
2453 }
2454
2455 fpOptFile = VSIFOpen( papszArgv[iArg+1], "rb" );
2456
2457 if( fpOptFile == NULL )
2458 {
2459 CPLError( CE_Failure, CPLE_AppDefined,
2460 "Unable to open optfile '%s'.\n%s",
2461 papszArgv[iArg+1], VSIStrerror( errno ) );
2462 CSLDestroy( papszReturn );
2463 return -1;
2464 }
2465
2466 while( (pszLine = CPLReadLine( fpOptFile )) != NULL )
2467 {
2468 char **papszTokens;
2469 int i;
2470
2471 if( pszLine[0] == '#' || strlen(pszLine) == 0 )
2472 continue;
2473
2474 papszTokens = CSLTokenizeString( pszLine );
2475 for( i = 0; papszTokens != NULL && papszTokens[i] != NULL; i++)
2476 papszReturn = CSLAddString( papszReturn, papszTokens[i] );
2477 CSLDestroy( papszTokens );
2478 }
2479
2480 VSIFClose( fpOptFile );
2481
2482 iArg += 1;
2483 }
2484
2485 /* -------------------------------------------------------------------- */
2486 /* --formats */
2487 /* -------------------------------------------------------------------- */
2488 else if( EQUAL(papszArgv[iArg], "--formats") )
2489 {
2490 int iDr;
2491
2492 if( nOptions == 0 )
2493 nOptions = GDAL_OF_RASTER;
2494
2495 printf( "Supported Formats:\n" );
2496 for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
2497 {
2498 GDALDriverH hDriver = GDALGetDriver(iDr);
2499
2500 const char *pszRFlag = "", *pszWFlag, *pszVirtualIO, *pszSubdatasets, *pszKind;
2501 char** papszMD = GDALGetMetadata( hDriver, NULL );
2502
2503 if( nOptions == GDAL_OF_RASTER &&
2504 !CSLFetchBoolean( papszMD, GDAL_DCAP_RASTER, FALSE ) )
2505 continue;
2506 if( nOptions == GDAL_OF_VECTOR &&
2507 !CSLFetchBoolean( papszMD, GDAL_DCAP_VECTOR, FALSE ) )
2508 continue;
2509
2510 if( CSLFetchBoolean( papszMD, GDAL_DCAP_OPEN, FALSE ) )
2511 pszRFlag = "r";
2512
2513 if( CSLFetchBoolean( papszMD, GDAL_DCAP_CREATE, FALSE ) )
2514 pszWFlag = "w+";
2515 else if( CSLFetchBoolean( papszMD, GDAL_DCAP_CREATECOPY, FALSE ) )
2516 pszWFlag = "w";
2517 else
2518 pszWFlag = "o";
2519
2520 if( CSLFetchBoolean( papszMD, GDAL_DCAP_VIRTUALIO, FALSE ) )
2521 pszVirtualIO = "v";
2522 else
2523 pszVirtualIO = "";
2524
2525 if( CSLFetchBoolean( papszMD, GDAL_DMD_SUBDATASETS, FALSE ) )
2526 pszSubdatasets = "s";
2527 else
2528 pszSubdatasets = "";
2529
2530 if( CSLFetchBoolean( papszMD, GDAL_DCAP_RASTER, FALSE ) &&
2531 CSLFetchBoolean( papszMD, GDAL_DCAP_VECTOR, FALSE ))
2532 pszKind = "raster,vector";
2533 else if( CSLFetchBoolean( papszMD, GDAL_DCAP_RASTER, FALSE ) )
2534 pszKind = "raster";
2535 else if( CSLFetchBoolean( papszMD, GDAL_DCAP_VECTOR, FALSE ) )
2536 pszKind = "vector";
2537 else
2538 pszKind = "unknown kind";
2539
2540 printf( " %s -%s- (%s%s%s%s): %s\n",
2541 GDALGetDriverShortName( hDriver ),
2542 pszKind,
2543 pszRFlag, pszWFlag, pszVirtualIO, pszSubdatasets,
2544 GDALGetDriverLongName( hDriver ) );
2545 }
2546
2547 CSLDestroy( papszReturn );
2548 return 0;
2549 }
2550
2551 /* -------------------------------------------------------------------- */
2552 /* --format */
2553 /* -------------------------------------------------------------------- */
2554 else if( EQUAL(papszArgv[iArg], "--format") )
2555 {
2556 GDALDriverH hDriver;
2557 char **papszMD;
2558
2559 CSLDestroy( papszReturn );
2560
2561 if( iArg + 1 >= nArgc )
2562 {
2563 CPLError( CE_Failure, CPLE_AppDefined,
2564 "--format option given without a format code." );
2565 return -1;
2566 }
2567
2568 hDriver = GDALGetDriverByName( papszArgv[iArg+1] );
2569 if( hDriver == NULL )
2570 {
2571 CPLError( CE_Failure, CPLE_AppDefined,
2572 "--format option given with format '%s', but that format not\n"
2573 "recognised. Use the --formats option to get a list of available formats,\n"
2574 "and use the short code (ie. GTiff or HFA) as the format identifier.\n",
2575 papszArgv[iArg+1] );
2576 return -1;
2577 }
2578
2579 printf( "Format Details:\n" );
2580 printf( " Short Name: %s\n", GDALGetDriverShortName( hDriver ) );
2581 printf( " Long Name: %s\n", GDALGetDriverLongName( hDriver ) );
2582
2583 papszMD = GDALGetMetadata( hDriver, NULL );
2584 if( CSLFetchBoolean( papszMD, GDAL_DCAP_RASTER, FALSE ) )
2585 printf( " Supports: Raster\n" );
2586 if( CSLFetchBoolean( papszMD, GDAL_DCAP_VECTOR, FALSE ) )
2587 printf( " Supports: Vector\n" );
2588
2589 const char* pszExt = CSLFetchNameValue( papszMD, GDAL_DMD_EXTENSIONS );
2590 if( pszExt != NULL )
2591 printf( " Extension%s: %s\n", (strchr(pszExt, ' ') ? "s" : ""),
2592 pszExt );
2593
2594 if( CSLFetchNameValue( papszMD, GDAL_DMD_MIMETYPE ) )
2595 printf( " Mime Type: %s\n",
2596 CSLFetchNameValue( papszMD, GDAL_DMD_MIMETYPE ) );
2597 if( CSLFetchNameValue( papszMD, GDAL_DMD_HELPTOPIC ) )
2598 printf( " Help Topic: %s\n",
2599 CSLFetchNameValue( papszMD, GDAL_DMD_HELPTOPIC ) );
2600
2601 if( CSLFetchBoolean( papszMD, GDAL_DMD_SUBDATASETS, FALSE ) )
2602 printf( " Supports: Subdatasets\n" );
2603 if( CSLFetchBoolean( papszMD, GDAL_DCAP_OPEN, FALSE ) )
2604 printf( " Supports: Open() - Open existing dataset.\n" );
2605 if( CSLFetchBoolean( papszMD, GDAL_DCAP_CREATE, FALSE ) )
2606 printf( " Supports: Create() - Create writeable dataset.\n" );
2607 if( CSLFetchBoolean( papszMD, GDAL_DCAP_CREATECOPY, FALSE ) )
2608 printf( " Supports: CreateCopy() - Create dataset by copying another.\n" );
2609 if( CSLFetchBoolean( papszMD, GDAL_DCAP_VIRTUALIO, FALSE ) )
2610 printf( " Supports: Virtual IO - eg. /vsimem/\n" );
2611 if( CSLFetchNameValue( papszMD, GDAL_DMD_CREATIONDATATYPES ) )
2612 printf( " Creation Datatypes: %s\n",
2613 CSLFetchNameValue( papszMD, GDAL_DMD_CREATIONDATATYPES ) );
2614 if( CSLFetchNameValue( papszMD, GDAL_DMD_CREATIONFIELDDATATYPES ) )
2615 printf( " Creation Field Datatypes: %s\n",
2616 CSLFetchNameValue( papszMD, GDAL_DMD_CREATIONFIELDDATATYPES ) );
2617 if( CSLFetchBoolean( papszMD, GDAL_DCAP_NOTNULL_FIELDS, FALSE ) )
2618 printf( " Supports: Creating fields with NOT NULL constraint.\n" );
2619 if( CSLFetchBoolean( papszMD, GDAL_DCAP_DEFAULT_FIELDS, FALSE ) )
2620 printf( " Supports: Creating fields with DEFAULT values.\n" );
2621 if( CSLFetchBoolean( papszMD, GDAL_DCAP_NOTNULL_GEOMFIELDS, FALSE ) )
2622 printf( " Supports: Creating geometry fields with NOT NULL constraint.\n" );
2623 if( CSLFetchNameValue( papszMD, GDAL_DMD_CREATIONOPTIONLIST ) )
2624 {
2625 CPLXMLNode *psCOL =
2626 CPLParseXMLString(
2627 CSLFetchNameValue( papszMD,
2628 GDAL_DMD_CREATIONOPTIONLIST ) );
2629 char *pszFormattedXML =
2630 CPLSerializeXMLTree( psCOL );
2631
2632 CPLDestroyXMLNode( psCOL );
2633
2634 printf( "\n%s\n", pszFormattedXML );
2635 CPLFree( pszFormattedXML );
2636 }
2637 if( CSLFetchNameValue( papszMD, GDAL_DS_LAYER_CREATIONOPTIONLIST ) )
2638 {
2639 CPLXMLNode *psCOL =
2640 CPLParseXMLString(
2641 CSLFetchNameValue( papszMD,
2642 GDAL_DS_LAYER_CREATIONOPTIONLIST ) );
2643 char *pszFormattedXML =
2644 CPLSerializeXMLTree( psCOL );
2645
2646 CPLDestroyXMLNode( psCOL );
2647
2648 printf( "\n%s\n", pszFormattedXML );
2649 CPLFree( pszFormattedXML );
2650 }
2651
2652 if( CSLFetchNameValue( papszMD, GDAL_DMD_CONNECTION_PREFIX ) )
2653 printf( " Connection prefix: %s\n",
2654 CSLFetchNameValue( papszMD, GDAL_DMD_CONNECTION_PREFIX ) );
2655
2656 if( CSLFetchNameValue( papszMD, GDAL_DMD_OPENOPTIONLIST ) )
2657 {
2658 CPLXMLNode *psCOL =
2659 CPLParseXMLString(
2660 CSLFetchNameValue( papszMD,
2661 GDAL_DMD_OPENOPTIONLIST ) );
2662 char *pszFormattedXML =
2663 CPLSerializeXMLTree( psCOL );
2664
2665 CPLDestroyXMLNode( psCOL );
2666
2667 printf( "%s\n", pszFormattedXML );
2668 CPLFree( pszFormattedXML );
2669 }
2670 return 0;
2671 }
2672
2673 /* -------------------------------------------------------------------- */
2674 /* --help-general */
2675 /* -------------------------------------------------------------------- */
2676 else if( EQUAL(papszArgv[iArg],"--help-general") )
2677 {
2678 printf( "Generic GDAL utility command options:\n" );
2679 printf( " --version: report version of GDAL in use.\n" );
2680 printf( " --license: report GDAL license info.\n" );
2681 printf( " --formats: report all configured format drivers.\n" );
2682 printf( " --format [format]: details of one format.\n" );
2683 printf( " --optfile filename: expand an option file into the argument list.\n" );
2684 printf( " --config key value: set system configuration option.\n" );
2685 printf( " --debug [on/off/value]: set debug level.\n" );
2686 printf( " --pause: wait for user input, time to attach debugger\n" );
2687 printf( " --locale [locale]: install locale for debugging (ie. en_US.UTF-8)\n" );
2688 printf( " --help-general: report detailed help on general options.\n" );
2689 CSLDestroy( papszReturn );
2690 return 0;
2691 }
2692
2693 /* -------------------------------------------------------------------- */
2694 /* --locale */
2695 /* -------------------------------------------------------------------- */
2696 else if( EQUAL(papszArgv[iArg],"--locale") && iArg < nArgc-1 )
2697 {
2698 CPLsetlocale( LC_ALL, papszArgv[++iArg] );
2699 }
2700
2701 /* -------------------------------------------------------------------- */
2702 /* --pause */
2703 /* -------------------------------------------------------------------- */
2704 else if( EQUAL(papszArgv[iArg],"--pause") )
2705 {
2706 printf( "Hit <ENTER> to Continue.\n" );
2707 CPLReadLine( stdin );
2708 }
2709
2710 /* -------------------------------------------------------------------- */
2711 /* carry through unrecognised options. */
2712 /* -------------------------------------------------------------------- */
2713 else
2714 {
2715 papszReturn = CSLAddString( papszReturn, papszArgv[iArg] );
2716 }
2717 }
2718
2719 *ppapszArgv = papszReturn;
2720
2721 return CSLCount( *ppapszArgv );
2722 }
2723
2724
2725 /************************************************************************/
2726 /* _FetchDblFromMD() */
2727 /************************************************************************/
2728
_FetchDblFromMD(char ** papszMD,const char * pszKey,double * padfTarget,int nCount,double dfDefault)2729 static int _FetchDblFromMD( char **papszMD, const char *pszKey,
2730 double *padfTarget, int nCount, double dfDefault )
2731
2732 {
2733 char szFullKey[200];
2734
2735 sprintf( szFullKey, "%s", pszKey );
2736
2737 const char *pszValue = CSLFetchNameValue( papszMD, szFullKey );
2738 int i;
2739
2740 for( i = 0; i < nCount; i++ )
2741 padfTarget[i] = dfDefault;
2742
2743 if( pszValue == NULL )
2744 return FALSE;
2745
2746 if( nCount == 1 )
2747 {
2748 *padfTarget = CPLAtofM( pszValue );
2749 return TRUE;
2750 }
2751
2752 char **papszTokens = CSLTokenizeStringComplex( pszValue, " ,",
2753 FALSE, FALSE );
2754
2755 if( CSLCount( papszTokens ) != nCount )
2756 {
2757 CSLDestroy( papszTokens );
2758 return FALSE;
2759 }
2760
2761 for( i = 0; i < nCount; i++ )
2762 padfTarget[i] = CPLAtofM(papszTokens[i]);
2763
2764 CSLDestroy( papszTokens );
2765
2766 return TRUE;
2767 }
2768
2769 /************************************************************************/
2770 /* GDALExtractRPCInfo() */
2771 /* */
2772 /* Extract RPC info from metadata, and apply to an RPCInfo */
2773 /* structure. The inverse of this function is RPCInfoToMD() in */
2774 /* alg/gdal_rpc.cpp (should it be needed). */
2775 /************************************************************************/
2776
GDALExtractRPCInfo(char ** papszMD,GDALRPCInfo * psRPC)2777 int CPL_STDCALL GDALExtractRPCInfo( char **papszMD, GDALRPCInfo *psRPC )
2778
2779 {
2780 if( CSLFetchNameValue( papszMD, RPC_LINE_NUM_COEFF ) == NULL )
2781 return FALSE;
2782
2783 if( CSLFetchNameValue( papszMD, RPC_LINE_NUM_COEFF ) == NULL
2784 || CSLFetchNameValue( papszMD, RPC_LINE_DEN_COEFF ) == NULL
2785 || CSLFetchNameValue( papszMD, RPC_SAMP_NUM_COEFF ) == NULL
2786 || CSLFetchNameValue( papszMD, RPC_SAMP_DEN_COEFF ) == NULL )
2787 {
2788 CPLError( CE_Failure, CPLE_AppDefined,
2789 "Some required RPC metadata missing in GDALExtractRPCInfo()");
2790 return FALSE;
2791 }
2792
2793 _FetchDblFromMD( papszMD, RPC_LINE_OFF, &(psRPC->dfLINE_OFF), 1, 0.0 );
2794 _FetchDblFromMD( papszMD, RPC_LINE_SCALE, &(psRPC->dfLINE_SCALE), 1, 1.0 );
2795 _FetchDblFromMD( papszMD, RPC_SAMP_OFF, &(psRPC->dfSAMP_OFF), 1, 0.0 );
2796 _FetchDblFromMD( papszMD, RPC_SAMP_SCALE, &(psRPC->dfSAMP_SCALE), 1, 1.0 );
2797 _FetchDblFromMD( papszMD, RPC_HEIGHT_OFF, &(psRPC->dfHEIGHT_OFF), 1, 0.0 );
2798 _FetchDblFromMD( papszMD, RPC_HEIGHT_SCALE, &(psRPC->dfHEIGHT_SCALE),1, 1.0);
2799 _FetchDblFromMD( papszMD, RPC_LAT_OFF, &(psRPC->dfLAT_OFF), 1, 0.0 );
2800 _FetchDblFromMD( papszMD, RPC_LAT_SCALE, &(psRPC->dfLAT_SCALE), 1, 1.0 );
2801 _FetchDblFromMD( papszMD, RPC_LONG_OFF, &(psRPC->dfLONG_OFF), 1, 0.0 );
2802 _FetchDblFromMD( papszMD, RPC_LONG_SCALE, &(psRPC->dfLONG_SCALE), 1, 1.0 );
2803
2804 _FetchDblFromMD( papszMD, RPC_LINE_NUM_COEFF, psRPC->adfLINE_NUM_COEFF,
2805 20, 0.0 );
2806 _FetchDblFromMD( papszMD, RPC_LINE_DEN_COEFF, psRPC->adfLINE_DEN_COEFF,
2807 20, 0.0 );
2808 _FetchDblFromMD( papszMD, RPC_SAMP_NUM_COEFF, psRPC->adfSAMP_NUM_COEFF,
2809 20, 0.0 );
2810 _FetchDblFromMD( papszMD, RPC_SAMP_DEN_COEFF, psRPC->adfSAMP_DEN_COEFF,
2811 20, 0.0 );
2812
2813 _FetchDblFromMD( papszMD, "MIN_LONG", &(psRPC->dfMIN_LONG), 1, -180.0 );
2814 _FetchDblFromMD( papszMD, "MIN_LAT", &(psRPC->dfMIN_LAT), 1, -90.0 );
2815 _FetchDblFromMD( papszMD, "MAX_LONG", &(psRPC->dfMAX_LONG), 1, 180.0 );
2816 _FetchDblFromMD( papszMD, "MAX_LAT", &(psRPC->dfMAX_LAT), 1, 90.0 );
2817
2818 return TRUE;
2819 }
2820
2821 /************************************************************************/
2822 /* GDALFindAssociatedAuxFile() */
2823 /************************************************************************/
2824
GDALFindAssociatedAuxFile(const char * pszBasename,GDALAccess eAccess,GDALDataset * poDependentDS)2825 GDALDataset *GDALFindAssociatedAuxFile( const char *pszBasename,
2826 GDALAccess eAccess,
2827 GDALDataset *poDependentDS )
2828
2829 {
2830 const char *pszAuxSuffixLC = "aux";
2831 const char *pszAuxSuffixUC = "AUX";
2832
2833 if( EQUAL(CPLGetExtension(pszBasename), pszAuxSuffixLC) )
2834 return NULL;
2835
2836 /* -------------------------------------------------------------------- */
2837 /* Don't even try to look for an .aux file if we don't have a */
2838 /* path of any kind. */
2839 /* -------------------------------------------------------------------- */
2840 if( strlen(pszBasename) == 0 )
2841 return NULL;
2842
2843 /* -------------------------------------------------------------------- */
2844 /* We didn't find that, so try and find a corresponding aux */
2845 /* file. Check that we are the dependent file of the aux */
2846 /* file, or if we aren't verify that the dependent file does */
2847 /* not exist, likely mean it is us but some sort of renaming */
2848 /* has occured. */
2849 /* -------------------------------------------------------------------- */
2850 CPLString osJustFile = CPLGetFilename(pszBasename); // without dir
2851 CPLString osAuxFilename = CPLResetExtension(pszBasename, pszAuxSuffixLC);
2852 GDALDataset *poODS = NULL;
2853 GByte abyHeader[32];
2854 VSILFILE *fp;
2855
2856 fp = VSIFOpenL( osAuxFilename, "rb" );
2857
2858
2859 if ( fp == NULL && VSIIsCaseSensitiveFS(osAuxFilename))
2860 {
2861 // Can't found file with lower case suffix. Try the upper case one.
2862 osAuxFilename = CPLResetExtension(pszBasename, pszAuxSuffixUC);
2863 fp = VSIFOpenL( osAuxFilename, "rb" );
2864 }
2865
2866 if( fp != NULL )
2867 {
2868 if( VSIFReadL( abyHeader, 1, 32, fp ) == 32 &&
2869 EQUALN((char *) abyHeader,"EHFA_HEADER_TAG",15) )
2870 {
2871 /* Avoid causing failure in opening of main file from SWIG bindings */
2872 /* when auxiliary file cannot be opened (#3269) */
2873 CPLTurnFailureIntoWarning(TRUE);
2874 if( poDependentDS != NULL && poDependentDS->GetShared() )
2875 poODS = (GDALDataset *) GDALOpenShared( osAuxFilename, eAccess );
2876 else
2877 poODS = (GDALDataset *) GDALOpen( osAuxFilename, eAccess );
2878 CPLTurnFailureIntoWarning(FALSE);
2879 }
2880 VSIFCloseL( fp );
2881 }
2882
2883 /* -------------------------------------------------------------------- */
2884 /* Try replacing extension with .aux */
2885 /* -------------------------------------------------------------------- */
2886 if( poODS != NULL )
2887 {
2888 const char *pszDep
2889 = poODS->GetMetadataItem( "HFA_DEPENDENT_FILE", "HFA" );
2890 if( pszDep == NULL )
2891 {
2892 CPLDebug( "AUX",
2893 "Found %s but it has no dependent file, ignoring.",
2894 osAuxFilename.c_str() );
2895 GDALClose( poODS );
2896 poODS = NULL;
2897 }
2898 else if( !EQUAL(pszDep,osJustFile) )
2899 {
2900 VSIStatBufL sStatBuf;
2901
2902 if( VSIStatExL( pszDep, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
2903 {
2904 CPLDebug( "AUX", "%s is for file %s, not %s, ignoring.",
2905 osAuxFilename.c_str(),
2906 pszDep, osJustFile.c_str() );
2907 GDALClose( poODS );
2908 poODS = NULL;
2909 }
2910 else
2911 {
2912 CPLDebug( "AUX", "%s is for file %s, not %s, but since\n"
2913 "%s does not exist, we will use .aux file as our own.",
2914 osAuxFilename.c_str(),
2915 pszDep, osJustFile.c_str(),
2916 pszDep );
2917 }
2918 }
2919
2920 /* -------------------------------------------------------------------- */
2921 /* Confirm that the aux file matches the configuration of the */
2922 /* dependent dataset. */
2923 /* -------------------------------------------------------------------- */
2924 if( poODS != NULL && poDependentDS != NULL
2925 && (poODS->GetRasterCount() != poDependentDS->GetRasterCount()
2926 || poODS->GetRasterXSize() != poDependentDS->GetRasterXSize()
2927 || poODS->GetRasterYSize() != poDependentDS->GetRasterYSize()) )
2928 {
2929 CPLDebug( "AUX",
2930 "Ignoring aux file %s as its raster configuration\n"
2931 "(%dP x %dL x %dB) does not match master file (%dP x %dL x %dB)",
2932 osAuxFilename.c_str(),
2933 poODS->GetRasterXSize(),
2934 poODS->GetRasterYSize(),
2935 poODS->GetRasterCount(),
2936 poDependentDS->GetRasterXSize(),
2937 poDependentDS->GetRasterYSize(),
2938 poDependentDS->GetRasterCount() );
2939
2940 GDALClose( poODS );
2941 poODS = NULL;
2942 }
2943 }
2944
2945 /* -------------------------------------------------------------------- */
2946 /* Try appending .aux to the end of the filename. */
2947 /* -------------------------------------------------------------------- */
2948 if( poODS == NULL )
2949 {
2950 osAuxFilename = pszBasename;
2951 osAuxFilename += ".";
2952 osAuxFilename += pszAuxSuffixLC;
2953 fp = VSIFOpenL( osAuxFilename, "rb" );
2954 if ( fp == NULL && VSIIsCaseSensitiveFS(osAuxFilename) )
2955 {
2956 // Can't found file with lower case suffix. Try the upper case one.
2957 osAuxFilename = pszBasename;
2958 osAuxFilename += ".";
2959 osAuxFilename += pszAuxSuffixUC;
2960 fp = VSIFOpenL( osAuxFilename, "rb" );
2961 }
2962
2963 if( fp != NULL )
2964 {
2965 if( VSIFReadL( abyHeader, 1, 32, fp ) == 32 &&
2966 EQUALN((char *) abyHeader,"EHFA_HEADER_TAG",15) )
2967 {
2968 /* Avoid causing failure in opening of main file from SWIG bindings */
2969 /* when auxiliary file cannot be opened (#3269) */
2970 CPLTurnFailureIntoWarning(TRUE);
2971 if( poDependentDS != NULL && poDependentDS->GetShared() )
2972 poODS = (GDALDataset *) GDALOpenShared( osAuxFilename, eAccess );
2973 else
2974 poODS = (GDALDataset *) GDALOpen( osAuxFilename, eAccess );
2975 CPLTurnFailureIntoWarning(FALSE);
2976 }
2977 VSIFCloseL( fp );
2978 }
2979
2980 if( poODS != NULL )
2981 {
2982 const char *pszDep
2983 = poODS->GetMetadataItem( "HFA_DEPENDENT_FILE", "HFA" );
2984 if( pszDep == NULL )
2985 {
2986 CPLDebug( "AUX",
2987 "Found %s but it has no dependent file, ignoring.",
2988 osAuxFilename.c_str() );
2989 GDALClose( poODS );
2990 poODS = NULL;
2991 }
2992 else if( !EQUAL(pszDep,osJustFile) )
2993 {
2994 VSIStatBufL sStatBuf;
2995
2996 if( VSIStatExL( pszDep, &sStatBuf, VSI_STAT_EXISTS_FLAG ) == 0 )
2997 {
2998 CPLDebug( "AUX", "%s is for file %s, not %s, ignoring.",
2999 osAuxFilename.c_str(),
3000 pszDep, osJustFile.c_str() );
3001 GDALClose( poODS );
3002 poODS = NULL;
3003 }
3004 else
3005 {
3006 CPLDebug( "AUX", "%s is for file %s, not %s, but since\n"
3007 "%s does not exist, we will use .aux file as our own.",
3008 osAuxFilename.c_str(),
3009 pszDep, osJustFile.c_str(),
3010 pszDep );
3011 }
3012 }
3013 }
3014 }
3015
3016 /* -------------------------------------------------------------------- */
3017 /* Confirm that the aux file matches the configuration of the */
3018 /* dependent dataset. */
3019 /* -------------------------------------------------------------------- */
3020 if( poODS != NULL && poDependentDS != NULL
3021 && (poODS->GetRasterCount() != poDependentDS->GetRasterCount()
3022 || poODS->GetRasterXSize() != poDependentDS->GetRasterXSize()
3023 || poODS->GetRasterYSize() != poDependentDS->GetRasterYSize()) )
3024 {
3025 CPLDebug( "AUX",
3026 "Ignoring aux file %s as its raster configuration\n"
3027 "(%dP x %dL x %dB) does not match master file (%dP x %dL x %dB)",
3028 osAuxFilename.c_str(),
3029 poODS->GetRasterXSize(),
3030 poODS->GetRasterYSize(),
3031 poODS->GetRasterCount(),
3032 poDependentDS->GetRasterXSize(),
3033 poDependentDS->GetRasterYSize(),
3034 poDependentDS->GetRasterCount() );
3035
3036 GDALClose( poODS );
3037 poODS = NULL;
3038 }
3039
3040 return poODS;
3041 }
3042
3043 /************************************************************************/
3044 /* -------------------------------------------------------------------- */
3045 /* The following stubs are present to ensure that older GDAL */
3046 /* bridges don't fail with newer libraries. */
3047 /* -------------------------------------------------------------------- */
3048 /************************************************************************/
3049
3050 CPL_C_START
3051
GDALCreateProjDef(const char *)3052 void * CPL_STDCALL GDALCreateProjDef( const char * )
3053 {
3054 CPLDebug( "GDAL", "GDALCreateProjDef no longer supported." );
3055 return NULL;
3056 }
3057
GDALReprojectToLongLat(void *,double *,double *)3058 CPLErr CPL_STDCALL GDALReprojectToLongLat( void *, double *, double * )
3059 {
3060 CPLDebug( "GDAL", "GDALReprojectToLatLong no longer supported." );
3061 return CE_Failure;
3062 }
3063
GDALReprojectFromLongLat(void *,double *,double *)3064 CPLErr CPL_STDCALL GDALReprojectFromLongLat( void *, double *, double * )
3065 {
3066 CPLDebug( "GDAL", "GDALReprojectFromLatLong no longer supported." );
3067 return CE_Failure;
3068 }
3069
GDALDestroyProjDef(void *)3070 void CPL_STDCALL GDALDestroyProjDef( void * )
3071
3072 {
3073 CPLDebug( "GDAL", "GDALDestroyProjDef no longer supported." );
3074 }
3075
3076 CPL_C_END
3077
3078 /************************************************************************/
3079 /* Infrastructure to check that dataset characteristics are valid */
3080 /************************************************************************/
3081
3082 CPL_C_START
3083
3084 /**
3085 * \brief Return TRUE if the dataset dimensions are valid.
3086 *
3087 * @param nXSize raster width
3088 * @param nYSize raster height
3089 *
3090 * @since GDAL 1.7.0
3091 */
GDALCheckDatasetDimensions(int nXSize,int nYSize)3092 int GDALCheckDatasetDimensions( int nXSize, int nYSize )
3093 {
3094 if (nXSize <= 0 || nYSize <= 0)
3095 {
3096 CPLError(CE_Failure, CPLE_AppDefined,
3097 "Invalid dataset dimensions : %d x %d", nXSize, nYSize);
3098 return FALSE;
3099 }
3100 return TRUE;
3101 }
3102
3103 /**
3104 * \brief Return TRUE if the band count is valid.
3105 *
3106 * If the configuration option GDAL_MAX_BAND_COUNT is defined,
3107 * the band count will be compared to the maximum number of band allowed.
3108 *
3109 * @param nBands the band count
3110 * @param bIsZeroAllowed TRUE if band count == 0 is allowed
3111 *
3112 * @since GDAL 1.7.0
3113 */
3114
GDALCheckBandCount(int nBands,int bIsZeroAllowed)3115 int GDALCheckBandCount( int nBands, int bIsZeroAllowed )
3116 {
3117 int nMaxBands = -1;
3118 const char* pszMaxBandCount = CPLGetConfigOption("GDAL_MAX_BAND_COUNT", NULL);
3119 if (pszMaxBandCount != NULL)
3120 {
3121 nMaxBands = atoi(pszMaxBandCount);
3122 }
3123 if (nBands < 0 || (!bIsZeroAllowed && nBands == 0) ||
3124 (nMaxBands >= 0 && nBands > nMaxBands) )
3125 {
3126 CPLError(CE_Failure, CPLE_AppDefined,
3127 "Invalid band count : %d", nBands);
3128 return FALSE;
3129 }
3130 return TRUE;
3131 }
3132
3133 CPL_C_END
3134
3135
3136 /************************************************************************/
3137 /* GDALSerializeGCPListToXML() */
3138 /************************************************************************/
3139
GDALSerializeGCPListToXML(CPLXMLNode * psParentNode,GDAL_GCP * pasGCPList,int nGCPCount,const char * pszGCPProjection)3140 void GDALSerializeGCPListToXML( CPLXMLNode* psParentNode,
3141 GDAL_GCP* pasGCPList,
3142 int nGCPCount,
3143 const char* pszGCPProjection )
3144 {
3145 CPLString oFmt;
3146
3147 CPLXMLNode *psPamGCPList = CPLCreateXMLNode( psParentNode, CXT_Element,
3148 "GCPList" );
3149
3150 CPLXMLNode* psLastChild = NULL;
3151
3152 if( pszGCPProjection != NULL
3153 && strlen(pszGCPProjection) > 0 )
3154 {
3155 CPLSetXMLValue( psPamGCPList, "#Projection",
3156 pszGCPProjection );
3157 psLastChild = psPamGCPList->psChild;
3158 }
3159
3160 for( int iGCP = 0; iGCP < nGCPCount; iGCP++ )
3161 {
3162 CPLXMLNode *psXMLGCP;
3163 GDAL_GCP *psGCP = pasGCPList + iGCP;
3164
3165 psXMLGCP = CPLCreateXMLNode( NULL, CXT_Element, "GCP" );
3166
3167 if( psLastChild == NULL )
3168 psPamGCPList->psChild = psXMLGCP;
3169 else
3170 psLastChild->psNext = psXMLGCP;
3171 psLastChild = psXMLGCP;
3172
3173 CPLSetXMLValue( psXMLGCP, "#Id", psGCP->pszId );
3174
3175 if( psGCP->pszInfo != NULL && strlen(psGCP->pszInfo) > 0 )
3176 CPLSetXMLValue( psXMLGCP, "Info", psGCP->pszInfo );
3177
3178 CPLSetXMLValue( psXMLGCP, "#Pixel",
3179 oFmt.Printf( "%.4f", psGCP->dfGCPPixel ) );
3180
3181 CPLSetXMLValue( psXMLGCP, "#Line",
3182 oFmt.Printf( "%.4f", psGCP->dfGCPLine ) );
3183
3184 CPLSetXMLValue( psXMLGCP, "#X",
3185 oFmt.Printf( "%.12E", psGCP->dfGCPX ) );
3186
3187 CPLSetXMLValue( psXMLGCP, "#Y",
3188 oFmt.Printf( "%.12E", psGCP->dfGCPY ) );
3189
3190 /* Note: GDAL 1.10.1 and older generated #GCPZ, but could not read it back */
3191 if( psGCP->dfGCPZ != 0.0 )
3192 CPLSetXMLValue( psXMLGCP, "#Z",
3193 oFmt.Printf( "%.12E", psGCP->dfGCPZ ) );
3194 }
3195 }
3196
3197 /************************************************************************/
3198 /* GDALDeserializeGCPListFromXML() */
3199 /************************************************************************/
3200
GDALDeserializeGCPListFromXML(CPLXMLNode * psGCPList,GDAL_GCP ** ppasGCPList,int * pnGCPCount,char ** ppszGCPProjection)3201 void GDALDeserializeGCPListFromXML( CPLXMLNode* psGCPList,
3202 GDAL_GCP** ppasGCPList,
3203 int* pnGCPCount,
3204 char** ppszGCPProjection )
3205 {
3206 CPLXMLNode *psXMLGCP;
3207 OGRSpatialReference oSRS;
3208
3209 if( ppszGCPProjection )
3210 {
3211 const char *pszRawProj = CPLGetXMLValue(psGCPList, "Projection", "");
3212
3213 if( strlen(pszRawProj) > 0
3214 && oSRS.SetFromUserInput( pszRawProj ) == OGRERR_NONE )
3215 oSRS.exportToWkt( ppszGCPProjection );
3216 else
3217 *ppszGCPProjection = CPLStrdup("");
3218 }
3219
3220 // Count GCPs.
3221 int nGCPMax = 0;
3222
3223 for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
3224 psXMLGCP = psXMLGCP->psNext )
3225 nGCPMax++;
3226
3227 *ppasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),nGCPMax);
3228 *pnGCPCount = 0;
3229
3230 for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
3231 psXMLGCP = psXMLGCP->psNext )
3232 {
3233 GDAL_GCP *psGCP = *ppasGCPList + *pnGCPCount;
3234
3235 if( !EQUAL(psXMLGCP->pszValue,"GCP") ||
3236 psXMLGCP->eType != CXT_Element )
3237 continue;
3238
3239 GDALInitGCPs( 1, psGCP );
3240
3241 CPLFree( psGCP->pszId );
3242 psGCP->pszId = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Id",""));
3243
3244 CPLFree( psGCP->pszInfo );
3245 psGCP->pszInfo = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Info",""));
3246
3247 psGCP->dfGCPPixel = CPLAtof(CPLGetXMLValue(psXMLGCP,"Pixel","0.0"));
3248 psGCP->dfGCPLine = CPLAtof(CPLGetXMLValue(psXMLGCP,"Line","0.0"));
3249
3250 psGCP->dfGCPX = CPLAtof(CPLGetXMLValue(psXMLGCP,"X","0.0"));
3251 psGCP->dfGCPY = CPLAtof(CPLGetXMLValue(psXMLGCP,"Y","0.0"));
3252 const char* pszZ = CPLGetXMLValue(psXMLGCP,"Z",NULL);
3253 if( pszZ == NULL )
3254 {
3255 /* Note: GDAL 1.10.1 and older generated #GCPZ, but could not read it back */
3256 pszZ = CPLGetXMLValue(psXMLGCP,"GCPZ","0.0");
3257 }
3258 psGCP->dfGCPZ = CPLAtof(pszZ);
3259
3260 (*pnGCPCount) ++;
3261 }
3262 }
3263
3264 /************************************************************************/
3265 /* GDALSerializeOpenOptionsToXML() */
3266 /************************************************************************/
3267
GDALSerializeOpenOptionsToXML(CPLXMLNode * psParentNode,char ** papszOpenOptions)3268 void GDALSerializeOpenOptionsToXML( CPLXMLNode* psParentNode, char** papszOpenOptions)
3269 {
3270 if( papszOpenOptions != NULL )
3271 {
3272 CPLXMLNode* psOpenOptions = CPLCreateXMLNode( psParentNode, CXT_Element, "OpenOptions" );
3273 CPLXMLNode* psLastChild = NULL;
3274
3275 for(char** papszIter = papszOpenOptions; *papszIter != NULL; papszIter ++ )
3276 {
3277 const char *pszRawValue;
3278 char *pszKey = NULL;
3279 CPLXMLNode *psOOI;
3280
3281 pszRawValue = CPLParseNameValue( *papszIter, &pszKey );
3282
3283 psOOI = CPLCreateXMLNode( NULL, CXT_Element, "OOI" );
3284 if( psLastChild == NULL )
3285 psOpenOptions->psChild = psOOI;
3286 else
3287 psLastChild->psNext = psOOI;
3288 psLastChild = psOOI;
3289
3290 CPLSetXMLValue( psOOI, "#key", pszKey );
3291 CPLCreateXMLNode( psOOI, CXT_Text, pszRawValue );
3292
3293 CPLFree( pszKey );
3294 }
3295 }
3296 }
3297
3298 /************************************************************************/
3299 /* GDALDeserializeOpenOptionsFromXML() */
3300 /************************************************************************/
3301
GDALDeserializeOpenOptionsFromXML(CPLXMLNode * psParentNode)3302 char** GDALDeserializeOpenOptionsFromXML( CPLXMLNode* psParentNode )
3303 {
3304 char** papszOpenOptions = NULL;
3305 CPLXMLNode* psOpenOptions = CPLGetXMLNode(psParentNode, "OpenOptions");
3306 if( psOpenOptions != NULL )
3307 {
3308 CPLXMLNode* psOOI;
3309 for( psOOI = psOpenOptions->psChild; psOOI != NULL;
3310 psOOI = psOOI->psNext )
3311 {
3312 if( !EQUAL(psOOI->pszValue,"OOI")
3313 || psOOI->eType != CXT_Element
3314 || psOOI->psChild == NULL
3315 || psOOI->psChild->psNext == NULL
3316 || psOOI->psChild->eType != CXT_Attribute
3317 || psOOI->psChild->psChild == NULL )
3318 continue;
3319
3320 char* pszName = psOOI->psChild->psChild->pszValue;
3321 char* pszValue = psOOI->psChild->psNext->pszValue;
3322 if( pszName != NULL && pszValue != NULL )
3323 papszOpenOptions = CSLSetNameValue( papszOpenOptions, pszName, pszValue );
3324 }
3325 }
3326 return papszOpenOptions;
3327 }
3328
3329 /************************************************************************/
3330 /* GDALRasterIOGetResampleAlg() */
3331 /************************************************************************/
3332
GDALRasterIOGetResampleAlg(const char * pszResampling)3333 GDALRIOResampleAlg GDALRasterIOGetResampleAlg(const char* pszResampling)
3334 {
3335 GDALRIOResampleAlg eResampleAlg = GRIORA_NearestNeighbour;
3336 if( EQUALN(pszResampling, "NEAR", 4) )
3337 eResampleAlg = GRIORA_NearestNeighbour;
3338 else if( EQUAL(pszResampling, "BILINEAR") )
3339 eResampleAlg = GRIORA_Bilinear;
3340 else if( EQUAL(pszResampling, "CUBIC") )
3341 eResampleAlg = GRIORA_Cubic;
3342 else if( EQUAL(pszResampling, "CUBICSPLINE") )
3343 eResampleAlg = GRIORA_CubicSpline;
3344 else if( EQUAL(pszResampling, "LANCZOS") )
3345 eResampleAlg = GRIORA_Lanczos;
3346 else if( EQUAL(pszResampling, "AVERAGE") )
3347 eResampleAlg = GRIORA_Average;
3348 else if( EQUAL(pszResampling, "MODE") )
3349 eResampleAlg = GRIORA_Mode;
3350 else if( EQUAL(pszResampling, "GAUSS") )
3351 eResampleAlg = GRIORA_Gauss;
3352 else
3353 CPLError(CE_Warning, CPLE_NotSupported,
3354 "GDAL_RASTERIO_RESAMPLING = %s not supported", pszResampling);
3355 return eResampleAlg;
3356 }
3357
3358 /************************************************************************/
3359 /* GDALRasterIOExtraArgSetResampleAlg() */
3360 /************************************************************************/
3361
GDALRasterIOExtraArgSetResampleAlg(GDALRasterIOExtraArg * psExtraArg,int nXSize,int nYSize,int nBufXSize,int nBufYSize)3362 void GDALRasterIOExtraArgSetResampleAlg(GDALRasterIOExtraArg* psExtraArg,
3363 int nXSize, int nYSize,
3364 int nBufXSize, int nBufYSize)
3365 {
3366 if( (nBufXSize != nXSize || nBufYSize != nYSize) &&
3367 psExtraArg->eResampleAlg == GRIORA_NearestNeighbour )
3368 {
3369 const char* pszResampling = CPLGetConfigOption("GDAL_RASTERIO_RESAMPLING", NULL);
3370 if( pszResampling != NULL )
3371 {
3372 psExtraArg->eResampleAlg = GDALRasterIOGetResampleAlg(pszResampling);
3373 }
3374 }
3375 }
3376
3377 /************************************************************************/
3378 /* GDALCanFileAcceptSidecarFile() */
3379 /************************************************************************/
3380
GDALCanFileAcceptSidecarFile(const char * pszFilename)3381 int GDALCanFileAcceptSidecarFile(const char* pszFilename)
3382 {
3383 if( strstr(pszFilename, "/vsicurl/") && strchr(pszFilename, '?') )
3384 return FALSE;
3385 return TRUE;
3386 }
3387