??本文介紹基于C++語言的GDAL模塊，按照給定的像元行數(shù)與列數(shù)，批量裁剪大量多波段柵格遙感影像文件，并將所得到的裁剪后新的多波段遙感影像文件保存在指定路徑中的方法。

??在之前的文章中，我們多次介紹了在不同平臺(tái)，或基于不同代碼語言，對(duì)柵格遙感影像加以裁剪、批量裁剪的方法，主要包括ArcPy：用矢量批量裁剪柵格圖像、Google Earth Engine谷歌地球引擎GEE矢量數(shù)據(jù)裁剪柵格影像數(shù)據(jù)、基于ENVI實(shí)現(xiàn)柵格圖像按像元行列號(hào)與像元數(shù)量劃定矩形裁剪區(qū)域等；而本文，我們就介紹一下基于C++語言的GDAL模塊，實(shí)現(xiàn)批量裁剪需求的方法。

??首先，來看一下本文的需求。現(xiàn)在有一個(gè)文件夾，如下圖所示，其中具有多個(gè).tiff格式的多波段遙感影像文件（為了方便，我們這里文件夾內(nèi)就只有2個(gè)文件，但實(shí)際上一般我們批量處理的需求肯定遠(yuǎn)遠(yuǎn)大于這個(gè)數(shù)量）。

??我們希望實(shí)現(xiàn)的，就是基于這個(gè)文件夾內(nèi)每一景遙感影像，將其左上角100 * 100像元的這一部分給裁剪下來（如下圖所示），并分別保存為新的遙感影像文件（其中，新的文件名稱就在原有文件名稱后加一個(gè)_C后綴即可），并存放在另一個(gè)指定的結(jié)果文件夾中。我們希望裁剪后的遙感影像，和原有的遙感影像對(duì)比起來，呈現(xiàn)如下圖所示的情況。

??本文所用代碼如下。

#include <iostream>
#include <string>
#include <gdal/gdal.h>
#include <gdal/gdal_priv.h>
using namespace std;

int main()
{
    GDALAllRegister();

    string inputFolder = "/home/cppGDAL/TIF/WFV1_2021_STB/test";
    string outputFolder = "/home/cppGDAL/TIF/WFV1_2021_STB_C";

    CPLStringList fileList;
    fileList = VSIReadDir(inputFolder.c_str());

    for (int i = 0; i < fileList.size(); i++)
    {
        string inputImagePath = fileList[i];
        if (!EQUAL(CPLGetExtension(inputImagePath.c_str()), "tiff"))
        {
            continue;
        }

        string full_path = inputFolder + "/" + inputImagePath;
        GDALDataset *poDataset = (GDALDataset *)GDALOpen(full_path.c_str(), GA_ReadOnly);

        int width = poDataset->GetRasterXSize();
        int height = poDataset->GetRasterYSize();

        int xOffset = 0;
        int yOffset = 0;
        int xSize = 100;
        int ySize = 100;
        cout << full_path;
        string outputImageName = (outputFolder + "/" + inputImagePath.substr(0, inputImagePath.size() - 5) + "_C.tiff");
        cout << outputImageName;
        GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
        GDALDataset *poOutputDataset = poDriver->Create(outputImageName.c_str(), xSize, ySize, 4, GDT_Float32, NULL);

        poOutputDataset->SetProjection(poDataset->GetProjectionRef());
        double adfGeoTransform[6];
        poDataset->GetGeoTransform(adfGeoTransform);
        // adfGeoTransform[1] = adfGeoTransform[1] / (width / xSize);
        // adfGeoTransform[5] = adfGeoTransform[5] / (height / ySize);
        poOutputDataset->SetGeoTransform(adfGeoTransform);

        for (int bandIndex = 1; bandIndex <= 4; bandIndex++)
        {
            float *buffer = new float[xSize * ySize];
            GDALRasterBand *poBand = poDataset->GetRasterBand(bandIndex);
            GDALRasterBand *poOutputBand = poOutputDataset->GetRasterBand(bandIndex);
            poBand->RasterIO(GF_Read, xOffset, yOffset, xSize, ySize, buffer, xSize, ySize, GDT_Float32, 0, 0, NULL);
            poOutputBand->RasterIO(GF_Write, 0, 0, xSize, ySize, buffer, xSize, ySize, GDT_Float32, 0, 0, NULL);
            delete[] buffer;
        }

        GDALClose(poOutputDataset);
        GDALClose(poDataset);
    }

    GDALDestroyDriverManager();

    return 0;
}

