找出长时序遥感影像的缺失日期并用像素均为0的栅格填充缺失日期的文件

· 浏览次数 : 0

小编点评

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

正文

  本文介绍基于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_sizey_size);我们后期的操作需要用到这个行列数,并且会将这个实际存在的栅格文件作为生成新的栅格文件的模板。

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

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

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

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

  至此,大功告成。

与找出长时序遥感影像的缺失日期并用像素均为0的栅格填充缺失日期的文件相似的内容:

找出长时序遥感影像的缺失日期并用像素均为0的栅格填充缺失日期的文件

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

[转帖]懒人有懒招:expect

https://www.cnblogs.com/amazingjxq/archive/2010/10/17/1852889.html 公司的开发机不知道怎么回事,无法用key登录ssh来免去输入密码的繁琐。昨天花了很长的时间来研究,也没找到问题所在。sshd_config里面的设置没有什么问题,反反

[转帖]优化工作中的取舍和大势判断

最近这段时间在参与一个优化项目,所以思考优化的事情比较多。因此最近的文章中,优化相关的内容会比较多,如果有朋友想看些其他方面的题材,可以留言告知。实际上我每天写的东西大多数都是前一天的思想内容,凭空找个话题来写上一篇,时间长了肯定会越来越空洞无物,我也没有那个本事,脑子里的积累写上几年还游刃有余。关

在 Net7.0环境下通过反射创建对象和调用方法

一、介绍 最近没事干,就用闲暇时间写点东西,也记录一下温习历程。老人说的好,好记性,不如烂笔头。时间一长,当时记忆的再清楚,都会变得模糊,索性就写博客记录下来,如果下次需要,直接打开博客就找到了,不用去网上乱找了。 今天我要写一些有关反射创建类型实例和调用方法的文章,泛型的类型的创建和泛型方法的调用

为医生打造专属数字分身!华为云联合万木健康打造医疗医学科普和患者教育数字人引擎

摘要:如今,医生出镜的视频已经成为喜闻乐见的医学科普和患者教育手段,但医生难以抽出时间拍摄、拍摄时间较长、成本较高等制作痛点也日益凸显。对此,国内首个医生AI数字人运营服务商——成都万木健康科技有限公司找到了破局之法。 本文分享自华为云社区《为医生打造专属数字分身!华为云联合万木健康打造医疗医学科普

性能分析: 快速定位SQL问题

在数据库性能调优的实践中,SQL性能分析是至关重要的一环。一个执行效率低下的SQL语句可能会导致整个系统的性能瓶颈。 为了快速定位并解决这些问题,我们需要对SQL进行性能分析。本文将介绍一些常用的方法和技术,帮助大家快速定位SQL问题。 1、找出执行时间最长的SQL 首先,我们需要找到执行时间最长的

文心一言 VS 讯飞星火 VS chatgpt (102)-- 算法导论9.3 8题

八、用go语言,设 X[1..n]和 Y[1..n]为两个数组,每个都包含n个有序的元素。请设计一个 O(lgn)时间的算法来找出数组 X和Y中所有 2n 个元素的中位数。 文心一言: 要在 O(lgn) 时间内找出两个有序数组 X 和 Y 中所有元素的中位数,可以使用二分查找算法。以下是用 Go

#Powerbi 季度时间进度计算

上一篇我们学习了月度时间进度的计算方法,今天我们学习季度时间进度的测算。 思路:找出目前共计消耗了多少天(季度),目前日期所在的季度共有多少天,两者相除即是季度的时间进度 首先列出DAX函数: 本季度第一天 = STARTOFQUARTER(TREATAS({TODAY()},'日期表'[日期]))

[转帖]耗时几个月,终于找到了JVM停顿十几秒的原因

https://www.cnblogs.com/codelogs/p/16060792.html 简介# 最近我们系统出现了一些奇怪的现象,系统每隔几个星期会在大半夜重启一次,分析过程花费了很长时间,令人印象深刻,故在此记录一下。 第一次排查# 由于重启后,进程现场信息都丢失了,所以这个问题非常难以

[转帖]耗时几个月,终于找到了JVM停顿十几秒的原因

https://www.cnblogs.com/codelogs/p/16060792.html 原创:打码日记(微信公众号ID:codelogs),欢迎分享,转载请保留出处。 简介# 最近我们系统出现了一些奇怪的现象,系统每隔几个星期会在大半夜重启一次,分析过程花费了很长时间,令人印象深刻,故在此