利用gdal读写DEM数据（.tif格式）

目标：

目标：读取DEM数据（.tif格式），控制台打印每一个像素值的高程值，同时输出该DEM数据（输出名为：out.tif）
后期目标：需要对读入的dem进行像素值的一些统计分析，并修改部分像素值，再输出，这里只进行读写操作，不修改像素值。

学习文件：

gdal官方说明网址（API）：https://gdal.org/api/index.html
学习API参数含义；重要的API：RasterIO、GetGeoTransform

一、步骤

1.主函数

main代码：

int main()
{
	tiffread("dem.tif");
	system("pause");
	return 0;
}

2.读数据和写函数

代码如下 ：

int tiffread(const char* file_path_name) //读数据并输出保存
{
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");	// 支持中文路径

	GDALAllRegister();  //注册所有的驱动 

	GDALDataset *poDataset;   //GDAL数据集 

	poDataset = (GDALDataset *)GDALOpen(file_path_name, GA_ReadOnly);
	if (poDataset == NULL)
	{
		cout << "fail in open files!!!" << endl;
		return 0;
	}

	//获取图像的尺寸 
	int nImgSizeX = poDataset->GetRasterXSize();
	int nImgSizeY = poDataset->GetRasterYSize();
	cout << nImgSizeX << "  *  " << nImgSizeY << endl;


	//获取坐标变换系数 
	double trans[6];
	CPLErr aaa = poDataset->GetGeoTransform(trans);
	cout << "坐标变换系数:";
	for (auto tran : trans)
		cout << tran << " ";
	cout << endl;


	//获取图像波段 
	GDALRasterBand *poBand;//*poBand1, *poBand2
	int bandcount;
	bandcount = poDataset->GetRasterCount();	// 获取波段数:DEM为1个波段
	cout << "波段数：" << bandcount << endl;
    
	
	BYTE *pafScanline;          //DEM数据为单波段数据 该代码只适用于单通道数据读取和写
	//pafScanline = new  BYTE [bandcount];
	//pafScanline = new BYTE[nImgSizeX*nImgSizeY];
	//for (int bandind = 0; bandind < bandcount; bandind++)
	//{
		pafScanline = new BYTE[nImgSizeX*nImgSizeY];     //pafScanline为int型数组指针 用于存储dem数据
	//}
	//for(int bandind = 1; bandind <= bandcount; bandind++)
	//{
		poBand = poDataset->GetRasterBand(bandcount);		// 获取对应波段

		//读取图像高程数据 
		int num_iamge_size = 0;		//数据计数
		
		poBand->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, pafScanline, nImgSizeX, nImgSizeY, /*GDALDataType(poBand->GetRasterDataType())*/GDT_Float32, 0, 0);
		for (int i = 0; i < nImgSizeX; i++)
		{
			for (int j = 0; j < nImgSizeY; j++)
			{
				num_iamge_size++;
				cout << pafScanline[i*nImgSizeY + j] <<" ";
			}
			cout << endl;
			//cout << i << "---" << pafScanline[i*nImgSizeY] << endl;
		}
		
//以下为输出dem数据
    char text[]= "out.tif";
	char *outFilePath = text;

	GDALDataType type = poBand->GetRasterDataType();
	GDALDataset *ods = createRaster(outFilePath, nImgSizeX, nImgSizeY, type);//创建空的栅格数据集
	ods->SetProjection(poDataset->GetProjectionRef());//给它设置投影
	ods->SetGeoTransform(trans);//给设置空间转换的六参数

	GDALRasterBand *oBand = ods->GetRasterBand(1);
	oBand->RasterIO(      
		GF_Write,
		0, 0,
		nImgSizeX,
		nImgSizeY,
		pafScanline,        //写入数据
		nImgSizeX,
		nImgSizeY,
		type,
		0,
		0
	);

3.创建栅格数据函数

代码如下 ：

GDALDataset * createRaster(char *filename, int nRow, int nCol, GDALDataType type)
{
	GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset *ds = driver->Create(
		filename,
		nRow,
		nCol,
		1,//波段数
		type,
		NULL);
	return ds;  //传回一个空的栅格数据集，但设置了它的一些基本参数
}

总结

结果：生成的out.tif文件，数值和输入dem一致，具有空间参考：

体会：代码主要是综合了下B站那位up主写的和参考网站1的代码，还是对用gdal库读取数据，以及该库的API的有了比较深刻的印象，后续继续对DEM数据进行一些统计分析，完善该程序，最好能可视化出该dem。

参考文章：

参考文章或网址
1、link.
2、link.
3、link.
4.bilibi：https://www.bilibili.com/video/BV1o7411i7pe?from=search&seid=12940848749693557054
原文链接