%% 清理变量
% clear; 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 = 'E:\2024.12 中国电厂数据';
Root_Path2 = 'E:\2024.12 中国电厂结果';
% Root_Path1 = 'F:\2024.7 美国原始数据';
% Root_Path2 = 'F:\2024.7 美国核查结果';
datafile = 'data/Togtoh.xlsx'; sheetName = '定量';

% datafile = 'data/Xinjiang-multisource.xlsx'; sheetName = 'Sheet1';
% datafile = 'data/美国碳源反演.xlsx'; sheetName = '2024.12';
%% 读取数据
GasType = 'CH4';
datafile2 = 'data/Band_CO2_v1.xlsx';
Tabledata = readtable(datafile2, 'VariableNamingRule', 'preserve', 'TreatAsEmpty', '');
fpath_k = 'data/CO2_k/CO2_k_v25.1.2.mat'; % 吸收系数查找表路径
fpath_wv = 'data/WaterVaporTable.mat';
load(fpath_wv, 'WaterVaporTable');
T = readtable(datafile, 'Sheet', sheetName, 'VariableNamingRule', 'preserve');
%% 筛选有排放的景号
% removed_img = cell2mat(table2cell(T(:, '排放'))) ~= 1; % 只保留1
% T(removed_img, :) = [];
T_cell = table2cell(T);
satellite = T{:, 'satellite'};
fprintf("当前选择的景数为：%d\n", size(T, 1));
%% 初始化
len = size(T, 1);
Emission = T{:,{'反演结果','不确定度'}};
windspeed = zeros(len, 1); % 实际风速
extract = T{:, '提取'};
IfPlot = 'on';
skip = false;
%% 加载气体吸收系数
load(fpath_k, 'k_GF5');load(fpath_k, 'k_GF5A');load(fpath_k, 'k_GF5B');load(fpath_k, 'k_ZY1E');load(fpath_k, 'k_ZY1F');
k_unit.('GF5') = k_GF5;k_unit.('GF5A') = k_GF5A;k_unit.('GF5B') = k_GF5B;k_unit.('ZY1E') = k_ZY1E;k_unit.('ZY1F') = k_ZY1F;

Altitude = 1067 / 1e3;
PARM_INPUT.HEIGHT = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0];
[~, idx_h] = min(abs(Altitude - PARM_INPUT.HEIGHT));
col = idx_h + 2;

