这里做个备份，原文链接
遥感影像小则几百兆，大则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 200void 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;}
