C++中GDAL：用像素均为0的栅格填补时序遥感数据中缺少的时相

  本文介绍基于C++语言的GDAL库，基于一个存储大量遥感影像的文件夹，依据每一景遥感影像的文件名中表示日期的那个字段，找出这些遥感影像中缺失的成像日期，并新生成多个像元值全部为0的栅格文件，作为这些缺失日期当日的遥感影像文件的方法。

  首先，我们来看一下本文需要实现的需求。现在有一个文件夹，存储了从2018年第001天到2022年第361天的全部遥感影像，其中每一景图像的像元个数、空间参考信息、NoData值等都是一致的。对于这些遥感影像，原本应该是每10天就有1景；但是由于遥感影像数据有缺失，因此部分日期没有对应的遥感影像。如下图所示，可以看到比如2018年的061这一天，它就没有对应的遥感影像。

  但是，由于后期处理的需要，我们现在希望对这些缺失日期的遥感影像文件加以填补——具体的需求是，我们新建若干个像元值全部为0的栅格文件，作为每一个缺失日期当日的遥感影像文件；这些填补的、新的遥感影像文件的各项信息（比如像元个数、空间参考信息等）都和原本的文件一致即可，只要保证全部的像元都是0就行。

  知道了需求，我们就可以开始代码的撰写。本文用到的代码具体如下所示。其中，关于C++语言配置GDAL库的方法，大家可以参考文章Visual Studio配置C++的GDAL库！。

1#include 

2#include 

3#include 

4#include "gdal_priv.h"

5#include "cpl_conv.h"

6

7using namespace std;

8

9void create_missing_raster(string path);

10

11int main() {

12    string file_path = R"(E:\02_Project\TIFF\TEST)";

13    create_missing_raster(file_path);

14    return 0;

15}

16

17void create_missing_raster(string path)

18{

19    vector all_file_path;

20    for (int year = 2018; year 

21    {

22        for (int day = 1; day 

23        {

24            ostringstream oss;

25            oss 

26            all_file_path.push_back(oss.str());

27        }

28    }

29

30    GDALAllRegister();

31    int x_size, y_size;

32    string one_actual_path;

33    for (string& one_file_path : all_file_path)

34    {

35        GDALDataset* poDataset_actual;

36        poDataset_actual = (GDALDataset*)GDALOpen(one_file_path.c_str(), GA_ReadOnly);

37        if (poDataset_actual != NULL)

38        {

39            one_actual_path = one_file_path;

40            x_size = poDataset_actual->GetRasterXSize();

41            y_size = poDataset_actual->GetRasterYSize();

42            GDALClose((GDALDatasetH)poDataset_actual);

43            break;

44        }

45    }

46

47    for (string& one_file_path : all_file_path)

48    {

49        if (CPLCheckForFile((char*)one_file_path.c_str(), NULL) == FALSE)

50        {

51            GDALDataset* poDataset;

52            GDALDriver* poDriver;

53            poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

54            GDALDataset* poDataset_actual;

55            poDataset_actual = (GDALDataset*)GDALOpen(one_actual_path.c_str(), GA_ReadOnly);

56            poDataset = poDriver->CreateCopy(one_file_path.c_str(), poDataset_actual, FALSE, NULL, NULL, NULL);

57            double* pafScanline = new double[x_size * y_size];

58            for (int i = 0; i 

59            {

60                pafScanline[i] = 0.0;

61            }

62            GDALRasterBand* poBand;

63            poBand = poDataset->GetRasterBand(1);

64            poBand->RasterIO(GF_Write, 0, 0, x_size, y_size, pafScanline, x_size, y_size, GDT_Float64, 0, 0);

65            delete[] pafScanline;

66            GDALClose((GDALDatasetH)poDataset);

67            GDALClose((GDALDatasetH)poDataset_actual);

68            cout 

69        }

70    }

71    GDALDestroyDriverManager();

72}

  上述代码主要都是在create_missing_raster(string path)函数内实现具体功能的，我们主要就对这一函数加以讲解。

  首先，我们需要基于文件夹中遥感影像文件的文件名称特征，遍历生成文件名列表。在这里，我们使用两个嵌套的for循环，生成所有可能的栅格图像文件名，并将这些文件名保存在all_file_path向量中。其中，栅格图像的文件名根据年份和天数生成，并通过setfill('0')与setw(3)这两个函数保证我们生成的日期满足YYYYDDD这种格式。

  随后，基于GDALAllRegister这一GDAL库的初始化函数，用于注册所有支持的数据格式驱动程序。接下来，我们使用GDALOpen函数，从2018001这一天开始，通过循环打开对应名字的文件，直到找到文件夹中第一个实际存在的栅格图像文件（poDataset_actual），并获取其栅格图像的行列数（x_size和y_size）；我们后期的操作需要用到这个行列数，并且会将这个实际存在的栅格文件作为生成新的栅格文件的模板。

  接下来，我们遍历文件名列表all_file_path，对每个文件名进行处理。对于不存在的栅格图像文件，使用GDALDriver创建一个新的数据集（poDataset），并将其中的像元值设置为0。如果栅格图像文件已经存在，则跳过不处理。其中，在对缺失的栅格图像加以生成时，我们首先使用GetGDALDriverManager()->GetDriverByName函数获取GDAL驱动程序对象，然后使用CreateCopy函数创建新的栅格图像；其中，我们就是以前期找到的文件夹中第一个实际存在的栅格图像文件one_actual_path为模板。随后，我们用0填充新创建的栅格图像，并使用RasterIO函数对栅格图像的像元进行写入操作。

  最后，在上述处理完成后，使用GDALClose函数关闭数据集，并输出新创建的栅格图像的文件名。随后，我们使用GDALDestroyDriverManager销毁GDAL驱动程序管理器，释放资源。

  随后，我们运行代码，可以看到每一个新生成的栅格图像文件（也就是原本当日没有成像的日期对应的遥感影像）都会打印出来。

  随后，我们打开文件夹，可以看到之前没有遥感影像的日期，目前也都存在一景遥感影像与其对应了。比如2018年的061这一天，目前已经有了一景遥感影像。

  至此，大功告成。