GDAL工具箱详解之gdalwarp.exe

xiaoxiao2022-06-22 34

GDAL自带的工具包含了众多的功能，可以看做是利用该开源库进行的典型开发案例。本节主要介绍其中基本的一项工具——gdalwarp.exe。

该工具箱位于GDAL文件夹下的apps文件夹中，属于命令行可执行程序。

使用方法

GDAL帮助文件中栅格工具箱里关于gdalwarp.exe的描述如下：

gdalwarp [--help-general] [--formats] [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"] [-order n | -tps | -rpc | -geoloc] [-et err_threshold] [-refine_gcps tolerance [minimum_gcps]] [-te xmin ymin xmax ymax] [-te_srs srs_def] [-tr xres yres] [-tap] [-ts width height] [-ovr level|AUTO|AUTO-n|NONE] [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16] [-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha [-r resampling_method] [-wm memory_in_mb] [-multi] [-q] [-cutline datasource] [-cl layer] [-cwhere expression] [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline] [-of format] [-co "NAME=VALUE"]* [-overwrite] [-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]* [-doo NAME=VALUE]* srcfile* dstfile 参数意义样例-s_srs原坐标系统+proj=longlat +datum=WGS84 +no_defs-t_srs目标坐标系统+proj=tmerc +lat_0=0 +lon_0=111 +k=1 +x_0=500000 +y_0=0 +a=6378140 +b=6356755.288157528 +units=m +no_defs-to变换选项暂时未知-order n变换多项式顺序1-3-tps强制基于已知控制点使用样条函数进行变换只需输入设定即可-rpc强制使用rpc模型只需输入设定即可-geoloc强制使用元数据中的地理位置数组详细格式请参照：https://trac.osgeo.org/gdal/wiki/rfc4_geolocate-et err_threshold变换允许的偏差值像素为单位-refine_gcps tolerance minimum_gcps精化后控制点的最少数量-te xmin ymin xmax ymax输出文件地理范围最大最小x、y值-te_srs srs_def特殊指定SRS输出文件的地理参照GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n or a file containing the WKT-tr xres yres输出图像的分辨率xres yres-tap输出图像范围的坐标系校正到图像分辨率一致只需输入设定即可-ts width height设置输出图像的宽和高像素数-ovr levelAUTOAUTO-n-wo “NAME=VALUE”这里可以特别指定重投影选项重投影选项都可以-ot type输出数据类型type-wt type处理时使用的数据类型type-r resampling_method设定重采样算法resampling_method-srcnodata value [value…]源图像无数据处的值value-dstnodata value [value…]输出图像无数据处的值value-dstalpha创建输出数据集中的透明波段以辨识无数据值只需输入设定即可-wm memory_in_mb设定该应用程序内存memory_in_mb-multi是否使用多线程只需输入设定即可-q是否在控制台输出运行步骤提示只需输入设定即可-of format输出格式format-co “NAME=VALUE”格式驱动所支持的选项-cutline datasource激活OGR支持的矢量作为数据分割线datasource-cl layername指定特定名称的图层layername-cwhere expression基于属性查询分割线expression-csql query使用sql语句选择在-cl指定图层上的分割线要素-cblend distance设置混合距离按像素-crop_to_cutline按照分割线裁剪数据集只需输入设定即可-overwrite是否覆写目标数据集只需输入设定即可-nomd不复制元数据只需输入设定即可-cvmd meta_conflict_value冲突元数据信息值“”-setci参照输入数据集的色彩索引表建立输出数据集波段的色彩索引表只需输入设定即可-oo NAME=VALUE打开选项驱动打开选项-doo NAME=VALUE输出选项驱动输出选项srcfile输入文件raw_spot.tifdstfile输出文件utm11.tif

关于-oo及doo的叙述在上一节gdalinfo中有所涉及，这里就不赘述了。借助于QGIS中集成的的GDAL可视化界面我们可以简易地生成命令，下图是该工具的GUI 在使用该界面设置完成一些参数以后可以在后面的文本框中得到命令行参数。

gdalwarp -ot Float32 -s_srs EPSG:4326 -t_srs EPSG:2382 -r bilinear -of GTiff -te 116.526854182 40.0601709856 116.532848182 40.0661649856 -te_srs EPSG:4326 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -wo OPTIMIZE_SIZE=TRUE E:\GDALhomework\000000.tif E://warped.tif

启用一些高级设置，得到：

gdalwarp -ot UInt16 -s_srs EPSG:4326 -t_srs EPSG:2382 -r bilinear -of GTiff -te 116.526854182 40.0481829856 116.532848182 40.0541769856 -te_srs EPSG:4326 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -co TILED=YES -wo OPTIMIZE_SIZE=TRUE E:\GDALhomework\000002.tif E://warped1.tif

在QGIS中打开变换后的图像，对比如图下层真彩色的图像为原图像，上层发黄的是变换后的影像。可以看出两幅影像相同的地理位置处的影像是相同的，也就是变换不会损坏图像的信息，但从旁边的黑色部分可以看出这一地理位置的信息是没有的。同时，查看两幅图像的信息，发现转换后的图像经过了

经过了压缩行列数目发生了变化图像顶点及范围变化像元大小发生变化参考系在参考椭球和高斯平面坐标间发生了变化（最重要）

下面，在控制台中直接使用原生程序gdalwarp.exe：进入控制台，并键入

D:\gdal-2.1.0\apps>gdalwarp -ot UInt16 -s_srs EPSG:4326 -t_srs EPSG:2382 -r bilinear -of GTiff -te 116.526854182 40.0481829856 116.532848182 40.0541769856 - te_srs EPSG:4326 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -co TILED=YES -wo OPTIMIZE_SIZE=TRUE E:\GDALhomework\000002.tif E:\\warpedcmd.tif 0...10...20...30...40...50...60...70...80...90...100 - done. bit length overflow code 6 bits 6->7 bit length overflow code 14 bits 6->7 code 13 bits 6->7

这样转换的图像和在QGIS中操作的结果是一样的。这里只是提供一种间接学习使用GDAL工具箱的方法。下面看程序源码。

源码分析

应用程序主函数体：

/************************************************************************/ /* 主函数() /************************************************************************/ int main( int argc, char ** argv ) { GDALDatasetH *pahSrcDS = NULL; int nSrcCount = 0; EarlySetConfigOptions(argc, argv); /* -------------------------------------------------------------------- */ /* 注册所有驱动,并且处理gdla常规命令选项。 /* -------------------------------------------------------------------- */ GDALAllRegister(); argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 ); if( argc < 1 ) GDALExit( -argc ); for( int i = 0; argv != NULL && argv[i] != NULL; i++ ) { if( EQUAL(argv[i], "--utility_version") ) { printf("%s was compiled against GDAL %s and is running against GDAL %s\n", argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME")); CSLDestroy( argv ); return 0; } else if( EQUAL(argv[i],"--help") ) { Usage(NULL); } } /* -------------------------------------------------------------------- */ /* 进行必要设置以应对超长型输入。 /* 将最大数据集池的大小设置为450的原因是:常规的Linux进程仅允许1024个 /* 文件描述词但我们还需要留出一部分空间给共享库，就剩了900个 /* 一些数据集需要两个文件描述词，所以稳妥起见，除以2得到了450。 /* -------------------------------------------------------------------- */ if( CPLGetConfigOption("GDAL_MAX_DATASET_POOL_SIZE", NULL) == NULL ) { CPLSetConfigOption("GDAL_MAX_DATASET_POOL_SIZE", "450"); } GDALWarpAppOptionsForBinary* psOptionsForBinary = GDALWarpAppOptionsForBinaryNew();/*二值化选项申请内存*/ GDALWarpAppOptions *psOptions = GDALWarpAppOptionsNew(argv + 1, psOptionsForBinary);/*根据命令行输入参数填充参数*/ CSLDestroy( argv );/*释放字符串数组*/ if( psOptions == NULL ) { Usage(NULL); } if( psOptionsForBinary->pszDstFilename == NULL ) { Usage("No target filename specified."); } if ( CSLCount(psOptionsForBinary->papszSrcFiles) == 1 && strcmp(psOptionsForBinary->papszSrcFiles[0], psOptionsForBinary->pszDstFilename) == 0 && psOptionsForBinary->bOverwrite) { CPLError(CE_Failure, CPLE_IllegalArg, "Source and destination datasets must be different.\n"); GDALExit(1);/*源文件和目标文件相同且唯一时又指定要求覆盖，将会弹出错误消息并退出。 } /* -------------------------------------------------------------------- */ /* 打开源文件. /* -------------------------------------------------------------------- */ for(int i = 0; psOptionsForBinary->papszSrcFiles[i] != NULL; i++) { nSrcCount++; pahSrcDS = (GDALDatasetH *) CPLRealloc(pahSrcDS, sizeof(GDALDatasetH) * nSrcCount); pahSrcDS[nSrcCount-1] = GDALOpenEx( psOptionsForBinary->papszSrcFiles[i], GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, NULL, (const char* const* )psOptionsForBinary->papszOpenOptions, NULL );/*将打开选项这个字符串数组的指针及数组内容全部固定，以适应函数输入参数的类型*/ if( pahSrcDS[nSrcCount-1] == NULL ) GDALExit(2); } /* -------------------------------------------------------------------- */ /* 输出数据集是否已经存在? /* -------------------------------------------------------------------- */ /* FIXME ? source filename=target filename and -overwrite is definitely */ /* an error. But I can't imagine of a valid case (without -overwrite), */ /* where it would make sense. In doubt, let's keep that dubious possibility... */ //这是作者的一段修正标识 int bOutStreaming = FALSE; if( strcmp(psOptionsForBinary->pszDstFilename, "/vsistdout/") == 0 ) { psOptionsForBinary->bQuiet = TRUE; bOutStreaming = TRUE; } #ifdef S_ISFIFO else { VSIStatBufL sStat; if( VSIStatExL(psOptionsForBinary->pszDstFilename, &sStat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG) == 0 && S_ISFIFO(sStat.st_mode) ) { bOutStreaming = TRUE; } } #endif GDALDatasetH hDstDS = NULL; if( bOutStreaming ) { GDALWarpAppOptionsSetWarpOption(psOptions, "STREAMABLE_OUTPUT", "YES"); } else { CPLPushErrorHandler( CPLQuietErrorHandler ); hDstDS = GDALOpenEx( psOptionsForBinary->pszDstFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR | GDAL_OF_UPDATE, NULL, psOptionsForBinary->papszDestOpenOptions, NULL ); CPLPopErrorHandler(); } if( hDstDS != NULL && psOptionsForBinary->bOverwrite ) { GDALClose(hDstDS); hDstDS = NULL; } if( hDstDS != NULL && psOptionsForBinary->bCreateOutput ) { CPLError( CE_Failure, CPLE_AppDefined, "Output dataset %s exists,\n" "but some command line options were provided indicating a new dataset\n" "should be created. Please delete existing dataset and run again.\n", psOptionsForBinary->pszDstFilename ); GDALExit(1); } /* 避免覆盖一个已经存在但即便是用提升模式也打不开的新的GTiff类型的目标文件 */ if ( !bOutStreaming && hDstDS == NULL && !psOptionsForBinary->bOverwrite ) { CPLPushErrorHandler( CPLQuietErrorHandler ); hDstDS = GDALOpen( psOptionsForBinary->pszDstFilename, GA_ReadOnly ); CPLPopErrorHandler(); if (hDstDS) { CPLError( CE_Failure, CPLE_AppDefined, "Output dataset %s exists, but cannot be opened in update mode\n", psOptionsForBinary->pszDstFilename ); GDALClose(hDstDS); GDALExit(1); } } if( !(psOptionsForBinary->bQuiet) ) { GDALWarpAppOptionsSetProgress(psOptions, GDALTermProgress, NULL); } if (hDstDS == NULL && !psOptionsForBinary->bQuiet && !psOptionsForBinary->bFormatExplicitlySet) CheckExtensionConsistency(psOptionsForBinary->pszDstFilename, psOptionsForBinary->pszFormat); int bUsageError = FALSE; GDALDatasetH hOutDS = GDALWarp(psOptionsForBinary->pszDstFilename, hDstDS, nSrcCount, pahSrcDS, psOptions, &bUsageError);//核心函数 if( bUsageError ) Usage();//错误处理 int nRetCode = (hOutDS) ? 0 : 1; GDALWarpAppOptionsFree(psOptions);/*释放选项*/ GDALWarpAppOptionsForBinaryFree(psOptionsForBinary);/*释放二值化选项*/ for(int i = 0; i < nSrcCount; i++) { GDALClose(pahSrcDS[i]); } CPLFree(pahSrcDS); GDALClose( hOutDS ? hOutDS : hDstDS ); GDALDumpOpenDatasets( stderr );/*关闭打开的数据集*/ GDALDestroyDriverManager();/*释放驱动管理*/ OGRCleanupAll(); return nRetCode; }

从上篇文章的经验可以看出，处理过程是相似的，这里附上该应用程序程序的流程图。所以这里就直接进入核心函数体进行学习。

变换库核心函数

函数体及其注释：

/** * Image reprojection and warping function. * 图像重投影及变换函数 * This is the equivalent of the <a href="gdal_warp.html">gdalwarp</a> utility. * * GDALWarpAppOptions* must be allocated and freed with GDALWarpAppOptionsNew() * and GDALWarpAppOptionsFree() respectively. * pszDest and hDstDS cannot be used at the same time. * * @param pszDest the destination dataset path or NULL. * @param hDstDS the destination dataset or NULL. * @param nSrcCount the number of input datasets. * @param pahSrcDS the list of input datasets. * @param psOptionsIn the options struct returned by GDALWarpAppOptionsNew() or NULL. * @param pbUsageError the pointer to int variable to determine any usage error has occurred. * @return the output dataset (new dataset that must be closed using GDALClose(), or hDstDS is not NULL) or NULL in case of error. * * @since GDAL 2.1 */ GDALDatasetH GDALWarp( const char *pszDest, GDALDatasetH hDstDS, int nSrcCount, GDALDatasetH *pahSrcDS, const GDALWarpAppOptions *psOptionsIn, int *pbUsageError ) { CPLErrorReset(); if( pszDest == NULL && hDstDS == NULL ) { CPLError( CE_Failure, CPLE_AppDefined, "pszDest == NULL && hDstDS == NULL");/*输出文件缺失处理*/ if(pbUsageError) *pbUsageError = TRUE; return NULL; } if( pszDest == NULL ) pszDest = GDALGetDescription(hDstDS);/*获取驱动描述*/ int bMustCloseDstDSInCaseOfError = (hDstDS == NULL); GDALTransformerFunc pfnTransformer = NULL; void *hTransformArg = NULL; int bHasGotErr = FALSE; int bVRT = FALSE; void *hCutline = NULL; GDALWarpAppOptions* psOptions = (psOptionsIn) ? GDALWarpAppOptionsClone(psOptionsIn) : GDALWarpAppOptionsNew(NULL, NULL);/*如果没有传入选项则重新初始化一个新的变换选项*/ if( EQUAL(psOptions->pszFormat,"VRT") ) { if( hDstDS != NULL ) { CPLError(CE_Warning, CPLE_NotSupported, "VRT output not compatible with existing dataset."); GDALWarpAppOptionsFree(psOptions); return NULL; } bVRT = TRUE; if( nSrcCount > 1 ) { CPLError(CE_Warning, CPLE_AppDefined, "gdalwarp -of VRT just takes into account " "the first source dataset.\nIf all source datasets " "are in the same projection, try making a mosaic of\n" "them with gdalbuildvrt, and use the resulting " "VRT file as the input of\ngdalwarp -of VRT."); } } /* -------------------------------------------------------------------- */ /* 检查互不兼容的选项 /* -------------------------------------------------------------------- */ if ((psOptions->nForcePixels != 0 || psOptions->nForceLines != 0) && (psOptions->dfXRes != 0 && psOptions->dfYRes != 0)) { CPLError(CE_Failure, CPLE_IllegalArg, "-tr and -ts options cannot be used at the same time.");//-tr和-ts选项不能兼容 if(pbUsageError) *pbUsageError = TRUE; GDALWarpAppOptionsFree(psOptions); return NULL; } if (psOptions->bTargetAlignedPixels && psOptions->dfXRes == 0 && psOptions->dfYRes == 0) { CPLError(CE_Failure, CPLE_IllegalArg, "-tap option cannot be used without using -tr.");//-tap不能被使用除非使用了-tr if(pbUsageError) *pbUsageError = TRUE; GDALWarpAppOptionsFree(psOptions); return NULL; } if( !psOptions->bQuiet && !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 && psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0) ) { if( psOptions->dfMinX >= psOptions->dfMaxX ) CPLError(CE_Warning, CPLE_AppDefined, "-ts values have minx >= maxx. This will result in a horizontally flipped image.");//警告：-ts中minx大于maxx if( psOptions->dfMinY >= psOptions->dfMaxY ) CPLError(CE_Warning, CPLE_AppDefined, "-ts values have miny >= maxy. This will result in a vertically flipped image.");//警告：-ts中miny大于maxy } if( psOptions->dfErrorThreshold < 0 ) { // 默认的, 使用默认的误差容许值除非指定了RPC_DEM if( CSLFetchNameValue(psOptions->papszWarpOptions, "RPC_DEM") != NULL ) psOptions->dfErrorThreshold = 0.0; else psOptions->dfErrorThreshold = 0.125; } /* -------------------------------------------------------------------- */ /* SRS目标区域选项 /* -------------------------------------------------------------------- */ if( psOptions->pszTE_SRS != NULL ) { if( psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 && psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0 ) { CPLError( CE_Warning, CPLE_AppDefined, "-te_srs ignored since -te is not specified."); } else { OGRSpatialReference oSRSIn; oSRSIn.SetFromUserInput(psOptions->pszTE_SRS);//根据传入选项构建空间参照类 OGRSpatialReference oSRSDS; int bOK = FALSE; if( CSLFetchNameValue( psOptions->papszTO, "DST_SRS" ) != NULL ) { oSRSDS.SetFromUserInput( CSLFetchNameValue( psOptions->papszTO, "DST_SRS" ) );//根据选项构建输出空间参照类 bOK = TRUE; } else if( CSLFetchNameValue( psOptions->papszTO, "SRC_SRS" ) != NULL ) { oSRSDS.SetFromUserInput( CSLFetchNameValue( psOptions->papszTO, "SRC_SRS" ) );//如果没有指定输出SRS参照，根据源SRS设置输出SRS参照 bOK = TRUE; } else { if( nSrcCount && pahSrcDS[0] && GDALGetProjectionRef(pahSrcDS[0]) && GDALGetProjectionRef(pahSrcDS[0])[0] )//如果没有源SRS，则取数据集中的投影参考系作为输出参考系 { oSRSDS.SetFromUserInput( GDALGetProjectionRef(pahSrcDS[0]) ); bOK = TRUE; } } if( !bOK ) { CPLError( CE_Failure, CPLE_AppDefined, "-te_srs ignored since none of -t_srs, -s_srs is specified or the input dataset has no projection.");//如果都没有设置成功，则弹出错误 GDALWarpAppOptionsFree(psOptions); return NULL; } if( !oSRSIn.IsSame(&oSRSDS) )//如果输入的SRS参照与输出SRS参照不同 { OGRCoordinateTransformation* poCT = OGRCreateCoordinateTransformation(&oSRSIn, &oSRSDS);//计算转换参数，由于对SRS不甚熟悉，所以内部究竟是如何计算的这里先不深入 if( !(poCT && poCT->Transform(1, &psOptions->dfMinX, &psOptions->dfMinY) && poCT->Transform(1, &psOptions->dfMaxX, &psOptions->dfMaxY)) ) { OGRCoordinateTransformation::DestroyCT(poCT); CPLError( CE_Failure, CPLE_AppDefined, "-te_srs ignored since coordinate transformation failed.");//转换失败就退出 GDALWarpAppOptionsFree(psOptions); return NULL; } delete poCT; } } } /* -------------------------------------------------------------------- */ /* 如果我们有数据源分割线，读取它并把它和变换选项捆绑在一起 /* -------------------------------------------------------------------- */ if( psOptions->pszCutlineDSName != NULL ) { CPLErr eError; eError = LoadCutline( psOptions->pszCutlineDSName, psOptions->pszCLayer, psOptions->pszCWHERE, psOptions->pszCSQL, &hCutline );//读取分割线 if(eError == CE_Failure) { GDALWarpAppOptionsFree(psOptions); return NULL; } } if ( psOptions->bCropToCutline && hCutline != NULL ) { CPLErr eError; eError = CropToCutline( hCutline, psOptions->papszTO, nSrcCount, pahSrcDS, psOptions->dfMinX, psOptions->dfMinY, psOptions->dfMaxX, psOptions->dfMaxY );// if(eError == CE_Failure) { GDALWarpAppOptionsFree(psOptions); OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); return NULL; } } /* -------------------------------------------------------------------- */ /* 如果没有，我们需要创建它 /* -------------------------------------------------------------------- */ void* hUniqueTransformArg = NULL; int bInitDestSetByUser = ( CSLFetchNameValue( psOptions->papszWarpOptions, "INIT_DEST" ) != NULL ); const char* pszWarpThreads = CSLFetchNameValue(psOptions->papszWarpOptions, "NUM_THREADS"); if( pszWarpThreads != NULL ) { /*在使用TPS变换参数的时候启用并行计算来加快矩阵计算 */ psOptions->papszTO = CSLSetNameValue(psOptions->papszTO, "NUM_THREADS", pszWarpThreads); } if( hDstDS == NULL ) { hDstDS = GDALWarpCreateOutput( nSrcCount, pahSrcDS, pszDest,psOptions->pszFormat, psOptions->papszTO, &psOptions->papszCreateOptions, psOptions->eOutputType, &hUniqueTransformArg, psOptions->bSetColorInterpretation, psOptions);//创建变换输出 if(hDstDS == NULL) { if( hUniqueTransformArg ) GDALDestroyTransformer( hUniqueTransformArg ); GDALWarpAppOptionsFree(psOptions); if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); return NULL; } psOptions->bCreateOutput = TRUE; if( !bInitDestSetByUser ) { if ( psOptions->pszDstNodata == NULL ) { psOptions->papszWarpOptions = CSLSetNameValue(psOptions->papszWarpOptions, "INIT_DEST", "0"); } else { psOptions->papszWarpOptions = CSLSetNameValue(psOptions->papszWarpOptions, "INIT_DEST", "NO_DATA" ); } } CSLDestroy( psOptions->papszCreateOptions ); psOptions->papszCreateOptions = NULL; } if( hDstDS == NULL ) { GDALWarpAppOptionsFree(psOptions); return NULL; } /* -------------------------------------------------------------------- */ /* 循环依次处理每个数据源 /* -------------------------------------------------------------------- */ int iSrc; for( iSrc = 0; iSrc < nSrcCount; iSrc++ ) { GDALDatasetH hSrcDS; /* -------------------------------------------------------------------- */ /* 打开文件. /* -------------------------------------------------------------------- */ hSrcDS = pahSrcDS[iSrc]; if( hSrcDS == NULL ) { GDALWarpAppOptionsFree(psOptions); if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); if( bMustCloseDstDSInCaseOfError ) GDALClose(hDstDS); return NULL; } /* -------------------------------------------------------------------- */ /* 检查数据中应至少存在一个波段. /* -------------------------------------------------------------------- */ if ( GDALGetRasterCount(hSrcDS) == 0 ) { CPLError(CE_Failure, CPLE_AppDefined, "Input file %s has no raster bands.", GDALGetDescription(hSrcDS) ); GDALWarpAppOptionsFree(psOptions); if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); if( bMustCloseDstDSInCaseOfError ) GDALClose(hDstDS); return NULL; } if( !psOptions->bQuiet ) printf( "Processing input file %s.\n", GDALGetDescription(hSrcDS) ); /* -------------------------------------------------------------------- */ /* 获取第一个数据集中的元数据并把它拷贝到目标数据集中. /* 仅当源波段和目标波段的数量相等的时候，复制波段级的元数据和其他信息. /* 在两者之间任何冲突的值都将被设置为pszMDConflictValue. /* -------------------------------------------------------------------- */ if ( psOptions->bCopyMetadata ) { char **papszMetadata = NULL; const char *pszSrcInfo = NULL; const char *pszDstInfo = NULL; GDALRasterBandH hSrcBand = NULL; GDALRasterBandH hDstBand = NULL; /* 从第一个数据集中复制元数据信息 */ if ( iSrc == 0 ) { CPLDebug("WARP", "Copying metadata from first source to destination dataset"); /* 复制数据集级别的元数据 */ papszMetadata = GDALGetMetadata( hSrcDS, NULL ); char** papszMetadataNew = NULL; for( int i = 0; papszMetadata != NULL && papszMetadata[i] != NULL; i++ ) { // 当输出包含透明波段的时候，不用提前设置 NODATA_VALUES if( psOptions->bEnableDstAlpha && STARTS_WITH_CI(papszMetadata[i], "NODATA_VALUES=") ) { continue; } papszMetadataNew = CSLAddString(papszMetadataNew, papszMetadata[i]); } if ( CSLCount(papszMetadataNew) > 0 ) { if ( GDALSetMetadata( hDstDS, papszMetadataNew, NULL ) != CE_None ) CPLError( CE_Warning, CPLE_AppDefined, "error copying metadata to destination dataset." ); } CSLDestroy(papszMetadataNew); /* 复制波段级别的元数据信息 */ if ( GDALGetRasterCount( hSrcDS ) == GDALGetRasterCount( hDstDS ) ) { for ( int iBand = 0; iBand < GDALGetRasterCount( hSrcDS ); iBand++ ) { hSrcBand = GDALGetRasterBand( hSrcDS, iBand + 1 ); hDstBand = GDALGetRasterBand( hDstDS, iBand + 1 ); /* 复制元数据except stats (#5319) */ papszMetadata = GDALGetMetadata( hSrcBand, NULL); if ( CSLCount(papszMetadata) > 0 ) { //GDALSetMetadata( hDstBand, papszMetadata, NULL ); papszMetadataNew = NULL; for( int i = 0; papszMetadata != NULL && papszMetadata[i] != NULL; i++ ) { if (!STARTS_WITH(papszMetadata[i], "STATISTICS_")) papszMetadataNew = CSLAddString(papszMetadataNew, papszMetadata[i]); } GDALSetMetadata( hDstBand, papszMetadataNew, NULL ); CSLDestroy(papszMetadataNew); } /* 复制其他信息 (描述, 数据类型等) */ if ( psOptions->bCopyBandInfo ) { pszSrcInfo = GDALGetDescription( hSrcBand ); if( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 ) GDALSetDescription( hDstBand, pszSrcInfo ); pszSrcInfo = GDALGetRasterUnitType( hSrcBand ); if( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 ) GDALSetRasterUnitType( hDstBand, pszSrcInfo ); } } } } /* 移除两个数据集间相互冲突的元数据 */ else { CPLDebug("WARP", "Removing conflicting metadata from destination dataset (source #%d)", iSrc ); /* 移除数据集级别的相互冲突的元数据信息 */ RemoveConflictingMetadata( hDstDS, GDALGetMetadata( hSrcDS, NULL ), psOptions->pszMDConflictValue ); /* 移除波段级别的相互冲突的元数据及其他信息 */ if ( GDALGetRasterCount( hSrcDS ) == GDALGetRasterCount( hDstDS ) ) { for ( int iBand = 0; iBand < GDALGetRasterCount( hSrcDS ); iBand++ ) { hSrcBand = GDALGetRasterBand( hSrcDS, iBand + 1 ); hDstBand = GDALGetRasterBand( hDstDS, iBand + 1 ); /* 移除冲突的元数据 */ RemoveConflictingMetadata( hDstBand, GDALGetMetadata( hSrcBand, NULL ), psOptions->pszMDConflictValue ); /* 移除冲突的其他信息 */ if ( psOptions->bCopyBandInfo ) { pszSrcInfo = GDALGetDescription( hSrcBand ); pszDstInfo = GDALGetDescription( hDstBand ); if( ! ( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 && pszDstInfo != NULL && strlen(pszDstInfo) > 0 && EQUAL( pszSrcInfo, pszDstInfo ) ) ) GDALSetDescription( hDstBand, "" );//若相同将输出波段描述置空 pszSrcInfo = GDALGetRasterUnitType( hSrcBand ); pszDstInfo = GDALGetRasterUnitType( hDstBand ); if( ! ( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 && pszDstInfo != NULL && strlen(pszDstInfo) > 0 && EQUAL( pszSrcInfo, pszDstInfo ) ) ) GDALSetRasterUnitType( hDstBand, "" );//若相同将输出波段数据类型置空 } } } } } /* -------------------------------------------------------------------- */ /* 警告如果输入文件有颜色表,将会询问是否使用最邻近法重采样 /* -------------------------------------------------------------------- */ if ( psOptions->eResampleAlg != GRA_NearestNeighbour && psOptions->eResampleAlg != GRA_Mode && GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL) { if( !psOptions->bQuiet ) CPLError( CE_Warning, CPLE_AppDefined, "Input file %s has a color table, which will likely lead to " "bad results when using a resampling method other than " "nearest neighbour or mode. Converting the dataset prior to 24/32 bit " "is advised.", GDALGetDescription(hSrcDS) ); } /* -------------------------------------------------------------------- */ /* 检查是否存在透明通道? /* -------------------------------------------------------------------- */ int bEnableSrcAlpha = psOptions->bEnableSrcAlpha; if( GDALGetRasterColorInterpretation( GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) ) == GCI_AlphaBand && !bEnableSrcAlpha ) { bEnableSrcAlpha = TRUE; if( !psOptions->bQuiet ) printf( "Using band %d of source image as alpha.\n", GDALGetRasterCount(hSrcDS) ); } /* -------------------------------------------------------------------- */ /* 根据源坐标系与目标坐标系创造一个变换对象 /* -------------------------------------------------------------------- */ if (hUniqueTransformArg) hTransformArg = hUniqueTransformArg; else hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, psOptions->papszTO );//创造常规变换参数 if( hTransformArg == NULL ) { GDALWarpAppOptionsFree(psOptions); if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); if( bMustCloseDstDSInCaseOfError ) GDALClose(hDstDS); return NULL; } pfnTransformer = GDALGenImgProjTransform; /* -------------------------------------------------------------------- */ /* 决定我们是处理全分辨率数据源或者一层缩略级别的数据源 /* -------------------------------------------------------------------- */ GDALDataset* poSrcDS = (GDALDataset*) hSrcDS; GDALDataset* poSrcOvrDS = NULL; int nOvCount = poSrcDS->GetRasterBand(1)->GetOverviewCount(); if( psOptions->nOvLevel <= -2 && nOvCount > 0 ) { double adfSuggestedGeoTransform[6]; double adfExtent[4]; int nPixels, nLines; /* 根据输入数据集和变换函数计算输出图像的分辨率(像素)和图像大小 */ if( GDALSuggestedWarpOutput2(hSrcDS, pfnTransformer, hTransformArg, adfSuggestedGeoTransform, &nPixels, &nLines, adfExtent, 0) == CE_None) { double dfTargetRatio = 1.0 / adfSuggestedGeoTransform[1]; if( dfTargetRatio > 1.0 ) { int iOvr; for( iOvr = -1; iOvr < nOvCount-1; iOvr++ ) { double dfOvrRatio = (iOvr < 0) ? 1.0 : (double)poSrcDS->GetRasterXSize() / poSrcDS->GetRasterBand(1)->GetOverview(iOvr)->GetXSize(); double dfNextOvrRatio = (double)poSrcDS->GetRasterXSize() / poSrcDS->GetRasterBand(1)->GetOverview(iOvr+1)->GetXSize(); if( dfOvrRatio < dfTargetRatio && dfNextOvrRatio > dfTargetRatio ) break; if( fabs(dfOvrRatio - dfTargetRatio) < 1e-1 ) break; } iOvr += (psOptions->nOvLevel+2); if( iOvr >= 0 ) { CPLDebug("WARP", "Selecting overview level %d for %s", iOvr, GDALGetDescription(hSrcDS)); poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE, FALSE ); } } } } else if( psOptions->nOvLevel >= 0 ) { poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, psOptions->nOvLevel, TRUE, FALSE ); if( poSrcOvrDS == NULL ) { if( !psOptions->bQuiet ) { CPLError(CE_Warning, CPLE_AppDefined, "cannot get overview level %d for " "dataset %s. Defaulting to level %d", psOptions->nOvLevel, GDALGetDescription(hSrcDS), nOvCount - 1); } if( nOvCount > 0 ) poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, nOvCount - 1, FALSE, FALSE ); } else { CPLDebug("WARP", "Selecting overview level %d for %s", psOptions->nOvLevel, GDALGetDescription(hSrcDS)); } } GDALDatasetH hWrkSrcDS = (poSrcOvrDS) ? (GDALDatasetH)poSrcOvrDS : hSrcDS; /* 当处理一个鸟瞰图，我们需要重新构建变换 */ if( poSrcOvrDS != NULL ) { GDALDestroyGenImgProjTransformer( hTransformArg ); hTransformArg = GDALCreateGenImgProjTransformer2( hWrkSrcDS, hDstDS, psOptions->papszTO );//变换参数构建函数 } /* -------------------------------------------------------------------- */ /* 使用线性近似器变形变换系数，除非允许误差为零. /* -------------------------------------------------------------------- */ if( psOptions->dfErrorThreshold != 0.0 ) { hTransformArg = GDALCreateApproxTransformer( GDALGenImgProjTransform, hTransformArg, psOptions->dfErrorThreshold); pfnTransformer = GDALApproxTransform; GDALApproxTransformerOwnsSubtransformer(hTransformArg, TRUE); } /* -------------------------------------------------------------------- */ /* 在处理第一张图像后清空选项类中的输出数据集初始化值. /* -------------------------------------------------------------------- */ if( psOptions->bCreateOutput && iSrc == 1 ) psOptions->papszWarpOptions = CSLSetNameValue( psOptions->papszWarpOptions, "INIT_DEST", NULL ); /* -------------------------------------------------------------------- */ /* 装载变换选项. /* -------------------------------------------------------------------- */ GDALWarpOptions *psWO = GDALCreateWarpOptions(); psWO->papszWarpOptions = CSLDuplicate(psOptions->papszWarpOptions); psWO->eWorkingDataType = psOptions->eWorkingType; psWO->eResampleAlg = psOptions->eResampleAlg; psWO->hSrcDS = hWrkSrcDS; psWO->hDstDS = hDstDS; psWO->pfnTransformer = pfnTransformer;//上面生成的变换函数 psWO->pTransformerArg = hTransformArg;//上面生成的变换系数 psWO->pfnProgress = psOptions->pfnProgress; psWO->pProgressArg = psOptions->pProgressData; if( psOptions->dfWarpMemoryLimit != 0.0 ) psWO->dfWarpMemoryLimit = psOptions->dfWarpMemoryLimit; /* -------------------------------------------------------------------- */ /* 装载有色波段. /* -------------------------------------------------------------------- */ if( bEnableSrcAlpha ) psWO->nBandCount = GDALGetRasterCount(hWrkSrcDS) - 1; else psWO->nBandCount = GDALGetRasterCount(hWrkSrcDS); psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int)); psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int)); for( int i = 0; i < psWO->nBandCount; i++ ) { psWO->panSrcBands[i] = i+1; psWO->panDstBands[i] = i+1; } /* -------------------------------------------------------------------- */ /* 如果存在使用的透明波段就装载它. /* -------------------------------------------------------------------- */ if( bEnableSrcAlpha ) psWO->nSrcAlphaBand = GDALGetRasterCount(hWrkSrcDS); int bEnableDstAlpha = psOptions->bEnableDstAlpha; if( !bEnableDstAlpha && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1 && GDALGetRasterColorInterpretation( GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS))) == GCI_AlphaBand ) { if( !psOptions->bQuiet ) printf( "Using band %d of destination image as alpha.\n", GDALGetRasterCount(hDstDS) ); bEnableDstAlpha = TRUE; } if( bEnableDstAlpha ) psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS); /* -------------------------------------------------------------------- */ /* 装载NODATA选项. /* -------------------------------------------------------------------- */ if( psOptions->pszSrcNodata != NULL && !EQUAL(psOptions->pszSrcNodata,"none") ) { char **papszTokens = CSLTokenizeString( psOptions->pszSrcNodata ); int nTokenCount = CSLCount(papszTokens); psWO->padfSrcNoDataReal = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); psWO->padfSrcNoDataImag = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); for( int i = 0; i < psWO->nBandCount; i++ ) { if( i < nTokenCount ) { CPLStringToComplex( papszTokens[i], psWO->padfSrcNoDataReal + i, psWO->padfSrcNoDataImag + i ); } else { psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i-1]; psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i-1]; } } CSLDestroy( papszTokens ); psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "UNIFIED_SRC_NODATA", "YES" ); } /* -------------------------------------------------------------------- */ /* 如果-srcnodata没有被指定, 而数据中有无数据值，那就使用它. /* -------------------------------------------------------------------- */ if( psOptions->pszSrcNodata == NULL ) { int bHaveNodata = FALSE; double dfReal = 0.0; for( int i = 0; !bHaveNodata && i < psWO->nBandCount; i++ ) { GDALRasterBandH hBand = GDALGetRasterBand( hWrkSrcDS, i+1 ); dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata ); } if( bHaveNodata ) { if( !psOptions->bQuiet ) { if (CPLIsNan(dfReal)) printf( "Using internal nodata values (e.g. nan) for image %s.\n", GDALGetDescription(hSrcDS) ); else printf( "Using internal nodata values (e.g. %g) for image %s.\n", dfReal, GDALGetDescription(hSrcDS) ); } psWO->padfSrcNoDataReal = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); psWO->padfSrcNoDataImag = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); for( int i = 0; i < psWO->nBandCount; i++ ) { GDALRasterBandH hBand = GDALGetRasterBand( hWrkSrcDS, i+1 ); dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata ); if( bHaveNodata ) { psWO->padfSrcNoDataReal[i] = dfReal; psWO->padfSrcNoDataImag[i] = 0.0; } else { psWO->padfSrcNoDataReal[i] = -123456.789; psWO->padfSrcNoDataImag[i] = 0.0; } } } } /* -------------------------------------------------------------------- */ /* 如果创建了输出数据集,而有想要指定的无数据值循环给每一波段标明该值. /* -------------------------------------------------------------------- */ if( psOptions->pszDstNodata != NULL && !EQUAL(psOptions->pszDstNodata,"none") ) { char **papszTokens = CSLTokenizeString( psOptions->pszDstNodata ); int nTokenCount = CSLCount(papszTokens); int bDstNoDataNone = TRUE; psWO->padfDstNoDataReal = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); psWO->padfDstNoDataImag = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); for( int i = 0; i < psWO->nBandCount; i++ ) { psWO->padfDstNoDataReal[i] = -1.1e20; psWO->padfDstNoDataImag[i] = 0.0; if( i < nTokenCount ) { if ( papszTokens[i] != NULL && EQUAL(papszTokens[i],"none") ) { CPLDebug( "WARP", "dstnodata of band %d not set", i ); bDstNoDataNone = TRUE; continue; } else if ( papszTokens[i] == NULL ) // this should not happen, but just in case { CPLError( CE_Failure, CPLE_AppDefined, "Error parsing dstnodata arg #%d", i ); bDstNoDataNone = TRUE; continue; } CPLStringToComplex( papszTokens[i], psWO->padfDstNoDataReal + i, psWO->padfDstNoDataImag + i ); bDstNoDataNone = FALSE; CPLDebug( "WARP", "dstnodata of band %d set to %f", i, psWO->padfDstNoDataReal[i] ); } else { if ( ! bDstNoDataNone ) { psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i-1]; psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i-1]; CPLDebug( "WARP", "dstnodata of band %d set from previous band", i ); } else { CPLDebug( "WARP", "dstnodata value of band %d not set", i ); continue; } } GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 ); int bClamped = FALSE, bRounded = FALSE; psWO->padfDstNoDataReal[i] = GDALAdjustValueToDataType( GDALGetRasterDataType(hBand), psWO->padfDstNoDataReal[i], &bClamped, &bRounded ); if (bClamped) { CPLError(CE_Warning, CPLE_AppDefined, "for band %d, destination nodata value has been clamped " "to %.0f, the original value being out of range.", i + 1, psWO->padfDstNoDataReal[i]); } else if(bRounded) { CPLError(CE_Warning, CPLE_AppDefined, "for band %d, destination nodata value has been rounded " "to %.0f, %s being an integer datatype.", i + 1, psWO->padfDstNoDataReal[i], GDALGetDataTypeName(GDALGetRasterDataType(hBand))); } if( psOptions->bCreateOutput && iSrc == 0 ) { GDALSetRasterNoDataValue( GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ), psWO->padfDstNoDataReal[i] ); } } CSLDestroy( papszTokens ); } /* 检查输出数据集中是否已经存在无数据值 */ if ( psOptions->pszDstNodata == NULL ) { int bHaveNodataAll = TRUE; for( int i = 0; i < psWO->nBandCount; i++ ) { GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 ); int bHaveNodata = FALSE; GDALGetRasterNoDataValue( hBand, &bHaveNodata ); bHaveNodataAll &= bHaveNodata; } if( bHaveNodataAll ) { psWO->padfDstNoDataReal = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); psWO->padfDstNoDataImag = (double *) CPLCalloc(psWO->nBandCount, sizeof(double)); for( int i = 0; i < psWO->nBandCount; i++ ) { GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 ); int bHaveNodata = FALSE; psWO->padfDstNoDataReal[i] = GDALGetRasterNoDataValue( hBand, &bHaveNodata ); CPLDebug("WARP", "band=%d dstNoData=%f", i, psWO->padfDstNoDataReal[i] ); } } } /* 否则使用资源图像中波段的无数据值 */ if ( psOptions->pszDstNodata == NULL && psWO->padfSrcNoDataReal != NULL && psWO->padfDstNoDataReal == NULL ) { psWO->padfDstNoDataReal = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); psWO->padfDstNoDataImag = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); if( !psOptions->bQuiet ) printf( "Copying nodata values from source %s to destination %s.\n", GDALGetDescription(hSrcDS), pszDest ); for( int i = 0; i < psWO->nBandCount; i++ ) { psWO->padfDstNoDataReal[i] = psWO->padfSrcNoDataReal[i]; psWO->padfDstNoDataImag[i] = psWO->padfSrcNoDataImag[i]; CPLDebug("WARP", "srcNoData=%f dstNoData=%f", psWO->padfSrcNoDataReal[i], psWO->padfDstNoDataReal[i] ); if( psOptions->bCreateOutput && iSrc == 0 ) { CPLDebug("WARP", "calling GDALSetRasterNoDataValue() for band#%d", i ); GDALSetRasterNoDataValue( GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ), psWO->padfDstNoDataReal[i] ); } } if( psOptions->bCreateOutput && !bInitDestSetByUser && iSrc == 0 ) { /*由于开始时不一定指定了无数据值，需要检查*/ psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA" ); /*我们已经初始化INIT_DEST=0. 现在使用NO_DATA赋写它*/ } } /* -------------------------------------------------------------------- */ /* 如果存在数据分割线, 把它转换到资源图像的像素坐标系并且将其插入变换选项. /* -------------------------------------------------------------------- */ if( hCutline != NULL ) { CPLErr eError; eError = TransformCutlineToSource( hWrkSrcDS, hCutline, &(psWO->papszWarpOptions), psOptions->papszTO ); if(eError == CE_Failure) { if( hTransformArg != NULL ) GDALDestroyTransformer( hTransformArg ); GDALDestroyWarpOptions( psWO ); GDALWarpAppOptionsFree(psOptions); if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); if( bMustCloseDstDSInCaseOfError ) GDALClose(hDstDS); if( poSrcOvrDS ) delete poSrcOvrDS; return NULL; } } /* -------------------------------------------------------------------- */ /* 如果我们生成了输出数据集的虚拟栅格的输出, 接下来使用变换选项初始化它 /* 并且现在写给虚拟栅格而不是使用处理运行类. /* -------------------------------------------------------------------- */ if( bVRT ) { GDALSetMetadataItem(hDstDS, "SrcOvrLevel", CPLSPrintf("%d", psOptions->nOvLevel), NULL); if( GDALInitializeWarpedVRT( hDstDS, psWO ) != CE_None ) { if( hTransformArg != NULL ) GDALDestroyTransformer( hTransformArg ); GDALDestroyWarpOptions( psWO ); GDALWarpAppOptionsFree(psOptions); if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); if( bMustCloseDstDSInCaseOfError ) GDALClose(hDstDS); if( poSrcOvrDS ) delete poSrcOvrDS; return NULL; } // 在释放poSrcOvrDS和warping options前,我们需要关闭它. CPLString osDstFilename(GDALGetDescription(hDstDS)); char** papszContent = NULL; if( osDstFilename.size() == 0 ) papszContent = CSLDuplicate(GDALGetMetadata(hDstDS, "xml:VRT")); GDALClose(hDstDS); if( poSrcOvrDS ) delete poSrcOvrDS; GDALDestroyWarpOptions( psWO ); GDALWarpAppOptionsFree(psOptions); if( papszContent == NULL ) return GDALOpen(osDstFilename, GA_Update); else { GDALDatasetH hDSRet = GDALOpen(papszContent[0], GA_ReadOnly); CSLDestroy(papszContent); return hDSRet; } } /* -------------------------------------------------------------------- */ /* 初始化并且执行变换. /* -------------------------------------------------------------------- */ GDALWarpOperation oWO; if( oWO.Initialize( psWO ) == CE_None ) { CPLErr eErr; if( psOptions->bMulti ) eErr = oWO.ChunkAndWarpMulti( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) );//如果选项中开启了多进程处理，就多线程处理. else eErr = oWO.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) );//常规处理 if (eErr != CE_None) bHasGotErr = TRUE; } /* -------------------------------------------------------------------- */ /* 清空变换参数和选项. /* -------------------------------------------------------------------- */ if( hTransformArg != NULL ) GDALDestroyTransformer( hTransformArg ); GDALDestroyWarpOptions( psWO ); if( poSrcOvrDS ) delete poSrcOvrDS; } /* -------------------------------------------------------------------- */ /* 最后清空. /* -------------------------------------------------------------------- */ CPLErr eErrBefore = CPLGetLastErrorType(); GDALFlushCache( hDstDS ); if (eErrBefore == CE_None && CPLGetLastErrorType() != CE_None) { bHasGotErr = TRUE; } if( hCutline != NULL ) OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); GDALWarpAppOptionsFree(psOptions); if( bHasGotErr && bMustCloseDstDSInCaseOfError ) GDALClose(hDstDS); return (bHasGotErr) ? NULL : hDstDS; }

分析这段代码后得到的流程图如下：该图中的几个函数都是生成的变换参数的。整个函数所做的事情基本就是读入源图像信息并根据输入项构建变换项，计算变换参数，初始化输出数据集，图像变换填充，释放内存。而这个变换的核心，是变换运行类——GDALWarpOperation、变换参数hTransformArg以及变换函数指针pfnTransformer。关于这些我摘录了一部分GDAL帮助内的说明并翻译如下：

专利

最新回复(0)

GDAL工具箱详解之gdalwarp.exe

目录

使用方法

源码分析

应用程序主函数体：

变换库核心函数

相关内容研究

GDALWarpOperation类

GDALWarpKernel类

专利