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