c++实现使用GDAL实现大幅影像的快速读取

        遥感影像小则几百兆，大则5,6GB，所以在使用GDAL进行图像读取时面临读写速度较慢的问题，我们可以深入研究gdal中RasterIO函数的机制，发现该函数是通过一行一行读取影像来实现影像读入内存的，所以我们在分块读取的时候也按照几行几行读取这样会加快读取速度，而行数可以通过我们设定的内存大小，即下面代码中的RAM_SIZE=200M来计算得到行数，这样读取速度大概可以提高19倍之多(对于5.5个GB的影像)，这样的算法可以加快计算机速度，提升我们工作效率，下面是快速读取影像的代码，请读者仔细分析其中的参数的意思，加深自己的理解，学习是个不断磨炼的过程，希望大家可以耐得住寂寞，写的代码也不是很完美，大家多多指正！！！

#include "stdafx.h"

#include <stdlib.h> 
#include "gdal_priv.h"
#include <algorithm>
#include <Windows.h>
#include <stdio.h> 
#include <iostream>
using namespace std;


#define RAM_RIZE 200


void main()
{
GDALAllRegister();
GDALDataset *poDataset2;
GDALDriver*poDriver2;
//申明数据格式
const char *pszFormat2 = "Gtiff";
//驱动获取数据格式
poDriver2 = GetGDALDriverManager()->GetDriverByName(pszFormat2);
//申明写出图像路径
const char*imgPath3 = ("F:\\gdal\\data\\DOM5.tif");
//数据集poDataset2初始化，用来定义数据格式
const char *strImg = ("F:\\gdal\\data\\DOM\\ddd.tif");
GDALDataset *ImgBef = (GDALDataset*)GDALOpen(strImg, GA_ReadOnly);
if (ImgBef == NULL)
{
printf("Open Img Failed:\n%s\n", strImg);

}
int nCols = ImgBef->GetRasterXSize();         //获取影像信息
int nRows = ImgBef->GetRasterYSize();
int nBands = ImgBef->GetRasterCount();
GDALDataType gBand = ImgBef->GetRasterBand(1)->GetRasterDataType();
int nBits = GDALGetDataTypeSize(gBand);


double GeoTransform[6];                      //获取坐标信息
ImgBef->GetGeoTransform(GeoTransform);
const char *sProRef = ImgBef->GetProjectionRef(); //获取投影信息


int nStepSize = (RAM_RIZE * 1024 * 1024) / (nCols*nBands);
int nStepNum = nRows / nStepSize; if (nRows%nStepSize) nStepNum++;
int *pBand = new int[nBands]; for (int gi = 0; gi<nBands; gi++) { pBand[gi] = gi + 1; }
int isize = GDALGetDataTypeSize(GDT_UInt16) / 8;
double nodata;
nodata = ImgBef->GetRasterBand(1)->GetNoDataValue();
poDataset2 = poDriver2->Create(imgPath3, nCols, nRows, nBands, GDT_UInt16, NULL);
for (int k = 0; k < nStepNum; k++)
{
int ybeg = max(0, min(nStepSize*k, nRows - 1));
int yend = max(0, min(nStepSize*(k + 1), nRows));
WORD *pImg = new WORD[(yend - ybeg)*nCols*nBands];       //存储开操作时输入影像
memset(pImg, 0, (yend - ybeg)*nCols*nBands*sizeof(WORD));


ImgBef->RasterIO(GF_Read, 0, ybeg, nCols, (yend - ybeg), pImg, nCols, (yend - ybeg),GDT_UInt16, nBands, pBand, isize*nBands, isize*nBands*nCols, isize);
poDataset2->RasterIO(GF_Write, 0, ybeg, nCols, (yend - ybeg), pImg, nCols, (yend - ybeg), GDT_UInt16, nBands, pBand, isize*nBands, isize*nBands*nCols, isize);
delete[]pImg; pImg = NULL;
}
poDataset2->GetRasterBand(1)->SetNoDataValue(nodata);
//设置图像坐标系
poDataset2->SetGeoTransform(GeoTransform);
//设置投影方式
poDataset2->SetProjection(sProRef);
cout << "分块成功" << endl;
delete ImgBef;
delete poDataset2;

}