C++ GDAL开发实现分区统计功能

代码仓库：ZonalStatistics: GDAL分区统计实现 (gitee.com)

背景

实现这个算法的想法萌生于去年，我的导师指导我去协助参与重庆农业项目，任务是更新重庆市所有地块一年的NDVI数据，这个任务看似简单但工作量极大，除去数据准备外，该任务最核心的步骤就是使用分区统计将栅格格式的NDVI数据赋值到以多边形为组织形式的地块数据中。然而就是这样简单的分区统计功能却耗费了我们大量时间，ArcGIS和QGIS软件提供的分区统计功能并不足以支持我们高效精准的完成这项任务，只能退而求其次选择了中心点赋值法完成了NDVI赋值的任务，但这样不利于数据不确定性的控制，进而影响后续计算与预测，造成误差的传递和放大。

现有GIS软件分区统计功能的不足

Arcgis：因为我只有arcgis10.8的版本，不太清楚新的Arcgis Pro 3的算法实现怎么样，但10.8版本中分区统计算法太慢了，比QGIS慢数十倍，输出结果是一个表格，还要再进行字段连接，不太行。

QGIS：QGIS处理倒挺快的，问题是不能写入原文件，只能输出一个新矢量文件，那我一次运算12个月的栅格图像，没法批处理运行了。

思来想去还是自己重写一个

使用GDALWarp来实现

一开始倒是想用现成的API来实现的，GDALWarp中有图像裁切的API可以拿来用。分区统计本质上就是统计每一个多边形内部像元点的值，其实也就是要对图像进行裁切。只不过我们不需要将裁切后的图像保存为新文件，只用把像元值存储的二维数组输入到内存里在对数组进行运算就好了。遗憾的是GDAL学习资料太少，李民录先生写的《GDAL源码剖析与开发指南》介绍了Warp的基本用法，但是并没说怎么输入到内存，只说了输入到内存是用WarpRegionToBuffer这个执行语句。我也是进行了一番尝试，这个Warp注意事项主要有两点，一是输入的AOI范围需要是多边形或多多边形，其坐标应该是像素行列号而不是地理坐标，且要求是WKTGeo格式，好在OGRGeometry类是有exportToWkt的函数。第二点就是使用Warp做裁切还需要设置坐标转换的参数，书中给的示例是用的GDALCreateGenImgProjTransformer2这个函数，输入输入和输出图像对应的GDALDataset就行了，但我们不需要输出到文件，自然也没有创建输出栅格图像的GDALDataset，如果创建势必会极大影响处理速度，于是我去查了其他几个GDALCreateGenImgProjTransformer函数，发现可以这样替代。

//书中示例代码
void *hTransformArg=GDALCreateGenImgProjTransformer2((GDALDatasetH)pSrcDS, (GDALDatasetH)pDstDS,NULL);
//只创建输出图像的地理转换六参数，调用两次原图像的OGRSpatialReference
void* hTransformArg = GDALCreateGenImgProjTransformer4((OGRSpatialReferenceH*)MyRasterdata->GetSpatialRef(), pSrcGeoTransform, (OGRSpatialReferenceH*)MyRasterdata->GetSpatialRef(), pDstGeoTransform, NULL);

调用GDALWarp进行裁剪的完整函数代码

//调用GDALWarp进行裁剪（结果输出至内存，相当于掩膜提取）
template <typename T>
T* ZonalStatistics<T>::imageCutByAoi(const OGRGeometry& poGeometry, int& size) {
	GDALAllRegister();

	//Geometry转行列号WKTPolygon
	OGRGeometry* pixelGeometry = NULL;
	//OGRPolygon* Pixelpolygon = (OGRPolygon*)OGRGeometryFactory::createGeometry(wkbPolygon);
	MyParcelManager->PLtoXY(poGeometry, &pixelGeometry);

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

	//获取图像尺寸
	int iBandCount = MyRasterdata->GetRasterCount();
	int iSrcWidth = MyRasterdata->GetRasterXSize();
	int iSrcHeight = MyRasterdata->GetRasterYSize();

	//OGRGeometry转WKT
	std::string temp = pixelGeometry->exportToWkt();
	char* WKTGeo = new char[temp.size() + 1];
	strcpy(WKTGeo, temp.c_str());
	//std::cout << WKTGeo << std::endl;
	OGREnvelope poRect;
	pixelGeometry->getEnvelope(&poRect);
	int iDstWidth = static_cast<int>(poRect.MaxX - poRect.MinX);
	int iDstHeight = static_cast<int>(poRect.MaxY - poRect.MinY);

	// 设置输出图像的左上角坐标
	GDALApplyGeoTransform(pSrcGeoTransform, poRect.MinX, poRect.MinY, (&pDstGeoTransform[0]), &(pDstGeoTransform[3]));

	//构造坐标转换关系
	void* hTransformArg = GDALCreateGenImgProjTransformer4((OGRSpatialReferenceH*)MyRasterdata->GetSpatialRef(), pSrcGeoTransform, (OGRSpatialReferenceH*)MyRasterdata->GetSpatialRef(), pDstGeoTransform, NULL);
	GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;

	T* pBuf = (T*)CPLMalloc(sizeof(float) * iDstWidth * iDstHeight);
	size = iDstWidth * iDstHeight;

	//构造GDALWarp的变换选项
	GDALWarpOptions* psWO = GDALCreateWarpOptions();
	psWO->papszWarpOptions = CSLDuplicate(NULL);
	psWO->eWorkingDataType = uDT;
	psWO->eResampleAlg = GRA_NearestNeighbour;
	//测试发现不能正常设置nodata值，还未解决
	double* nodata = new double{ static_cast<double>(NoDataValue) };
	psWO->padfDstNoDataReal = nodata;

	psWO->hSrcDS = (GDALDatasetH)(MyRasterdata.get());
	psWO->pfnTransformer = pfnTransformer;
	psWO->pTransformerArg = hTransformArg;
	psWO->nBandCount = 1;
	psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "CUTLINE", WKTGeo);
	psWO->panSrcBands = (int*)CPLMalloc(psWO->nBandCount * sizeof(int));
	psWO->panSrcBands[0] = 1;
	//创建GDALWarp执行对象，并使用GDALWarpOptions
	GDALWarpOperation oWO;
	oWO.Initialize(psWO);


	//执行处理(输出到内存)
	//oWO.ChunkAndWarpImage(0, 0, iDstWidth, iDstHeight);
	if (oWO.WarpRegionToBuffer(0, 0, iDstWidth, iDstHeight, pBuf, uDT)) {
		std::cout << "cut failed" << std::endl;
		delete[] WKTGeo;
		CPLFree(pBuf);
		// 释放资源和关闭文件
		GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
		GDALDestroyWarpOptions(psWO);
		exit(0);
	}
	else {
		delete[] WKTGeo;
		// 释放资源和关闭文件
		GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
		GDALDestroyWarpOptions(psWO);
	}
	return pBuf;
}

