Gdal总结

珞珈Lena

已于 2023-03-18 18:10:17 修改

阅读量2.3k

点赞数 5

分类专栏： Gdal库文章标签： c++ 开发语言

于 2022-10-08 09:43:34 首次发布

本文链接：https://blog.csdn.net/Zm1366172479/article/details/127200434

版权

Gdal库专栏收录该内容

2 篇文章 0 订阅

订阅专栏

1、头文件

#include<gdal.h>
#include<gdal_priv.h>
#include<gdalwarper.h>

2、注册驱动+设置中文路径+加载数据

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


    //设置支持中文路径和文件名
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); 
	


    //加载tif数据
	string file_path_name = "test.tif";
	//std::cout << "请输入图片路径：" << std::endl;
	//std::cin >> file_path_name;


    //将数据读入poDataset数据集
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(file_path_name.c_str(),GA_ReadOnly);

	if (poDataset == NULL)
	{
		std::cout << "指定的文件不能打开!" << std::endl;
		return 0;
	}

注意：GDALOpen中参数分为GA_Update（写数据），GA_ReadOnly（仅仅读数据）

3、获取图像的尺寸

  int nImgSizeX = poDataset->GetRasterXSize();
  int nImgSizeY = poDataset->GetRasterYSize();

注意：此处得Xsize指的是有多少列（即宽度），Ysize指的是有多少行（即高度）。如下图

4、获取图像的通道数

    int bandCount = poDataset->GetRasterCount();

5、获取特定波段

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

注意：此处得波段数是从1计数，不是从0；

6、获取数据类型

    GDALDataType g_type = GDALDataType(poBand1->GetRasterDataType());

注意：Gdal总共有12中数据类型，具体如下

GDT_Unknown : 未知数据类型
GDT_Byte : 8bit正整型 (C++中对应unsigned char)
GDT_UInt16 : 16bit正整型 (C++中对应 unsigned short)
GDT_Int16 : 16bit整型 (C++中对应 short 或 short int)
GDT_UInt32 : 32bit 正整型 (C++中对应unsigned long)
GDT_Int32 : 32bit整型 (C++中对应int 或 long 或 long int)
GDT_Float32 : 32bit 浮点型 (C++中对应float)
GDT_Float64 : 64bit 浮点型 (C++中对应double)
GDT_CInt16 : 16bit复整型 (?)
GDT_CInt32 : 32bit复整型 (?)
GDT_CFloat32 : 32bit复浮点型 (?)
GDT_CFloat64 : 64bit复浮点型 (?)

7、读写GDAL数据

（1）读取Gdal数据

poDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, data, nImgSizeX, nImgSizeY, g_type, 0, 0);

参数解释：

1、GF_Read：代表是读取tif文件数据。

2、3：代表初始位置（a，b），注意第二个参数a指的是X方向，即多少列。第三个参数b指的是Y方向，即多少行。原点在左上角。如下图所示。

4，5指的是偏移量（X,Y）。X代表距离参数（a，b）的水平平移量和竖直偏移量。即宽和高。

如下图所示

6、中data用来将读取到的影像数据进行存储。想到于读到缓存中。data是一个数组指针类型，数据类型需要和GDAL通过GetRasterDataType获取的数据类型一致。创建格式如下

数组的大小要能存放下计划读入缓存的数据。如果是读取图像，那么读取出来的图像像素值就存储在这个Data中，如果是写入图像，那么这个Data中的数据会被写入到图像上指定的位置中去

unsigned short* data_ref = new unsigned short[100 * nRows];

7、8 表示data数据的大小。比如4、5设置读取图像的范围是300×200，如果我要读取图像中的原始数据，那么这两个参数分别应该设置为nImgSizeX=300和nImgSizeY = 200。如果设置比这个小或者大就会进行重采样生成放大或者缩小的图像。

9、gtype标记data 的数据类型。此时用的是GDAL格式的数据类型。

10、11默认0即可。具体不再阐述。

（2）写数据

写数据需要在加载数据的时候改为GA_Update，然后用以下代码写入

