基于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