对于特定的日期,进行制图

某些日期可能需要检查具体条带的数据,可视化制图查看更加方便

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值