如何使用GDAL进行AOI裁剪
来源：互联网 发布：代步滑板哪种好 知乎 编辑：程序博客网 时间：2024/04/30 09:14
在工作中，会经常使用一个行政区的矢量图去裁剪一个遥感影像图，得到该行政区的影像图，一般的商业遥感软件都具有这个功能。今天就是用GDAL来实现这一个很实用的功能。首先用到的是GDAL中的gdalwarp，又是warp，呵呵，上一篇就是使用warp进行重采样的。
首先需要用到gdal源码目录里面的app文件夹下的gdalwarp.cpp文件中的几个函数，大概行数是1651行，直到文件结尾，代码如下：
/************************************************************************//* GeoTransform_Transformer() *//* *//* Convert points from georef coordinates to pixel/line based *//* on a geotransform. *//************************************************************************/class CutlineTransformer : public OGRCoordinateTransformation{public: void *hSrcImageTransformer; virtual OGRSpatialReference *GetSourceCS() { return NULL; } virtual OGRSpatialReference *GetTargetCS() { return NULL; } 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; } virtual int TransformEx( int nCount, double *x, double *y, double *z = NULL, int *pabSuccess = NULL ) { return GDALGenImgProjTransform( hSrcImageTransformer, TRUE, nCount, x, y, z, pabSuccess ); }};/************************************************************************//* LoadCutline() *//* *//* Load blend cutline from OGR datasource. *//************************************************************************/static voidLoadCutline( const char *pszCutlineDSName, const char *pszCLayer, const char *pszCWHERE, const char *pszCSQL, void **phCutlineRet ){#ifndef OGR_ENABLED CPLError( CE_Failure, CPLE_AppDefined, "Request to load a cutline failed, this build does not support OGR features./n" ); exit( 1 );#else // def OGR_ENABLED OGRRegisterAll();/* -------------------------------------------------------------------- *//* Open source vector dataset. *//* -------------------------------------------------------------------- */ OGRDataSourceH hSrcDS; hSrcDS = OGROpen( pszCutlineDSName, FALSE, NULL ); if( hSrcDS == NULL ) exit( 1 );/* -------------------------------------------------------------------- *//* 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 ) { fprintf( stderr, "Failed to identify source layer from datasource./n" ); exit( 1 ); }/* -------------------------------------------------------------------- *//* 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 ) { fprintf( stderr, "ERROR: Cutline feature without a geometry./n" ); exit( 1 ); } 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 { fprintf( stderr, "ERROR: Cutline not of polygon type./n" ); exit( 1 ); } OGR_F_Destroy( hFeat ); } if( OGR_G_GetGeometryCount( hMultiPolygon ) == 0 ) { fprintf( stderr, "ERROR: Did not get any cutline features./n" ); exit( 1 ); }/* -------------------------------------------------------------------- *//* 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 );#endif}/************************************************************************//* TransformCutlineToSource() *//* *//* Transform cutline from its SRS to source pixel/line coordinates.*//************************************************************************/static voidTransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline, char ***ppapszWarpOptions, char **papszTO_In ){#ifdef OGR_ENABLED 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 hCutlineSRS = OGR_G_GetSpatialReference( hMultiPolygon ); if( hRasterSRS != NULL && hCutlineSRS != NULL ) { /* ok, we will reproject */ } else if( hRasterSRS != NULL && hCutlineSRS == NULL ) { fprintf(stderr, "Warning : the source raster dataset has a SRS, but the cutline features/n" "not. We assume that the cutline coordinates are expressed in the destination SRS./n" "If not, cutline results may be incorrect./n"); } else if( hRasterSRS == NULL && hCutlineSRS != 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( hCutlineSRS != NULL ) { char *pszCutlineSRS_WKT = NULL; OSRExportToWkt( hCutlineSRS, &pszCutlineSRS_WKT ); papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszCutlineSRS_WKT ); CPLFree( pszCutlineSRS_WKT ); }/* -------------------------------------------------------------------- *//* Transform the geometry to pixel/line coordinates. *//* -------------------------------------------------------------------- */ CutlineTransformer oTransformer; /* The cutline transformer will *invert* the hSrcImageTransformer */ /* so it will convert from the cutline SRS to the source pixel/line */ /* coordinates */ 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 );#endif}需要将上面的代码稍微进行改造，也就是把exit(1)之类的去掉，换成返回错误代码，具体修改过的代码就不贴了，就那么几个，还有就是记得把函数的返回值void改成int，用于获取错误代码。下面就是我写的函数代码：
/*** @brief AOI截图(GDAL)* @param pszInFile 输入文件的路径* @param pszOutFile 截图后输出文件的路径* @param pszAOIFile AOI文件路径* @param pszSQL 指定AOI文件中的属性字段值来裁剪* @param pBandInex 指定裁剪的波段，默认为NULL，表示裁剪所有波段* @param piBandCount 指定裁剪波段的个数，同上一个参数同时使用* @param iBuffer 指定AOI文件外扩范围，默认为0，表示不外扩* @param pszFormat 截图后输出文件的格式* @param pProgress 进度条指针* @return 成功返回*/ int CutImageByAOIGDAL(const char* pszInFile, const char* pszOutFile, const char* pszAOIFile, const char* pszSQL, int *pBandInex, int *piBandCount, int iBuffer, const char* pszFormat, LT_Progress *pProgress){ if(pProgress != NULL) { pProgress->SetProgressCaption("AOI裁剪"); pProgress->SetProgressTip("开始裁剪图像..."); } GDALAllRegister(); void *hCutLine = NULL; int iRev = LoadCutline( pszAOIFile, "", "", pszSQL, &hCutline ); GDALDataset * pSrcDS = (GDALDataset*) GDALOpen(pszInFile, GA_ReadOnly); if (pSrcDS == NULL) { if (pProgress != NULL) pProgress->SetProgressTip("输入的栅格文件不能打开，请检查文件是否存在！"); return RE_NOFILE; } int iBandCount = pSrcDS->GetRasterCount(); const char* pszWkt = pSrcDS->GetProjectionRef(); GDALDataType dT = pSrcDS->GetRasterBand(1)->GetRasterDataType(); double adfGeoTransform[6] = {0}; double newGeoTransform[6] = {0}; pSrcDS->GetGeoTransform(adfGeoTransform); memcpy(newGeoTransform, adfGeoTransform, sizeof(double)*6); int *pSrcBand = NULL; int iNewBandCount = iBandCount; if (pBandInex != NULL && piBandCount != NULL) { int iMaxBandIndex = pBandInex[0]; for (int i=1; i<*piBandCount; i++) { if (iMaxBandIndex < pBandInex[i]) iMaxBandIndex = pBandInex[i]; } if (iMaxBandIndex > iBandCount) { if (pProgress != NULL) pProgress->SetProgressTip("输入的波段序号没有指定的波段数！"); GDALClose((GDALDatasetH) pSrcDS); return RE_PARAMERROR; } iNewBandCount = *piBandCount; pSrcBand = new int[iNewBandCount]; for (int i=0; i<iNewBandCount; i++) pSrcBand[i] = pBandInex[i]; } else { pSrcBand = new int[iNewBandCount]; for (int i=0; i<iNewBandCount; i++) pSrcBand[i] = i+1; } OGRGeometryH hGeometry = (OGRGeometryH) hCutline; OGRGeometryH hMultiPoly = NULL; if (iBuffer > 0) { double dDistance = ABS(adfGeoTransform[1]*iBuffer); hMultiPoly = OGR_G_Buffer(hGeometry, dDistance, 30); OGR_G_AssignSpatialReference(hMultiPoly, OGR_G_GetSpatialReference(hGeometry)); OGR_G_DestroyGeometry(hGeometry); } else hMultiPoly = hGeometry; OGREnvelope eRect; OGR_G_GetEnvelope(hMultiPoly, &eRect); int iNewWidth = 0, iNewHeigh = 0; int iBeginRow = 0, iBeginCol = 0; newGeoTransform[0] = eRect.MinX; newGeoTransform[3] = eRect.MaxY; iNewWidth = static_cast<int>((eRect.MaxX -eRect.MinX) / ABS(adfGeoTransform[1])); iNewHeigh = static_cast<int>((eRect.MaxY -eRect.MinY) / ABS(adfGeoTransform[5])); Projection2ImageRowCol(adfGeoTransform, newGeoTransform[0], newGeoTransform[3], iBeginCol, iBeginRow); ImageRowCol2Projection(adfGeoTransform, iBeginCol, iBeginRow, newGeoTransform[0], newGeoTransform[3]); GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); if (pDriver == NULL) { if (pProgress != NULL) pProgress->SetProgressTip("不能创建指定类型的文件！"); GDALClose((GDALDatasetH) pSrcDS); return RE_NOFILE; } GDALDataset* pDstDS = pDriver->Create(pszOutFile, iNewWidth, iNewHeigh, iNewBandCount, dT, NULL); if (pDstDS == NULL) { if (pProgress != NULL) pProgress->SetProgressTip("创建文件失败！"); GDALClose((GDALDatasetH) pSrcDS); return RE_CREATEFILE; } pDstDS->SetGeoTransform(newGeoTransform); pDstDS->SetProjection(pszWkt); void *hTransformArg, *hGenImgProjArg=NULL; char **papszTO = NULL; /* -------------------------------------------------------------------- */ /* Create a transformation object from the source to */ /* destination coordinate system. */ /* -------------------------------------------------------------------- */ hTransformArg = hGenImgProjArg = GDALCreateGenImgProjTransformer2( hSrcDS, (GDALDatasetH)pDstDS, papszTO ); if( hTransformArg == NULL ) { if (pProgress != NULL) pProgress->SetProgressTip("创建投影失败，请检查输入参数！"); GDALClose((GDALDatasetH) pSrcDS); GDALClose((GDALDatasetH) pDstDS); RELEASE(pSrcBand); return RE_PARAMERROR; } GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform; GDALWarpOptions *psWO = GDALCreateWarpOptions(); psWO->papszWarpOptions = CSLDuplicate(NULL); psWO->eWorkingDataType = dT; psWO->eResampleAlg = GRA_Bilinear ; psWO->hSrcDS = (GDALDatasetH) pSrcDS; psWO->hDstDS = (GDALDatasetH) pDstDS; psWO->pfnTransformer = pfnTransformer; psWO->pTransformerArg = hTransformArg; psWO->pfnProgress = GDALProgress; psWO->pProgressArg = pProgress; psWO->nBandCount = iNewBandCount; psWO->panSrcBands = (int *) CPLMalloc(iNewBandCount*sizeof(int)); psWO->panDstBands = (int *) CPLMalloc(iNewBandCount*sizeof(int)); for (int i=0; i<iNewBandCount; i++) { psWO->panSrcBands[i] = pSrcBand[i]; psWO->panDstBands[i] = i+1; } RELEASE(pSrcBand); psWO->hCutline = (void*) hMultiPoly; TransformCutlineToSource((GDALDatasetH) pSrcDS, (void*)hMultiPoly, &(psWO->papszWarpOptions), papszTO ); GDALWarpOperation oWO; if (oWO.Initialize(psWO) != CE_None) { if(pProgress != NULL) pProgress->SetProgressTip("转换参数错误！"); GDALClose((GDALDatasetH) pSrcDS); GDALClose((GDALDatasetH) pDstDS); return RE_PARAMERROR; } oWO.ChunkAndWarpImage(0, 0, iNewWidth, iNewHeigh); GDALDestroyGenImgProjTransformer(psWO->pTransformerArg); GDALDestroyWarpOptions( psWO ); GDALClose((GDALDatasetH) pSrcDS); GDALClose((GDALDatasetH) pDstDS); if(pProgress != NULL) pProgress->SetProgressTip("图像裁剪完成！"); return RE_SUCCESS;}裁剪后的效果图就不贴了，希望对大家有用。
………………………………………………分割线 下面更新于2013-5-2…………………………………………………………
上面的代码中有两处粘贴错误，已经修改，同时有人问代码中有两个坐标转换的函数，就是行列号和投影坐标相互转换的函数，其实很简单，就是根据六参数计算各仿射变换就OK，下面把代码贴出来。
bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow){try{double dTemp = adfGeoTransform[1]*adfGeoTransform[5] - adfGeoTransform[2]*adfGeoTransform[4];double dCol = 0.0, dRow = 0.0;dCol = (adfGeoTransform[5]*(dProjX - adfGeoTransform[0]) - adfGeoTransform[2]*(dProjY - adfGeoTransform[3])) / dTemp + 0.5;dRow = (adfGeoTransform[1]*(dProjY - adfGeoTransform[3]) - adfGeoTransform[4]*(dProjX - adfGeoTransform[0])) / dTemp + 0.5;iCol = static_cast<int>(dCol);iRow = static_cast<int>(dRow);return true;}catch(...){return false;}}bool ImageRowCol2Projection(double *adfGeoTransform, int iCol, int iRow, double &dProjX, double &dProjY){//adfGeoTransform[6] 数组adfGeoTransform保存的是仿射变换中的一些参数，分别含义见下//adfGeoTransform[0] 左上角x坐标 //adfGeoTransform[1] 东西方向分辨率//adfGeoTransform[2] 旋转角度, 0表示图像 "北方朝上"//adfGeoTransform[3] 左上角y坐标 //adfGeoTransform[4] 旋转角度, 0表示图像 "北方朝上"//adfGeoTransform[5] 南北方向分辨率try{dProjX = adfGeoTransform[0] + adfGeoTransform[1] * iCol + adfGeoTransform[2] * iRow;dProjY = adfGeoTransform[3] + adfGeoTransform[4] * iCol + adfGeoTransform[5] * iRow;return true;}catch(...){return false;}}
- 如何使用GDAL进行AOI裁剪
- 如何使用GDAL进行图像镶嵌
- 使用IDirectDrawClipper进行裁剪
- 使用IDirectDrawClipper进行裁剪
- 如何进行研发过程裁剪
- 使用GDAL进行栅格数据读取示例
- 使用GDAL对HDF数据进行校正
- 使用GDAL对DEM进行彩色渲染
- 使用GDAL进行RPC坐标转换
- 如何使用GDAL重采样图像
- Ubantu 下如何安装使用GDAL库
- 如何使用GDAL重采样图像 .
- 如何进行研发项目中的过程裁剪
- 如何对PCB界面进行裁剪
- Linux下使用GDAL进行开发（automake使用）
- 使用 Java 进行图像处理 - 图像裁剪
- 使用裁剪测试进行多窗口渲染
- 使用ThinkPHP进行图片批量裁剪
- 与JDK无关的吐槽
- QSplitter
- 一个打印菱形的小例子
- linux、windows中实现gettickcount()
- 云计算与分布式计算的区别？
- 如何使用GDAL进行AOI裁剪
- 最新翻译书籍《物理数据库设计—索引、视图和存储技术》出版
- 华中电网7.1停电事故教训
- 相对布局简解
- 刻录mac安装盘
- jquery + php 删除数据库中多条记录
- （转）游戏程序员养成计划
- JAVA中的向上转型与向下转型
- 修正RO6.0.43中Oracle取字段长度不正确的问题.
