GDAL：利用Shape文件对遥感影像进行裁剪

支持单个或多个多边形要素的shp文件，对于没有仿射信息的影像裁剪不成功。
直接上代码：

//ImageCut.h
#include "gdal_priv.h"
#include <iostream>
#include "ogr_spatialref.h"
#include "ogrsf_frmts.h"
#include "gdalwarper.h"

using namespace std;

class ImageCut{
public:
	ImageCut();
	~ImageCut();
	/*设置影像输入路径*/
	void setSrcImageFilename(const char* pszSrcFile);
	/*设置裁剪后影像输出的路径*/
	void setDstImageFilename(const char* pszDstFile);
	/*设置shape文件的路径*/
	void setShapeFilename(const char*pszShapeName);
	/*设置输出影像的格式*/
	void setImageFormat(const char* pszFormat);
	/*开始影像AOI裁剪*/
	int ImageCutByAOI();
private:
	const char* pszSrcFile;
	const char* pszDstFile;
	const char*	pszShapeName;
	const char* pszFormat;
	/*输入影像的指针*/
	GDALDataset *ptrSrcImage;
	/*投影坐标转行列号*/
	bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow);
	/*行列号转投影坐标*/
	bool ImageRowCol2Projection(double *adfGeoTransform, int iCol, int iRow, double &dProjX, double &dProjY);
	/*根据shape文件和输入的影像输出多边形WKT文件*/
	int outputPolygonWkt(string &pszAOIWKT);
	/*影像裁剪的代码*/
	int ImageCutting(const char* pszAOIWKT);

};

//ImageCut.cpp
#include "ImageCut.h"
ImageCut::ImageCut(){
	GDALAllRegister();
	OGRRegisterAll();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	CPLSetConfigOption("SHAPE_ENCODING", "");
	pszDstFile = "";
	pszFormat = "";
	pszShapeName = "";
	pszSrcFile = "";
	ptrSrcImage = NULL;
}
ImageCut::~ImageCut(){
	GDALClose(ptrSrcImage);//关闭源影像
}
void ImageCut::setSrcImageFilename(const char* pszSrcFile){
	this->pszSrcFile = pszSrcFile;
}
void ImageCut::setDstImageFilename(const char* pszDstFile){
	this->pszDstFile = pszDstFile;
}
void ImageCut::setShapeFilename(const char*pszShapeName){
	this->pszShapeName = pszShapeName;
}
void ImageCut::setImageFormat(const char* pszFormat){
	this->pszFormat = pszFormat;
}
int ImageCut::ImageCutByAOI(){
	string pszAOIWKT = "";
	int isoutput = outputPolygonWkt(pszAOIWKT);
	if (isoutput == 0)
	{
		cout << "output polygon wkt false" << endl;
		return 0;
	}
	int iscut = ImageCutting(pszAOIWKT.c_str());
	if (iscut != 1)
	{
		cout << "Image Cut false" << endl;
		return 0;
	}
	return 1;
}