%% 循环处理
for i = 1:len
    close all % 清除图窗
    tic;
    %% 读取当前排放设施的坐标以及时间信息
    Item = GetFacilityDataNew(T, i);
    PrintBasicInfo(Item);
    %% 数据路径
    origin_path = fullfile(Root_Path1, Item.region, Item.filename);
    process_path = fullfile(Root_Path2, Item.foldername); mymkdir(process_path) % 若不存在就新建文件夹
    result_path = fullfile(process_path, Item.facility); mymkdir(result_path)
    tempdata_path = fullfile(process_path, 'INPUT'); mymkdir(tempdata_path) % 数据缓存路径
    %% 读取文件名
    flist = ReadFiles(origin_path); % 读取数据
    xmlStruct = ReadXmlData(flist.xml); % 读取xml文件
    SatelliteID = xmlStruct.SatelliteID;
    soa = str2double(xmlStruct.SolarAzimuth);
    soz = str2double(xmlStruct.SolarZenith);
    %%
    outputname = [GasType, '_', Item.foldername, '.tif'];
    tiffname_RMF = fullfile(process_path, [GasType, '_', Item.foldername, '_RMF.tif']);
    % tiffname_ILMF = fullfile(process_path, [GasType, '_', Item.foldername, '_ILMF.tif']);
    tiffname_RGB = fullfile(process_path, 'RGB.tif');
    % a=readgeoraster(tiffname_RGB);
    % figure('Position', [100 0 1000 800])
    % imagesc(a);
    % axis square
    % disp(T{i,'反演结果'});
    % continue;
    %% 数据读取
    if(~exist(tiffname_RMF, 'file') || ~exist(tiffname_RGB, 'file'))
        FWHM_SWIR = Read_Raw(flist.swir_fwhm, 2); % FWHM 包含central wavelength 中心波长 和 full width at half maxima 半宽
        bandindex = GetBandIndex(Tabledata, SatelliteID, FWHM_SWIR);
        swir_inv = ReadSWIR(tempdata_path, flist, bandindex);
        cw_inv = FWHM_SWIR(bandindex, 1);
        %% RGB
        FWHM_VNIR = Read_Raw(flist.vnir_fwhm, 2);
        WV = ReadDataWaterVapor(tempdata_path, flist, FWHM_VNIR); % 读取水汽透过率
        WaterVapor = GetWaterVapor(WV, WaterVaporTable);
        rgb = ReadDataRGB(tempdata_path, flist, FWHM_VNIR); % 读取RGB
        %% 气体吸收系数
        % col = ChooseKtableCol(xmlStruct, WaterVapor); % 选择数据

        k_ = k_unit.(SatelliteID);
        k_ = k_(bandindex, col);
        amf = double(1 + 1 / cos(soz / 180 * pi));
        k_ = k_ * amf / 2;

        %% 匹配滤波算法
        IsAC = false; % 不使用albedo factor矫正
        swir_inv = PreProcessSWIR(swir_inv);
        swir_unit = ConvertUnit(swir_inv, cw_inv); % 单位转换为 W/(cm2×Sr×cm-1)
        result_RMF = RMF(swir_unit, k_, IsAC);
        % result_ILMF = ILMF(swir_unit, k_, IsAC);
        % mergedata = MergeData(rgb, result_ILMF);
        % TIFF_1 = ProcessTIFFdata(mergedata, Item, flist, 150);
        %         figure('Position', [ 20 50 600 600])
        %         ShowImage(TIFF_1.cutrgb, TIFF_1.cutdata, '结果');
        %         colorbar
        %         caxis([0 50])
        %% TIFF 数据
        [rgb, result_RMF] = CheckSize(rgb, result_RMF); % 检查尺寸是否一致
        % if(~exist(tiffname_RGB, 'file'))
        rgb = rgb / 256; rgb_s = ImgStretch(rgb, 0.5, 99.5); % 拉伸RGB色彩
        rgb_t = color_AHSI(rgb_s);
        ExportResultTIFF(rgb_t, flist.SWIR_rpb, tiffname_RGB);
        % end
        ExportResultTIFF(result_RMF, flist.SWIR_rpb, tiffname_RMF);
        % ExportResultTIFF(result_ILMF, flist.SWIR_rpb, tiffname_ILMF);
    end
    % continue;
    %% 获取TIFF
    range = 300;
    percent = 2; % 阈值提取的分位数
    thresh = 50; % 保留的相连像素数
    % TIFF_1 = GetTIFFData(tiffname_ILMF, tiffname_RGB, Item, range);    % 获取导出的TIFF反演数据
    TIFF_2 = GetTIFFData(tiffname_RMF, tiffname_RGB, Item, range);    % 获取导出的TIFF反演数据
    TIFF_2 = ThreshTIFF(TIFF_2, percent, thresh, range); % 羽流自动提取
    [TIFF_3, flag] = MaskPlume2(TIFF_2);
    if(~flag)
        % continue; % 可选
        fprintf("未探测到羽流！请在图中确认羽流点源位置\n");
        TIFF_2 = SelectPlumePosition(TIFF_2);
        [TIFF_3, flag] = MaskPlume2(TIFF_2);
    end
    if(strcmp(GasType, 'CO2'))
        mydisp(['反演不确定度为',num2str(round(TIFF_3.finaldata_u,0)),'ppm']);
    elseif(strcmp(GasType, 'CH4'))
        mydisp(['反演不确定度为',num2str(round(TIFF_3.finaldata_u * 1000,0)),'ppb']);
    end
    if(strcmp(IfPlot, 'on'))
        visualizeTIFF(TIFF_3)
    end
    % TIFF_3.cutdata = TIFF_1.cutdata;
    % TIFF_3.finaldata = TIFF_3.finalmask .* TIFF_1.cutdata;
    %% 风速
    [VMA, Wind] = MainFunctionIME(result_path, Item, flist, GasType); % IME方法所需的相关数据
    windspeed(i) = Wind.speed; % 风速
    %% 核算自动提取结果
    IME = calc_IME(TIFF_3, VMA, GasType);    % 二氧化碳整体质量增量 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));
    [Q_manual, TIFF_4, IME] = CutTIFFByHand(TIFF_3, VMA, Wind, result_path, skip, IfPlot, GasType);
    if(~skip) % 如果跳过
        extract(i) = input(sprintf("请选择提取结果：(当前值：%d)\n[1]：自动  [-1]：手动\n", extract(i)));
    end
    if(extract(i) == 1) % 用自动提取结果
        Q = Q_auto;
        TIFF_5 = TIFF_3;
    elseif(extract(i) == -1) % 用手动提取结果
        Q = Q_manual;
        TIFF_5 = TIFF_4;
    end
    fprintf("提取结果: ");cprintf('*black',sprintf('%.2f±%.2f(t/h)\n', Q.ton_h, Q.ton_h_u));
    %% 不确定度
    fprintf("风速不确定度：%.2f\n", Q.wind_u);
    fprintf("反演不确定度：%.2f\n", Q.inversion_u);
    Emission_1 = Q.ton_h; % 排放强度
    Emission_2 = Q.ton_h_u; % 排放不确定度
    Emission(i, :) = [Emission_1, Emission_2];
    %% 高斯后处理
    sigma = 0.5;
    TIFF_5.gaussiandata = imgaussfilt(TIFF_5.finaldata, sigma);
    %% 绘图
    fig = figure('Visible', IfPlot);
    showtitle = strrep(Item.foldername, '_', '\_');
    ShowImage(TIFF_5.cutrgb, TIFF_5.gaussiandata, showtitle);
    colorbar
    figname = [result_path,'\', Item.foldername, '.png'];
    saveas(fig, figname);
    %% 绘图
    fig = figure('visible', IfPlot);
    % showtitle = strrep(Item.foldername, '_', '\_');
    ShowImage(TIFF_5.cutrgb, TIFF_5.gaussiandata, '');
    cmap = mycolorbar(5);
    colormap(gca, cmap);  % 使用颜色映射
    % figname = ['C:\Users\huayi\Desktop\大唐托克托全年发电量\pic\', Item.foldername, '.png'];
    % % % print(fig, figname, '-dpng', '-r600');
    % saveas(fig, figname);
    %% 保存裁剪后的羽流
    outtiffname = [result_path,'\', Item.foldername, '_', Item.facility, '.tif'];
    cropped_tiff = SaveCroppedTiffData(TIFF_5, 'gaussiandata', outtiffname);
end
%% 写入表格
currentdata = [Emission, windspeed, extract];
columnname = {'反演结果', '不确定度', '风速', '提取'};
olddata = cell2mat(table2cell(T(:, columnname)));
zeroIndices = currentdata == 0; % 创建逻辑索引，指示 currentdata 中是否为零
newData = currentdata; % 选择非零值
newData(zeroIndices) = olddata(zeroIndices); % 对应位置的旧数据
T{:, columnname} = newData;
writetable(T, datafile, 'Sheet', sheetName);