??我們來介紹一下上述代碼的具體內(nèi)容。

??首先，我們需要通過GDALAllRegister();，來注冊(cè)所有的GDAL驅(qū)動(dòng)器。同時(shí)，我們定義了輸入和輸出文件夾路徑——inputFolder就是存儲(chǔ)輸入遙感影像（待裁剪的遙感影像）的文件夾路徑，outputFolder則是存儲(chǔ)結(jié)果遙感影像的文件夾路徑。

??其次，我們通過CPLStringList fileList;定義一個(gè)字符串列表，用于存儲(chǔ)文件夾中的文件列表；并使用VSIReadDir函數(shù)讀取輸入文件夾中的所有文件，并將結(jié)果存儲(chǔ)在fileList中。

??接下來，我們使用循環(huán)迭代處理每個(gè)文件。首先，通過string inputImagePath = fileList[i];獲取當(dāng)前文件的路徑；如果文件的擴(kuò)展名不是tiff，則跳過該文件。接下來，對(duì)于文件的擴(kuò)展名是tiff的，我們構(gòu)建完整的輸入文件路徑，并使用GDALOpen函數(shù)打開輸入文件，返回一個(gè)GDALDataset對(duì)象，存儲(chǔ)在poDataset中。

??接下來，我們即可獲取輸入文件的寬度和高度，并定義裁剪區(qū)域的偏移量（左上角像元的位置）、寬度和高度。前面提到了，我這里就是需要在原本遙感影像的最左上角（所以偏移量均為0），裁剪下來100 * 100像元的這一部分。其次，構(gòu)建輸出文件的路徑，并使用GetGDALDriverManager()->GetDriverByName函數(shù)獲取GTiff驅(qū)動(dòng)器對(duì)象，存儲(chǔ)在poDriver中。隨后，我們使用poDriver->Create函數(shù)創(chuàng)建輸出文件，返回一個(gè)GDALDataset對(duì)象，存儲(chǔ)在poOutputDataset中。

??接下來這個(gè)部分需要稍微注意一下。首先，我們使用poOutputDataset->SetProjection設(shè)置輸出文件的投影信息，即與輸入文件相同的投影；其次，使用poDataset->GetGeoTransform獲取輸入文件的地理變換參數(shù)，存儲(chǔ)在adfGeoTransform數(shù)組中。由于在我這里，裁剪后遙感影像的像元大小（即單個(gè)像元的長度與寬度）沒有改變，且裁剪前后柵格遙感影像的左上角像元沒有發(fā)生變化，所以新的柵格遙感影像的地理變換參數(shù)和老的柵格遙感影像比起來，無需有任何改變；但是如果大家的裁剪需求不是這樣的話（比如像元大小發(fā)生變化了，或者是裁剪并不是從左上角像元開始的），那么就需要調(diào)整這6個(gè)地理變換參數(shù)——至于怎么變，這就比較復(fù)雜了，我也還沒完全搞清楚，大家就結(jié)合自己的實(shí)際需求，到GDAL官網(wǎng)查閱即可。最后，我們使用poOutputDataset->SetGeoTransform，設(shè)置輸出文件的地理變換參數(shù)，在我這里就是與輸入文件完全相同的地理變換參數(shù)。

??代碼最后，我們使用循環(huán)迭代處理每個(gè)波段（我這里每一個(gè)遙感影像都是共4個(gè)波段）。首先，創(chuàng)建一個(gè)大小為xSize * ySize的浮點(diǎn)型緩沖區(qū)，并使用poBand->RasterIO從輸入文件中讀取對(duì)應(yīng)波段的像元數(shù)據(jù)到緩沖區(qū)；接下來，使用poOutputBand->RasterIO將緩沖區(qū)中的數(shù)據(jù)寫入到輸出文件對(duì)應(yīng)波段中。隨后，即可釋放緩沖區(qū)內(nèi)存，并關(guān)閉輸出文件和輸入文件。

??運(yùn)行上述代碼，我們即可在結(jié)果文件夾中看到已經(jīng)裁剪好的遙感影像文件，且新的文件的文件名稱也符合我們的要求；如下圖所示。

??至此，大功告成。