之前写的博客《GDAL源码剖析（十三）之GDAL网格插值说明》里面大致说明了GDAL插值的一些方法和原理，今天使用一部分示例数据进行说明。
首先准备插值用的点数据，这里使用的数据大概是126个点组成的，排列按照X、Y、Z顺序。内容如下：
53414.28,31421.88,39.555 53387.8,31425.02,36.8774 53359.06,31426.62,31.225 53348.04,31425.53,27.416 53344.57,31440.31,27.7945 53352.89,31454.84,28.4999 53402.88,31442.45,37.951 53393.47,31393.86,32.5395 53358.85,31387.57,29.426 53358.59,31376.62,29.223 53348.66,31364.21,28.2538 53362.8,31340.89,26.8212 53335.73,31347.62,26.2299 53331.84,31362.69,26.6612 53351.82,31402.35,28.4848 53335.09,31399.61,26.6922 53331.15,31333.34,24.6894 53344.1,31322.26,24.3684 53326.8,31381.66,26.7581 53396.59,31331.42,28.7137 53372.7,31317.25,25.8215 53404.54,31313.74,26.9055 53416.04,31349.16,31.7509 53424.7,31367.77,34.8919 53414.85,31383.5,37.4818 53399.59,31370.21,32.5866 53386.89,31353.32,30.5459 53383.23,31336.82,29.2504 53421.13,31322.59,27.7593 53468.29,31316.53,27.0276 53434.93,31313.32,26.5662 53456.71,31324.05,28.8742 53491.54,31372.23,33.0459 53470.83,31363.75,32.6194 53414.07,31397.87,39.5041 53446.84,31368.77,34.5865 53438,31337.25,30.1398 53456.08,31344.36,30.1871 53472.76,31401.02,38.5963 53479.13,31432.82,42.3405 53499.63,31422.7,40.7577 53492.35,31402.93,37.9286 53489.32,31390.16,36.34 53449.7,31383.25,36.9367 53444.56,31406.78,41.6945 53464.5,31427.23,43.5075 53503.60,31439.93,41.0365 53515.85,31437.89,39.9929 53505.89,31493.01,33.8673 53485.63,31490.18,36.4479 53493.83,31472.49,39.8801 53482.52,31450.67,42.3327 53508.10,31461.98,40.2172 53513.7,31482.32,36.919 53532.14,31466.79,37.7726 53547.81,31446.42,38.3385 53555.09,31429.29,37.5484 53550.59,31466.8,37.233 53542.39,31490.08,33.8351 53562.32,31481.02,34.4659 53582.49,31461.23,34.1418 53585.43,31474.09,32.5303 53596.28,31464.84,31.7542 53573.76,31494.1,31.3335 53567.67,31505.61,29.5076 53435.85,31465.51,40.6101 53446.92,31438.86,43.5085 53451.21,31462.16,42.1787 53425.14,31442.84,42.2289 53465.83,31471.64,41.4353 53444.54,31478.39,39.4434 53436.92,31489.73,35.15 53446.87,31502.75,32.0204 53424.8,31499.59,31.1902 53415.55,31477.23,36.1728 53405.41,31459.37,36.2749 53407.98,31487.67,31.2808 53410.21,31511.39,29.09 53392.51,31495.92,29.1024 53434.96,31524.73,28.8913 53415.92,31523.57,28.6408 53381.47,31504.81,28.3432 53365.02,31495.85,27.6838 53352.99,31500.58,26.1047 53348.82,31486.82,26.2623 53350.47,31471.08,27.5493 53357.64,31521.57,25.8516 53387.14,31526.14,25.9322 53418.4,31538.09,25.2847 53448.49,31533.77,25.8802 53465.16,31494.46,34.603 53458.79,31508.1,31.0874 53503.63,31519.64,28.9581 53504.38,31505.74,29.4945 53487.85,31513.90,29.3574 53473.73,31522.38,28.2401 53474.65,31534.46,28.1753 53501.33,31537.72,27.601 53526.89,31530.08,29.1993 53538.63,31519.63,29.2767 53327.38,31431.95,26.6129 53326.22,31419.39,26.2965 53591.41,31400.14,31.3133 53580.18,31375.51,30.2236 53603.18,31381.36,28.8233 53529.15,31341.71,27.0601 53583.86,31408.9,33.6867 53565.07,31419.84,36.2957 53544.34,31402.71,36.3476 53568.31,31403.82,34.7975 53517.51,31388.7,35.1603 53521.11,31366.48,31.672 53539.64,31349.28,29.138 53506.62,31335.71,26.9185 53505.15,31350.97,29.6483 53494.79,31379.74,33.5259 53511.44,31374.64,33.1928 53493.14,31354.75,30.8649 53549.66,31508.67,30.3146 53369.18,31500.28,27.7 53442.8,31431.06,43.90 53455.78,31444.08,43.460 53461.02,31374.57,35.5 53437.30,31388.9,38.2 53399.61,31407.54,37 53514.12,31450.64,40.1
下面就开始读取坐标点数据进行插值处理，大致流程是：
首先读取点坐标存入数组中，并统计点的最大最小值，用于计算输出图像的大小；
其次是构造插值算法的参数，并调用GDAL插值函数进行插值；
最后将插值的结果写入图像，释放资源等。
代码如下：
void GDALGridTest()
{
const char* filename = "E:\\Datum\\GDALGrid\\离散点数据文件.dat";
const char* outputfullname = "E:\\Datum\\GDALGrid\\插值结果.tif";
//计算最大最小值
double dfXMin;
double dfXMax;
double dfYMin;
double dfYMax;
double dfZMin;
double dfZMax;
GUInt32 nPoints = 126;
double * padfX = new double[nPoints];
double * padfY = new double[nPoints];
double * padfZ = new double[nPoints];
//将文本文件读入数组并统计出最大最小值
SET_LOCAL;
ifstream ifile(filename, ios::in);
string strBuf;
int i = 0;
while(getline(ifile, strBuf))
{
vector<string> vStr;
// 使用boost库的split拆分字符串
split(vStr, strBuf, is_any_of(","), token_compress_on);
double tmpX = atof(vStr[0].c_str());
double tmpY = atof(vStr[1].c_str());
double tmpZ = atof(vStr[2].c_str());
padfX[i] = tmpX;
padfY[i] = tmpY;
padfZ[i] = tmpZ;
if(i = = 0)
{
dfXMin = tmpX;
dfXMax = tmpX;
dfYMin = tmpY;
dfYMax = tmpY;
dfZMin = tmpZ;
dfZMax = tmpZ;
}
dfXMin = (tmpX<dfXMin ) ? tmpX : dfXMin;
dfXMax = (tmpX>dfXMax ) ? tmpX : dfXMax;
dfYMin = (tmpY<dfYMin ) ? tmpY : dfYMin;
dfYMax = (tmpY>dfYMax ) ? tmpY : dfYMax;
dfZMin = (tmpZ<dfZMin ) ? tmpZ : dfZMin;
dfZMax = (tmpZ>dfZMax ) ? tmpZ : dfZMax;
i++;
}
ifile.close();
REVERT_LOCAL;
//计算图像大小
double pixResoultion = 0.5; //设置分辨率为0.5
GUInt32 nXSize = (dfXMax-dfXMin) / pixResoultion;
GUInt32 nYSize = (dfYMax-dfYMin) / pixResoultion;
GDALAllRegister();
// 离散点内插方法，使用反距离权重插值法
GDALGridInverseDistanceToAPowerOptions *poOptions = new GDALGridInverseDistanceToAPowerOptions();
poOptions->dfPower = 2;
poOptions->dfRadius1 = 20;
poOptions->dfRadius2 = 15;
char *pData = new char[nXSize*nYSize];
// 离散点内插方法，使用反距离权重插值法。使用其他的插值算法，这里换成其他的，还有下面的GDALGridCreate函数的对应参数
GDALGridCreate(GGA_InverseDistanceToAPower, poOptions, nPoints, padfX, padfY, padfZ,
dfXMin, dfXMax, dfYMin, dfYMax, nXSize, nYSize, GDT_Byte, pData, NULL, NULL);
// 创建输出数据集，格式为GeoTiff格式
GDALDriver * pDriver = NULL;
pDriver = GetGDALDriverManager()->GetDriverByName("Gtiff");
GDALDataset *poDataset = pDriver->Create(outputfullname, nXSize,nYSize, 1, GDT_Byte, NULL);
// 设置六参数
double adfGeoTransform[6] = {dfXMin, pixResoultion, 0 , dfYMax, 0, -pixResoultion};
poDataset->SetGeoTransform(adfGeoTransform );
// 写入影像
poDataset->RasterIO(GF_Write, 0, 0, nXSize, nYSize, pData, nXSize, nYSize, GDT_Byte, 1, 0, 0, 0, 0);
// 释放资源 关闭图像
delete poOptions;
delete []pData;
GDALClose(poDataset);
poDataset = NULL;
}顺便说一下，上面代码中有两个宏定义代码，其作用是使用防止ifstream打开中文路径的文件错误，具体的定义如下：
/**
* @brief 将全局区域设为操作系统默认区域
*/
#define SET_LOCAL { locale::global(locale("")); setlocale(LC_ALL,"Chinese-simplified"); }
/**
* @brief 还原全局区域设定
*/
#define REVERT_LOCAL locale::global(locale("C"))
最后放一张插值后的图片：
原始图像
使用ArcMap渲染后的图像