如何使用GDAL进行AOI裁剪

    在工作中，会经常使用一个行政区的矢量图去裁剪一个遥感影像图，得到该行政区的影像图，一般的商业遥感软件都具有这个功能。今天就是用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 void
LoadCutline( 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 void
TransformCutlineToSource( 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;
	}
}