如何使用GDAL进行图像镶嵌 |
您所在的位置：网站首页 › gdal影像拼接 › 如何使用GDAL进行图像镶嵌 |
|
相信看过和用过gdalwarp工具的同学都会对gdalwarp这个工具的强悍功能所震惊。今天主要就是用gdalwarp来进行图像镶嵌，当然这个镶嵌是比较简单的镶嵌，就是上层的图像会覆盖下层的图像。对于gdalwarp的介绍，我还要搬过来一下，见下面： The gdalwarp utility is an imagemosaicing, reprojection and warping utility. The program can reproject to any supported projection, and can also apply GCPs stored with the image if the image is "raw" with control information. 总觉得自己翻译的不好，就不翻译了，关于这篇博客的主要内容，就是上面的红色字体，mosaicing。gdalwarp使用的镶嵌算法很简单，就是根据图像的地理位置进行拼接，对于重叠去的处理，就是后面处理的图像会覆盖前面处理的图像，下面的代码都是改造gdalwarp.cpp中的代码。这个镶嵌可以支持镶嵌块文件镶嵌，是一个面状的shp文件，通过shp的属性表字段和对应的影像进行关联。废话不多说，上代码（代码比较粗糙，板砖轻拍：））。 先是几个子函数，第一个获取文件名称的，这个代码是很久之前写的，用的比较笨的办法，其实可以用boost库的filesysterm中的path类，一个函数就搞定了。反正都一样，呵呵。 /** * @brief 获取文件名称，除去后缀名 */ inline string GetFileName(const char* pszFile) { string temp = pszFile; size_t a = temp.find_last_of('\\'); size_t aa = temp.find_last_of('/'); if (a != string::npos && aa != string::npos) a = MAX(a,aa); else a = MIN(a,aa); size_t b = temp.find_last_of('.'); string strLayerName = temp.substr(a+1); return strLayerName; } 接下来是创建输出图像的代码，这个函数比较长，简单的思想就是，遍历所有的输入文件，然后计算所有的输入文件的空间范围的并集，然后根据图像的分辨率，计算结果图像的大小投影等信息，这个函数里面其实可以直接进行投影转换，等一系列的坐标转换处理，这里不进行深入讨论，以免偏离本文重点，这个代码是从gdalwarp.cpp中搬过来的，中间我加了一些注释，修改了一些返回值之类的，好，看代码： /** * @brief 创建输出文件 * @param papszSrcFiles 输入文件列表 * @param pszFilename 输出文件路径 * @param pszFormat 输出文件格式 * @param papszTO 转换选项 * @param ppapszCreateOptions 创建文件选项 * @param eDT 创建文件数据类型 * @return 返回输出文件的句柄 */ static GDALDatasetH GDALWarpCreateOutput( char **papszSrcFiles, const char *pszFilename, const char *pszFormat, char **papszTO, char ***ppapszCreateOptions, GDALDataType eDT, int &iRst ) { GDALDriverH hDriver; GDALDatasetH hDstDS; void *hTransformArg; GDALColorTableH hCT = NULL; double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0; double dfWrkResX=0, dfWrkResY=0; int nDstBandCount = 0; iRst = RE_SUCCESS; /* -------------------------------------------------------------------- */ /* 获取输出文件类型驱动 */ /* -------------------------------------------------------------------- */ hDriver = GDALGetDriverByName( pszFormat ); if( hDriver == NULL || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL ) { iRst = RE_FILENOTSUPPORT; return NULL; } /* -------------------------------------------------------------------- */ /* 循环所有的输入文件，计算图像范围 */ /* -------------------------------------------------------------------- */ char *pszThisTargetSRS = (char*)CSLFetchNameValue( papszTO, "DST_SRS" ); if( pszThisTargetSRS != NULL ) pszThisTargetSRS = CPLStrdup( pszThisTargetSRS ); for(int iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ ) { GDALDatasetH hSrcDS; const char *pszThisSourceSRS = CSLFetchNameValue(papszTO,"SRC_SRS"); hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly ); if( hSrcDS == NULL ) { iRst = RE_FILENOTEXIST; return NULL; } /* -------------------------------------------------------------------- */ /* 检查当前文件是否存在波段 */ /* -------------------------------------------------------------------- */ if ( GDALGetRasterCount(hSrcDS) == 0 ) { iRst = RE_FILETYPEERROR; return NULL; } if( eDT == GDT_Unknown ) eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)); /* -------------------------------------------------------------------- */ /* If we are processing the first file, and it has a color */ /* table, then we will copy it to the destination file. */ /* -------------------------------------------------------------------- */ if( iSrc == 0 ) { nDstBandCount = GDALGetRasterCount(hSrcDS); hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) ); if( hCT != NULL ) hCT = GDALCloneColorTable( hCT ); } /* -------------------------------------------------------------------- */ /* Get the sourcesrs from the dataset, if not set already. */ /* -------------------------------------------------------------------- */ if( pszThisSourceSRS == NULL ) { const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" ); if( GDALGetProjectionRef( hSrcDS ) != NULL && strlen(GDALGetProjectionRef( hSrcDS )) > 0 && (pszMethod == NULL || EQUAL(pszMethod,"GEOTRANSFORM")) ) pszThisSourceSRS = GDALGetProjectionRef( hSrcDS ); else if( GDALGetGCPProjection( hSrcDS ) != NULL && strlen(GDALGetGCPProjection(hSrcDS)) > 0 && GDALGetGCPCount( hSrcDS ) > 1 && (pszMethod == NULL || EQUALN(pszMethod,"GCP_",4)) ) pszThisSourceSRS = GDALGetGCPProjection( hSrcDS ); else if( pszMethod != NULL && EQUAL(pszMethod,"RPC") ) pszThisSourceSRS = SRS_WKT_WGS84; else pszThisSourceSRS = ""; } if( pszThisTargetSRS == NULL ) pszThisTargetSRS = CPLStrdup( pszThisSourceSRS ); /* -------------------------------------------------------------------- */ /* Create a transformation object from the source to */ /* destination coordinate system. */ /* -------------------------------------------------------------------- */ hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO ); if( hTransformArg == NULL ) { CPLFree( pszThisTargetSRS ); GDALClose( hSrcDS ); return NULL; } /* -------------------------------------------------------------------- */ /* Get approximate output definition. */ /* -------------------------------------------------------------------- */ double adfThisGeoTransform[6]; double adfExtent[4]; int nThisPixels, nThisLines; if( GDALSuggestedWarpOutput2( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfThisGeoTransform, &nThisPixels, &nThisLines, adfExtent, 0 ) != CE_None ) { CPLFree( pszThisTargetSRS ); GDALClose( hSrcDS ); return NULL; } if (CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", NULL ) == NULL) { double MinX = adfExtent[0]; double MaxX = adfExtent[2]; double MaxY = adfExtent[3]; double MinY = adfExtent[1]; int bSuccess = TRUE; /* Check that the the edges of the target image are in the validity area */ /* of the target projection */ #define N_STEPS 20 int i,j; for(i=0;i (MaxY - MinY) / nThisLines) bSuccess = FALSE; } } /* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces ogrct.cpp */ /* to check the consistency of each requested projection result with the */ /* invert projection */ if (!bSuccess) { CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" ); CPLDebug("WARP", "Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE"); GDALDestroyGenImgProjTransformer(hTransformArg); hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO ); if( GDALSuggestedWarpOutput2( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfThisGeoTransform, &nThisPixels, &nThisLines, adfExtent, 0 ) != CE_None ) { CPLFree( pszThisTargetSRS ); GDALClose( hSrcDS ); return NULL; } } } /* -------------------------------------------------------------------- */ /* Expand the working bounds to include this region, ensure the */ /* working resolution is no more than this resolution. */ /* -------------------------------------------------------------------- */ if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 ) { dfWrkMinX = adfExtent[0]; dfWrkMaxX = adfExtent[2]; dfWrkMaxY = adfExtent[3]; dfWrkMinY = adfExtent[1]; dfWrkResX = adfThisGeoTransform[1]; dfWrkResY = ABS(adfThisGeoTransform[5]); } else { dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]); dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]); dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]); dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]); dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]); dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5])); } GDALDestroyGenImgProjTransformer( hTransformArg ); GDALClose( hSrcDS ); } /* -------------------------------------------------------------------- */ /* Did we have any usable sources? */ /* -------------------------------------------------------------------- */ if( nDstBandCount == 0 ) { CPLError( CE_Failure, CPLE_AppDefined, "No usable source images." ); CPLFree( pszThisTargetSRS ); return NULL; } /* -------------------------------------------------------------------- */ /* Turn the suggested region into a geotransform and suggested */ /* number of pixels and lines. */ /* -------------------------------------------------------------------- */ double adfDstGeoTransform[6]; int nPixels, nLines; adfDstGeoTransform[0] = dfWrkMinX; adfDstGeoTransform[1] = dfWrkResX; adfDstGeoTransform[2] = 0.0; adfDstGeoTransform[3] = dfWrkMaxY; adfDstGeoTransform[4] = 0.0; adfDstGeoTransform[5] = -1 * dfWrkResY; nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5); nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5); /* -------------------------------------------------------------------- */ /* Did the user override some parameters? */ /* -------------------------------------------------------------------- */ if( dfXRes != 0.0 && dfYRes != 0.0 ) { if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 ) { dfMinX = adfDstGeoTransform[0]; dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels; dfMaxY = adfDstGeoTransform[3]; dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines; } nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes); nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes); adfDstGeoTransform[0] = dfMinX; adfDstGeoTransform[3] = dfMaxY; adfDstGeoTransform[1] = dfXRes; adfDstGeoTransform[5] = -dfYRes; } else if( nForcePixels != 0 && nForceLines != 0 ) { if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 ) { dfMinX = dfWrkMinX; dfMaxX = dfWrkMaxX; dfMaxY = dfWrkMaxY; dfMinY = dfWrkMinY; } dfXRes = (dfMaxX - dfMinX) / nForcePixels; dfYRes = (dfMaxY - dfMinY) / nForceLines; adfDstGeoTransform[0] = dfMinX; adfDstGeoTransform[3] = dfMaxY; adfDstGeoTransform[1] = dfXRes; adfDstGeoTransform[5] = -dfYRes; nPixels = nForcePixels; nLines = nForceLines; } else if( nForcePixels != 0 ) { if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 ) { dfMinX = dfWrkMinX; dfMaxX = dfWrkMaxX; dfMaxY = dfWrkMaxY; dfMinY = dfWrkMinY; } dfXRes = (dfMaxX - dfMinX) / nForcePixels; dfYRes = dfXRes; adfDstGeoTransform[0] = dfMinX; adfDstGeoTransform[3] = dfMaxY; adfDstGeoTransform[1] = dfXRes; adfDstGeoTransform[5] = -dfYRes; nPixels = nForcePixels; nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes); } else if( nForceLines != 0 ) { if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 ) { dfMinX = dfWrkMinX; dfMaxX = dfWrkMaxX; dfMaxY = dfWrkMaxY; dfMinY = dfWrkMinY; } dfYRes = (dfMaxY - dfMinY) / nForceLines; dfXRes = dfYRes; adfDstGeoTransform[0] = dfMinX; adfDstGeoTransform[3] = dfMaxY; adfDstGeoTransform[1] = dfXRes; adfDstGeoTransform[5] = -dfYRes; nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes); nLines = nForceLines; } else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 ) { dfXRes = adfDstGeoTransform[1]; dfYRes = fabs(adfDstGeoTransform[5]); nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes); nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes); dfXRes = (dfMaxX - dfMinX) / nPixels; dfYRes = (dfMaxY - dfMinY) / nLines; adfDstGeoTransform[0] = dfMinX; adfDstGeoTransform[3] = dfMaxY; adfDstGeoTransform[1] = dfXRes; adfDstGeoTransform[5] = -dfYRes; } /* -------------------------------------------------------------------- */ /* Do we want to generate an alpha band in the output file? */ /* -------------------------------------------------------------------- */ if( bEnableSrcAlpha ) nDstBandCount--; if( bEnableDstAlpha ) nDstBandCount++; /* -------------------------------------------------------------------- */ /* Create the output file. */ /* -------------------------------------------------------------------- */ hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines, nDstBandCount, eDT, *ppapszCreateOptions ); if( hDstDS == NULL ) { CPLFree( pszThisTargetSRS ); return NULL; } /* -------------------------------------------------------------------- */ /* Write out the projection definition. */ /* -------------------------------------------------------------------- */ GDALSetProjection( hDstDS, pszThisTargetSRS ); GDALSetGeoTransform( hDstDS, adfDstGeoTransform ); /* -------------------------------------------------------------------- */ /* Try to set color interpretation of output file alpha band. */ /* TODO: We should likely try to copy the other bands too. */ /* -------------------------------------------------------------------- */ if( bEnableDstAlpha ) { GDALSetRasterColorInterpretation( GDALGetRasterBand( hDstDS, nDstBandCount ), GCI_AlphaBand ); } /* -------------------------------------------------------------------- */ /* Copy the color table, if required. */ /* -------------------------------------------------------------------- */ if( hCT != NULL ) { GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT ); GDALDestroyColorTable( hCT ); } CPLFree( pszThisTargetSRS ); return hDstDS; } 接下来是投影坐标到图像行列号坐标转换的一个类，也是从gdalwarp.cpp中搬过来的： /** * @brief 转换投影坐标到图像行列号坐标 * * 用于图像裁切中，按照AOI来裁切时转换AOI坐标到图像行列号坐标 */ class CutlineTransformer : public OGRCoordinateTransformation { public: /** * @brief 源图像转换参数 */ void *hSrcImageTransformer; /** * @brief 获取源数据空间参考坐标系 * @return 返回源数据空间参考坐标系 */ virtual OGRSpatialReference *GetSourceCS() { return NULL; } /** * @brief 获取目标数据空间参考坐标系 * @return 返回目标数据空间参考坐标系 */ virtual OGRSpatialReference *GetTargetCS() { return NULL; } /** * @brief 批量转换投影坐标到行列号 * @param nCount 坐标点个数 * @param x x坐标串 * @param y y坐标串 * @param z z坐标串 * @return 返回是否执行正确 */ virtual int Transform( int nCount, double *x, double *y, double *z = NULL ) { int nResult; int *pabSuccess = (int *) CPLCalloc(sizeof(int),nCount); nResult = TransformEx( nCount, x, y, z, pabSuccess ); CPLFree( pabSuccess ); return nResult; } /** * @brief 批量转换投影坐标到行列号 * @param nCount 坐标点个数 * @param x x坐标串 * @param y y坐标串 * @param z z坐标串 * @param pabSuccess 保存执行成功的数组 * @return 返回是否执行正确 */ virtual int TransformEx( int nCount, double *x, double *y, double *z = NULL, int *pabSuccess = NULL ) { return GDALGenImgProjTransform( hSrcImageTransformer, TRUE, nCount, x, y, z, pabSuccess ); } }; 接下来这个函数比较重要，主要作用是读取矢量文件中的镶嵌块文件，构造一个镶嵌线，在镶嵌的时候使用： /** * @brief 加载镶嵌块 * @param pszCutlineDSName 镶嵌块文件 * @param pszCLayer 镶嵌块图层名称，可以为NULL * @param pszCWHERE 镶嵌块过滤字段 * @param pszCSQL 镶嵌块SQL过滤字段 * @param phCutlineRet 返回的镶嵌块数据指针 * @return 返回值0表示该正确 */ static int LoadCutline( const char *pszCutlineDSName, const char *pszCLayer, const char *pszCWHERE, const char *pszCSQL, void **phCutlineRet ) { OGRRegisterAll(); /* -------------------------------------------------------------------- */ /* Open source vector dataset. */ /* -------------------------------------------------------------------- */ OGRDataSourceH hSrcDS; hSrcDS = OGROpen( pszCutlineDSName, FALSE, NULL ); if( hSrcDS == NULL ) { if(m_pProcess != NULL) m_pProcess->SetMessage("指定的镶嵌块文件不存在..."); return RE_FILENOTEXIST; } /* -------------------------------------------------------------------- */ /* Get the source layer */ /* -------------------------------------------------------------------- */ OGRLayerH hLayer = NULL; if( pszCSQL != NULL ) hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszCSQL, NULL, NULL ); else if( pszCLayer != NULL ) hLayer = OGR_DS_GetLayerByName( hSrcDS, pszCLayer ); else hLayer = OGR_DS_GetLayer( hSrcDS, 0 ); if( hLayer == NULL ) { if(m_pProcess != NULL) m_pProcess->SetMessage("指定的矢量文件中没有找到图层..."); OGR_DS_Destroy( hSrcDS ); return RE_FAILED; } /* -------------------------------------------------------------------- */ /* Apply WHERE clause if there is one. */ /* -------------------------------------------------------------------- */ if( pszCWHERE != NULL ) OGR_L_SetAttributeFilter( hLayer, pszCWHERE ); /* -------------------------------------------------------------------- */ /* Collect the geometries from this layer, and build list of */ /* burn values. */ /* -------------------------------------------------------------------- */ OGRFeatureH hFeat; OGRGeometryH hMultiPolygon = OGR_G_CreateGeometry( wkbMultiPolygon ); OGR_L_ResetReading( hLayer ); while( (hFeat = OGR_L_GetNextFeature( hLayer )) != NULL ) { OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat); if( hGeom == NULL ) { if(m_pProcess != NULL) m_pProcess->SetMessage("指定的矢量文件图层中不存在图形..."); OGR_DS_Destroy( hSrcDS ); return RE_FAILED; } OGRwkbGeometryType eType = wkbFlatten(OGR_G_GetGeometryType( hGeom )); if( eType == wkbPolygon ) OGR_G_AddGeometry( hMultiPolygon, hGeom ); else if( eType == wkbMultiPolygon ) { int iGeom; for( iGeom = 0; iGeom < OGR_G_GetGeometryCount( hGeom ); iGeom++ ) { OGR_G_AddGeometry( hMultiPolygon, OGR_G_GetGeometryRef(hGeom,iGeom) ); } } else { if(m_pProcess != NULL) m_pProcess->SetMessage("指定的镶嵌块文件不是多边形矢量..."); OGR_DS_Destroy( hSrcDS ); return RE_FAILED; } OGR_F_Destroy( hFeat ); } if( OGR_G_GetGeometryCount( hMultiPolygon ) == 0 ) { if(m_pProcess != NULL) m_pProcess->SetMessage("没有从指定的镶嵌块文件图层中读取镶嵌块..."); return RE_FAILED; } /* -------------------------------------------------------------------- */ /* Ensure the coordinate system gets set on the geometry. */ /* -------------------------------------------------------------------- */ OGR_G_AssignSpatialReference( hMultiPolygon, OGR_L_GetSpatialRef(hLayer) ); *phCutlineRet = (void *) hMultiPolygon; /* -------------------------------------------------------------------- */ /* Cleanup */ /* -------------------------------------------------------------------- */ if( pszCSQL != NULL ) OGR_DS_ReleaseResultSet( hSrcDS, hLayer ); OGR_DS_Destroy( hSrcDS ); return RE_SUCCESS; } 下面的函数是用来将镶嵌块的投影坐标转换到图像的行列号坐标中去： /** * @brief 转换AOI到源文件之间的行列号 * @param hSrcDS 输入文件GDAL数据集句柄 * @param hCutline 裁切的几何形状 * @param ppapszWarpOptions 转换选项，用于配置裁切参数 * @param papszTO_In 转换选项 * @return 返回值0表示该正确 */ static void TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline, char ***ppapszWarpOptions, char **papszTO_In ) { OGRGeometryH hMultiPolygon = OGR_G_Clone( (OGRGeometryH) hCutline ); char **papszTO = CSLDuplicate( papszTO_In ); /* -------------------------------------------------------------------- */ /* Checkout that SRS are the same. */ /* -------------------------------------------------------------------- */ OGRSpatialReferenceH hRasterSRS = NULL; const char *pszProjection = NULL; if( GDALGetProjectionRef( hSrcDS ) != NULL && strlen(GDALGetProjectionRef( hSrcDS )) > 0 ) pszProjection = GDALGetProjectionRef( hSrcDS ); else if( GDALGetGCPProjection( hSrcDS ) != NULL ) pszProjection = GDALGetGCPProjection( hSrcDS ); if( pszProjection != NULL ) { hRasterSRS = OSRNewSpatialReference(NULL); if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None ) { OSRDestroySpatialReference(hRasterSRS); hRasterSRS = NULL; } } OGRSpatialReferenceH hSrcSRS = OGR_G_GetSpatialReference( hMultiPolygon ); if( hRasterSRS != NULL && hSrcSRS != NULL ) { /* ok, we will reproject */ } else if( hRasterSRS != NULL && hSrcSRS == NULL ) { fprintf(stderr, "Warning : the source raster dataset has a SRS, but the input vector layer\n" "not. Cutline results may be incorrect.\n"); } else if( hRasterSRS == NULL && hSrcSRS != NULL ) { fprintf(stderr, "Warning : the input vector layer has a SRS, but the source raster dataset does not.\n" "Cutline results may be incorrect.\n"); } if( hRasterSRS != NULL ) OSRDestroySpatialReference(hRasterSRS); /* -------------------------------------------------------------------- */ /* Extract the cutline SRS WKT. */ /* -------------------------------------------------------------------- */ if( hSrcSRS != NULL ) { char *pszCutlineSRS_WKT = NULL; OSRExportToWkt( hSrcSRS, &pszCutlineSRS_WKT ); papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszCutlineSRS_WKT ); CPLFree( pszCutlineSRS_WKT ); } else { int iDstSRS = CSLFindString( papszTO, "DST_SRS" ); if( iDstSRS >= 0 ) papszTO = CSLRemoveStrings( papszTO, iDstSRS, 1, NULL ); } /* -------------------------------------------------------------------- */ /* Transform the geometry to pixel/line coordinates. */ /* -------------------------------------------------------------------- */ CutlineTransformer oTransformer; oTransformer.hSrcImageTransformer = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO ); CSLDestroy( papszTO ); if( oTransformer.hSrcImageTransformer == NULL ) exit( 1 ); OGR_G_Transform( hMultiPolygon, (OGRCoordinateTransformationH) &oTransformer ); GDALDestroyGenImgProjTransformer( oTransformer.hSrcImageTransformer ); /* -------------------------------------------------------------------- */ /* Convert aggregate geometry into WKT. */ /* -------------------------------------------------------------------- */ char *pszWKT = NULL; OGR_G_ExportToWkt( hMultiPolygon, &pszWKT ); OGR_G_DestroyGeometry( hMultiPolygon ); *ppapszWarpOptions = CSLSetNameValue( *ppapszWarpOptions, "CUTLINE", pszWKT ); CPLFree( pszWKT ); } 接下来就是最重要的函数了，镶嵌函数的主体函数，注意两点，函数里面有关于进度条的信息，是我自己定义的，你可以参考我之前的博客，如果不想要，删掉即可；此外返回值是我自己定义的int类型，用来标记函数为什么返回了。下面的代码也是改造gdalwarp.cpp而来。 /** * @brief 图像镶嵌处理 * 注意事项，调用此函数之前，请先对影像进行几何纠正到同一投影坐标系统下，分辨率可以不同， * 但是投影信息以及输入的各个数据的波段个数必须一致，否则会出现不能正常完成镶嵌操作。 * 图像的分辨率会按照输入影像的第一个影像获取，包括投影等信息 * @param vStrSrcFiles 输入文件以及镶嵌线路径数组，默认第一个为参考影像，第一个图像在最下层，后面的依次向上 * @param pszCutLineFile 输入镶嵌块文件文件路径数组 * @param pszOutFile 输出文件路径 * @param eResampleMethod 重采样方式 * @param pszFormat 输出文件格式，详细参考GDAL支持数据类型 * @param pProcess 进度条指针，这个是我自己定义的类型，你可以自己写一个，或者删除 * @return 返回值，表示计算过程中出现的各种错误信息，返回值是我自己定义的int类型 */ int ImageMosaicing(vector vStrSrcFiles, const char* pszCutLineFile, const char* pszOutFile, GDALResampleAlg eResampleMethod, const char *pszFormat, CBatProcessBase* pProcess ) { CProcessBase* pAllProcess = NULL; CProcessBase* pSubProcess = NULL; if(pProcess != NULL) { pAllProcess = pProcess->m_pAllProcess; pSubProcess = pProcess->m_pSubProcess; pProcess->SetMessage("正在执行镶嵌..."); pAllProcess->SetStepCount(vStrSrcFiles.size()); } //输入输出文件 GDALDatasetH hDstDS; char **papszSrcFiles = NULL; const char *pszDstFilename = NULL; int bCreateOutput = FALSE; //坐标转换对象 void *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL; char **papszWarpOptions = NULL; double dfErrorThreshold = 0.125; GDALTransformerFunc pfnTransformer = NULL; //输出文件选项 char **papszCreateOptions = NULL; GDALDataType eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown; GDALResampleAlg eResampleAlg = (GDALResampleAlg)eResampleMethod; //NODATA设置 const char *pszSrcNodata = NULL; const char *pszDstNodata = NULL; int bMulti = FALSE; char **papszTO = NULL; //镶嵌线文件 const char *pszCutlineDSName = NULL; if(vStrSrcFiles.empty()) return RE_PARAMERROR; GDALAllRegister(); OGRRegisterAll(); pszSrcNodata = "0 0 0"; pszDstNodata = "0 0 0"; dfXRes = 0.0; dfYRes = 0.0; bCreateOutput = TRUE; eResampleAlg = GRA_Bilinear; pszCutlineDSName = pszCutLineFile; for(size_t i=0; ipapszWarpOptions = CSLDuplicate(papszWarpOptions); psWO->eWorkingDataType = eWorkingType; psWO->eResampleAlg = eResampleAlg; psWO->hSrcDS = hSrcDS; psWO->hDstDS = hDstDS; psWO->pfnTransformer = pfnTransformer; psWO->pTransformerArg = hTransformArg; psWO->pfnProgress = ALGTermProgress; psWO->pProgressArg = pSubProcess; psWO->dfWarpMemoryLimit = 52428800; //使用50M的内存 /* -------------------------------------------------------------------- */ /* 创建波段映射关系 */ /* -------------------------------------------------------------------- */ if( bEnableSrcAlpha ) psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1; else psWO->nBandCount = GDALGetRasterCount(hSrcDS); 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; } /* -------------------------------------------------------------------- */ /* 构建alpha波段 */ /* -------------------------------------------------------------------- */ if( bEnableSrcAlpha ) psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS); if( !bEnableDstAlpha && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1 && GDALGetRasterColorInterpretation( GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS))) == GCI_AlphaBand ) { printf( "Using band %d of destination image as alpha.\n", GDALGetRasterCount(hDstDS) ); bEnableDstAlpha = TRUE; } if( bEnableDstAlpha ) psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS); /* -------------------------------------------------------------------- */ /* 创建 NODATA 选项 */ /* -------------------------------------------------------------------- */ if( pszSrcNodata != NULL && !EQUALN(pszSrcNodata,"n",1) ) { char **papszTokens = CSLTokenizeString( 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" ); } /* -------------------------------------------------------------------- */ /* 没有指定NODATA值，从源文件中读取 */ /* -------------------------------------------------------------------- */ if( pszSrcNodata == NULL ) { int bHaveNodata = FALSE; double dfReal = 0.0; for(int i = 0; !bHaveNodata && i < psWO->nBandCount; i++ ) { GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 ); dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata ); } if( bHaveNodata ) { printf( "Using internal nodata values (eg. %g) for image %s.\n", dfReal, papszSrcFiles[iSrc] ); 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( hSrcDS, 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; } } } } /* -------------------------------------------------------------------- */ /* 设置输出文件的NODATA为最大值 */ /* -------------------------------------------------------------------- */ if( pszDstNodata != NULL ) { char **papszTokens = CSLTokenizeString( pszDstNodata ); int nTokenCount = CSLCount(papszTokens); psWO->padfDstNoDataReal = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); psWO->padfDstNoDataImag = (double *) CPLMalloc(psWO->nBandCount*sizeof(double)); for(int i = 0; i < psWO->nBandCount; i++ ) { if( i < nTokenCount ) { CPLStringToComplex( papszTokens[i], psWO->padfDstNoDataReal + i, psWO->padfDstNoDataImag + i ); } else { psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i-1]; psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i-1]; } GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 ); int bClamped = FALSE, bRounded = FALSE; #define CLAMP(val,type,minval,maxval) \ do { if (val < minval) { bClamped = TRUE; val = minval; } \ else if (val > maxval) { bClamped = TRUE; val = maxval; } \ else if (val != (type)val) { bRounded = TRUE; val = (type)(val + 0.5); } } \ while(0) switch(GDALGetRasterDataType(hBand)) { case GDT_Byte: CLAMP(psWO->padfDstNoDataReal[i], GByte, 0.0, 255.0); break; case GDT_Int16: CLAMP(psWO->padfDstNoDataReal[i], GInt16, -32768.0, 32767.0); break; case GDT_UInt16: CLAMP(psWO->padfDstNoDataReal[i], GUInt16, 0.0, 65535.0); break; case GDT_Int32: CLAMP(psWO->padfDstNoDataReal[i], GInt32, -2147483648.0, 2147483647.0); break; case GDT_UInt32: CLAMP(psWO->padfDstNoDataReal[i], GUInt32, 0.0, 4294967295.0); break; default: break; } if (bClamped) { printf( "for band %d, destination nodata value has been clamped " "to %.0f, the original value being out of range.\n", i + 1, psWO->padfDstNoDataReal[i]); } else if(bRounded) { printf("for band %d, destination nodata value has been rounded " "to %.0f, %s being an integer datatype.\n", i + 1, psWO->padfDstNoDataReal[i], GDALGetDataTypeName(GDALGetRasterDataType(hBand))); } if( bCreateOutput ) { GDALSetRasterNoDataValue( GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ), psWO->padfDstNoDataReal[i] ); } } CSLDestroy( papszTokens ); } /* -------------------------------------------------------------------- */ /* 读取镶嵌线 */ /* -------------------------------------------------------------------- */ void *hCutline = NULL; if( pszCutlineDSName != NULL ) { string strFileName = GetFileName(papszSrcFiles[iSrc]); string strWhere = "影像路径=\"" + strFileName + "\""; LoadCutline( pszCutlineDSName, NULL, strWhere.c_str(), NULL, &hCutline ); } if( hCutline != NULL ) { TransformCutlineToSource( hSrcDS, hCutline, &(psWO->papszWarpOptions), papszTO ); } /* -------------------------------------------------------------------- */ /* 初始化执行warp */ /* -------------------------------------------------------------------- */ GDALWarpOperation oWO; CPLErr CE = CE_None; if( oWO.Initialize( psWO ) == CE_None ) { if( bMulti ) CE = oWO.ChunkAndWarpMulti( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) ); else CE = oWO.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) ); } /* -------------------------------------------------------------------- */ /* 清理资源 */ /* -------------------------------------------------------------------- */ if( hApproxArg != NULL ) GDALDestroyApproxTransformer( hApproxArg ); if( hGenImgProjArg != NULL ) GDALDestroyGenImgProjTransformer( hGenImgProjArg ); GDALDestroyWarpOptions( psWO ); if( hCutline != NULL ) //释放镶嵌块资源 { OGR_G_DestroyGeometry( (OGRGeometryH) hCutline ); hCutline = NULL; } GDALClose( hSrcDS ); if (CE != CE_None) { GDALClose( hDstDS ); CSLDestroy( papszSrcFiles ); CSLDestroy( papszWarpOptions ); CSLDestroy( papszTO ); GDALDumpOpenDatasets( stderr ); OGRCleanupAll(); RasterDelete(pszOutFile); if(pProcess != NULL) { if (!pAllProcess->m_bIsContinue) { pProcess->SetMessage("取消计算！"); return RE_USERCANCEL; } else { pProcess->SetMessage("执行失败！"); return RE_FAILED; } } return RE_FAILED; } if(pAllProcess != NULL) pAllProcess->StepIt(); } /* -------------------------------------------------------------------- */ /* 善后工作 */ /* -------------------------------------------------------------------- */ CLEAN: GDALClose( hDstDS ); CSLDestroy( papszSrcFiles ); CSLDestroy( papszWarpOptions ); CSLDestroy( papszTO ); GDALDumpOpenDatasets( stderr ); OGRCleanupAll(); if(pProcess != NULL) pProcess->SetMessage("镶嵌完成..."); return RE_SUCCESS; } 这里读取镶嵌块文件的时候，镶嵌块文件保存了参与镶嵌的文件的镶嵌块，同时有一个属性字段是“影像路径”，其值是镶嵌块对应的影像的绝对路径，当然你可以不用给镶嵌块文件，也是可以进行处理的，只不过效果可能就……你懂得。有任何疑问，可以看gdalwarp.cpp中的代码，基本上全是那里面的代码。 PS：这里面需要定义几个全局变量，要不然可能有问题，全局变量定义如下： //全局变量 图像范围 static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0; //全局变量 图像分辨率 static double dfXRes=0.0, dfYRes=0.0; //全局变量 图像宽高 static int nForcePixels=0, nForceLines=0; //全局变量 图像是否具有Alpha通道 static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE; //全局变量 进度条指针 CProcessBase* m_pProcess = NULL;忽然，发现代码好多，唉，我也不想搞这么长的代码啊，没办法，您就凑合着看吧，欢迎大家指导。 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |