IT小白

这个屌丝很懒，什么也没留下！

热门标签

C++下GDAL的详细使用案例（含项目配置、tif读取为cv::Mat、Mat保存为tif）_c++ gdal

作者：IT小白 | 2024-03-28 11:56:22

踩

c++ gdal

1、GDAL的安装与配置

1.1 GDAL的下载安装

直接到下列链接下载即可，按照说明，将bin目录添加的系统环境变量中即可

windows下GDAL322的库-深度学习文档类资源-CSDN下载

1.2 vs中GDAL的配置

库目录中设置lib的路径

附加依赖项中设置gdal_i.lib

2、GDAL读取数据

GDAL读取数钱需要注册一下驱动（用于编码解码图像的驱动），同时可以设置一下支持中文路径。加载数据时需要注意，GA_Update和GA_ReadOnly两种模式。


	GDALAllRegister();//注册所有的驱动
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //设置支持中文路径和文件名
	//1、加载tif数据
	string file_path_name = "test.tif";
	//std::cout << "请输入图片路径：" << std::endl;
	//std::cin >> file_path_name;
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(file_path_name.c_str(), GA_ReadOnly);//GA_Update和GA_ReadOnly两种模式
	if (poDataset == NULL)
	{
		std::cout << "指定的文件不能打开!" << std::endl;
		return 0;
	}

2.1 获取图像的尺寸

这里所有的波段的size都是一样的


    int nImgSizeX = poDataset->GetRasterXSize();
    int nImgSizeY = poDataset->GetRasterYSize();
    std::cout << "ImageX = " << nImgSizeX << ",	ImageY = " << nImgSizeY << std::endl;

2.2 获取图像的通道数


    int bandCount = poDataset->GetRasterCount();
    std::cout << "bandCount = " << bandCount << std::endl;

2.3 获取特定波段

需要注意的是在GDAL中波段数的起始数是1，而非0

    GDALRasterBand* poBand1 = poDataset->GetRasterBand(1);

2.4 获取数据类型

根据波段获取数据类型，而数据中有多个波段。由此可知，每一个波段都可以有不同的数据类型。


    GDALDataType g_type = GDALDataType(poBand1->GetRasterDataType());
	std::cout << "g_type = " << g_type << std::endl;

GDAL中共有12种数据类型，具体如下所示


	typedef enum {
		GDT_Unknown = 0,
		GDT_Byte = 1,
		GDT_UInt16 = 2,
		GDT_Int16 = 3,
		GDT_UInt32 = 4,
		GDT_Int32 = 5,GDT_UInt64,GDT_Int64
		GDT_Float32 = 6,
		GDT_Float64 = 7,
		GDT_CInt16 = 8,
		GDT_CInt32 = 9,GDT_CInt64
		GDT_CFloat32 = 10,
		GDT_CFloat64 = 11,
		GDT_TypeCount = 12
	} GDALDataType;

3、坐标系与投影信息处理

3.1 获取及设置坐标变换系数

坐标变换系数的具体格式为左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率。对应同一区域不同级别的的图像，只有水平分辨率（trans[1]）与竖直分辨率（trans[5]）的值不一样。其默认值为{ 0,1,0,0,0,1 }，即x、y分辨率为1，其他信息为0 。


    double trans[6] = { 0,1,0,0,0,1 };//定义为默认值，即x、y分辨率为1，其他信息为0 
	//具体格式为左上角x坐标， 水平分辨率，旋转参数， 左上角y坐标，旋转参数，竖直分辨率。对应同一区域不同级别的的图像，只有水平分辨率（trans[1]）与竖直分辨率（trans[5]）的值不一样
	CPLErr aaa = poDataset->GetGeoTransform(trans);
	trans[2] = 0.3;//修改x的旋转参数信息
	trans[4] = 0.1;//修改y的旋转参数信息
	//poDataset->SetGeoTransform(trans);//设置坐标变换系数
	std::cout << "trans = " << trans[0] << "," << trans[1] << "," << trans[2] << "," << trans[3] << "," << trans[4] << "," << trans[5] << std::endl;

3.2 像素坐标与投影坐标的相互转换

下面实现了经纬度坐标与像素坐标的相互转换


	double dProjX, dProjY;
	int iCol, iRow;
	iCol = 111;
	iRow = 111;
	ImageRowCol2Projection(trans, iCol, iRow, dProjX, dProjY);
	std::cout << "在trans中，像素坐标=》经纬度：" << iCol << "，" << iRow << "====》" << dProjX << "，" << dProjY << std::endl;
	Projection2ImageRowCol(trans, dProjX, dProjY, iCol, iRow);
	std::cout << "在trans中，经纬度=》像素坐标：" << dProjX << "，" << dProjY << "====》" << iCol << "，" << iRow << std::endl;

两个转换函数的定义如下所示，参考了关于GDAL计算图像坐标的几个问题_IvanLJF的博客-CSDN博客


//计算trans中图片xy点的经纬度信息
//adfGeoTransform的6个参数分别为左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率,一般来说，旋转参数都为0
bool Projection2ImageRowCol(double* adfGeoTransform, double dProjX, double dProjY, int& iCol, int& iRow)
{
	try
	{
		double dTemp = adfGeoTransform[1] * adfGeoTransform[5] - adfGeoTransform[2] * adfGeoTransform[4];
		double dCol = 0.0, dRow = 0.0;
		dCol = (adfGeoTransform[5] * (dProjX - adfGeoTransform[0]) -
			adfGeoTransform[2] * (dProjY - adfGeoTransform[3])) / dTemp + 0.5;
		dRow = (adfGeoTransform[1] * (dProjY - adfGeoTransform[3]) -
			adfGeoTransform[4] * (dProjX - adfGeoTransform[0])) / dTemp + 0.5;
 
		iCol = int(dCol);
		iRow = int(dRow);
		return true;
	}
	catch (...)
	{
		return false;
	}
}
 
bool ImageRowCol2Projection(double* adfGeoTransform, int iCol, int iRow, double& dProjX, double& dProjY)
{
	try
	{
		dProjX = adfGeoTransform[0] + adfGeoTransform[1] * iCol + adfGeoTransform[2] * iRow;
		dProjY = adfGeoTransform[3] + adfGeoTransform[4] * iCol + adfGeoTransform[5] * iRow;
		return true;
	}
	catch (...)
	{
		return false;
	}
}

3.3 获取与设置投影坐标系信息

在GDAL数据坐标系中有WGS84坐标系、CGCS2000坐标系、GCJ02坐标系、BD09坐标系等。具体可分为地心坐标系、投影坐标系、原始坐标系、加密坐标系4类。更多坐标系信息可以参考https://blog.csdn.net/sinat_41310868/article/details/115551276


    std::string projs = poDataset->GetProjectionRef();
	//设置地理坐标系信息
	//poDataset->SetProjection(projs.c_str());
	std::cout << "projs = " << projs << std::endl;

4、读写GDAL数据

4.1 将数据读取到Mat中

需要注意的是cv::Mat创建时，是(height,width)的格式，与GDAL的(width,height)刚好相反。同时GADL的起始通道数是1，不是0.


	cv::Mat gdal_mat1(nImgSizeY, nImgSizeX, CV_8UC1, Scalar(0));
	cv::Mat gdal_mat2(nImgSizeY, nImgSizeX, CV_8UC1, Scalar(0));
	cv::Mat gdal_mat3(nImgSizeY, nImgSizeX, CV_8UC1, Scalar(0));
	poDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, gdal_mat1.data, nImgSizeX, nImgSizeY, g_type, 0, 0);
	poDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, gdal_mat2.data, nImgSizeX, nImgSizeY, g_type, 0, 0);
	poDataset->GetRasterBand(3)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, gdal_mat3.data, nImgSizeX, nImgSizeY, g_type, 0, 0);
	cv::Mat mg;
	cv::merge(vector<cv::Mat>{ gdal_mat3, gdal_mat2, gdal_mat1,}, mg);
	cv::imwrite("read_save.jpg", mg);

RasterIO函数的参数列表如下所示。从参数列表中是可以看到，GDAL是支持将数据分块读入的内存中的。RasterIO参数列表的详细说明可以参考 https://blog.51cto.com/u_15469043/4903358


    CPLErr GDALRasterBand::RasterIO (   GDALRWFlag eRWFlag,
										int     nXOff,//x的起始点
										int     nYOff,//y的起始点
										int     nXSize,//读取窗口的宽
										int     nYSize,//读取窗口的高
										void * pData,
										int     nBufXSize,//与nXSize相同
										int     nBufYSize,//与nYSize相同
										GDALDataType    eBufType,
										int     nPixelSpace,//通常默认为0
										int     nLineSpace //通常默认为0
										)

除了将数据读取到mat中外，我们还可以将数据读取到指针或者数组中


    //void * malloc(size_t n)：给指针分配相应的内存，并返回内存空间的首地址。当内存不再使用的时候，应使用free()函数将内存块释放掉。
	uint8_t* srcData = (uint8_t*)malloc(sizeof(uint8_t) * nImgSizeX * nImgSizeY);
	//void * memset (void * p,int c,size_t n)：将p中的n个字节都赋值为c
	memset(srcData, 0, sizeof(uint8_t) * 1 * nImgSizeX * nImgSizeY);//为空间赋默认值0
	poDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, srcData, nImgSizeX, nImgSizeY, g_type, 0, 0);