poDataset->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, nImgSizeX, nImgSizeY, data, nImgSizeX, nImgSizeY, g_type, 0, 0);

指的是将data中的数据写入打开的tiff文件指定位置

注意：data中的数据是一行，和影像对应的关系如下图

（3）GDAL与OPENCV联动

注意：cv::Mat创建时，是(height,width)的格式，与GDAL的(width,height)刚好相反。

如以下代码所示

pDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, ncols, nrows, data, ncols, nrows, dataType, 0, 0);
Mat temp= (Mat(nrows, ncols, dataType, data)).clone();

也可以直接读取数据到mat中

cv::Mat gdal_mat1(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);

8、Gdal复制影像

//加载原始影像数据
		GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(file_path_name.c_str(), GA_ReadOnly);
		//创建要存储的数据集
		GDALDataset* pDatasetc;
		//注册create的驱动
		GDALDriver* poDriver;
		//设置驱动格式
		const char* pszFormat = "GTiff";
		//接下来的两步是完成复制，file_path_nameout是输出位置+文件名，poSrcDS是原始影像数据集
		poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
		pDatasetc = poDriver->CreateCopy(file_path_nameout.c_str(), poSrcDS, FALSE, NULL, NULL, NULL);
		//判断复制后的影像是不是空的，是不是都关闭数据集。
		if (pDatasetc != NULL) {
			GDALClose((GDALDatasetH)pDatasetc);
		}
		GDALClose((GDALDatasetH)poSrcDS);

注意：这段程序只需要修改输入文件：file_path_name，和输出位置文件：file_path_nameout，即可直接复制运行。

9、Gdal地理信息

影像左上角横坐标（经度）：geoTransform[0]
影像左上角纵坐标（纬度）：geoTransform[3]

遥感图像的水平空间分辨率为geoTransform[1]
遥感图像的垂直空间分辨率为geoTransform[5]
通常geoTransform[5] 与 geoTransform[1]相等

double dfTransform[6] = { 0 };//获取影像六参数
readDataset->GetGeoTransform(dfTransform);
cout << dfTransform[0] << endl;//获取起始经度
cout << dfTransform[3] << endl;//获取起始纬度

10、获取投影

	char* Project = NULL;
	Project = const_cast<char*>(readDataset->GetProjectionRef());

11、获取影像转换参数

	double dfTransform[6] = { 0 };
	readDataset->GetGeoTransform(dfTransform);

12、创建数据集并打开

	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* pDDst = pDriver->Create(filepathChong.c_str(), nNewWidth, nNewHeight, nBandCount,g_type, NULL);
	if (pDDst == NULL)
	{
		GDALClose((GDALDatasetH)readDataset);
		cout << "创建重采样文件失败" << endl;
		return -2;
	}

13、对创建的影像设置投影和六参数

	pDDst->SetProjection(Project);//设置投影
	pDDst->SetGeoTransform(dfTransformnew);//设置6参数

14、遥感影像重采样

#include<iostream>
#include<string>
using namespace std;
#include<gdal.h>
#include<gdal_priv.h>
#include<gdalwarper.h>
#define PI 3.1415926535897932384626433832795


