Main Page Class Hierarchy Compound List File List Compound Members File Members Related Pages

# gdal_misc.cpp

    00001 /******************************************************************************
    00002  * $Id: gdal_misc_cpp-source.html,v 1.13 2002/12/21 19:13:12 warmerda Exp $
    00003  *
    00004  * Project:  GDAL Core
    00005  * Purpose:  Free standing functions for GDAL.
    00006  * Author:   Frank Warmerdam, warmerda@home.com
    00007  *
    00008  ******************************************************************************
    00009  * Copyright (c) 1999, Frank Warmerdam
    00010  *
    00011  * Permission is hereby granted, free of charge, to any person obtaining a
    00012  * copy of this software and associated documentation files (the "Software"),
    00013  * to deal in the Software without restriction, including without limitation
    00014  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
    00015  * and/or sell copies of the Software, and to permit persons to whom the
    00016  * Software is furnished to do so, subject to the following conditions:
    00017  *
    00018  * The above copyright notice and this permission notice shall be included
    00019  * in all copies or substantial portions of the Software.
    00020  *
    00021  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
    00022  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    00023  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
    00024  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    00025  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
    00026  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
    00027  * DEALINGS IN THE SOFTWARE.
    00028  ******************************************************************************
    00029  *
    00030  * $Log: gdal_misc_cpp-source.html,v $
    00030  * Revision 1.13  2002/12/21 19:13:12  warmerda
    00030  * updated
    00030  *
    00031  * Revision 1.40  2002/12/11 21:21:46  warmerda
    00032  * fixed debug format problem
    00033  *
    00034  * Revision 1.39  2002/12/09 20:05:31  warmerda
    00035  * fixed return flag from GDALReadTabFile
    00036  *
    00037  * Revision 1.38  2002/12/09 18:53:25  warmerda
    00038  * GDALDecToDMS() now calls CPLDecToDMS()
    00039  *
    00040  * Revision 1.37  2002/12/05 17:55:30  warmerda
    00041  * gdalreadtabfile should not be static
    00042  *
    00043  * Revision 1.36  2002/12/05 15:46:38  warmerda
    00044  * added GDALReadTabFile()
    00045  *
    00046  * Revision 1.35  2002/07/09 20:33:12  warmerda
    00047  * expand tabs
    00048  *
    00049  * Revision 1.34  2002/05/06 21:37:29  warmerda
    00050  * added GDALGCPsToGeoTransform
    00051  *
    00052  * Revision 1.33  2002/04/25 16:18:41  warmerda
    00053  * added extra checking
    00054  *
    00055  * Revision 1.32  2002/04/24 19:21:26  warmerda
    00056  * Include <ctype.h> for toupper(), tolower().
    00057  *
    00058  * Revision 1.31  2002/04/19 12:22:05  dron
    00059  * added GDALWriteWorldFile()
    00060  *
    00061  * Revision 1.30  2002/04/16 13:59:33  warmerda
    00062  * added GDALVersionInfo
    00063  *
    00064  * Revision 1.29  2001/12/07 20:04:21  warmerda
    00065  * fixed serious bug in random sampler
    00066  *
    00067  * Revision 1.28  2001/11/30 03:41:26  warmerda
    00068  * Fixed bug with the block sampling rate being too low to satisfy large
    00069  * sample count values.  Fixed bug with tiled images including some uninitialized
    00070  * or zero data in the sample set on partial edge tiles.
    00071  *
    00072  * Revision 1.27  2001/11/26 20:14:01  warmerda
    00073  * added GDALProjDef stubs for old 'bridges'
    00074  *
    00075  * Revision 1.26  2001/11/19 16:03:16  warmerda
    00076  * moved GDALDectoDMS here
    00077  *
    00078  * Revision 1.25  2001/08/15 15:05:44  warmerda
    00079  * return magnitude for complex samples in random sampler
    00080  *
    00081  * Revision 1.24  2001/07/18 04:04:30  warmerda
    00082  * added CPL_CVSID
    00083  *
    00084  * Revision 1.23  2001/05/01 18:09:25  warmerda
    00085  * added GDALReadWorldFile()
    00086  *
    00087  * Revision 1.22  2000/12/04 20:45:14  warmerda
    00088  * removed unused variable.
    00089  *
    00090  * Revision 1.21  2000/10/06 15:22:49  warmerda
    00091  * added GDALDataTypeUnion
    00092  *
    00093  * Revision 1.20  2000/08/18 15:24:48  warmerda
    00094  * added GDALTermProgress
    00095  *
    00096  * Revision 1.19  2000/08/09 16:25:42  warmerda
    00097  * don't crash if block is null
    00098  *
    00099  * Revision 1.18  2000/07/11 14:35:43  warmerda
    00100  * added documentation
    00101  *
    00102  * Revision 1.17  2000/07/05 17:53:33  warmerda
    00103  * Removed unused code related to nXCheck.
    00104  *
    00105  * Revision 1.16  2000/06/27 17:21:26  warmerda
    00106  * added GDALGetRasterSampleOverview
    00107  *
    00108  * Revision 1.15  2000/06/26 22:17:49  warmerda
    00109  * added scaled progress support
    00110  *
    00111  * Revision 1.14  2000/06/05 17:24:05  warmerda
    00112  * added real complex support
    00113  *
    00114  * Revision 1.13  2000/04/21 21:55:32  warmerda
    00115  * made more robust if block read fails
    00116  *
    00117  * Revision 1.12  2000/04/17 20:59:40  warmerda
    00118  * Removed printf.
    00119  *
    00120  * Revision 1.11  2000/04/17 20:59:14  warmerda
    00121  * fixed sampling bug
    00122  *
    00123  * Revision 1.10  2000/03/31 13:41:45  warmerda
    00124  * added gcp support functions
    00125  *
    00126  * Revision 1.9  2000/03/24 00:09:19  warmerda
    00127  * added sort-of random sampling
    00128  *
    00129  * Revision 1.8  2000/03/23 16:53:21  warmerda
    00130  * use overviews for approximate min/max
    00131  *
    00132  * Revision 1.7  2000/03/09 23:21:44  warmerda
    00133  * added GDALDummyProgress
    00134  *
    00135  * Revision 1.6  2000/03/06 21:59:44  warmerda
    00136  * added min/max calculate
    00137  *
    00138  * Revision 1.5  2000/03/06 02:20:15  warmerda
    00139  * added getname functions for colour interpretations
    00140  *
    00141  * Revision 1.4  1999/07/23 19:35:47  warmerda
    00142  * added GDALGetDataTypeName
    00143  *
    00144  * Revision 1.3  1999/05/17 02:00:45  vgough
    00145  * made pure_virtual C linkage
    00146  *
    00147  * Revision 1.2  1999/05/16 19:32:13  warmerda
    00148  * Added __pure_virtual.
    00149  *
    00150  * Revision 1.1  1998/12/06 02:50:16  warmerda
    00151  * New
    00152  *
    00153  */
    00154
    00155 #include "gdal_priv.h"
    00156 #include "cpl_string.h"
    00157 #include <ctype.h>
    00158
    00159 CPL_CVSID("$Id: gdal_misc_cpp-source.html,v 1.13 2002/12/21 19:13:12 warmerda Exp $");
    00160
    00161 #ifdef HAVE_MITAB
    00162 #include "ogr_spatialref.h"
    00163 // from mitab component.
    00164 OGRSpatialReference * MITABCoordSys2SpatialRef( const char * pszCoordSys );
    00165 #endif
    00166
    00167 /************************************************************************/
    00168 /*                           __pure_virtual()                           */
    00169 /*                                                                      */
    00170 /*      The following is a gross hack to remove the last remaining      */
    00171 /*      dependency on the GNU C++ standard library.                     */
    00172 /************************************************************************/
    00173
    00174 #ifdef __GNUC__
    00175
    00176 extern "C"
    00177 void __pure_virtual()
    00178
    00179 {
    00180 }
    00181
    00182 #endif
    00183
    00184 /************************************************************************/
    00185 /*                         GDALDataTypeUnion()                          */
    00186 /************************************************************************/
    00187
    00198 GDALDataType GDALDataTypeUnion( GDALDataType eType1, GDALDataType eType2 )
    00199
    00200 {
    00201     int         bFloating, bComplex, nBits, bSigned;
    00202
    00203     bComplex = GDALDataTypeIsComplex(eType1) | GDALDataTypeIsComplex(eType2);
    00204
    00205     switch( eType1 )
    00206     {
    00207       case GDT_Byte:
    00208         nBits = 8;
    00209         bSigned = FALSE;
    00210         bFloating = FALSE;
    00211         break;
    00212
    00213       case GDT_Int16:
    00214       case GDT_CInt16:
    00215         nBits = 16;
    00216         bSigned = TRUE;
    00217         bFloating = FALSE;
    00218         break;
    00219
    00220       case GDT_UInt16:
    00221         nBits = 16;
    00222         bSigned = FALSE;
    00223         bFloating = FALSE;
    00224         break;
    00225
    00226       case GDT_Int32:
    00227       case GDT_CInt32:
    00228         nBits = 32;
    00229         bSigned = TRUE;
    00230         bFloating = FALSE;
    00231         break;
    00232
    00233       case GDT_UInt32:
    00234         nBits = 32;
    00235         bSigned = FALSE;
    00236         bFloating = FALSE;
    00237         break;
    00238
    00239       case GDT_Float32:
    00240       case GDT_CFloat32:
    00241         nBits = 32;
    00242         bSigned = TRUE;
    00243         bFloating = TRUE;
    00244         break;
    00245
    00246       case GDT_Float64:
    00247       case GDT_CFloat64:
    00248         nBits = 64;
    00249         bSigned = TRUE;
    00250         bFloating = TRUE;
    00251         break;
    00252
    00253       default:
    00254         CPLAssert( FALSE );
    00255         return GDT_Unknown;
    00256     }
    00257
    00258     switch( eType2 )
    00259     {
    00260       case GDT_Byte:
    00261         break;
    00262
    00263       case GDT_Int16:
    00264         nBits = MAX(nBits,16);
    00265         bSigned = TRUE;
    00266         break;
    00267
    00268       case GDT_UInt16:
    00269         nBits = MAX(nBits,16);
    00270         break;
    00271
    00272       case GDT_Int32:
    00273       case GDT_CInt32:
    00274         nBits = MAX(nBits,32);
    00275         bSigned = TRUE;
    00276         break;
    00277
    00278       case GDT_UInt32:
    00279         nBits = MAX(nBits,32);
    00280         break;
    00281
    00282       case GDT_Float32:
    00283       case GDT_CFloat32:
    00284         nBits = MAX(nBits,32);
    00285         bSigned = TRUE;
    00286         bFloating = TRUE;
    00287         break;
    00288
    00289       case GDT_Float64:
    00290       case GDT_CFloat64:
    00291         nBits = MAX(nBits,64);
    00292         bSigned = TRUE;
    00293         bFloating = TRUE;
    00294         break;
    00295
    00296       default:
    00297         CPLAssert( FALSE );
    00298         return GDT_Unknown;
    00299     }
    00300
    00301     if( nBits == 8 )
    00302         return GDT_Byte;
    00303     else if( nBits == 16 && bComplex )
    00304         return GDT_CInt16;
    00305     else if( nBits == 16 && bSigned )
    00306         return GDT_Int16;
    00307     else if( nBits == 16 && !bSigned )
    00308         return GDT_UInt16;
    00309     else if( nBits == 32 && bFloating && bComplex )
    00310         return GDT_CFloat32;
    00311     else if( nBits == 32 && bFloating )
    00312         return GDT_Float32;
    00313     else if( nBits == 32 && bComplex )
    00314         return GDT_CInt32;
    00315     else if( nBits == 32 && bSigned )
    00316         return GDT_Int32;
    00317     else if( nBits == 32 && !bSigned )
    00318         return GDT_UInt32;
    00319     else if( nBits == 64 && bComplex )
    00320         return GDT_CFloat64;
    00321     else
    00322         return GDT_Float64;
    00323 }
    00324
    00325
    00326 /************************************************************************/
    00327 /*                        GDALGetDataTypeSize()                         */
    00328 /************************************************************************/
    00329 int GDALGetDataTypeSize( GDALDataType eDataType )
    00330
    00331 {
    00332     switch( eDataType )
    00333     {
    00334       case GDT_Byte:
    00335         return 8;
    00336
    00337       case GDT_UInt16:
    00338       case GDT_Int16:
    00339         return 16;
    00340
    00341       case GDT_UInt32:
    00342       case GDT_Int32:
    00343       case GDT_Float32:
    00344       case GDT_CInt16:
    00345         return 32;
    00346
    00347       case GDT_Float64:
    00348       case GDT_CInt32:
    00349       case GDT_CFloat32:
    00350         return 64;
    00351
    00352       case GDT_CFloat64:
    00353         return 128;
    00354
    00355       default:
    00356         CPLAssert( FALSE );
    00357         return 0;
    00358     }
    00359 }
    00360
    00361 /************************************************************************/
    00362 /*                       GDALDataTypeIsComplex()                        */
    00363 /************************************************************************/
    00364
    00365 int GDALDataTypeIsComplex( GDALDataType eDataType )
    00366
    00367 {
    00368     switch( eDataType )
    00369     {
    00370       case GDT_CInt16:
    00371       case GDT_CInt32:
    00372       case GDT_CFloat32:
    00373       case GDT_CFloat64:
    00374         return TRUE;
    00375
    00376       default:
    00377         return FALSE;
    00378     }
    00379 }
    00380
    00381 /************************************************************************/
    00382 /*                        GDALGetDataTypeName()                         */
    00383 /************************************************************************/
    00384
    00385 const char *GDALGetDataTypeName( GDALDataType eDataType )
    00386
    00387 {
    00388     switch( eDataType )
    00389     {
    00390       case GDT_Byte:
    00391         return "Byte";
    00392
    00393       case GDT_UInt16:
    00394         return "UInt16";
    00395
    00396       case GDT_Int16:
    00397         return "Int16";
    00398
    00399       case GDT_UInt32:
    00400         return "UInt32";
    00401
    00402       case GDT_Int32:
    00403         return "Int32";
    00404
    00405       case GDT_Float32:
    00406         return "Float32";
    00407
    00408       case GDT_Float64:
    00409         return "Float64";
    00410
    00411       case GDT_CInt16:
    00412         return "CInt16";
    00413
    00414       case GDT_CInt32:
    00415         return "CInt32";
    00416
    00417       case GDT_CFloat32:
    00418         return "CFloat32";
    00419
    00420       case GDT_CFloat64:
    00421         return "CFloat64";
    00422
    00423       default:
    00424         return NULL;
    00425     }
    00426 }
    00427
    00428 /************************************************************************/
    00429 /*                  GDALGetPaletteInterpretationName()                  */
    00430 /************************************************************************/
    00431
    00432 const char *GDALGetPaletteInterpretationName( GDALPaletteInterp eInterp )
    00433
    00434 {
    00435     switch( eInterp )
    00436     {
    00437       case GPI_Gray:
    00438         return "Gray";
    00439
    00440       case GPI_RGB:
    00441         return "RGB";
    00442
    00443       case GPI_CMYK:
    00444         return "CMYK";
    00445
    00446       case GPI_HLS:
    00447         return "HLS";
    00448
    00449       default:
    00450         return "Unknown";
    00451     }
    00452 }
    00453
    00454 /************************************************************************/
    00455 /*                   GDALGetColorInterpretationName()                   */
    00456 /************************************************************************/
    00457
    00458 const char *GDALGetColorInterpretationName( GDALColorInterp eInterp )
    00459
    00460 {
    00461     switch( eInterp )
    00462     {
    00463       case GCI_Undefined:
    00464         return "Undefined";
    00465
    00466       case GCI_GrayIndex:
    00467         return "Gray";
    00468
    00469       case GCI_PaletteIndex:
    00470         return "Palette";
    00471
    00472       case GCI_RedBand:
    00473         return "Red";
    00474
    00475       case GCI_GreenBand:
    00476         return "Green";
    00477
    00478       case GCI_BlueBand:
    00479         return "Blue";
    00480
    00481       case GCI_AlphaBand:
    00482         return "Alpha";
    00483
    00484       case GCI_HueBand:
    00485         return "Hue";
    00486
    00487       case GCI_SaturationBand:
    00488         return "Saturation";
    00489
    00490       case GCI_LightnessBand:
    00491         return "Lightness";
    00492
    00493       case GCI_CyanBand:
    00494         return "Cyan";
    00495
    00496       case GCI_MagentaBand:
    00497         return "Magenta";
    00498
    00499       case GCI_YellowBand:
    00500         return "Yellow";
    00501
    00502       case GCI_BlackBand:
    00503         return "Black";
    00504
    00505       default:
    00506         return "Unknown";
    00507     }
    00508 }
    00509
    00510 /************************************************************************/
    00511 /*                      GDALComputeRasterMinMax()                       */
    00512 /************************************************************************/
    00513
    00532 void GDALComputeRasterMinMax( GDALRasterBandH hBand, int bApproxOK,
    00533                               double adfMinMax[2] )
    00534
    00535 {
    00536     double       dfMin=0.0, dfMax=0.0;
    00537     GDALRasterBand *poBand;
    00538
    00539 /* -------------------------------------------------------------------- */
    00540 /*      Does the driver already know the min/max?                       */
    00541 /* -------------------------------------------------------------------- */
    00542     if( bApproxOK )
    00543     {
    00544         int          bSuccessMin, bSuccessMax;
    00545
    00546         dfMin = GDALGetRasterMinimum( hBand, &bSuccessMin );
    00547         dfMax = GDALGetRasterMaximum( hBand, &bSuccessMax );
    00548
    00549         if( bSuccessMin && bSuccessMax )
    00550         {
    00551             adfMinMax[0] = dfMin;
    00552             adfMinMax[1] = dfMax;
    00553             return;
    00554         }
    00555     }
    00556
    00557 /* -------------------------------------------------------------------- */
    00558 /*      If we have overview bands, use them for min/max.                */
    00559 /* -------------------------------------------------------------------- */
    00560     if( bApproxOK )
    00561         poBand = (GDALRasterBand *) GDALGetRasterSampleOverview( hBand, 2500 );
    00562     else
    00563         poBand = (GDALRasterBand *) hBand;
    00564
    00565 /* -------------------------------------------------------------------- */
    00566 /*      Figure out the ratio of blocks we will read to get an           */
    00567 /*      approximate value.                                              */
    00568 /* -------------------------------------------------------------------- */
    00569     int         nBlockXSize, nBlockYSize;
    00570     int         nBlocksPerRow, nBlocksPerColumn;
    00571     int         nSampleRate;
    00572     int         bGotNoDataValue, bFirstValue = TRUE;
    00573     double      dfNoDataValue;
    00574
    00575     dfNoDataValue = poBand->GetNoDataValue( &bGotNoDataValue );
    00576
    00577     poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
    00578     nBlocksPerRow = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
    00579     nBlocksPerColumn = (poBand->GetYSize() + nBlockYSize - 1) / nBlockYSize;
    00580
    00581     if( bApproxOK )
    00582         nSampleRate =
    00583             (int) MAX(1,sqrt((double) nBlocksPerRow * nBlocksPerColumn));
    00584     else
    00585         nSampleRate = 1;
    00586
    00587     for( int iSampleBlock = 0;
    00588          iSampleBlock < nBlocksPerRow * nBlocksPerColumn;
    00589          iSampleBlock += nSampleRate )
    00590     {
    00591         double dfValue = 0.0;
    00592         int  iXBlock, iYBlock, nXCheck, nYCheck;
    00593         GDALRasterBlock *poBlock;
    00594
    00595         iYBlock = iSampleBlock / nBlocksPerRow;
    00596         iXBlock = iSampleBlock - nBlocksPerRow * iYBlock;
    00597
    00598         poBlock = poBand->GetBlockRef( iXBlock, iYBlock );
    00599         if( poBlock == NULL )
    00600             continue;
    00601
    00602         if( (iXBlock+1) * nBlockXSize > poBand->GetXSize() )
    00603             nXCheck = poBand->GetXSize() - iXBlock * nBlockXSize;
    00604         else
    00605             nXCheck = nBlockXSize;
    00606
    00607         if( (iYBlock+1) * nBlockYSize > poBand->GetYSize() )
    00608             nYCheck = poBand->GetYSize() - iYBlock * nBlockYSize;
    00609         else
    00610             nYCheck = nBlockYSize;
    00611
    00612         /* this isn't the fastest way to do this, but is easier for now */
    00613         for( int iY = 0; iY < nYCheck; iY++ )
    00614         {
    00615             for( int iX = 0; iX < nXCheck; iX++ )
    00616             {
    00617                 int    iOffset = iX + iY * nBlockXSize;
    00618
    00619                 switch( poBlock->GetDataType() )
    00620                 {
    00621                   case GDT_Byte:
    00622                     dfValue = ((GByte *) poBlock->GetDataRef())[iOffset];
    00623                     break;
    00624                   case GDT_UInt16:
    00625                     dfValue = ((GUInt16 *) poBlock->GetDataRef())[iOffset];
    00626                     break;
    00627                   case GDT_Int16:
    00628                     dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset];
    00629                     break;
    00630                   case GDT_UInt32:
    00631                     dfValue = ((GUInt32 *) poBlock->GetDataRef())[iOffset];
    00632                     break;
    00633                   case GDT_Int32:
    00634                     dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset];
    00635                     break;
    00636                   case GDT_Float32:
    00637                     dfValue = ((float *) poBlock->GetDataRef())[iOffset];
    00638                     break;
    00639                   case GDT_Float64:
    00640                     dfValue = ((double *) poBlock->GetDataRef())[iOffset];
    00641                     break;
    00642                   case GDT_CInt16:
    00643                     dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset*2];
    00644                     break;
    00645                   case GDT_CInt32:
    00646                     dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset*2];
    00647                     break;
    00648                   case GDT_CFloat32:
    00649                     dfValue = ((float *) poBlock->GetDataRef())[iOffset*2];
    00650                     break;
    00651                   case GDT_CFloat64:
    00652                     dfValue = ((double *) poBlock->GetDataRef())[iOffset*2];
    00653                     break;
    00654                   default:
    00655                     CPLAssert( FALSE );
    00656                 }
    00657
    00658                 if( bGotNoDataValue && dfValue == dfNoDataValue )
    00659                     continue;
    00660
    00661                 if( bFirstValue )
    00662                 {
    00663                     dfMin = dfMax = dfValue;
    00664                     bFirstValue = FALSE;
    00665                 }
    00666                 else
    00667                 {
    00668                     dfMin = MIN(dfMin,dfValue);
    00669                     dfMax = MAX(dfMax,dfValue);
    00670                 }
    00671             }
    00672         }
    00673     }
    00674
    00675     adfMinMax[0] = dfMin;
    00676     adfMinMax[1] = dfMax;
    00677 }
    00678
    00679 /************************************************************************/
    00680 /*                         GDALDummyProgress()                          */
    00681 /************************************************************************/
    00682
    00747 int GDALDummyProgress( double, const char *, void * )
    00748
    00749 {
    00750     return TRUE;
    00751 }
    00752
    00753 /************************************************************************/
    00754 /*                         GDALScaledProgress()                         */
    00755 /************************************************************************/
    00756 typedef struct {
    00757     GDALProgressFunc pfnProgress;
    00758     void *pData;
    00759     double dfMin;
    00760     double dfMax;
    00761 } GDALScaledProgressInfo;
    00762
    00763 int GDALScaledProgress( double dfComplete, const char *pszMessage,
    00764                         void *pData )
    00765
    00766 {
    00767     GDALScaledProgressInfo *psInfo = (GDALScaledProgressInfo *) pData;
    00768
    00769     return psInfo->pfnProgress( dfComplete * (psInfo->dfMax - psInfo->dfMin)
    00770                                 + psInfo->dfMin,
    00771                                 pszMessage, psInfo->pData );
    00772 }
    00773
    00774 /************************************************************************/
    00775 /*                      GDALCreateScaledProgress()                      */
    00776 /************************************************************************/
    00777
    00778 void *GDALCreateScaledProgress( double dfMin, double dfMax,
    00779                                 GDALProgressFunc pfnProgress,
    00780                                 void * pData )
    00781
    00782 {
    00783     GDALScaledProgressInfo *psInfo;
    00784
    00785     psInfo = (GDALScaledProgressInfo *)
    00786         CPLCalloc(sizeof(GDALScaledProgressInfo),1);
    00787
    00788     if( ABS(dfMin-dfMax) < 0.0000001 )
    00789         dfMax = dfMin + 0.01;
    00790
    00791     psInfo->pData = pData;
    00792     psInfo->pfnProgress = pfnProgress;
    00793     psInfo->dfMin = dfMin;
    00794     psInfo->dfMax = dfMax;
    00795
    00796     return (void *) psInfo;
    00797 }
    00798
    00799 /************************************************************************/
    00800 /*                     GDALDestroyScaledProgress()                      */
    00801 /************************************************************************/
    00802
    00803 void GDALDestroyScaledProgress( void * pData )
    00804
    00805 {
    00806     CPLFree( pData );
    00807 }
    00808
    00809 /************************************************************************/
    00810 /*                          GDALTermProgress()                          */
    00811 /************************************************************************/
    00812
    00813 int GDALTermProgress( double dfComplete, const char *pszMessage, void * )
    00814
    00815 {
    00816     static double dfLastComplete = -1.0;
    00817
    00818     if( dfLastComplete > dfComplete )
    00819     {
    00820         if( dfLastComplete > 1.0 )
    00821             dfLastComplete = -1.0;
    00822         else
    00823             dfLastComplete = dfComplete;
    00824     }
    00825
    00826     if( floor(dfLastComplete*10) != floor(dfComplete*10) )
    00827     {
    00828         int    nPercent = (int) floor(dfComplete*100);
    00829
    00830         if( nPercent == 0 && pszMessage != NULL )
    00831             fprintf( stdout, "%s:", pszMessage );
    00832
    00833         if( nPercent == 100 )
    00834             fprintf( stdout, "%d - done.\n", (int) floor(dfComplete*100) );
    00835         else
    00836         {
    00837             fprintf( stdout, "%d.", (int) floor(dfComplete*100) );
    00838             fflush( stdout );
    00839         }
    00840     }
    00841     else if( floor(dfLastComplete*30) != floor(dfComplete*30) )
    00842     {
    00843         fprintf( stdout, "." );
    00844         fflush( stdout );
    00845     }
    00846
    00847     dfLastComplete = dfComplete;
    00848
    00849     return TRUE;
    00850 }
    00851
    00852 /************************************************************************/
    00853 /*                    GDALGetRasterSampleOverview()                     */
    00854 /************************************************************************/
    00855
    00871 GDALRasterBandH GDALGetRasterSampleOverview( GDALRasterBandH hBand,
    00872                                              int nDesiredSamples )
    00873
    00874 {
    00875     int     nBestSamples;
    00876     GDALRasterBandH hBestBand = hBand;
    00877
    00878     nBestSamples = GDALGetRasterBandXSize(hBand)
    00879         * GDALGetRasterBandYSize(hBand);
    00880
    00881     for( int iOverview = 0;
    00882          iOverview < GDALGetOverviewCount( hBand );
    00883          iOverview++ )
    00884     {
    00885         GDALRasterBandH hOBand = GDALGetOverview( hBand, iOverview );
    00886         int    nOSamples;
    00887
    00888         nOSamples = GDALGetRasterBandXSize(hOBand)
    00889             * GDALGetRasterBandYSize(hOBand);
    00890
    00891         if( nOSamples < nBestSamples && nOSamples > nDesiredSamples )
    00892         {
    00893             nBestSamples = nOSamples;
    00894             hBestBand = hOBand;
    00895         }
    00896     }
    00897
    00898     return hBestBand;
    00899 }
    00900
    00901 /************************************************************************/
    00902 /*                     GDALGetRandomRasterSample()                      */
    00903 /************************************************************************/
    00904
    00905 int GDALGetRandomRasterSample( GDALRasterBandH hBand, int nSamples,
    00906                                float *pafSampleBuf )
    00907
    00908 {
    00909     GDALRasterBand *poBand;
    00910
    00911     poBand = (GDALRasterBand *) GDALGetRasterSampleOverview( hBand, nSamples );
    00912
    00913 /* -------------------------------------------------------------------- */
    00914 /*      Figure out the ratio of blocks we will read to get an           */
    00915 /*      approximate value.                                              */
    00916 /* -------------------------------------------------------------------- */
    00917     int         nBlockXSize, nBlockYSize;
    00918     int         nBlocksPerRow, nBlocksPerColumn;
    00919     int         nSampleRate;
    00920     int         bGotNoDataValue;
    00921     double      dfNoDataValue;
    00922     int         nActualSamples = 0;
    00923     int         nBlockSampleRate;
    00924     int         nBlockPixels, nBlockCount;
    00925
    00926     dfNoDataValue = poBand->GetNoDataValue( &bGotNoDataValue );
    00927
    00928     poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
    00929
    00930     nBlocksPerRow = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
    00931     nBlocksPerColumn = (poBand->GetYSize() + nBlockYSize - 1) / nBlockYSize;
    00932
    00933     nBlockPixels = nBlockXSize * nBlockYSize;
    00934     nBlockCount = nBlocksPerRow * nBlocksPerColumn;
    00935
    00936     if( nBlocksPerRow == 0 || nBlocksPerColumn == 0 || nBlockPixels == 0
    00937         || nBlockCount == 0 )
    00938     {
    00939         CPLError( CE_Failure, CPLE_AppDefined,
    00940                   "GDALGetRandomSample(): returning because band"
    00941                   " appears degenerate." );
    00942
    00943         return FALSE;
    00944     }
    00945
    00946     nSampleRate = (int) MAX(1,sqrt((double) nBlockCount)-2.0);
    00947
    00948     if( nSampleRate == nBlocksPerRow && nSampleRate > 1 )
    00949         nSampleRate--;
    00950
    00951     while( nSampleRate > 1
    00952            && ((nBlockCount-1) / nSampleRate + 1) * nBlockPixels < nSamples )
    00953         nSampleRate--;
    00954
    00955     nBlockSampleRate =
    00956         MAX(1,nBlockPixels / (nSamples / ((nBlockCount-1) / nSampleRate + 1)));
    00957
    00958     for( int iSampleBlock = 0;
    00959          iSampleBlock < nBlockCount;
    00960          iSampleBlock += nSampleRate )
    00961     {
    00962         double dfValue = 0.0, dfReal, dfImag;
    00963         int  iXBlock, iYBlock, iX, iY, iXValid, iYValid, iRemainder = 0;
    00964         GDALRasterBlock *poBlock;
    00965
    00966         iYBlock = iSampleBlock / nBlocksPerRow;
    00967         iXBlock = iSampleBlock - nBlocksPerRow * iYBlock;
    00968
    00969         poBlock = poBand->GetBlockRef( iXBlock, iYBlock );
    00970         if( poBlock == NULL )
    00971             continue;
    00972
    00973         if( iXBlock * nBlockXSize > poBand->GetXSize() )
    00974             iXValid = poBand->GetXSize() - iXBlock * nBlockXSize;
    00975         else
    00976             iXValid = nBlockXSize;
    00977
    00978         if( iYBlock * nBlockYSize > poBand->GetYSize() )
    00979             iYValid = poBand->GetYSize() - iYBlock * nBlockYSize;
    00980         else
    00981             iYValid = nBlockYSize;
    00982
    00983         for( iY = 0; iY < iYValid; iY++ )
    00984         {
    00985             for( iX = iRemainder; iX < iXValid; iX += nBlockSampleRate )
    00986             {
    00987                 int     iOffset;
    00988
    00989                 iOffset = iX + iY * nBlockXSize;
    00990                 switch( poBlock->GetDataType() )
    00991                 {
    00992                   case GDT_Byte:
    00993                     dfValue = ((GByte *) poBlock->GetDataRef())[iOffset];
    00994                     break;
    00995                   case GDT_UInt16:
    00996                     dfValue = ((GUInt16 *) poBlock->GetDataRef())[iOffset];
    00997                     break;
    00998                   case GDT_Int16:
    00999                     dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset];
    01000                     break;
    01001                   case GDT_UInt32:
    01002                     dfValue = ((GUInt32 *) poBlock->GetDataRef())[iOffset];
    01003                     break;
    01004                   case GDT_Int32:
    01005                     dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset];
    01006                     break;
    01007                   case GDT_Float32:
    01008                     dfValue = ((float *) poBlock->GetDataRef())[iOffset];
    01009                     break;
    01010                   case GDT_Float64:
    01011                     dfValue = ((double *) poBlock->GetDataRef())[iOffset];
    01012                     break;
    01013                   case GDT_CInt16:
    01014                     dfReal = ((GInt16 *) poBlock->GetDataRef())[iOffset*2];
    01015                     dfImag = ((GInt16 *) poBlock->GetDataRef())[iOffset*2+1];
    01016                     dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
    01017                     break;
    01018                   case GDT_CInt32:
    01019                     dfReal = ((GInt32 *) poBlock->GetDataRef())[iOffset*2];
    01020                     dfImag = ((GInt32 *) poBlock->GetDataRef())[iOffset*2+1];
    01021                     dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
    01022                     break;
    01023                   case GDT_CFloat32:
    01024                     dfReal = ((float *) poBlock->GetDataRef())[iOffset*2];
    01025                     dfImag = ((float *) poBlock->GetDataRef())[iOffset*2+1];
    01026                     dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
    01027                     break;
    01028                   case GDT_CFloat64:
    01029                     dfReal = ((double *) poBlock->GetDataRef())[iOffset*2];
    01030                     dfImag = ((double *) poBlock->GetDataRef())[iOffset*2+1];
    01031                     dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
    01032                     break;
    01033                   default:
    01034                     CPLAssert( FALSE );
    01035                 }
    01036
    01037                 if( bGotNoDataValue && dfValue == dfNoDataValue )
    01038                     continue;
    01039
    01040                 if( nActualSamples < nSamples )
    01041                     pafSampleBuf[nActualSamples++] = dfValue;
    01042             }
    01043
    01044             iRemainder = iX - iXValid;
    01045         }
    01046     }
    01047
    01048     return nActualSamples;
    01049 }
    01050
    01051 /************************************************************************/
    01052 /*                            GDALInitGCPs()                            */
    01053 /************************************************************************/
    01054
    01055 void GDALInitGCPs( int nCount, GDAL_GCP * psGCP )
    01056
    01057 {
    01058     for( int iGCP = 0; iGCP < nCount; iGCP++ )
    01059     {
    01060         memset( psGCP, 0, sizeof(GDAL_GCP) );
    01061         psGCP->pszId = CPLStrdup("");
    01062         psGCP->pszInfo = CPLStrdup("");
    01063         psGCP++;
    01064     }
    01065 }
    01066
    01067 /************************************************************************/
    01068 /*                           GDALDeinitGCPs()                           */
    01069 /************************************************************************/
    01070
    01071 void GDALDeinitGCPs( int nCount, GDAL_GCP * psGCP )
    01072
    01073 {
    01074     for( int iGCP = 0; iGCP < nCount; iGCP++ )
    01075     {
    01076         CPLFree( psGCP->pszId );
    01077         CPLFree( psGCP->pszInfo );
    01078         psGCP++;
    01079     }
    01080 }
    01081
    01082 /************************************************************************/
    01083 /*                         GDALDuplicateGCPs()                          */
    01084 /************************************************************************/
    01085
    01086 GDAL_GCP *GDALDuplicateGCPs( int nCount, const GDAL_GCP *pasGCPList )
    01087
    01088 {
    01089     GDAL_GCP    *pasReturn;
    01090
    01091     pasReturn = (GDAL_GCP *) CPLMalloc(sizeof(GDAL_GCP) * nCount);
    01092     GDALInitGCPs( nCount, pasReturn );
    01093
    01094     for( int iGCP = 0; iGCP < nCount; iGCP++ )
    01095     {
    01096         CPLFree( pasReturn[iGCP].pszId );
    01097         pasReturn[iGCP].pszId = CPLStrdup( pasGCPList[iGCP].pszId );
    01098
    01099         CPLFree( pasReturn[iGCP].pszInfo );
    01100         pasReturn[iGCP].pszInfo = CPLStrdup( pasGCPList[iGCP].pszInfo );
    01101
    01102         pasReturn[iGCP].dfGCPPixel = pasGCPList[iGCP].dfGCPPixel;
    01103         pasReturn[iGCP].dfGCPLine = pasGCPList[iGCP].dfGCPLine;
    01104         pasReturn[iGCP].dfGCPX = pasGCPList[iGCP].dfGCPX;
    01105         pasReturn[iGCP].dfGCPY = pasGCPList[iGCP].dfGCPY;
    01106         pasReturn[iGCP].dfGCPZ = pasGCPList[iGCP].dfGCPZ;
    01107     }
    01108
    01109     return pasReturn;
    01110 }
    01111
    01112 /************************************************************************/
    01113 /*                         GDALReadTabFile()                            */
    01114 /*                                                                      */
    01115 /*      Helper function for translator implementators wanting           */
    01116 /*      support for MapInfo .tab-files.                                 */
    01117 /************************************************************************/
    01118
    01119 #define MAX_GCP 256
    01120
    01121 int GDALReadTabFile( const char * pszBaseFilename,
    01122                      double *padfGeoTransform, char **ppszWKT,
    01123                      int *pnGCPCount, GDAL_GCP **ppasGCPs )
    01124
    01125
    01126 {
    01127     const char  *pszTAB;
    01128     FILE        *fpTAB;
    01129     char        **papszLines;
    01130     char    **papszTok=NULL;
    01131     int         bTypeRasterFound = FALSE;
    01132     int         bInsideTableDef = FALSE;
    01133     int         iLine, numLines=0;
    01134     int         nCoordinateCount = 0;
    01135     GDAL_GCP    asGCPs[MAX_GCP];
    01136
    01137
    01138 /* -------------------------------------------------------------------- */
    01139 /*      Try lower case, then upper case.                                */
    01140 /* -------------------------------------------------------------------- */
    01141     pszTAB = CPLResetExtension( pszBaseFilename, "tab" );
    01142
    01143     fpTAB = VSIFOpen( pszTAB, "rt" );
    01144
    01145 #ifndef WIN32
    01146     if( fpTAB == NULL )
    01147     {
    01148         pszTAB = CPLResetExtension( pszBaseFilename, "TAB" );
    01149         fpTAB = VSIFOpen( pszTAB, "rt" );
    01150     }
    01151 #endif
    01152
    01153     if( fpTAB == NULL )
    01154         return FALSE;
    01155
    01156     VSIFClose( fpTAB );
    01157
    01158 /* -------------------------------------------------------------------- */
    01159 /*      We found the file, now load and parse it.                       */
    01160 /* -------------------------------------------------------------------- */
    01161     papszLines = CSLLoad( pszTAB );
    01162
    01163     numLines = CSLCount(papszLines);
    01164
    01165     // Iterate all lines in the TAB-file
    01166     for(iLine=0; iLine<numLines; iLine++)
    01167     {
    01168         CSLDestroy(papszTok);
    01169         papszTok = CSLTokenizeStringComplex(papszLines[iLine], " \t(),;",
    01170                                             TRUE, FALSE);
    01171
    01172         if (CSLCount(papszTok) < 2)
    01173             continue;
    01174
    01175         // Did we find table definition
    01176         if (EQUAL(papszTok[0], "Definition") && EQUAL(papszTok[1], "Table") )
    01177         {
    01178             bInsideTableDef = TRUE;
    01179         }
    01180         else if (bInsideTableDef && (EQUAL(papszTok[0], "Type")) )
    01181         {
    01182             // Only RASTER-type will be handled
    01183             if (EQUAL(papszTok[1], "RASTER"))
    01184             {
    01185                 bTypeRasterFound = TRUE;
    01186             }
    01187             else
    01188             {
    01189                 CSLDestroy(papszTok);
    01190                 CSLDestroy(papszLines);
    01191                 return FALSE;
    01192             }
    01193         }
    01194         else if (bTypeRasterFound && bInsideTableDef
    01195                  && CSLCount(papszTok) > 5
    01196                  && EQUAL(papszTok[4], "Label")
    01197                  && nCoordinateCount < MAX_GCP )
    01198         {
    01199             GDALInitGCPs( 1, asGCPs + nCoordinateCount );
    01200
    01201             asGCPs[nCoordinateCount].dfGCPPixel = atof(papszTok[2]);
    01202             asGCPs[nCoordinateCount].dfGCPLine = atof(papszTok[3]);
    01203             asGCPs[nCoordinateCount].dfGCPX = atof(papszTok[0]);
    01204             asGCPs[nCoordinateCount].dfGCPY = atof(papszTok[1]);
    01205             CPLFree( asGCPs[nCoordinateCount].pszId );
    01206             asGCPs[nCoordinateCount].pszId = CPLStrdup(papszTok[5]);
    01207
    01208             nCoordinateCount++;
    01209         }
    01210         else if( bTypeRasterFound && bInsideTableDef
    01211                  && EQUAL(papszTok[0],"CoordSys")
    01212                  && ppszWKT != NULL )
    01213         {
    01214 #ifdef HAVE_MITAB
    01215             OGRSpatialReference *poSRS = NULL;
    01216
    01217             poSRS = MITABCoordSys2SpatialRef( papszLines[iLine] );
    01218             if( poSRS != NULL )
    01219             {
    01220                 poSRS->exportToWkt( ppszWKT );
    01221                 delete poSRS;
    01222             }
    01223 #else
    01224             CPLDebug( "GDAL", "GDALReadTabFile(): Found `%s',\n"
    01225                  "but GDALReadTabFile() not configured with MITAB callout.",
    01226                       papszLines[iLine] );
    01227 #endif
    01228         }
    01229     }
    01230
    01231 /* -------------------------------------------------------------------- */
    01232 /*      Try to convert the GCPs into a geotransform definition, if      */
    01233 /*      possible.  Otherwise we will need to use them as GCPs.          */
    01234 /* -------------------------------------------------------------------- */
    01235     if( !GDALGCPsToGeoTransform( nCoordinateCount, asGCPs, padfGeoTransform,
    01236                                  FALSE ) )
    01237     {
    01238         CPLDebug( "GDAL",
    01239                   "GDALReadTabFile(%s) found file, wasn't able to derive a\n"
    01240                   "first order geotransform.  Using points as GCPs.",
    01241                   pszTAB );
    01242
    01243         *ppasGCPs = (GDAL_GCP *)
    01244             CPLCalloc(sizeof(GDAL_GCP),nCoordinateCount);
    01245         memcpy( *ppasGCPs, asGCPs, sizeof(GDAL_GCP) * nCoordinateCount );
    01246         *pnGCPCount = nCoordinateCount;
    01247     }
    01248     else
    01249     {
    01250         GDALDeinitGCPs( nCoordinateCount, asGCPs );
    01251     }
    01252
    01253     CSLDestroy(papszTok);
    01254     CSLDestroy(papszLines);
    01255     return *pnGCPCount == 0;
    01256 }
    01257
    01258 /************************************************************************/
    01259 /*                         GDALReadWorldFile()                          */
    01260 /*                                                                      */
    01261 /*      Helper function for translator implementators wanting           */
    01262 /*      support for ESRI world files.                                   */
    01263 /************************************************************************/
    01264
    01265 int GDALReadWorldFile( const char * pszBaseFilename, const char *pszExtension,
    01266                        double *padfGeoTransform )
    01267
    01268 {
    01269     const char  *pszTFW;
    01270     char        szExtUpper[32], szExtLower[32];
    01271     int         i;
    01272     FILE        *fpTFW;
    01273     char        **papszLines;
    01274
    01275     if( *pszExtension == '.' )
    01276         pszExtension++;
    01277
    01278 /* -------------------------------------------------------------------- */
    01279 /*      Generate upper and lower case versions of the extension.        */
    01280 /* -------------------------------------------------------------------- */
    01281     strncpy( szExtUpper, pszExtension, 32 );
    01282     strncpy( szExtLower, pszExtension, 32 );
    01283
    01284     for( i = 0; szExtUpper[i] != '\0' && i < 32; i++ )
    01285     {
    01286         szExtUpper[i] = toupper(szExtUpper[i]);
    01287         szExtLower[i] = tolower(szExtLower[i]);
    01288     }
    01289
    01290 /* -------------------------------------------------------------------- */
    01291 /*      Try lower case, then upper case.                                */
    01292 /* -------------------------------------------------------------------- */
    01293     pszTFW = CPLResetExtension( pszBaseFilename, szExtLower );
    01294
    01295     fpTFW = VSIFOpen( pszTFW, "rt" );
    01296
    01297 #ifndef WIN32
    01298     if( fpTFW == NULL )
    01299     {
    01300         pszTFW = CPLResetExtension( pszBaseFilename, szExtUpper );
    01301         fpTFW = VSIFOpen( pszTFW, "rt" );
    01302     }
    01303 #endif
    01304
    01305     if( fpTFW == NULL )
    01306         return FALSE;
    01307
    01308     VSIFClose( fpTFW );
    01309
    01310 /* -------------------------------------------------------------------- */
    01311 /*      We found the file, now load and parse it.                       */
    01312 /* -------------------------------------------------------------------- */
    01313     papszLines = CSLLoad( pszTFW );
    01314     if( CSLCount(papszLines) >= 6
    01315         && atof(papszLines[0]) != 0.0
    01316         && atof(papszLines[3]) != 0.0 )
    01317     {
    01318         padfGeoTransform[0] = atof(papszLines[4]);
    01319         padfGeoTransform[1] = atof(papszLines[0]);
    01320         padfGeoTransform[2] = atof(papszLines[2]);
    01321         padfGeoTransform[3] = atof(papszLines[5]);
    01322         padfGeoTransform[4] = atof(papszLines[1]);
    01323         padfGeoTransform[5] = atof(papszLines[3]);
    01324
    01325         // correct for center of pixel vs. top left of pixel
    01326         padfGeoTransform[0] -= 0.5 * padfGeoTransform[1];
    01327         padfGeoTransform[0] -= 0.5 * padfGeoTransform[2];
    01328         padfGeoTransform[3] -= 0.5 * padfGeoTransform[4];
    01329         padfGeoTransform[3] -= 0.5 * padfGeoTransform[5];
    01330
    01331         CSLDestroy(papszLines);
    01332
    01333         return TRUE;
    01334     }
    01335     else
    01336     {
    01337         CPLDebug( "GDAL",
    01338                   "GDALReadWorldFile(%s) found file, but it was corrupt.",
    01339                   pszTFW );
    01340         CSLDestroy(papszLines);
    01341         return FALSE;
    01342     }
    01343 }
    01344
    01345 /************************************************************************/
    01346 /*                         GDALWriteWorldFile()                          */
    01347 /*                                                                      */
    01348 /*      Helper function for translator implementators wanting           */
    01349 /*      support for ESRI world files.                                   */
    01350 /************************************************************************/
    01351
    01352 int GDALWriteWorldFile( const char * pszBaseFilename, const char *pszExtension,
    01353                        double *padfGeoTransform )
    01354
    01355 {
    01356     const char  *pszTFW;
    01357         FILE    *fpTFW;
    01358
    01359         pszTFW = CPLResetExtension( pszBaseFilename, pszExtension );
    01360         fpTFW = VSIFOpen( pszTFW, "wt" );
    01361         if( fpTFW == NULL )
    01362         return FALSE;
    01363
    01364 /* -------------------------------------------------------------------- */
    01365 /*      We open the file, now fill it with the world data.                        */
    01366 /* -------------------------------------------------------------------- */
    01367         fprintf( fpTFW, "%.10f\n", padfGeoTransform[1] );
    01368         fprintf( fpTFW, "%.10f\n", padfGeoTransform[4] );
    01369         fprintf( fpTFW, "%.10f\n", padfGeoTransform[2] );
    01370         fprintf( fpTFW, "%.10f\n", padfGeoTransform[5] );
    01371         fprintf( fpTFW, "%.10f\n", padfGeoTransform[0]
    01372                  + 0.5 * padfGeoTransform[1]
    01373                  + 0.5 * padfGeoTransform[2] );
    01374         fprintf( fpTFW, "%.10f\n", padfGeoTransform[3]
    01375                  + 0.5 * padfGeoTransform[4]
    01376                  + 0.5 * padfGeoTransform[5] );
    01377
    01378         VSIFClose( fpTFW );
    01379         return TRUE;
    01380 }
    01381
    01382 /************************************************************************/
    01383 /*                          GDALVersionInfo()                           */
    01384 /************************************************************************/
    01385
    01404 const char *GDALVersionInfo( const char *pszRequest )
    01405
    01406 {
    01407     static char szResult[128];
    01408
    01409
    01410     if( pszRequest == NULL || EQUAL(pszRequest,"VERSION_NUM") )
    01411         sprintf( szResult, "%d", GDAL_VERSION_NUM );
    01412     else if( EQUAL(pszRequest,"RELEASE_DATE") )
    01413         sprintf( szResult, "%d", GDAL_RELEASE_DATE );
    01414     else if( EQUAL(pszRequest,"RELEASE_NAME") )
    01415         sprintf( szResult, "%s", GDAL_RELEASE_NAME );
    01416     else // --version
    01417         sprintf( szResult, "GDAL %s, released %d/%02d/%02d",
    01418                  GDAL_RELEASE_NAME,
    01419                  GDAL_RELEASE_DATE / 10000,
    01420                  (GDAL_RELEASE_DATE % 10000) / 100,
    01421                  GDAL_RELEASE_DATE % 100 );
    01422
    01423     return szResult;
    01424 }
    01425
    01426 /************************************************************************/
    01427 /*                            GDALDecToDMS()                            */
    01428 /*                                                                      */
    01429 /*      Translate a decimal degrees value to a DMS string with          */
    01430 /*      hemisphere.                                                     */
    01431 /************************************************************************/
    01432
    01433 const char *GDALDecToDMS( double dfAngle, const char * pszAxis,
    01434                           int nPrecision )
    01435
    01436 {
    01437     return CPLDecToDMS( dfAngle, pszAxis, nPrecision );
    01438 }
    01439
    01440 /************************************************************************/
    01441 /*                       GDALGCPsToGeoTransform()                       */
    01442 /************************************************************************/
    01443
    01460 int GDALGCPsToGeoTransform( int nGCPCount, const GDAL_GCP *pasGCPs,
    01461                             double *padfGeoTransform, int bApproxOK )
    01462
    01463 {
    01464     int   iAnchor=0, iPnt1, iPnt2, i;
    01465     double adfDPixel[2], adfDLine[2], adfDX[2], adfDY[2];
    01466
    01467 /* -------------------------------------------------------------------- */
    01468 /*      Recognise a few special cases.                                  */
    01469 /* -------------------------------------------------------------------- */
    01470     if( nGCPCount < 2 )
    01471         return FALSE;
    01472
    01473     if( nGCPCount == 2 )
    01474     {
    01475         if( pasGCPs[1].dfGCPPixel == pasGCPs[0].dfGCPPixel
    01476             || pasGCPs[1].dfGCPLine == pasGCPs[0].dfGCPLine )
    01477             return FALSE;
    01478
    01479         padfGeoTransform[1] = (pasGCPs[1].dfGCPX - pasGCPs[0].dfGCPX)
    01480             / (pasGCPs[1].dfGCPPixel - pasGCPs[0].dfGCPPixel);
    01481         padfGeoTransform[2] = 0.0;
    01482
    01483         padfGeoTransform[4] = 0.0;
    01484         padfGeoTransform[5] = (pasGCPs[1].dfGCPY - pasGCPs[0].dfGCPY)
    01485             / (pasGCPs[1].dfGCPLine - pasGCPs[0].dfGCPLine);
    01486
    01487         padfGeoTransform[0] = pasGCPs[0].dfGCPX
    01488             - pasGCPs[0].dfGCPPixel * padfGeoTransform[1]
    01489             - pasGCPs[0].dfGCPLine * padfGeoTransform[2];
    01490
    01491         padfGeoTransform[3] = pasGCPs[0].dfGCPY
    01492             - pasGCPs[0].dfGCPPixel * padfGeoTransform[4]
    01493             - pasGCPs[0].dfGCPLine * padfGeoTransform[5];
    01494
    01495         return TRUE;
    01496     }
    01497
    01498 /* -------------------------------------------------------------------- */
    01499 /*      We use the first point as our anchor.  Select two other         */
    01500 /*      points that don't have the same Pixel value to analyse.         */
    01501 /* -------------------------------------------------------------------- */
    01502     iPnt1 = -1;
    01503     iPnt2 = -1;
    01504     for( i = 1; (iPnt1 == -1 || iPnt2 == -1 ) && i < nGCPCount; i++ )
    01505     {
    01506         double dfDPixel = pasGCPs[i].dfGCPPixel-pasGCPs[iAnchor].dfGCPPixel;
    01507         double dfDLine = pasGCPs[i].dfGCPLine - pasGCPs[iAnchor].dfGCPLine;
    01508         double dfDX = pasGCPs[i].dfGCPX - pasGCPs[iAnchor].dfGCPX;
    01509         double dfDY = pasGCPs[i].dfGCPY - pasGCPs[iAnchor].dfGCPY;
    01510
    01511         if( iPnt1 == -1 && ABS(dfDPixel) > 0.001 )
    01512         {
    01513             iPnt1 = i;
    01514             adfDPixel[0] = dfDPixel;
    01515             adfDLine[0] = dfDLine;
    01516             adfDX[0] = dfDX;
    01517             adfDY[0] = dfDY;
    01518         }
    01519         else if( iPnt2 == -1 )
    01520         {
    01521             iPnt2 = i;
    01522             adfDPixel[1] = dfDPixel;
    01523             adfDLine[1] = dfDLine;
    01524             adfDX[1] = dfDX;
    01525             adfDY[1] = dfDY;
    01526         }
    01527     }
    01528
    01529 /* -------------------------------------------------------------------- */
    01530 /*      If necessary, scale one of the points to avoid divide by        */
    01531 /*      zeros.                                                          */
    01532 /* -------------------------------------------------------------------- */
    01533     if( ABS((adfDLine[0] / adfDPixel[0] - adfDLine[1])) < 0.0001 )
    01534     {
    01535         adfDX[1] *= 2;
    01536         adfDY[1] *= 2;
    01537         adfDPixel[1] *= 2;
    01538         adfDLine[1] *= 2;
    01539     }
    01540
    01541 /* -------------------------------------------------------------------- */
    01542 /*      Compute X related coefficients.                                 */
    01543 /* -------------------------------------------------------------------- */
    01544     padfGeoTransform[2] =
    01545         (adfDX[1] - (adfDPixel[1] * adfDX[0]) / adfDPixel[0])
    01546         / (adfDLine[1] - (adfDLine[0]*adfDPixel[1]) / adfDPixel[0]);
    01547
    01548     padfGeoTransform[1] = (adfDX[0] - adfDLine[0] * padfGeoTransform[2])
    01549         / adfDPixel[0];
    01550
    01551 /* -------------------------------------------------------------------- */
    01552 /*      Compute Y related coefficients.                                 */
    01553 /* -------------------------------------------------------------------- */
    01554     padfGeoTransform[5] =
    01555         (adfDY[1] - (adfDPixel[1] * adfDY[0]) / adfDPixel[0])
    01556         / (adfDLine[1] - (adfDLine[0]*adfDPixel[1]) / adfDPixel[0]);
    01557
    01558     padfGeoTransform[4] =
    01559         (adfDY[0] - adfDLine[0] * padfGeoTransform[5]) / adfDPixel[0];
    01560
    01561 /* -------------------------------------------------------------------- */
    01562 /*      Compute top/left origin.                                        */
    01563 /* -------------------------------------------------------------------- */
    01564
    01565     padfGeoTransform[0] = pasGCPs[0].dfGCPX
    01566         - pasGCPs[0].dfGCPPixel * padfGeoTransform[1]
    01567         - pasGCPs[0].dfGCPLine * padfGeoTransform[2];
    01568
    01569     padfGeoTransform[3] = pasGCPs[0].dfGCPY
    01570         - pasGCPs[0].dfGCPPixel * padfGeoTransform[4]
    01571         - pasGCPs[0].dfGCPLine * padfGeoTransform[5];
    01572
    01573 /* -------------------------------------------------------------------- */
    01574 /*      Now check if any of the input points fit this poorly.           */
    01575 /* -------------------------------------------------------------------- */
    01576     if( !bApproxOK )
    01577     {
    01578         double dfPixelSize = ABS(padfGeoTransform[1])
    01579             + ABS(padfGeoTransform[2])
    01580             + ABS(padfGeoTransform[4])
    01581             + ABS(padfGeoTransform[5]);
    01582
    01583         for( i = 0; i < nGCPCount; i++ )
    01584         {
    01585             double      dfErrorX, dfErrorY;
    01586
    01587             dfErrorX =
    01588                 (pasGCPs[i].dfGCPPixel * padfGeoTransform[1]
    01589                  + pasGCPs[i].dfGCPLine * padfGeoTransform[2]
    01590                  + padfGeoTransform[0])
    01591                 - pasGCPs[i].dfGCPX;
    01592             dfErrorY =
    01593                 (pasGCPs[i].dfGCPPixel * padfGeoTransform[4]
    01594                  + pasGCPs[i].dfGCPLine * padfGeoTransform[5]
    01595                  + padfGeoTransform[3])
    01596                 - pasGCPs[i].dfGCPY;
    01597
    01598             if( ABS(dfErrorX) > 0.25 * dfPixelSize
    01599                 || ABS(dfErrorY) > 0.25 * dfPixelSize )
    01600                 return FALSE;
    01601         }
    01602     }
    01603
    01604     return TRUE;
    01605 }
    01606
    01607 /************************************************************************/
    01608 /* -------------------------------------------------------------------- */
    01609 /*      The following stubs are present to ensure that older GDAL       */
    01610 /*      bridges don't fail with newer libraries.                        */
    01611 /* -------------------------------------------------------------------- */
    01612 /************************************************************************/
    01613
    01614 CPL_C_START
    01615
    01616 void *GDALCreateProjDef( const char * pszDef )
    01617 {
    01618     CPLDebug( "GDAL", "GDALCreateProjDef no longer supported." );
    01619     return NULL;
    01620 }
    01621
    01622 CPLErr GDALReprojectToLongLat( void *pDef, double *, double * )
    01623 {
    01624     CPLDebug( "GDAL", "GDALReprojectToLatLong no longer supported." );
    01625     return CE_Failure;
    01626 }
    01627
    01628 CPLErr GDALReprojectFromLongLat( void *pDef, double *, double * )
    01629 {
    01630     CPLDebug( "GDAL", "GDALReprojectFromLatLong no longer supported." );
    01631     return CE_Failure;
    01632 }
    01633
    01634 void GDALDestroyProjDef( void *pDef )
    01635
    01636 {
    01637     CPLDebug( "GDAL", "GDALDestroyProjDef no longer supported." );
    01638 }
    01639
    01640 CPL_C_END