分区统计算法实现（C++ GDAL）

Arcgis和QGIS都有分区统计工具，geopanda库也有相关函数，但无法满足我们课题组的需求。其中Arcgis只能输出一张表格，还需要进行连接，QGIS只能输出至新shp文件，geopandas计算速度相比慢了点。

我们课题组需要对上千万的地块以及TB级别的卫星遥感影像进行分区统计计算，所以我重造了这个轮子，处理速度与QGIS基本一致，且能输出至原始shp文件的属性表以便批量处理。

代码仓库：gitee.com/gislxz/zonal-statistics

一开始倒是想用现成的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来实现存在以下几点不足：

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

2不能精确捕获NODATA值

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

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

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

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

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

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

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

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

另外要特别注意矢量与栅格边缘相交的情况，即多边形只有一部分与栅格相交。

完整代码见代码仓库

template <typename T>
void ZonalStatistics<T>::MyStatistics()
{
	//检查并初始化
	checkAndInitial();
	//计算矢量图层计算范围
	MyParcelManager->getLayerExtent();
	//根据范围逐行读取栅格文件
	//防止矢量范围越界
	//lineNum_begin起始绝对行号
	int lineNum_begin = MyParcelManager->layerPixelExtent->top > 0 ? MyParcelManager->layerPixelExtent->top : 0;
	//leftColumn起始列的绝对列号
	int leftColumn = MyParcelManager->layerPixelExtent->left > 0 ? MyParcelManager->layerPixelExtent->left : 0;
	//lineWidth读取的长度
	int iSrcWidth = MyRasterdata->GetRasterXSize();
	int lineWidth = (MyParcelManager->layerPixelExtent->right < iSrcWidth ? MyParcelManager->layerPixelExtent->right : iSrcWidth - 1) - leftColumn + 1;
	//lineHeight读取的绝对高度
	int iSrcHeight = MyRasterdata->GetRasterYSize();
	int lineHeight = MyParcelManager->layerPixelExtent->bottom < iSrcHeight ? MyParcelManager->layerPixelExtent->bottom : iSrcHeight - 1;
	//循环所有要素，计算掩膜提取区域
	MyParcelManager->initialAllParcel(iSrcHeight);
	//读取波段
	auto* poBand = MyRasterdata->GetRasterBand(1);
	//读取图层
	auto poLayer = MyParcelManager->CRSFile->GetLayer(0);
	//分配内存
	T* pafScanline = (T*)CPLMalloc(sizeof(T) * lineWidth);
	for (int lineNum = lineNum_begin; lineNum <= lineHeight; lineNum++) {
		//读取一行
		poBand->RasterIO(GF_Read, leftColumn, lineNum, lineWidth, 1, pafScanline, lineWidth, 1, uDT, 0, 0);
		//遍历begin链表
		for (auto* beginParcel : *(MyParcelManager->beginList[lineNum - lineNum_begin])) {
			//parcelUnit初始化
			MyParcelManager->initialParcelMask(beginParcel, iSrcWidth, iSrcHeight);
			//加入处理链表
			MyParcelManager->processList->insert(beginParcel);
		}
		//对处理链表内的所有要素输入数据
		typename std::unordered_set<ParcelUnit<T>*>::iterator proParcel = MyParcelManager->processList->begin();
		while (proParcel != MyParcelManager->processList->end()) {
			//发送栅格数据
			MyParcelManager->assignByLine(pafScanline, *proParcel, lineNum, leftColumn, lineWidth, NoDataValue ,lineNum_begin);
			proParcel++;
		}
		//遍历release链表
		for (auto* reParcel : *(MyParcelManager->releaseList[lineNum - lineNum_begin])) {
			//从处理列表中释放，释放之前计算出结果
			//计算结果并写入属性列表
			MyParcelManager->finalStatistic(reParcel, myAttributeIndex, &attriNameIndex, poLayer);
			//释放要素
			MyParcelManager->releaseParcel(reParcel);
			//从处理链表中删除
			MyParcelManager->processList->erase(reParcel);
		}
	}
	CPLFree(pafScanline);
	std::cout << "statistics...finish!" << std::endl;
}

GDAL学习笔记及打包后的工具见个人网站：www.gislxz.com