基于gdal的格网插值

      格网插值就是使用离散的数据点创建一个栅格图像的过程。通常情况下，有一系列研究区域的离散点，如果我们想将这些点转换为规则的网格数据来进行进一步的处理，或者和其他网格数据进行合并 等处理，就需要使用格网插值算法。

        我们的应用是在海洋数据的处理的处理上,像海洋温度数据，在海洋与陆地的接触部分，数据会有缺失 ，对这些数据缺失的点就需要进行数据的插值处理，才能进行接下来其他的处理，所以插值处理尤为重要。

        在gdal库中有对应各个插值算法的算法参数结构体，我一直用的是反距离权重插值算法（GDALGridInverseDistanceToAPowerOptions），这个根据你应用的数据可以使用不同的算法，gdal中插值的算伐有四种。下面是数据插值处理的代码。

#include <QCoreApplication>

#include <gdal_priv.h>

#include <iostream>

#include <fstream>

#include <string>

#include <vector>

#include <gdalgrid.h>

#include <QDebug>

using namespace std;

using std::string;  //使用string对象

using std::vector; //使用vector

//自己定义的split函数，作用相当于boost库的split函数

void Split(const std::string& src, const std::string& separator, std::vector<std::string>& dest) //字符串分割到数组

{

        //参数1：要分割的字符串；参数2：作为分隔符的字符；参数3：存放分割后的字符串的vector向量

    string str = src;

    string substring;

    string::size_type start = 0, index;

    dest.clear();

    index = str.find_first_of(separator,start);

    do

    {

        if (index != string::npos)

        {

            substring = str.substr(start,index-start );

            dest.push_back(substring);

            start =index+separator.size();

            index = str.find(separator,start);

            if (start == string::npos) break;

        }

    }while(index != string::npos);

    //the last part

    substring = str.substr(start);

    dest.push_back(substring);

}

void DpdInterpolation(const char*pszPoints,const char*pszDst)

{

    //const char*pszPoints = "";//目标文件

    //const char*pszDst = "E:\\odp_workplace\\odp_data\\testdata\\interpolation\\Dpdinterpolation_test.grd";//生成的结果文件

    //设置输出图像的分辨率为0.5

    double pixResoultion = 0.5;

    //计算离散点XY坐标的最大最小值，用于确定输出图像的大小

    double dfXMin;

    double dfXMax;

    double dfYMin;

    double dfYMax;

    double dfZMin;

    double dfZMax;

    double tmpX;

    double tmpY;

    double tmpZ;

    int i = 0;

    int nPoints = 125;//设置离散点的个数

    double * padfX = new double[nPoints];

    double * padfY = new double[nPoints];

    double * padfZ = new double[nPoints];

    //将文本文件读入数组

    ifstream ifile(pszPoints,ios::in);

    string strBuf;

    while(getline(ifile,strBuf))

    {

        vector<string> vstr;

        Split(strBuf,",",vstr);

        tmpX = atof(vstr[0].c_str());

        tmpY = atof(vstr[1].c_str());

        tmpZ = atof(vstr[2].c_str());

        qDebug()<<tmpX<<tmpY<<tmpZ<<endl;

        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();

    //计算图像大小

    int nXSize = (dfXMax-dfXMin)/pixResoultion;

    int nYSize = (dfYMax-dfYMin)/pixResoultion;

    //构造插值算法参数并设置参数值

    GDALGridInverseDistanceToAPowerOptions *poOptions = new GDALGridInverseDistanceToAPowerOptions();

    poOptions->dfPower = 2;//设置权重指数p

    poOptions->dfRadius1 = 20;//设置搜索椭圆第一半径

    poOptions->dfRadius2 = 15;//设置搜索椭圆第二半径

    char *pData = new char[nXSize*nYSize];

    GDALGridCreate(GGA_InverseDistanceToAPower,poOptions,nPoints,padfX,padfY,padfZ,dfXMin,dfXMax,dfYMin,dfYMax,nXSize,nYSize,GDT_Float32,pData,NULL,NULL);

    delete  poOptions;

    //创建新的数据集

    GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GS7BG");//GSAG GSBG GS7BG

    GDALDataset *poDataset = pDriver->Create(pszDst,nXSize,nYSize,1,GDT_Float32,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_Float32,1,0,0,0,0);

    delete []pData;

    GDALClose(poDataset);

}

要使用这段代码就得配置好环境，我的编译环境为qt5.11.2+vs2017+gdal