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