int main(int argc, char* argv[])
{
	//读取设置参数
	GDALResampleAlg eResample; 
	int method;
	if (argc < 4)
	{
		method = 2;
	}
	else
	{
	sscanf(argv[3], "%d", &method);
	}

	
	
	switch (method)
	{
	case 1:
		eResample = GRA_NearestNeighbour;//最邻近采样 速度最快
		break;
	case 2:
		eResample = GRA_Bilinear; //双线性内插采样
		break;
	case 3:
		eResample = GRA_Cubic;//三次立方卷积采样
		break;
	}
	

	double newscale;
	sscanf(argv[2], "%lf", &newscale);

	FileInfo file(argv[1]);
	string filepath = file.path +"/"+ file.name +"."+ file.format;
	string filepathChong = file.path + "/" + file.name + "-Resampling"+"." + file.format;


	//注册所有的驱动
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //设置支持中文路径和文件名



	GDALDataset* readDataset = (GDALDataset*)GDALOpen(filepath.c_str(), GA_ReadOnly);//加载tif原图
	if (readDataset == NULL)
	{
		std::cout << "指定的文件不能打开!" << std::endl;
		return 0;
	}

	//获取影像分辨率
	double dfTransform[6] = { 0 };
	readDataset->GetGeoTransform(dfTransform);
	double resolution_lat = -dfTransform[5];
	double resolution_lon = dfTransform[1];
	double resolutiony = Latitude_difference_to_distance(resolution_lat);
	double resolutionx = Longitude_difference_to_distance(resolution_lon, dfTransform[3]);
	//double resolution = std::max(resolutionx, resolutiony);

	//获取影像大小
	int nImgSizeX = readDataset->GetRasterXSize();
	int nImgSizeY = readDataset->GetRasterYSize();

	//获取投影
	char* Project = NULL;
	Project = const_cast<char*>(readDataset->GetProjectionRef());

	//获取影像数据类型
	GDALDataType g_type = GDALDataType(readDataset->GetRasterBand(1)->GetRasterDataType());

	//获取通道数
	int nBandCount = readDataset->GetRasterCount();


	double dfTransformnew[6] = { 0 };
	dfTransformnew[0]= dfTransform[0];
	dfTransformnew[1]= distance_to_latitude_difference(newscale);
	dfTransformnew[2]= dfTransform[2];
	dfTransformnew[3]= dfTransform[3];
	dfTransformnew[4]= dfTransform[4];
	dfTransformnew[5]= -distance_to_longitude_difference(newscale, dfTransform[3]);
	
	double fResX = dfTransform[1] / dfTransformnew[1];//计算重采样前后的比例
	double fResY = dfTransform[5] / dfTransformnew[5];

	cout << dfTransform[0] << endl;
	int nNewWidth = static_cast<int>(nImgSizeX * fResX + 0.5);//计算重采样后的宽度
	int nNewHeight = static_cast<int>(nImgSizeY * fResY + 0.5);//计算重采样后的高度

	//创建结果数据集  
	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* pDDst = pDriver->Create(filepathChong.c_str(), nNewWidth, nNewHeight, nBandCount,g_type, NULL);
	if (pDDst == NULL)
	{
		GDALClose((GDALDatasetH)readDataset);
		cout << "创建重采样文件失败" << endl;
		return -2;
	}



	pDDst->SetProjection(Project);//设置投影
	pDDst->SetGeoTransform(dfTransformnew);//设置6参数


	//重采样
	GDALWarpOptions* psWo = GDALCreateWarpOptions();
	
	psWo->eResampleAlg = eResample;//设置重采样方式
	psWo->eWorkingDataType = g_type;//设置数据格式

	psWo->hSrcDS = (GDALDatasetH)readDataset;//源数据
	psWo->hDstDS = (GDALDatasetH)pDDst;//设置目标数据

	psWo->pfnTransformer = GDALGenImgProjTransform;
	psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH)readDataset, Project, (GDALDatasetH)pDDst, Project, FALSE, 0.0, 1);;

	psWo->nBandCount = nBandCount;
	psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));

	for (int i = 0; i < nBandCount; i++)
	{
		psWo->panSrcBands[i] = i + 1;
		psWo->panDstBands[i] = i + 1;
	}

	GDALWarpOperation oWo;
	if (oWo.Initialize(psWo) != CE_None)
	{
		GDALClose((GDALDatasetH)readDataset);
		GDALClose((GDALDatasetH)pDDst);
		cout << "处理出错" << endl;
		return -3;
	}

	oWo.ChunkAndWarpImage(0, 0, nNewWidth, nNewHeight);
	GDALFlushCache(pDDst);

	GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
	GDALDestroyWarpOptions(psWo);
	GDALClose((GDALDatasetH)readDataset);
	GDALClose((GDALDatasetH)pDDst);

	cout << "重采样完成,生成文件地址为：" << filepathChong <<endl;
	return 0;
}

15、如果一直循环IO读取数据会导致大量缓存，可通过以下代码清楚缓存

		ref.readDataset->FlushCache();
		src.readDataset->FlushCache();