某些日期可能需要检查具体条带的数据,可视化制图查看更加方便
clc;
clear;
% 设置文件夹路径和文件名 图像分辨率
data_folder = 'FY_3E_data'; % HDF5数据文件夹路径
grid_scale = 0.1; %画图的尺度大小,以度为单位
hdf5_files = dir(fullfile(data_folder, '*.hdf')); % 获取所有HDF5文件
yeardoy = zeros(size(hdf5_files)); % 存储每个文件的年日数据
% 提取年日数据
for i = 1:length(hdf5_files)
yeardoy(i,1) = double(string(hdf5_files(i).name(20:27)));
end
% 去除重复的年日数据,保留唯一值
yeardoy = unique(yeardoy);
% 个性化定义需要制图的日期列表
dates_C = {'20230901', '20230907', '20230913', '20230915', '20230918', '20230921'};
dates_Ku = {'20230908', '20230910', '20230911', '20230917', '20230927', '20231001'};
% 定义路径
longitude_data_path = '/10km/Geolocation/GridInfo/Longitude';
latitude_data_path = '/10km/Geolocation/GridInfo/Latitude';
Zenith_path = '/10km/Geolocation/VV/SensorZenith';
backscatter_values_path = '/10km/Data/VV/Sigma0';
dim1_global_longitude = -180:grid_scale:180;
dim1_global_latitude = -90:grid_scale:90;
% 创建进度条
%h = waitbar(0, 'all_days_Processing...');
for ii = 1:length(yeardoy)
% 处理C波段和Ku波段数据
yeardoy_str = num2str(yeardoy(ii));
yeardoy_cell = {yeardoy_str};
if ismember (yeardoy_cell,dates_C)
file_names = {hdf5_files.name};
%cellfun()对每个cell应用函数,contains()检查是否包括特定字母
cbandlist = file_names(cellfun(@(x) contains(x, 'WRADC') && contains(x, num2str(yeardoy(ii))) , file_names));
band = 'C';
fileList = cbandlist;
% 将文件名转换回结构体数组
fileStructList = struct('name', {});
for k = 1:length(fileList)
fileStructList(k).name = fileList{k};
end
[errorOccurred] = ribbon_draw(data_folder, fileStructList, ...
Zenith_path, backscatter_values_path, longitude_data_path, latitude_data_path);
if errorOccurred
fprintf('Skipping %d due to error in processing.\n', yeardoy(ii));
end
end
if ismember (yeardoy_cell,dates_Ku)
file_names = {hdf5_files.name};
kubandlist = file_names(cellfun(@(x) contains(x, 'WRADK') && contains(x, num2str(yeardoy(ii))) , file_names));
band = 'Ku';
fileList = kubandlist;
fileStructList = struct('name', {});
for k = 1:length(fileList)
fileStructList(k).name = fileList{k};
end
[errorOccurred] = ribbon_draw(data_folder, fileStructList, ...
Zenith_path, backscatter_values_path, longitude_data_path, latitude_data_path);
if errorOccurred
fprintf('Skipping %d due to error in processing.\n', yeardoy(ii));
end
end
end
条带成图函数如下
function [errorOccurred] = ribbon_draw(data_folder, bandlist, Zenith_path, backscatter_values_path, longitude_data_path, latitude_data_path)
%这个函数用于生成FY-3E WRAD条带图像, data_folder为数据文件夹,bandlist为当日C波段或Ku波段列表,
%backscatter_values_path, longitude_data_path, latitude_data_path分别为各自在hdf5文件中的路径,
% longitude_data_path = '/10km/Geolocation/GridInfo/Longitude';
% latitude_data_path = '/10km/Geolocation/GridInfo/Latitude';
% Zenith_path = '/10km/Geolocation/VV/SensorZenith';
% backscatter_values_path = '/10km/Data/VV/Sigma0';
% 初始化错误标志
errorOccurred = false;
% 打开错误日志文件,用于追加错误信息
logFile = 'error_log.txt'; % 错误日志文件名
fid = fopen(logFile, 'a'); % 打开文件用于追加,如果文件不存在则创建
if fid == -1
error('无法打开错误日志文件 %s', logFile);
end
% 循环处理每个波段文件
length_bandlisth = length(bandlist);
for i = 1:length_bandlisth
yeardoy = bandlist(i).name(20:32);
code = bandlist(i).name(10:15);
file = fullfile(data_folder, bandlist(i).name); % 确保bandlist是一个包含文件名的结构体数组
%disp(['Processing file: ', file]);
try
% 读取数据
Zenith = h5read(file, Zenith_path);
Zenith = permute(Zenith, [2, 3, 1]);
Zenith = double(Zenith) / 100;
backscatter_values = h5read(file, backscatter_values_path);
backscatter_values = permute(backscatter_values, [2, 3, 1]);
longitude_value = h5read(file, longitude_data_path);
latitude_value = h5read(file, latitude_data_path);
%逐条带检查图像 %这几个是为了减小运行压力,因此间隔取值
backscatter = backscatter_values(:,:,1);
latitude_value = latitude_value(1:3:end, 1:3:end);
longitude_value = longitude_value(1:3:end, 1:3:end);
backscatter = backscatter(1:3:end, 1:3:end);
%figure;
figure('Visible','off');
ax3 = geoaxes;
geoscatter(ax3, latitude_value(:), longitude_value(:), 2, backscatter(:), 'filled', 'square');
titleStr = sprintf('%s %s %drd out of %d', yeardoy, code, i, length_bandlisth);
title(titleStr);
geobasemap(ax3, 'grayterrain');
geolimits(ax3, 'auto');
caxis(ax3, [-50, 10]);
colorbar;
outputFilename = sprintf('%s_%s_%drd_out_of_%d',yeardoy,code, i, length_bandlisth);
fullPath = fullfile('D:\temporary\all_projects\FY-3E_WindRAD\figures_trial', outputFilename);
print(gcf, fullPath, '-dpng', '-r100');
close(gcf);
catch ME
% 捕获错误并输出错误信息
fprintf('Error processing file: %s\n', file);
fprintf('Error message: %s\n', ME.message);
% 将错误信息写入日志文件
fprintf(fid, 'Error processing file: %s\n', file);
fprintf(fid, 'Error message: %s\n', ME.message);
fprintf(fid, '----------------------------------------\n');
errorOccurred = true; % 设置错误标志
break; % 跳出当前文件的处理
end
end
% 关闭错误日志文件
fclose(fid);
end