4.2 将mat数据保存到GDAL中（tiff格式）

创建gdal对象一次性写入数据，只支持tiff数据。PEN、JPEG等驱动没有实现相应的Create方法。


    int nImgSizeX3 = gdal_mat1.cols;
	int nImgSizeY3 = gdal_mat1.rows;
	GDALDriver* pDriverMEM3 = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (!pDriverMEM3) {
		fprintf(stderr, "get driver by name failed\n");
		return -1;
	}
	int nBands3 = 1;
	GDALDataset* poDataset3 = pDriverMEM3->Create("saved3.tif", nImgSizeX3, nImgSizeY3, nBands3, g_type, NULL);
	if (!poDataset3) {
		fprintf(stderr, "Create GDALDataset failed\n");
		return -1;
	}
	poDataset3->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, nImgSizeX3, nImgSizeY3, gdal_mat1.data, nImgSizeX3, nImgSizeY3, GDT_Byte, 0, 0);
	std::cout << "tif 文件保存成功" << std::endl;

4.3 将mat数据保存到GDAL中（任意格式）

GDAL中可用的驱动格式还有：BMP、JPEG、PNG、GTiff、GIF、HFA、BT、ECW、FITS、HDF4、EHdr。这些驱动格式分别对应着不同的文件类型，MEM表示为内存对象，可以快速的分块追加写入数据。MEM文件大小是和你的系统内存大小有关系，并不会存储到磁盘中。


    int nImgSizeX2 = gdal_mat1.cols;
	int nImgSizeY2 = gdal_mat1.rows;
	//获取GDAL驱动
	GDALDriver* pDriverMEM = GetGDALDriverManager()->GetDriverByName("MEM");
	int nBands = 1;
	//创建GDAL对象，只保存原图的一个通道
	//Create(const char * pszName,int nXSize, int nYSize, int nBands, GDALDataType eType, char** papszOptions)
	GDALDataset* poDataset2 = pDriverMEM->Create("", nImgSizeX2, nImgSizeY2, nBands, g_type, NULL);
	//将mat数据写入到GDALDataset中
	poDataset2->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, nImgSizeX2, nImgSizeY2, gdal_mat1.data, nImgSizeX2, nImgSizeY2, GDT_Byte, 0, 0);
	//获取GDAL驱动，PNG表示为用png驱动保存数据
	GDALDriver* pDriverSave = GetGDALDriverManager()->GetDriverByName("PNG");
	pDriverSave->CreateCopy("saved.png", poDataset2, TRUE, 0, 0, 0); //创建png文件
	std::cout << "png 文件保存成功" << std::endl;

4.4 关闭GDAL对象，并注销所有驱动


	GDALClose(poDataset);
	GDALClose(poDataset2);
	GDALClose(poDataset3);
	GDALDestroyDriverManager();

5、代码及效果

5.1 全部代码


#include <iostream>
#include <memory>
#include <chrono>
#include <fstream>
#include <string>
#include <iomanip>
#include <opencv2/highgui.hpp>
#include <opencv2/opencv.hpp>
#include "gdal_priv.h"
#include <gdal_alg_priv.h>
#include <gdal.h>
using namespace std;
using namespace cv;
//参考https://blog.csdn.net/ivan_ljf/article/details/9226463
//计算trans中图片xy点的经纬度信息
//adfGeoTransform的6个参数分别为左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率,一般来说，旋转参数都为0
bool Projection2ImageRowCol(double* adfGeoTransform, double dProjX, double dProjY, int& iCol, int& iRow)
{
	try
	{
		double dTemp = adfGeoTransform[1] * adfGeoTransform[5] - adfGeoTransform[2] * adfGeoTransform[4];
		double dCol = 0.0, dRow = 0.0;
		dCol = (adfGeoTransform[5] * (dProjX - adfGeoTransform[0]) -
			adfGeoTransform[2] * (dProjY - adfGeoTransform[3])) / dTemp + 0.5;
		dRow = (adfGeoTransform[1] * (dProjY - adfGeoTransform[3]) -
			adfGeoTransform[4] * (dProjX - adfGeoTransform[0])) / dTemp + 0.5;
 
		iCol = int(dCol);
		iRow = int(dRow);
		return true;
	}
	catch (...)
	{
		return false;
	}
}
 
