GDAL切割重采样遥感图像（航拍影像、卫片）

一个小测试程序开发全过程实录，完全新手入门级的实例，如果你还在为处理大影像而发愁，来试试这个称手的工具吧。

Imagec 开发日记

2013-6-25

需求:

影像数据切割，重采样

数据切割的要求是简单的给予矩形的等分切割，并以2的幂次为分割单元，无需使用AOI裁切，

重采样需要实现多种采样模式，用户可以切换采样模式（下文中所提供的代码只是利用了RasterIO的一个特性使用了默认的最近邻重采样方法）

基本思路

考虑是否存在使实用的gdal接口，

自行设计，利用GDAL的读写接口完成数据输入输出工作

初步了解，GDAL并不提供现成的切割和重采样的接口。

/////////////////////////////////////////////////////////////////////

2013-6-26

确定使用自己的方法实现:

裁切

1、读入数据

2、设置剪切范围

3、完成裁切，将数据写入目标文件

2013-6-27

重采样

1.读入数据

2.借鉴Warp工具代码重写重采样代码

将数据写入目标文件

2013-6-28

重采样2

1.读入数据

2.使用RasterIO自带的重采样功能

3.数据写入目标文件

////////////////////////////////////////////////////////////////////////////

2013-6-27

实现算法

裁切：

读入文件

GDALDataset *poDataset;

GDALAllRegister();

//open a dataset

poDataset =(GDALDataset *)GDALOpen(fn,GA_ReadOnly);

if(poDataset == NULL)

设置相关参数为实现裁切做准备

char **papsOptions = NULL;

const int iXSize = poSrcDS->GetRasterXSize()/2;

const int iYSize = poSrcDS->GetRasterYSize()/2;

int iSize = GDALGetBandTypeSize(GDT_Byte)/8;

//生成一个用于存放数据的缓存空间

poDstDS = poDriver->Create(fnDst,iXSize,iYSize,3,GDT_Byte,papsOptions);

GByte *abyByte = new GByte[iXSize*iYSize*bandCount];

int bandMap[3] = {1,2,3};

const char *pszSRS_WKT = poSrcDS->GetProjectionRef();

poDstDS->SetProjection(pszSRS_WKT);

poDstDS->SetGeoTransform(adfGeoTransform);

//进行裁切工作，读取原始的数据，写入目标文件

poSrcDS->RasterIO(GF_Read,0,0,iXSize,iYSize,abyRater,

iXSize,iYSize,GDT_Byte,bandCount,bandMap,iSize*bandCount,iSize*iXSize*3,iSize);

poDstDS->RasterIO(GF_Write,0,0,iXSize,iYSize,abyRater,

iXSize,iYSize,GDT_Byte,bandCount,bandMap,iSize*bandCount,iSize*iXSize*3,iSize);

////////////////////////////////////////////////////////////////

//数据重采样实现：

打开数据

GDALDataset *SrcDS =  (GDALDataset *)GDALOpen(fn,GA_ReadOnly);

为RasterIO设置一些必要参数：

const char *strWkt =pSrcDS->GetProjectionRef();

int bandCount = pSrcDS->GetRasterCount();

int iSize = GDALGetDataTypeSize(GDT_Byte)/8;

int bandMap = {1,2,3};

double adfGeoTransform[6];

pSrcDS->GetGeoTransformation(adfGeoTransform);

double dProj[4] ={0};

ImageRowCol2Projection(adfGeoTransform,0,0,dProj[0],dProj[2]);//具体参考李明录重采样博文

ImageRowCol2Projection(adfGeoTransform,pSrcDS->GetRasterXSize(),pSrcDS->GetRasterYSize(),dProj[1],dProj[3]);

double maxX = dProj[0] >dProj[1]?dProj[0]:dProj[1];

double minX = dProj[0]<dProj[1]?dProj[1]:dProj[0];

double maxY = dProj[2]>dProj[3]?dProj[2]:dProj[3];

double minY = dProj[2]<dProj[3]?dProj[2]:dProj[3];

double fResX = 2;

double fResY = 2;

adfGeoTransform[0] = maxX;

adfGeoTransform[3] = maxY;

adfGeoTransform[1] = adfGeoTransform[1]/fResX;

adfGeoTransform[4] = adfGeoTransform[5]/fResY;

//得到输出图像长宽

int iNewWidth = static_cast<int>((maxX-minX)/ABS(adfGeoTransform[1]) + 0.5);

int iNewHeight = static_cast<int>((maxY-minY)/ABS(adfGeoTransform[5]) + 0.5);

GByte *abyraster2= new GByte[iNewWidth*iNewHeight*bandCount];

创建输出文件：

GDALDataset *pDstDS2 = poDriver->Create("good.img",

iNewWidth,iNewHeight,3,GDT_Byte,NULL);

设置输出数据的投影:

pDstDS2->SetGeoTransform(adfGeoTransform);

pDstDS2->SetProjection(strWkt);

使用RasterIO实现数据读写：（重采样全在这读写之间，参考第一篇文献的 第四种方式 一节）

pSrcDS->RasterIO(GF_Read,0,0,pSrcDS->GetRasterXSize(),pSrcDS->GetRasterYSize(),abyraster2,

iNewWidth,iNewHeight,GDT_Byte,bandCount,bandMap,iSize*bandCount,

iSize*iNewWidth*3,iSize);

pDstDS->RasterIO(GF_Write,0,0,iNewWidth,iNewHeight,GDT_Byte,bandCount,bandMap,iSize*bandCount,

iSize*iNewWidth*3,iSize);

完成读写后关闭打开的数据集，释放占用的内存

if(pDstDS != NULL)

{

GDALClose((GDALDatasetH)pDstDS2);

}

GDALClose((GDALDatasetH)pSrcDS);

 

报错实例

没有注册驱动导致数据读取失败

GDALAllRegister();

ERROR 4'../image/result.img' not recognized as a supported format.

投影问题

setGeoTransform 无法为输出图像计算一个投影

参考文献

GDAL源码剖析（七）之GDAL RasterIO使用说明

如何使用GDAL重采样图像

如何使用GDAL进行AOI裁剪