Matlab基于经纬度点并行提取指定日期的tiff栅格位置的值

前言

该 MATLAB 代码用于从 GeoTIFF 文件中提取基于特定地理位置（经纬度）和日期的某个点的相关数据。代码首先读取一个包含事件数据（日期、经纬度）的 Excel 文件，然后根据日期和位置尝试从存储滑坡风险数据的 GeoTIFF 文件中提取相应的数值。为了加快数据提取过程，代码使用了并行处理。

一、基本说明

数据加载：从 Excel 文件中读取包含日期、经纬度的事件数据。
路径设置：指定 GeoTIFF 文件存储的文件夹路径。
GeoTIFF 数据提取：针对每个事件的日期和位置：
构建文件名并检查文件是否存在。
若文件存在，计算指定经纬度对应的行列索引并提取相应数值。
若坐标超出边界，则赋值为 NaN。
并行处理：使用 parfor 来加速处理每个日期-位置对。
结果存储：将提取的数值存储在数组 gpmresult 或 mswepresult 中。
该代码使得从大规模数据集中高效提取滑坡风险数据成为可能，并且具备处理缺失或越界数据的鲁棒性。

二、代码

% 清空环境、命令窗口，并关闭所有打开的图形窗口
clear all;
clc;
close all;% 加载滑坡事件数据，包括日期、纬度和经度
data = readtable('E:\data.xlsx');
longitude = data{:, 3};  % 提取经度数据
latitude = data{:, 2};   % 提取纬度数据
dates = data{:, 1};      % 提取日期，格式为 'yyyy/MM/dd'% 设置 GeoTIFF 文件夹路径
tif_folder = 'E:\';% 初始化结果数组
gpmresult = [];% 使用并行处理遍历每个数据点
parfor i = 1:length(dates)% 检查日期是单元格格式还是字符数组if iscell(dates)date_str = dates{i};elsedate_str = dates(i, :);end% 将日期字符串转换为 datetime 对象date_obj = datetime(date_str, 'InputFormat', 'yyyy/MM/dd');year_folder = datestr(date_obj, 'yyyy'); % 从日期中提取年份作为文件夹名称% 获取当前点的经纬度lon = longitude(i);lat = latitude(i);% 使用 'yyyymmdd' 格式构建 GeoTIFF 文件路径tif_file = fullfile(tif_folder, year_folder, sprintf('%s.tif', datestr(date_obj, 'yyyymmdd')));% 检查 GeoTIFF 文件是否存在if isfile(tif_file)% 获取 GeoTIFF 文件的元数据并读取数据info = geotiffinfo(tif_file);[data, R] = readgeoraster(tif_file);% 计算指定纬度和经度的行列索引row = ceil((max(R.LatitudeLimits) - lat) / R.CellExtentInLatitude);col = ceil((lon - min(R.LongitudeLimits)) / R.CellExtentInLongitude);% 检查计算的行列是否越界if row <= 0 || col <= 0 || row > size(data, 1) || col > size(data, 2)value30 = NaN; % 如果越界则设置为 NaNelsevalue30 = data(row, col); % 如果在范围内则提取数据值end% 将提取的值添加到结果数组中gpmresult = [gpmresult; value30];else% 如果文件不存在则输出提示信息fprintf('File not found: %s\n', tif_file);end
end
disp('GPM well done');% 第
% 以下部分与第一个代码块类似，但用于另一个数据集
% 使用日期格式 'yyyydoy' (一年中的第几天)% % 设置 MSWEP GeoTIFF 文件夹路径
% tif_folder = 'H:\';
% mswepresult = [];
% 
% % 使用并行处理遍历每个点
% parfor i = 1:length(dates)
%     % 检查日期格式并转换
%     if iscell(dates)
%         date_str = dates{i};
%     else
%         date_str = dates(i, :);
%     end
%     
%     date_obj = datetime(date_str, 'InputFormat', 'yyyy/MM/dd');
%     year_folder = datestr(date_obj, 'yyyy');
%     day_of_year = day(date_obj, 'dayofyear'); % 获取一年中的第几天
%     
%     lon = longitude(i);
%     lat = latitude(i);
%     tif_file = fullfile(tif_folder, year_folder, sprintf('_%d%03d.tif', year(date_obj), day_of_year));
%     
%     if isfile(tif_file)
%         info = geotiffinfo(tif_file);
%         [data, R] = readgeoraster(tif_file);
%         
%         row = ceil((max(R.LatitudeLimits) - lat) / R.CellExtentInLatitude);
%         col = ceil((lon - min(R.LongitudeLimits)) / R.CellExtentInLongitude);
%         
%         if row <= 0 || col <= 0 || row > size(data, 1) || col > size(data, 2)
%             value30 = NaN; % 如果越界则设置为 NaN
%         else
%             value30 = data(row, col); % 如果在范围内则提取数据值
%         end
%         
%         mswepresult = [mswepresult; value30];
%     else
%         fprintf('File not found: %s\n', tif_file);
%     end
% end
% disp('MSWEP well done');