int ImageCut::outputPolygonWkt(string &pszAOIWKT){
	ptrSrcImage = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
	if (ptrSrcImage == NULL)
	{
		cout << "read image false" << endl;
		return 0;
	}
	double SrcGeoTransform[6] = { 0 };
	ptrSrcImage->GetGeoTransform(SrcGeoTransform);
	//读取shape数据--------------
	OGRDataSource *ptrShape = OGRSFDriverRegistrar::Open(pszShapeName);
	if (ptrShape == NULL)
	{
		cout << "read shape false" << endl;
		return 0;
	}
	int numLayer = ptrShape->GetLayerCount();//获取图层个数
	if (numLayer != 1)
	{
		cout << "the method only support one layer" << endl;
		return 0;
	}
	OGRLayer *ptrLayer = ptrShape->GetLayer(0);//这里只支持一个图层
	if (ptrLayer == NULL)
	{
		cout << "open the layer false" << endl;
		return 0;
	}
	ptrLayer->ResetReading();//将图层重置一下，标号从0开始
	//OGRFeatureDefn *ptrDefn = ptrLayer->GetLayerDefn();//获取当前图层的属性表结构,主要包涵一些字段信息

	OGRFeature *poFeature = NULL;//定义包含几何和属性的要素指针


	//统计一下图层内多边形个数
	int count_ = 0, cout_1(0);
	while ((poFeature = ptrLayer->GetNextFeature()) != NULL){
		count_++;
	}
	poFeature = NULL;
	ptrLayer->ResetReading();
	//-*-----------------------------------------------------------
	//string pszAOIWKT = "";//待输出的多边形wkt文件
	pszAOIWKT.clear();
	if (count_ == 1)
	{
		pszAOIWKT = "POLYGON ((";
	}
	else
	{
		pszAOIWKT = "MULTIPOLYGON (((";
	}
	while ((poFeature = ptrLayer->GetNextFeature()) != NULL)
	{
		OGRGeometry *ptrGeometry = poFeature->GetGeometryRef();//获得要素的几何属性
		if (wkbFlatten(ptrGeometry->getGeometryType()) != wkbPolygon)
		{
			cout << "Only polygon feature cutting images are supported(wkbPolygon)" << endl;
			return 0;
		}
		OGRPolygon *ptrPolygon = (OGRPolygon *)ptrGeometry->clone();//将几何要素指针强制转化为多边形要素指针
		OGRLinearRing *ptrOGRLinearRing = ptrPolygon->getExteriorRing();//获取多边形的外环点的指针
		int numExterRing = ptrOGRLinearRing->getNumPoints();//外环点数
		int numInteriorRing = ptrPolygon->getNumInteriorRings();//内环数(暂时不加入考虑)
		double x, y;//存储临时的xy值
		int row_, col_;//存储临时的行列值
		for (int i = 0; i < numExterRing; i++)
		{
			if (i != 0){
				pszAOIWKT += ",";
			}
			x = ptrOGRLinearRing->getX(i);
			y = ptrOGRLinearRing->getY(i);
			row_, col_;
			Projection2ImageRowCol(SrcGeoTransform, x, y, row_, col_);
			pszAOIWKT += (to_string(row_) + " " + to_string(col_));
		}
		cout_1++;
		if (count_ == 1)
		{
			pszAOIWKT += "))";
		}
		else if (cout_1 == count_)
		{
			pszAOIWKT += ")))";
		}
		else
		{
			pszAOIWKT += ")),((";
		}
		OGRFeature::DestroyFeature(poFeature);
	}

	//GDALClose(ptrSrcImage);//关闭影像，放在了析构函数，避免重复打开影像
	OGRDataSource::DestroyDataSource(ptrShape);
	OGRCleanupAll();
}
int ImageCut::ImageCutting(const char* pszAOIWKT){
	// 打开原始图像并计算图像信息
	//GDALDataset *pSrcDS = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
	GDALDataType eDT = ptrSrcImage->GetRasterBand(1)->GetRasterDataType();

	int iBandCount = ptrSrcImage->GetRasterCount();
	int iSrcWidth = ptrSrcImage->GetRasterXSize();
	int iSrcHeight = ptrSrcImage->GetRasterYSize();

	double pSrcGeoTransform[6] = { 0 };
	double pDstGeoTransform[6] = { 0 };
	ptrSrcImage->GetGeoTransform(pSrcGeoTransform);			//图像的仿射变换信息
	memcpy(pDstGeoTransform, pSrcGeoTransform, sizeof(double) * 6);

	// 将传入的AOI的WKT处理为一个OGRGeometry类型，用于后续处理
	char *pszWKT = (char*)pszAOIWKT;
	OGRGeometry* pAOIGeometry = NULL;
	//cout << *pszWKT << endl;
	if (*pszWKT == 'P')
	{
		pAOIGeometry = OGRGeometryFactory::createGeometry(wkbPolygon);
	}
	else
	{
		pAOIGeometry = OGRGeometryFactory::createGeometry(wkbMultiPolygon);
	}

	pAOIGeometry->importFromWkt(&pszWKT);

	OGREnvelope eRect;
	pAOIGeometry->getEnvelope(&eRect);

	// 设置输出图像的左上角坐标
	GDALApplyGeoTransform(pSrcGeoTransform, eRect.MinX, eRect.MinY, (&pDstGeoTransform[0]), &(pDstGeoTransform[3]));
	// 根据裁切范围确定裁切后的图像宽高
	int iDstWidth = static_cast<int>(eRect.MaxX - eRect.MinX);
	int iDstHeight = static_cast<int>(eRect.MaxY - eRect.MinY);

	// 创建输出图像
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
	GDALDataset *pDstDS = poDriver->Create(pszDstFile, iDstWidth, iDstHeight, iBandCount, eDT, NULL);
	pDstDS->SetGeoTransform(pDstGeoTransform);
	pDstDS->SetProjection(ptrSrcImage->GetProjectionRef());

	// 构造坐标转换关系
	void *hTransformArg = GDALCreateGenImgProjTransformer2((GDALDatasetH)ptrSrcImage, (GDALDatasetH)pDstDS, NULL);
	GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;

	// 构造GDALWarp的变换选项
	GDALWarpOptions *psWO = GDALCreateWarpOptions();

	psWO->papszWarpOptions = CSLDuplicate(NULL);
	psWO->eWorkingDataType = eDT;
	psWO->eResampleAlg = GRA_NearestNeighbour;

	psWO->hSrcDS = (GDALDatasetH)ptrSrcImage;
	psWO->hDstDS = (GDALDatasetH)pDstDS;

	psWO->pfnTransformer = pfnTransformer;
	psWO->pTransformerArg = hTransformArg;

	psWO->nBandCount = iBandCount;
	psWO->panSrcBands = (int *)CPLMalloc(psWO->nBandCount*sizeof(int));
	psWO->panDstBands = (int *)CPLMalloc(psWO->nBandCount*sizeof(int));
	for (int i = 0; i < iBandCount; i++)
	{
		psWO->panSrcBands[i] = i + 1;
		psWO->panDstBands[i] = i + 1;
	}

	// 设置裁切AOI，AOI中的坐标必须是图像的行列号坐标，否则不能进行裁切
	psWO->hCutline = (void*)pAOIGeometry;
	// 设置上面的hCutline的值和使用下面的CUTLINE配置项的效果一样，两个选择一个即可
	psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "CUTLINE", pszAOIWKT);

	// 创建GDALWarp执行对象，并使用GDALWarpOptions来进行初始化
	GDALWarpOperation oWO;
	oWO.Initialize(psWO);

	// 执行处理
	oWO.ChunkAndWarpImage(0, 0, iDstWidth, iDstHeight);

	// 释放资源和关闭文件
	GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
	GDALDestroyWarpOptions(psWO);

	GDALClose((GDALDatasetH)pDstDS);

	return 1;
}
bool ImageCut::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 ImageCut::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;
	}
}