测试后发现调用GDALWarp来实现存在以下几点不足：

1处理速度比QGIS还是慢四到五倍

2不能精确捕获NODATA值

3 面积过小的地块发生错误

于是我只能尝试从头实现，自己来写裁切功能。

从头开始底层实现

算法设计思路

在判断多边形与栅格图像相交范围中，我采用了中心点判断法，及像元中心点被多边形所包含，即认为多边形与栅格图像的相交范围包含该像元。

计算多边形与栅格图像相交范围步骤：

1通过栅格图像行列号转地理坐标的转换六参数反推地理坐标转行列号的六参数。（需注意转换六参数中的起始点为左上角像素点的左上角位置而非中心点，因此反算时需要将起始点的坐标分别加上对应方向像素宽度的一半）

2对任一多边形遍历其外环（ExteriorRing），将得到的点坐标转为行列号，通过前后两个点连成的直线计算跨越行号时与对应直线的交点。

3对每一行的交点进行排序，从小到大即从左到右，这样包含在俩点间的像素点必然都在多边形内。对于凹多边形，多多边形，内环等情况该算法都适用。

在获得多边形包含的像元点后即可读取栅格将对应像素点的数值赋给多边形，再进行统计计算。这里我注意到如果遍历多边形以其四至范围去读取栅格并赋值会造成极大的硬盘IO开销，因而将策略调整成逐行读取栅格图像，并对该行所有相交的多边形进行赋值处理。

回忆NDVI赋值的任务，每月对应的一幅NDVI图像的尺寸，角度和范围可以认为是相同的，这就意味着每个多边形计算的相交（掩膜）范围是可以重复利用的，因此我又增加了一个批处理加速的策略，当输入多幅图像时可以使用批处理加速，在处理完每个地块时不释放相交范围和记录像元值缓冲区而只是重置有效像元计数，用更多的内存占用换取处理时间上的减少。

算法结构设计

说明：

1整个代码使用Entrance函数作为入口，打开第一幅栅格图像以获取栅格数据类型从而创建相应ZonalStatistics（模板类），并把打开栅格图像所用的GDALDatasetUniquePtr转给（std::move）ZonalStatistics。

2 ZonalStatistics类中具有一个ParcelManager指针，在构造函数中会构造一个其独有的ParcelManager。ZonalStatistics类包含了栅格图像文件的指针，需要添加的属性信息，其成员函数void poiStatistics()和void MyStatistics()分别实现了借助GDALWrap裁剪图像和底层构建的分区统计算法。

3 ParcelManager作为管理ParcelUnit的类，其关键成员在于具有两个std::forward_list<ParcelUnit<T>*>**成员分别用于记录栅格图像某一行需要启动和结算的地块对象，将其加入处理队列或从处理队列中删除并进行结算（批处理加速中只重置计数）。处理队列使用std::unordered_set<ParcelUnit<T>*>*结构实现方便插入和删除。

4 ParcelUnit作为地块单元，使用自定义结构体MyExtent记录几何体范围，以便插入ParcelManager的启动和结算队列；使用std::vector<double>**记录相交（掩膜）范围；使用double数组记录包含的像元值。同时其具有一个中心点坐标和对应像元值（如果在栅格图像范围内）以便处理地块过小不包含任一像元中心点的情况。

结果对比

实验数据为

·重庆市2020年10月及2020年11月NDVI合成图像（单波段，尺寸54615×44970，单幅数据大小9.89GB）

·重庆市綦江区耕地地块数据，共761496个几何体，包含凹多边形，多多边形，内环等各种情况。

单位（秒/s）	分项	QGIS	分项	本算法
处理时间（单幅图像）	总计	107s	总计	116s
	计算得到临时图层	64s	向原文件属性表中添加列	67s
	写入硬盘	43s	计算并写入	49s
处理第二幅图象	总计	121s	总计（批处理加速）	97s