bool ImageRowCol2Projection(double* adfGeoTransform, int iCol, int iRow, double& dProjX, double& dProjY)
{
	try
	{
		dProjX = adfGeoTransform[0] + adfGeoTransform[1] * iCol + adfGeoTransform[2] * iRow;
		dProjY = adfGeoTransform[3] + adfGeoTransform[4] * iCol + adfGeoTransform[5] * iRow;
		return true;
	}
	catch (...)
	{
		return false;
	}
}
int main() {
	GDALAllRegister();//注册所有的驱动
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //设置支持中文路径和文件名
	//1、加载tif数据
	string file_path_name = "test.tif";
	//std::cout << "请输入图片路径：" << std::endl;
	//std::cin >> file_path_name;
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(file_path_name.c_str(), GA_ReadOnly);//GA_Update和GA_ReadOnly两种模式
	if (poDataset == NULL)
	{
		std::cout << "指定的文件不能打开!" << std::endl;
		return 0;
	}
 
	//获取图像的尺寸
	int nImgSizeX = poDataset->GetRasterXSize();
	int nImgSizeY = poDataset->GetRasterYSize();
	std::cout << "ImageX = " << nImgSizeX << ",	ImageY = " << nImgSizeY << std::endl;
 
	//获取图像的通道数（波段数量）
	int bandCount = poDataset->GetRasterCount();
	std::cout << "bandCount = " << bandCount << std::endl;
 
	//获取图像波段  在GDAL中波段数的起始数是1，而非0
	GDALRasterBand* poBand1 = poDataset->GetRasterBand(1);
 
	//GDAL中的数据类型 由此可知，每一个波段都可以有不同的数据类型
	/* 一共包含以下12种数据类型
	typedef enum {
		GDT_Unknown = 0,
		GDT_Byte = 1,
		GDT_UInt16 = 2,
		GDT_Int16 = 3,
		GDT_UInt32 = 4,
		GDT_Int32 = 5,GDT_UInt64,GDT_Int64
		GDT_Float32 = 6,
		GDT_Float64 = 7,
		GDT_CInt16 = 8,
		GDT_CInt32 = 9,GDT_CInt64
		GDT_CFloat32 = 10,
		GDT_CFloat64 = 11,
		GDT_TypeCount = 12
	} GDALDataType;
	*/
	GDALDataType g_type = GDALDataType(poBand1->GetRasterDataType());
	std::cout << "g_type = " << g_type << std::endl;
 
 
	//获取坐标变换系数
	double trans[6] = { 0,1,0,0,0,1 };//定义为默认值，即x、y分辨率为1，其他信息为0 
	CPLErr aaa = poDataset->GetGeoTransform(trans);
	trans[2] = 0.3;//修改x的旋转参数信息
	trans[4] = 0.1;//修改y的旋转参数信息
	//poDataset->SetGeoTransform(trans);//设置坐标变换系数
	std::cout << "trans = " << trans[0] << "," << trans[1] << "," << trans[2] << "," << trans[3] << "," << trans[4] << "," << trans[5] << std::endl;
 
	//像素坐标与投影坐标的换算
	double dProjX, dProjY;
	int iCol, iRow;
	iCol = 111;
	iRow = 111;
	ImageRowCol2Projection(trans, iCol, iRow, dProjX, dProjY);
	std::cout << "在trans中，像素坐标=》经纬度：" << iCol << "，" << iRow << "====》" << dProjX << "，" << dProjY << std::endl;
	Projection2ImageRowCol(trans, dProjX, dProjY, iCol, iRow);
	std::cout << "在trans中，经纬度=》像素坐标：" << dProjX << "，" << dProjY << "====》" << iCol << "，" << iRow << std::endl;
 
 
	//获取图像投影坐标系信息，
	std::string projs = poDataset->GetProjectionRef();
	//设置地理坐标系信息
	//poDataset->SetProjection(projs.c_str());
	std::cout << "projs = " << projs << std::endl;
 
 
	//读取gadl中第一个通道的数据到mat中 【通道数是从1开始的】
	//RasterIO参数列表的详细说明可以参考 https://blog.51cto.com/u_15469043/4903358
	/*从参数列表中是可以看到，GDAL是支持将数据分块读入的内存中的
	CPLErr GDALRasterBand::RasterIO (   GDALRWFlag eRWFlag,
										int     nXOff,//x的起始点
										int     nYOff,//y的起始点
										int     nXSize,//读取窗口的宽
										int     nYSize,//读取窗口的高
										void * pData,
										int     nBufXSize,//与nXSize相同
										int     nBufYSize,//与nYSize相同
										GDALDataType    eBufType,
										int     nPixelSpace,//通常默认为0
										int     nLineSpace //通常默认为0
										)
	*/
	//
	cv::Mat gdal_mat1(nImgSizeY, nImgSizeX, CV_8UC1, Scalar(0));
	cv::Mat gdal_mat2(nImgSizeY, nImgSizeX, CV_8UC1, Scalar(0));
	cv::Mat gdal_mat3(nImgSizeY, nImgSizeX, CV_8UC1, Scalar(0));
	poDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, gdal_mat1.data, nImgSizeX, nImgSizeY, g_type, 0, 0);
	poDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, gdal_mat2.data, nImgSizeX, nImgSizeY, g_type, 0, 0);
	poDataset->GetRasterBand(3)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, gdal_mat3.data, nImgSizeX, nImgSizeY, g_type, 0, 0);
	cv::Mat mg;
	cv::merge(vector<cv::Mat>{ gdal_mat3, gdal_mat2, gdal_mat1,}, mg);
	cv::imwrite("read_save.jpg", mg);
	/*
	//读取gadl中第一个通道的数据到指针中
	//void * malloc(size_t n)：给指针分配相应的内存，并返回内存空间的首地址。当内存不再使用的时候，应使用free()函数将内存块释放掉。
	uint8_t* srcData = (uint8_t*)malloc(sizeof(uint8_t) * nImgSizeX * nImgSizeY);
	//void * memset (void * p,int c,size_t n)：将p中的n个字节都赋值为c
	memset(srcData, 0, sizeof(uint8_t) * 1 * nImgSizeX * nImgSizeY);//为空间赋默认值0
	poDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, srcData, nImgSizeX, nImgSizeY, g_type, 0, 0);
	*/
 
	//-----------创建gdal对象   MEM追加，CreateCopy保存支持tif、png、jpg等格式-----
	int nImgSizeX2 = gdal_mat1.cols;
	int nImgSizeY2 = gdal_mat1.rows;
	//获取GDAL驱动，MEM表示为内存对象，可以快速的分块追加写入数据。MEM文件大小是和你的系统内存大小有关系，并不会存储到磁盘中。可用的驱动格式还有：BMP、JPEG、PNG、GTiff、GIF、HFA、BT、ECW、FITS、HDF4、EHdr。分别对应着不同的文件类型
	GDALDriver* pDriverMEM = GetGDALDriverManager()->GetDriverByName("MEM");
	int nBands = 1;
	//创建GDAL对象，只保存原图的一个通道
	//Create(const char * pszName,int nXSize, int nYSize, int nBands, GDALDataType eType, char** papszOptions)
	GDALDataset* poDataset2 = pDriverMEM->Create("", nImgSizeX2, nImgSizeY2, nBands, g_type, NULL);
	//将mat数据写入到GDALDataset中
	poDataset2->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, nImgSizeX2, nImgSizeY2, gdal_mat1.data, nImgSizeX2, nImgSizeY2, GDT_Byte, 0, 0);
	//获取GDAL驱动，PNG表示为用png驱动保存数据
	GDALDriver* pDriverSave = GetGDALDriverManager()->GetDriverByName("PNG");
	pDriverSave->CreateCopy("saved.png", poDataset2, TRUE, 0, 0, 0); //创建png文件
	std::cout << "png 文件保存成功" << std::endl;
 
	//-----------创建gdal对象   一次性写入，只支持tiff数据。PEN、JPEG等驱动没有实现相应的Create方法-----
	int nImgSizeX3 = gdal_mat1.cols;
	int nImgSizeY3 = gdal_mat1.rows;
	GDALDriver* pDriverMEM3 = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (!pDriverMEM3) {
		fprintf(stderr, "get driver by name failed\n");
		return -1;
	}
	int nBands3 = 1;
	GDALDataset* poDataset3 = pDriverMEM3->Create("saved3.tif", nImgSizeX3, nImgSizeY3, nBands3, g_type, NULL);
	if (!poDataset3) {
		fprintf(stderr, "Create GDALDataset failed\n");
		return -1;
	}
	poDataset3->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, nImgSizeX3, nImgSizeY3, gdal_mat1.data, nImgSizeX3, nImgSizeY3, GDT_Byte, 0, 0);
	std::cout << "tif 文件保存成功" << std::endl;
 
	//关闭GDAL对象，并注销所有驱动
	GDALClose(poDataset);
	GDALClose(poDataset2);
	GDALClose(poDataset3);
	GDALDestroyDriverManager();
 
	return -1;
}

5.2 运行结果

运行界面输出

运行过程中的测试数据及保存结果

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/IT小白/article/detail/329869