使用GDAL进行大幅图像的分块读写

本博客为博主原创，如果觉得有用的话，关注博主表示支持哦！！！

#include "stdafx.h"

#include "windows.h"
#include <iostream> 
#include "gdal.h"
#include"gdal_priv.h"
using namespace std;


#define block_size 2048


int main()
{
GDALAllRegister();
//申明文件路径
const char* imgPath = "F:\\gdal\\data\\DOM\\E98D9_N33D0_20170103_GF1_WFV1_DOM.tif";
//创建一个数据集poDataset
GDALDataset* poDataset;
//打开数据集
poDataset = (GDALDataset *)GDALOpen(imgPath, GA_ReadOnly);
//获取数据集的宽，高，波段数和投影信息
int nBandCount = poDataset->GetRasterCount();
int nImgSizeX = poDataset->GetRasterXSize();
int nImgSizeY = poDataset->GetRasterYSize();
const char *Pro = poDataset->GetProjectionRef();
//申请内存空间，使用WORD类型申明，因为数据是无符号整形16的格式
WORD *pBuf = new WORD[block_size*block_size*nBandCount];
//初始化内存使用memset语句
//memset(pBuf, 0, nImgSizeX*nImgSizeY*nBandCount*sizeof(WORD));
//申明数据读入方式按照1，2，3，4方式读入
int panbandmap[4] = { 1, 2, 3, 4 };
//新建数据集poDataset2，新建驱动
GDALDataset *poDataset2; //GDAL数据集 poDriver2
GDALDriver*poDriver2;
//申明数据格式
const char *pszFormat2 = "Gtiff";
//驱动获取数据格式
poDriver2 = GetGDALDriverManager()->GetDriverByName(pszFormat2);
//申明写出图像路径
const char*imgPath3 = ("F:\\gdal\\data\\DOM4.tif");
//数据集poDataset2初始化，用来定义数据格式
poDataset2 = poDriver2->Create(imgPath3, nImgSizeX, nImgSizeY, nBandCount, GDT_UInt16, NULL);
//poDataset2->RasterIO(GF_Write, 0, 0, nImgSizeX, nImgSizeY, pBuf, nImgSizeX, nImgSizeY, GDT_UInt16, nBandCount, panbandmap, 0, 0, 0);
//设置写出数据6参数
double padfTransform[6];
poDataset->GetGeoTransform(padfTransform);
//获取原始影像nodata值
double nodata;
nodata=poDataset->GetRasterBand(1)->GetNoDataValue();
//分块算法
for (int i = 0; i < nImgSizeY; i += block_size)
{
for (int j = 0; j < nImgSizeX; j += block_size)
{
//定义两个变量来保存分块大小
int nXBK = block_size;
int nYBK = block_size;
//如果最下面和最右边的块不够256，剩下多少读取多少
if (i + block_size > nImgSizeY)     //最下面的剩余块
nYBK = nImgSizeY - i;
if (j + block_size > nImgSizeX)     //最右侧的剩余块
nXBK = nImgSizeX - j;
//将图像数据集读到pBuf指针中
poDataset->RasterIO(GF_Read, j, i, nXBK, nYBK, pBuf, nXBK, nYBK, GDT_UInt16, nBandCount, panbandmap, 0, 0, 0);
//写到结果图像
poDataset2->RasterIO(GF_Write, j, i, nXBK, nYBK, pBuf, nXBK, nYBK, GDT_UInt16, nBandCount, panbandmap, 0, 0, 0);


}
}
//设置写出图像的nodata值
poDataset2->GetRasterBand(1)->SetNoDataValue(nodata);
//设置图像坐标系
poDataset2->SetGeoTransform(padfTransform);
//设置投影方式
poDataset2->SetProjection(Pro);
cout << "分块成功" << endl;
//释放内存
delete[]pBuf;
pBuf = nullptr;
delete poDataset;
delete poDataset2;
}