基于C++ GDAL依据像元个数批量裁剪栅格图像

疯狂学习GIS

发布于 2024-08-01 09:29:55

950

发布于 2024-08-01 09:29:55

文章被收录于专栏：疯狂学习GIS

本文介绍基于C++语言的GDAL模块，按照给定的像元行数与列数，批量裁剪大量多波段栅格遥感影像文件，并将所得到的裁剪后新的多波段遥感影像文件保存在指定路径中的方法。

在之前的文章中，我们多次介绍了在不同平台，或基于不同代码语言，对栅格遥感影像加以裁剪、批量裁剪的方法，主要包括ArcPy依据矢量要素裁剪多张栅格图像，以及ArcPy栅格裁剪：对齐多个栅格图像的范围、统一行数与列数，还有Google Earth Engine谷歌地球引擎矢量数据裁剪栅格数据等；而本文，我们就介绍一下基于C++语言的GDAL模块，实现批量裁剪需求的方法。

首先，来看一下本文的需求。现在有一个文件夹，如下图所示，其中具有多个.tiff格式的多波段遥感影像文件（为了方便，我们这里文件夹内就只有2个文件，但实际上一般我们批量处理的需求肯定远远大于这个数量）。

我们希望实现的，就是基于这个文件夹内每一景遥感影像，将其左上角100 * 100像元的这一部分给裁剪下来（如下图所示），并分别保存为新的遥感影像文件（其中，新的文件名称就在原有文件名称后加一个_C后缀即可），并存放在另一个指定的结果文件夹中。我们希望裁剪后的遥感影像，和原有的遥感影像对比起来，呈现如下图所示的情况。

本文所用代码如下。

#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;
}

我们来介绍一下上述代码的具体内容。

首先，我们需要通过GDALAllRegister();，来注册所有的GDAL驱动器。同时，我们定义了输入和输出文件夹路径——inputFolder就是存储输入遥感影像（待裁剪的遥感影像）的文件夹路径，outputFolder则是存储结果遥感影像的文件夹路径。

其次，我们通过CPLStringList fileList;定义一个字符串列表，用于存储文件夹中的文件列表；并使用VSIReadDir函数读取输入文件夹中的所有文件，并将结果存储在fileList中。

接下来，我们使用循环迭代处理每个文件。首先，通过string inputImagePath = fileList[i];获取当前文件的路径；如果文件的扩展名不是tiff，则跳过该文件。接下来，对于文件的扩展名是tiff的，我们构建完整的输入文件路径，并使用GDALOpen函数打开输入文件，返回一个GDALDataset对象，存储在poDataset中。

接下来，我们即可获取输入文件的宽度和高度，并定义裁剪区域的偏移量（左上角像元的位置）、宽度和高度。前面提到了，我这里就是需要在原本遥感影像的最左上角（所以偏移量均为0），裁剪下来100 * 100像元的这一部分。其次，构建输出文件的路径，并使用GetGDALDriverManager()->GetDriverByName函数获取GTiff驱动器对象，存储在poDriver中。随后，我们使用poDriver->Create函数创建输出文件，返回一个GDALDataset对象，存储在poOutputDataset中。

接下来这个部分需要稍微注意一下。首先，我们使用poOutputDataset->SetProjection设置输出文件的投影信息，即与输入文件相同的投影；其次，使用poDataset->GetGeoTransform获取输入文件的地理变换参数，存储在adfGeoTransform数组中。由于在我这里，裁剪后遥感影像的像元大小（即单个像元的长度与宽度）没有改变，且裁剪前后栅格遥感影像的左上角像元没有发生变化，所以新的栅格遥感影像的地理变换参数和老的栅格遥感影像比起来，无需有任何改变；但是如果大家的裁剪需求不是这样的话（比如像元大小发生变化了，或者是裁剪并不是从左上角像元开始的），那么就需要调整这6个地理变换参数——至于怎么变，这就比较复杂了，我也还没完全搞清楚，大家就结合自己的实际需求，到GDAL官网查阅即可。最后，我们使用poOutputDataset->SetGeoTransform，设置输出文件的地理变换参数，在我这里就是与输入文件完全相同的地理变换参数。

代码最后，我们使用循环迭代处理每个波段（我这里每一个遥感影像都是共4个波段）。首先，创建一个大小为xSize * ySize的浮点型缓冲区，并使用poBand->RasterIO从输入文件中读取对应波段的像元数据到缓冲区；接下来，使用poOutputBand->RasterIO将缓冲区中的数据写入到输出文件对应波段中。随后，即可释放缓冲区内存，并关闭输出文件和输入文件。

运行上述代码，我们即可在结果文件夹中看到已经裁剪好的遥感影像文件，且新的文件的文件名称也符合我们的要求；如下图所示。