C++矢量转栅格

本文主要参考了以下三篇博文资料：

1. gdal矢量栅格化接口GDALRasterizeLayers介绍。https://blog.csdn.net/clever101/article/details/10419911
2. gdal矢量转栅格，这篇文章里的代码有一些小问题，已在本文修正。http://www.voidcn.com/article/p-udriyjzi-om.html
3. GDAL2.0读写矢量。https://blog.csdn.net/xzhh19921019/article/details/53408099

主要利用gdal提供的接口GDALRasterizeLayers，构建接口需要的参数，调用接口转化即可。

先介绍接口GDALRasterizeLayers。函数原型如下：

CPLErr CPL_DLL
GDALRasterizeLayers( GDALDatasetH hDS,
                     int nBandCount, int *panBandList,
                     int nLayerCount, OGRLayerH *pahLayers,
                     GDALTransformerFunc pfnTransformer,
                     void *pTransformArg,
                     double *padfLayerBurnValues,
                     char **papszOptions,
                     GDALProgressFunc pfnProgress,
                     void *pProgressArg );

接口参数说明：

hDS：输出的栅格数据集

nBandCount：栅格要更新的波段数

panBandList：要更新的波段列表，列表长度为nBandCount

nLayerCount：矢量图层数

pahLayers：矢量图层数组，数组长度为nLayerCount

pfnTransformer：将几何图形转换为栅格图像的像素/行列的转换器，如果该参数为空，则在内部创建一个

pTransformArg：将几何图形转换为栅格图像的像素/行列的转换器所用到的参数

padfLayerBurnValues：指定各个输出波段的输出像素值，长度nBandCount

papszOptions：控制栅格化的一系列选项值，具体在下面会介绍。

pfnProgress：gdal进度函数

pProgressArg：传给gdal进度函数的参数

papszOptions选项介绍

"ATTRIBUTE":

指定属性字段中的字段值作为栅格值写入栅格文件中，该值将输出到所有输出波段中。假如该值指定了，padfLayerBurnValue参数的值将失效并且padfLayerBurnValues参数可以设置为NULL。

"CHUNKYSIZE":

指定该运行操作的块的高度。该值越大所需的计算时间越小。如果该值没有设置或者设置为0则由GDAL的缓存大小根据公式：缓存所占的字节数/扫描函数的字节数得到。所以该值不会超出缓存的大小。

"ALL_TOUCHED":

设置为TRUE表示所有的像素接触到矢量中线或多边形，否则只是多边形中心或被brezenhams line algorithm算法（注：brezenhams line algorithm算法为一有名的矢量栅格化算法）选中的部分。

"BURN_VALUE_FROM":

用于设置几何图形的Z值，就是高程值。

"MERGE_ALG":

设置重写或增加新值到栅格数据中。选择REPLACE为 重写，选择ADD为添加一个新值到已存在的栅格数据中。默认值为REPLACE

返回值：CE_None 表示处理成功，CE_Failure表示有错误发生。

下面介绍该函数的一个调用例子，用gdal201实现的，演示了该函数用法，具体有一些地方要根据实际情况做改动。代码如下：

//矢量栅格化
int Rasterization(const char* ShpFile, const char* TifFile)
{
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");//支持中文路径
	GDALAllRegister();
	OGRRegisterAll();
	//打开矢量图层
	GDALDataset* shpData = (GDALDataset*)GDALOpenEx(ShpFile, GDAL_OF_READONLY, NULL, NULL, NULL);
	if (shpData == NULL)
	{
		return FALSE;
	}
	OGRLayer* shpLayer = shpData->GetLayer(0);
	OGREnvelope env;
	shpLayer->GetExtent(&env);//获取图层的坐标范围到env指向的内存中
	int m_nImageWidth = 1440;
	int m_nImageHeight = 720;
	OGRSpatialReference* pOgrSRS = shpLayer->GetSpatialRef();//数据投影信息
	char* pPrj = NULL;
	if (pOgrSRS == NULL)
	{
		cout << "无投影信息...\n";
		m_nImageHeight = (int)env.MaxX;
		m_nImageWidth = (int)env.MaxY;
	}
	else
	{
		pOgrSRS->exportToWkt(&pPrj);
	}

	//获得驱动创建数据集
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* poNewDS = poDriver->Create(TifFile, m_nImageWidth, m_nImageHeight, 1, GDT_Float32, NULL);

	double adfGeoTransform[6];
	adfGeoTransform[0] = env.MinX;//左上角经度
	adfGeoTransform[1] = (env.MaxX - env.MinX) / m_nImageWidth;//像元宽度(影像在宽度上的分辨率)
	adfGeoTransform[2] = 0; //如果影像是指北的, 这个参数的值为0
	adfGeoTransform[3] = env.MaxY;//左上角纬度
	adfGeoTransform[4] = 0; //如果影像是指北的, 这个参数的值为0。
	adfGeoTransform[5] = (env.MinY - env.MaxY) / m_nImageHeight; //像元高度(影像在高度上的分辨率)
	GDALSetGeoTransform(poNewDS, adfGeoTransform);
	if (pOgrSRS != NULL)
	{
		poNewDS->SetProjection(pPrj);
	}

	int* pnbandlist = new int[1];
	pnbandlist[0] = 1;

	OGRLayerH* player;
	player = new OGRLayerH[1];
	player[0] = (OGRLayerH)shpLayer;

	char** papszOptions = NULL;
	papszOptions = CSLSetNameValue(papszOptions, "CHUNKSIZE", "1");
	papszOptions = CSLSetNameValue(papszOptions, "ATTRIBUTE", "DF");//矢量图层属性字段名,我这里是DF

	CPLErr err = GDALRasterizeLayers((GDALDatasetH)poNewDS, 1, pnbandlist, 1, player,
		NULL, NULL, NULL, papszOptions, NULL, NULL);//这里传了很多空值，只简单展示矢量转栅格的功能，详情以后介绍

	GDALClose(shpData);
	GDALClose(poNewDS);
	GDALDestroyDriverManager();
	delete[]player;
	delete[]pnbandlist;
	if (err != CE_None)
	{
		return 0;
	}
	return 1;
}

测试已通过，其中栅格行列数，矢量图层属性字段名，以及一些设置为NULL的参数，需要根据实际情况设置。