利用GDAL库实现点云转栅格（GeoTif格式）

一、主要借助于GDAL库中的GDALRasterBand类的RasterIO成员函数来操作

CPLErr GDALRasterBand::RasterIO (   
GDALRWFlag eRWFlag,
int     nXOff,
int     nYOff,
int     nXSize,
int     nYSize,
void * pData,
int     nBufXSize,
int     nBufYSize,
GDALDataType    eBufType,
int     nPixelSpace,
int     nLineSpace 
)

接口中各个参数意义详见CSDN博客：
https://blog.csdn.net/liminlu0314/article/details/7072224

（1）考虑到GeoTif图像一般存在多个波段，如下图所示，根据点云RGB生成的GeoTif就存在三个波段，分别对应于R、G、B。

（2）设定像元步长step，利用点云在XOY平面上投影点的范围计算图像的总行数与总列数：

rows=int((Ymax-Ymin)/step);
cols=int((Xmax-Xmin)/step);

（3）在C++中建议使用二级指针来存储图像的像素值，image[波段号][像元索引号]=像元值，注意语法。

int band=3;  //波段数
//二级指针申请内存空间用于存储像元值
uint8_t** image=new uint8_t[band];  
for (int k = 0; k < band; ++k)
{
	image[k] = new uint8_t[rows * cols];  
}

//二级指针内存空间释放方法
if(image)
{
for (int i = 0; i < channel; ++i)
{
	delete[] imageData[i];
}
delete[] image;
image = nullptr;
}

（4）图像像元值写入

for(int r=0;r<rows;++r)
{
	for(int c=0;c<cols;++c)
	{
		for(int b=0;b<band;++b)
		{
			pixelValue[b]=XXX;  //设置图像第b+1个波段中位于[r,c]位置处的像元值XXX
			image[b][r * cols + c] = pixelValue[b];//从图像左下角开始依次写入像元值
		}
	}
}

（5）调用GDAL库的RasterIO接口生成GeoTif图像

bool createGeoTif(std::string savePath, int width, int height, int bandNums, uint8_t** pImage, double xOff, double yOff,double step) 
{
GDALAllRegister();  //GDAL注册驱动
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //允许使用中文字符

GDALDriver *pMemDriver = NULL;
pMemDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (pMemDriver == NULL)
{
	return false;
}


GDALDataset * pMemDataSet = pMemDriver->Create(savePath.c_str(), width, height, bandNums, GDT_Byte, NULL);
if (pMemDataSet == NULL)
{
	return  false;
}

//设置仿射变换参数
//geoTransfoem[0]：图像原点X坐标
//geoTransfoem[1]：图像在东西方向上分辨率（点云X轴）
//geoTransfoem[2]：旋转角，0表示图像北方向朝上
//geoTransfoem[3]：图像原点Y坐标
//geoTransfoem[4]：旋转角，0表示图像北方向朝上
//geoTransfoem[5]：图像在南北方向上分辨率（点云Y轴）
double geoTransfoem[6] = { xOff,step,0,yOff,0,step};
pMemDataSet->SetGeoTransform(geoTransfoem);
GDALRasterBand *pBand = NULL;


for (int i = 1; i <= bandNums; i++)
{
	pBand = pMemDataSet->GetRasterBand(i);
	CPLErr err = pBand->RasterIO(GF_Write, 0, 0, width, height, pImage[i - 1], width, height, GDT_Byte, 0, 0);
	if (err == CE_Failure)
	{
		return false;
	}
}

//关闭驱动
GDALClose(pMemDataSet);
GetGDALDriverManager()->DeregisterDriver(pMemDriver);
return true;
}

二、在点云转栅格过程中遇到的几个问题：

问题1：当像元值为浮点型数据时，GDAL所生成的GeoTif图像的像元为纯黑色或者纯白色偶尔会出现灰色。
原因：GDAL库在生成浮点型geoTif时所能接受的像元值大小必须在[0,1]区间内。若像元值>1，则强制变成1；若像元值<0，则强制变成0；若像元值>0且<1，则不改变其值。当时是根据点云Z坐标来生的GeoTif，测区平均高程在120米左右，导致所有非空像元的值都被强制变成了1，得到大面积黑色图像。因此必须提前对二级指针image中存储的像元值做归一化处理，且必须归一化至[0,1]区间内。归一化公式如下：

image[b][r * cols + c] = (pixelValue-pixelMinValue)/(pixelMaxValue-pixelMinValue);

此外，当你是根据点云的Z值来生成栅格时，像元值必然是浮点型，RasterIO的第9个参数不再是GDT_Byte而应该改成GDT_Float32，如果坚持使用前者的话虽然不会报错但输出结果显然是错的，数据被截断。
GDT_Byte等价于unsigned char，取值范围是[0,255]
GDT_Float32等价于float，取值范围是[0,1]

问题2：将生成的geoTif图像与原始点云同时在CloudCompare中打开，发现两者无法重叠在一起，如下图所示，位于下方黄色包围盒内的是原始点云，位于上分的是生成的GeoTif。

原因：图像坐标系原点位于左上角，点云坐标系原点位于左下角，但是在生GeoTif图像时，图像坐标系原点必须确保与点云坐标系原点保持一致，否则GeoTif图像必然发生偏移导致无法与点云重叠。
此处一开始设置的GeoTif图像原点为（Xmin，Ymax）（黄色包围盒左上角点），然而点云坐标系原点为（Xmin，Ymin）（黄色包围盒左下角点），加载到CC里面后就产生了上述偏移现象。改变图像坐标原点之后两者即可重叠
在向二级指针image中写入像元值时，无论是从左上角点开始写入还是从左下角点开始写入甚至从右下角点开始写入其实都可以，只要像元值与像元的位置相匹配就不会出错！个人习惯还是从左下角点开始写如比较舒服，不容易出错。

三、代码验证
使用CloudCompare设置相同的参数也生成一张geoTif进行对比。
原始点云：

CloudCompare软件生成的GeoTif：

两者可以重叠在一起：

代码生成的GeoTif：

代码生成的GeoTif与原始点云可以重叠：

注意：直接用windows自带的图片查看器打开GeoTif图像时发现图像变成了关于原始点云的镜像，但是加载到CloudCompare里面是可以和原始点云重叠的。不用担心，这是正常现象，因为windows图片查看器定义的原点在左上角，它会把位于左下角（Xmin，Ymin）的像元转换到图像左上角，依次类推，所有像元的行列号都会发生上述转换。同样的，如果打开CloudCompare生成的GeoTif会发现也变成了关于原始点云的镜像图。因此，还是要以CloudCompare为准。

四、展望
（1）当点云区域太大时，new运算申请堆内存可能会失败，因此后续需要研究如何对点云进行分块处理，分块生成对应GeoTif然后做图像拼接生成完整的GeoTif。
（2）当某个像元内不存在点云投影时，应为空像元，但目前还未找到合适的值来填充到空像元中。