GISAlgorithm学习笔记第3篇：GDAL影像降采样

内容摘要：

本文介绍利用GDAL开源库对tif格式的影像进行降采样的操作，并介绍其原理涉及的重要函数。

往期回顾：

GISAlgorithm学习笔记第1篇：VS2017配置GDAL库

GISAlgorithm学习笔记第2篇：GDAL影像裁剪

目录

一、开发环境

二、实现原理

三、完整代码

---

一、开发环境

1.语言：C++

2.编程环境：Visual Studio2017

3.开源库：GDAL

二、实现原理

降采样原理

降采样的目的是减少数据量。对于图像，降采样意味着减少像素数量。例如，将一个图像的宽度和高度各减半，意味着每个维度的像素数量减少一半，图像总像素数量减少到原来的四分之一，但图像大小保持不变。

本例通过仿射变换矩阵调整实现降采样：

• 像素宽度（adfGeoTransform[1]）和高度（adfGeoTransform[5]）需要根据降采样比例进行调整。例如，如果宽度和高度都缩小一半，那么它们需要乘以2。

代码如下：

// 设置新的仿射变换矩阵
double adfGeoTransform[6];
if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None) {
    adfGeoTransform[1] *= 2; // 像素宽度增加一倍
    adfGeoTransform[5] *= 2; // 像素高度增加一倍
    poDstDS->SetGeoTransform(adfGeoTransform);
}

部分与裁剪相同的函数如

GDALOpen函数

RasterIO函数

GDALDriver类和相关函数

GetGeoTransform函数

SetGeoTransform函数

GetProjectionRef函数

SetProjectionRef函数等内容见上篇笔记讲解。

GISAlgorithm学习笔记第2篇：GDAL影像裁剪

三、完整代码

降采样结果在ArcMap或ENVI软件中显示更加明显。

注意：文件路径必须为\\双斜杠！！！ 

#include<gdal_priv.h>
#include<cpl_conv.h>
#include<iostream>
#include<gdal.h>
using namespace std;
int main()
{
	GDALAllRegister();//注册所有GDAL驱动
	const char* pSrcFilename = "D:\\xx.tif";//修改为你的源文件路径
	const char* pDstFilename = "D:\\xxafter.tif";//修改为你的目标文件路径,保证修改前和修改后文件名不一样


	//打开源栅格文件
	GDALDataset *poSrcDS = (GDALDataset *)GDALOpen(pSrcFilename, GA_ReadOnly);
	if (poSrcDS == NULL) {//若打开文件失败，输出错误信息
		cout << "Cann't open file" << pSrcFilename << endl;
	}
	//获取源文件中的第一个波段（波段从1开始计数）
	GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);

	// 获取源文件的宽度和高度
	int srcWidth = poBand->GetXSize();
	int srcHeight = poBand->GetYSize();

	// 计算新图像的宽度和高度（原分辨率的一半）
	int dstWidth = srcWidth / 2;
	int dstHeight = srcHeight / 2;

	// 计算新的仿射变换矩阵
	double adfGeoTransform[6];
	if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None) {

	//为扫描线分配内存
		float *pafScanline = (float *)CPLMalloc(sizeof(float) * dstWidth * dstHeight);

		// 从源波段中读取数据到分配的内存中
		// RasterIO函数用于执行实际的读取操作
		poBand->RasterIO(GF_Read, 0, 0, srcWidth, srcHeight, pafScanline, dstWidth, dstHeight, GDT_Float32, 0, 0);

		// 获取GTiff驱动，用于创建目标文件
		GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); //HFA for .img

		//若获取失败，退出程序
		if (poDriver == NULL) exit(1);

		// 使用GTiff驱动创建一个新的栅格文件
		GDALDataset *poDstDS = poDriver->Create(pDstFilename, dstWidth, dstHeight, 1, GDT_Float32, NULL);

		/*GDALDataset *poDstDS = poDriver->Create(pDstFilename, poSrcDS->GetRasterXSize() / 2, poSrcDS->GetRasterYSize() / 2, 1, GDT_Float32, NULL);
		if (poDstDS == NULL) {
			std::cerr << "Could not create dataset" << std::endl;
			GDALClose(poSrcDS);
		}*/

		// 设置新的仿射变换矩阵
		poDstDS->SetGeoTransform(adfGeoTransform);

		// 获取源文件的地理变换参数,投影信息
		double adfGeoTransform[6];
		if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None) {
			// 如果成功获取，设置目标文件的地理变换参数
			//poDstDS->SetGeoTransform(adfGeoTransform);
			adfGeoTransform[1] *= 2; // 像素宽度增加一倍
			adfGeoTransform[5] *= 2; // 像素高度增加一倍
			poDstDS->SetGeoTransform(adfGeoTransform);
		}

		// 如果成功获取，设置目标文件的地理变换参数
		const char* p_WKT = poSrcDS->GetProjectionRef();
		if (p_WKT != NULL) {

			// 如果成功获取，设置目标文件的地理变换参数
			poDstDS->SetProjection(p_WKT);
		}

		// 获取目标文件中的第一个波段
		GDALRasterBand *poDstBand = poDstDS->GetRasterBand(1);

		// 将读取的数据写入目标波段
		poDstBand->RasterIO(GF_Write, 0, 0, dstWidth, dstHeight, (void*)pafScanline, dstWidth, dstHeight, GDT_Float32, 0, 0);

		// 释放扫描线内存
		CPLFree(pafScanline);

		// 关闭目标文件和源文件
		GDALClose((GDALDatasetH)poDstDS);
		GDALClose((GDALDatasetH)poSrcDS);
		cout << "降采样完成"<<endl;
		return 0;
	}
}

---

欢迎交流🌹🌹