找出長時序遙感影像的缺失日期並用畫素均為0的柵格填充缺失日期的檔案

疯狂学习GIS發表於2024-05-31

  本文介紹基於C++語言的GDAL庫，基於一個儲存大量遙感影像的資料夾，依據每一景遙感影像的檔名中表示日期的那個欄位，找出這些遙感影像中缺失的成像日期，並新生成多個像元值全部為0的柵格檔案，作為這些缺失日期當日的遙感影像檔案的方法。

  首先，我們來看一下本文需要實現的需求。現在有一個資料夾，儲存了從2018年第001天到2022年第361天的全部遙感影像，其中每一景影像的像元個數、空間參考資訊、NoData值等都是一致的。對於這些遙感影像，原本應該是每10天就有1景；但是由於遙感影像資料有缺失，因此部分日期沒有對應的遙感影像。如下圖所示，可以看到比如2018年的061這一天，它就沒有對應的遙感影像。

image

  但是，由於後期處理的需要，我們現在希望對這些缺失日期的遙感影像檔案加以填補——具體的需求是，我們新建若干個像元值全部為0的柵格檔案，作為每一個缺失日期當日的遙感影像檔案；這些填補的、新的遙感影像檔案的各項資訊（比如像元個數、空間參考資訊等）都和原本的檔案一致即可，只要保證全部的像元都是0就行。

  知道了需求，我們就可以開始程式碼的撰寫。本文用到的程式碼具體如下所示。其中，關於C++語言配置GDAL庫的方法，大家可以參考文章在Visual Studio中部署GDAL庫的C++版本（包括SQLite、PROJ等依賴）。

#include <iostream>
#include <iomanip>
#include <sstream>
#include "gdal_priv.h"
#include "cpl_conv.h"

using namespace std;

void create_missing_raster(string path);

int main() {
    string file_path = R"(E:\02_Project\TIFF\TEST)";
    create_missing_raster(file_path);
    return 0;
}

void create_missing_raster(string path)
{
	vector<string> all_file_path;
	for (int year = 2018; year <= 2022; year++)
	{
		for (int day = 1; day <= 361; day += 10)
		{
			ostringstream oss;
			oss << path << "/Albers_MuSyQ.NDVI.16m." << to_string(year) << setfill('0') << setw(3) << to_string(day) << "000000.49SDB.001.h5NDVI-NDVI.tif";
			all_file_path.push_back(oss.str());
		}
	}

	GDALAllRegister();
	int x_size, y_size;
	string one_actual_path;
	for (string& one_file_path : all_file_path)
	{
		GDALDataset* poDataset_actual;
		poDataset_actual = (GDALDataset*)GDALOpen(one_file_path.c_str(), GA_ReadOnly);
		if (poDataset_actual != NULL)
		{
			one_actual_path = one_file_path;
			x_size = poDataset_actual->GetRasterXSize();
			y_size = poDataset_actual->GetRasterYSize();
			GDALClose((GDALDatasetH)poDataset_actual);
			break;
		}
	}

	for (string& one_file_path : all_file_path)
	{
		if (CPLCheckForFile((char*)one_file_path.c_str(), NULL) == FALSE)
		{
			GDALDataset* poDataset;
			GDALDriver* poDriver;
			poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
			GDALDataset* poDataset_actual;
			poDataset_actual = (GDALDataset*)GDALOpen(one_actual_path.c_str(), GA_ReadOnly);
			poDataset = poDriver->CreateCopy(one_file_path.c_str(), poDataset_actual, FALSE, NULL, NULL, NULL);
			double* pafScanline = new double[x_size * y_size];
			for (int i = 0; i < x_size * y_size; i++)
			{
				pafScanline[i] = 0.0;
			}
			GDALRasterBand* poBand;
			poBand = poDataset->GetRasterBand(1);
			poBand->RasterIO(GF_Write, 0, 0, x_size, y_size, pafScanline, x_size, y_size, GDT_Float64, 0, 0);
			delete[] pafScanline;
			GDALClose((GDALDatasetH)poDataset);
			GDALClose((GDALDatasetH)poDataset_actual);
			cout << "New file is :" << one_file_path << endl;
		}
	}
	GDALDestroyDriverManager();
}

  上述程式碼主要都是在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這一天，目前已經有了一景遙感影像。

  至此，大功告成。

相關文章