%% 清理变量|加载路径
clear; close all; close all; fun_path = genpath('functions'); addpath(fun_path)
%% 设置python环境
% pyenv('Version', 'C:\Users\huayi\.conda\envs\gdal_env\python.exe');
pyenv('Version', 'D:\Program Files\Python311\python.exe');
try
    py.importlib.import_module('osgeo.gdal');
    disp('GDAL模块测试成功');
catch ME
    disp('导入GDAL失败');
    disp(ME.message);
    return;
end
%% 设置原始数据及反演结果路径
Root_Path1 = 'F:\2024.7 美国原始数据'; % 原始数据根目录
Root_Path2 = 'F:\2024.7 美国核查结果'; % 核查结果根目录
% Root_Path1 = 'E:\2024.12 中国电厂数据';
% Root_Path2 = 'E:\2024.12 中国电厂结果';

%% 设置吸收系数、反演波段及数据库版本
% 气体吸收系数
version = '25.1.2'; % US
% version = '25.1.3'; % Togtoh
[k_unit, PARM_INPUT] = load_gas_absorption_coef(version);
% 反演波段
band_file = 'data/Band_CO2_v1.xlsx'; band_data = readtable(band_file, 'VariableNamingRule', 'preserve', 'TreatAsEmpty', '');
% 排放数据库
data_file = 'data/US_electricity.xlsx'; sheetName = '2025.1.12';
% data_file = 'data/Togtoh.xlsx'; sheetName = '2025.1.16';
% data_file = 'data/Xinjiang-multisource.xlsx'; sheetName = '2025.1.16';
T = readtable(data_file, 'Sheet', sheetName, 'VariableNamingRule', 'preserve');
%% 设置反演参数
force = false; % 是否强制反演
skip = false; % 是否跳过裁剪
IfPlot = 'on'; % 是否开启绘图
PARM_CROP.range = 350; % 裁剪范围
PARM_CROP.percent = 2; % 阈值提取的分位数
PARM_CROP.thresh = 50; % 保留的相连像素数
sigma = 0.5; % 后处理参数
gas = 'CO2'; % 反演气体
algorithm = 'RMF'; % 反演算法
prefix = '20250108_'; % 输出文件名前缀
%% 循环处理
len = size(T, 1); fprintf("当前景数为：%d\n", len);
for i = 1:len
    close all % 清除图窗
    %% 读取基本信息
    Item = GetFacilityData(T, i);
    Item.algorithm = algorithm;
    Item = prepare_basis(Root_Path1, Root_Path2, Item);
    PrintBasicInfo(Item);
    %% 数据反演
    InversionData(Item, band_data, PARM_INPUT, k_unit, force);
    manualCorrection(Item, force) % 手动校准导出结果
    % continue;
    %% 获取结果
    TIFF_2 = GetTIFFData(Item, PARM_CROP);    % 获取导出的TIFF反演数据
    TIFF_2 = ThreshTIFF(TIFF_2, PARM_CROP); % 羽流自动提取
    [TIFF_3, flag] = MaskPlume2(TIFF_2);
    %% 气象资料
    [VMA, Wind] = MainFunctionIME(Item); % IME方法所需的相关数据
    if(~flag) % 无反演结果
        %% 赋值
        T.wind_speed(i) = Wind.speed; % 风速
        T.wind_direction(i) = Wind.direction; % 风向
        T.wind_speed_std(i) = Wind.speed_std; % 风速不确定度
        T.wind_fit(i) = T.emission_epa(i) * Q.L / IME.data / 3.6; % 风速拟合
        T.inversion_std(i) = TIFF_5.finaldata_u; % 反演不确定度
        continue; % 可选
    else
        % 绘图
        if(strcmp(IfPlot, 'on'))
            visualizeTIFF(TIFF_3)
        end
        %% 核算排放强度结果
        IME = calc_IME(TIFF_3, VMA, gas);    % 二氧化碳整体质量增量 integrated mass enhancement(IME)
        Q_auto = calc_emission(TIFF_3, Wind, IME);    % 计算Emission的函数
        fprintf("自动提取: "); cprintf('*red',sprintf('%.2f±%.2f(t/h)\n', Q_auto.ton_h, Q_auto.ton_h_u));
        Err_auto = (Q_auto.ton_h - T.emission_epa(i)) / T.emission_epa(i);
        fprintf("自动提取误差: %.2f%%\n", Err_auto * 100);
        [Q_manual, TIFF_4, IME] = CutTIFFByHand(TIFF_3, VMA, Wind, Item.path.result_path, skip, IfPlot, gas);
        Err_mannual = (Q_manual.ton_h - T.emission_epa(i)) / T.emission_epa(i);
        fprintf("手动提取误差: %.2f%%\n", Err_mannual * 100);
        if(~skip) % 如果跳过
            enter_content = input(sprintf("请选择提取结果：(当前值：%d)\n[1]：自动  [-1]：手动\n", T.extraction(i)));
            if(~isempty(enter_content))
                T.extraction(i) = enter_content;
            end
        end
        if(T.extraction(i) == 1) % 用自动提取结果
            Q = Q_auto;
            TIFF_5 = TIFF_3;
            Item.extract_str = 'Extraction: Auto';
        elseif(T.extraction(i) == -1) % 用手动提取结果
            Q = Q_manual;
            TIFF_5 = TIFF_4;
            Item.extract_str = 'Extraction: Manual';
        end
        Item.emission_str = sprintf('Emission: %.2f±%.2f(t/h)', Q.ton_h, Q.ton_h_u);
        cprintf('*black', [Item.emission_str, '\n']);
        if(~ismember(Item.plume_quality, [1, 0.7])) % 如果不属于可定量
            Item.emission_str = 'Emission: Unquantified';
        end
        fprintf("风速不确定度：%.2f\n", Q.wind_u);
        fprintf("反演不确定度：%.2f\n", Q.inversion_u);

        %% 后处理
        TIFF_5.gaussiandata = imgaussfilt(TIFF_5.finaldata, sigma);
        %% 输出TIFF图像
        outtiffname = fullfile(Item.path.result_path, [prefix, Item.foldername, '_spa.tif']);
        [cropped_plume, plume_bounds] = SaveCroppedTiffData(TIFF_5, 'gaussiandata', outtiffname); % 保存裁剪后的羽流（空间域高斯）
        %% 计算分析
        % 计算图像统计信息
        pan_img = mean(TIFF_3.cutrgb, 3);
        pan_img(pan_img == 0) = nan;
        stat = calculateOutlierRate(pan_img);
        fprintf("背景离群度：%.2f%%\n", 100 * stat.ratio_outlier);
        % 计算羽流强度与不确定度之比
        plume_intensity = max(TIFF_5.finaldata(:)); % 计算提取出来的去噪羽流的最高强度
        ratio_plume_bkg = plume_intensity / TIFF_5.finaldata_u; % 估计比例
        fprintf("羽流比强度: %.2f\n", T.ratio_plume_bkg(i));
        % 计算羽流方向与实际风向之间的比例偏差
        [wind_direction_offset, wind_deviation] = calculateWindPlumeDeviation(TIFF_5, Wind);

        %%  绘制Profile图
        if(strcmp(IfPlot, 'on'))
            fig = PlotProfile(Item, TIFF_5, Wind, pan_img, stat); % 绘制profile
            figname = fullfile(Item.path.result_path, [prefix, Item.foldername, '_Profile.png']);
            print(fig, figname, '-dpng', '-r300')
        end

        %% 赋值
        T.plume_bounds{i} = plume_bounds; % 羽流边界
        T.emission(i) = Q.ton_h; % 排放强度
        T.emission_uncertainty(i) = Q.ton_h_u; % 排放不确定度
        T.quantification_error(i) = abs(Q.ton_h / T.emission_epa(i) - 1); % 与EPA清单排放量的绝对比例偏差
        T.wind_speed(i) = Wind.speed; % 风速
        T.wind_direction(i) = Wind.direction; % 风向
        T.wind_speed_std(i) = Wind.speed_std; % 风速不确定度
        T.wind_fit(i) = T.emission_epa(i) * Q.L / IME.data / 3.6; % 风速拟合
        T.plume_length(i) = Q.L;
        T.vertical_mass_atmosphere(i) = IME.data;
        T.inversion_std(i) = TIFF_5.finaldata_u; % 反演不确定度
        T.uncertainty_wind_part(i) = Q.wind_u; % 风速部分的不确定度
        T.uncertainty_ime_part(i) = Q.inversion_u; % 反演部分的不确定度
        T.ratio_rgb_outlier(i) = stat.ratio_outlier;
        T.ratio_plume_bkg(i) = ratio_plume_bkg;
        T.wind_direction_offset(i) = wind_direction_offset;
        T.wind_deviation(i) = wind_deviation;
    end
end
% mydisp(['自动提取景数为：',num2str(sum(T.extraction==1))]);
% %% 有效风速拟合
% FitUeffPlot(T); % 有效风速拟合并绘图
%% 绘制散点图
outputfolder = 'C:\Users\huayi\Desktop\NC论文\Figure1';
Plot_scatter_all(T, outputfolder)
%% 写入表格
writetable(T, data_file, 'Sheet', sheetName); % 写入数据到指定的Excel文件的工作表中
disp('Data has been updated in the specified sheet in the Excel file.'); % 显示确认信息
%% 检查反演结果
% ImageBrowser(T, len, Root_Path1, Root_Path2, prefix)
% %% 复制所有结果
% dstfolder = 'D:\huayi\Documents\研究资料\20250116 USplumes';
% for i = 1:len
%     figname = fullfile(Root_Path2, T.foldername{i}, T.facility{i}, [prefix, T.foldername{i}, '_Profile.png']);
%     % 检查文件是否存在
%     if exist(figname, 'file') == 2
%         % 构建目标文件的路径（可以保持原文件名或者根据需求重命名）
%         dstfile = fullfile(dstfolder, [prefix, T.foldername{i}, '_Profile.png']);
% 
%         % 复制文件到目标文件夹
%         copyfile(figname, dstfile);
%     else
%         disp(['File does not exist: ' figname]);
%     end
% end