准备资料
1. 一张wgs84投影的大tiff文件，建议初学者使用一张全球 2048 * 1024 / 4096 * 2048 的完整数据(有助于观察验证)
2. 准备C++开发环境，配置好gdal （笔者使用的环境是 vs2022 + gdal-2.3.0) c++ 开发环境
3. 建立一个测试工程（笔者建立的控制台)
4. 实现wgs84投影类/初期可以指建立wgs84,后续需要用到Mercator，则一次递进
5. 关键算法实现：如何从一个大图像中获取指定的部分
6. GDALAllRegister() 函数实现对GDAL的初始化
GDALAllRegister(); char szCWD[FE_PATH_LENGHT] = {0}; /// 获取当前exe的路径 std::string path = getcwd(szCWD,sizeof(szCWD)); std::string gdalData = path + "/gdal_data"; /// 设置gdal_data路径，否则gdal无法正确工作 CPLSetConfigOption( "GDAL_DATA", gdalData.c_str() );
其中 CPLSetConfigOption( "GDAL_DATA", gdalData.c_str() ); 告诉GDAL初始化数据从什么位置读取，很重要，否则会影响后续功能
7. 关键算法实现
1 bool createTile(double geoTrans[6],TileDataSet& result) 2 { 3 // The transformation option list 4 CPLStringList transformOptions; 5 6 // The source, sink and grid srs 7 const char* pszSrcWKT = GDALGetProjectionRef(_srcDS); 8 const char* pszGridWKT = pszSrcWKT; 9 if (pszSrcWKT == nullptr || !strlen(pszSrcWKT)) 10 return false; 11 12 // Populate the SRS WKT strings if we need to reproject 13 if (requiresReprojection()) 14 { 15 transformOptions.SetNameValue("SRC_SRS", pszSrcWKT); 16 transformOptions.SetNameValue("DST_SRS", pszGridWKT); 17 } 18 19 GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); 20 psWarpOptions->eResampleAlg = GRA_Average; 21 psWarpOptions->dfWarpMemoryLimit= 0; 22 psWarpOptions->hSrcDS = _srcDS; 23 psWarpOptions->nBandCount = _srcDS->GetRasterCount(); 24 psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); 25 psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); 26 27 for (int i = 0; i < psWarpOptions->nBandCount; ++i) 28 { 29 psWarpOptions->panDstBands[i] = i + 1; 30 psWarpOptions->panSrcBands[i] = i + 1; 31 } 32 33 void* transformerArg = GDALCreateGenImgProjTransformer2(_srcDS, NULL, transformOptions.List()); 34 if(transformerArg == NULL) 35 { 36 GDALDestroyWarpOptions(psWarpOptions); 37 return false; 38 } 39 40 auto hWrkSrcDS = psWarpOptions->hSrcDS = _srcDS; 41 GDALSetGenImgProjTransformerDstGeoTransform(transformerArg, geoTrans ); 42 43 if (false) 44 { 45 // approximate: wrap the transformer with a linear approximator 46 psWarpOptions->pTransformerArg = GDALCreateApproxTransformer(GDALGenImgProjTransform, transformerArg, 0.125000000); 47 48 if (psWarpOptions->pTransformerArg == nullptr) 49 { 50 GDALDestroyWarpOptions(psWarpOptions); 51 GDALDestroyGenImgProjTransformer(transformerArg); 52 return false; 53 } 54 psWarpOptions->pfnTransformer = GDALApproxTransform; 55 } 56 else 57 { 58 psWarpOptions->pTransformerArg = transformerArg; 59 psWarpOptions->pfnTransformer = GDALGenImgProjTransform; 60 } 61 62 // Specify a multi threaded warp operation using all CPU cores 63 CPLStringList warpOptions(psWarpOptions->papszWarpOptions, false); 64 warpOptions.SetNameValue("NUM_THREADS", "ALL_CPUS"); 65 psWarpOptions->papszWarpOptions = warpOptions.StealList(); 66 67 // The raster tile is represented as a VRT dataset 68 auto dstDS = GDALCreateWarpedVRT(hWrkSrcDS, _tileSize, _tileSize, geoTrans, psWarpOptions); 69 70 if (dstDS == nullptr) 71 { 72 GDALDestroyGenImgProjTransformer(transformerArg); 73 return false; 74 } 75 76 /// Set the projection information on the dataset. This will always be the grid 77 /// SRS. 78 if (GDALSetProjection( dstDS, pszGridWKT ) != CE_None) 79 { 80 GDALClose(dstDS); 81 if (transformerArg != nullptr) 82 { 83 GDALDestroyGenImgProjTransformer(transformerArg); 84 } 85 return false; 86 } 87 88 result._dataSet = dstDS; 89 result._transArg = transformerArg; 90 return true; 91 }
其中重点 ： GDALCreateWarpedVRT函数
8. 框架代码实现如下:
1 int main(int argc,char** args) 2 { 3 system("color 2"); 4 GDALAllRegister(); 5 char szCWD[FE_PATH_LENGHT] = {0}; 6 /// 获取当前exe的路径 7 std::string path = getcwd(szCWD,sizeof(szCWD)); 8 std::string gdalData = path + "/gdal_data"; 9 /// 设置gdal_data路径，否则gdal无法正确工作 10 CPLSetConfigOption( "GDAL_DATA", gdalData.c_str() ); 11 12 /// conv("d:/tms",int3(4,1,2)); 13 /// return 0; 14 std::cout << std::endl; 15 /// 版本号说明 1.0.0.0 16 /// 主版本号.副版本号.增加的功能个数.修复的bug数量 17 std::cout << "版权所有:成都经纬快图科技有限公司" << std::endl; 18 std::cout << "咨询微信:JingWeiKuaiTu" << std::endl; 19 std::cout << "当前版本:V2.0.0.0" << std::endl; 20 21 22 cmdline::parser a; 23 a.add<string>("in", 'i', "src file", true, ""); 24 a.add<string>("out", 'o', "dst dir", false, ""); 25 a.add<string>("lev", 'l', "levels", false, ""); 26 a.add<string>("ext", 'e', "jpg,png", false, "jpg"); 27 a.add<bool>("mercator", 'm', "mercator", false, false); 28 a.add("help", 0, "-i d:/src.tif -o d:/dst -l 1-4 -e jpg -m 0\n"); 29 30 a.parse_check(argc, args); 31 32 Levels levels; 33 34 std::string srcFile = a.get<string>("in"); 35 std::string outDir = a.get<string>("out"); 36 std::string levs = a.get<string>("lev"); 37 std::string ext = a.get<string>("ext"); 38 bool isMercator = a.get<bool>("mercator"); 39 if (!levs.empty()) 40 { 41 levels = toLevels(levs.c_str()); 42 } 43 if (outDir.empty()) 44 { 45 outDir = getFilePath(srcFile) + "/tile"; 46 } 47 createDir(outDir.c_str()); 48 49 Tiff tifObject(256,isMercator); 50 if (!tifObject.open(srcFile.c_str())) 51 { 52 return 0; 53 } 54 FETimestamp tms; 55 tifObject.exportTile(outDir,levels,ext.c_str()); 56 57 printf("\nall cost %lf second",tms.getElapsedSecond()); 58 /// 获取tif的基本信息 59 60 return 0; 61 }
代码已经上传