#include "ImageCut.h"

int main(){
	const char *ptrShapeFilename = "D:/shp/shp2.shp";
	const char *ptrImageFilename = "D:/software/practice/GF2_PMS2_E100.6_N26.5_20190104_L1A0003735651_AC - 副本.tiff";
	//const char *ptrImageFilename = "F:/测试数据/GF1/GF1_PMS2_E113.2_N23.5_20180113_L1A0002927191-MSS2.tiff";
	const char *ptrDstImageFilename = "D:/OUTPUT5.tiff";
	const char *pszFormat = "GTiff";
	ImageCut ptrImageCut;
	ptrImageCut.setSrcImageFilename(ptrImageFilename);
	ptrImageCut.setDstImageFilename(ptrDstImageFilename);
	ptrImageCut.setImageFormat(pszFormat);
	ptrImageCut.setShapeFilename(ptrShapeFilename);
	cout << "开始AOI裁剪" << endl;
	ptrImageCut.ImageCutByAOI();
	cout << "计算结束" << endl;
	system("pause");
	return 1;
}

参考李明录大哥的代码：
https://blog.csdn.net/liminlu0314/article/details/6136512?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522158717665319195239802256%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fblog.%2522%257D&request_id=158717665319195239802256&biz_id=0&utm_source=distribute.pc_search_result.none-task-blog-2_blogfirst_rank_v2~rank_v25-1