%% 清理变量
clear; close all;
%% 加载路径
fun_path = genpath('functions'); addpath(fun_path)
%% 设置python环境
% 需要先启动python gdal环境
% pyenv('Version', 'C:\Users\huayi\.conda\envs\gdal_env\python.exe');
% pyenv('Version', 'D:\Program Files\Python311\python.exe');
% % 测试gdal模块是否正常
try
    py.importlib.import_module('osgeo.gdal');
    disp('GDAL模块测试成功');
catch ME
    disp('导入GDAL失败');
    disp(ME.message);
    return;
end
%% 反演气体
GasType = 'CO2'; % or CO2
% cmax = 60;
% cmax = 1;
method = 'RMF'; % or ILMF
%% 读取数据
datafile = ['data/Band_', GasType, '_v1.xlsx'];
Tabledata = readtable(datafile, 'VariableNamingRule', 'preserve');
fpath_k = ['data/', GasType, '_k/', GasType, '_k_v25.1.2.mat']; % 吸收系数查找表路径
fpath_wv = 'data/WaterVaporTable.mat';
load(fpath_wv, 'WaterVaporTable');
%% 加载气体吸收系数
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;

%% 设置路径
% Root_Path1 = 'E:\2024.11 全球油田数据';
% Root_Path2 = 'E:\2024.11 全球油田结果';
% Region = 'VEN';
% Root_Path1 = 'E:/2024.10 甲烷比赛数据';
% Root_Path2 = 'E:/2024.10 甲烷比赛结果';
% Region = 'US';
% Root_Path1 = 'E:/2024.10 EMIT电厂数据';
% Root_Path2 = 'E:/2024.10 EMIT电厂结果';
% Region = 'Martin Lake';
Root_Path1 = 'E:/2024.12 山东东营数据';
Root_Path2 = 'E:/2024.12 山东东营结果';
Region = 'shengli';
% Root_Path1 = 'F:\2025.1 加州野火数据';
% Root_Path2 = 'F:\2025.1 加州野火结果';
% Region = 'CA';
%% 获取文件夹信息
folders = dir(fullfile(Root_Path1, Region)); folders = folders([folders.isdir]); folders = folders(~ismember({folders.name}, {'.', '..'}));
len = size(folders, 1);
%%
Process_WV = false;
Process_Cloud = false;
%% 循环处理
for i = 1:len
    mydisp(i);
    filename = folders(i).name;
    mydisp(filename);
    %% date
    underscoreIndices = strfind(filename, '_');
    date = filename(underscoreIndices(4) + 1 : underscoreIndices(4) + 8);
    %% foldername
    dataname = filename(underscoreIndices(end) + 1 : end);
    foldername = [Region, '_', date, '_', dataname];
    %% 数据路径
    origin_path = fullfile(Root_Path1, Region, filename);
    result_path = fullfile(Root_Path2, foldername); mymkdir(result_path) % 若不存在就新建文件夹
    datapath = fullfile(result_path, 'INPUT'); mymkdir(datapath) % 输入数据路径
    %% 读取文件名
    flist = ReadFiles(origin_path); % 读取数据
    xmlStruct = ReadXmlData(flist.xml); % 读取xml文件
    SatelliteID = xmlStruct.SatelliteID;
    %% 输出名称
    outputname = [GasType, '_', foldername, '_', method, '.tif'];
    outputname2 = [GasType, '_', foldername,'_', method, '_thresh', '.tif'];
    tiffname_RMF = fullfile(result_path, outputname);
    tiffname_RMF_thr = fullfile(result_path, outputname2);
    tiffname_RGB = fullfile(result_path, ['RGB_', foldername, '.tif']);
    %% 数据读取
    FWHM_SWIR = Read_Raw(flist.swir_fwhm, 2);
    FWHM_VNIR = Read_Raw(flist.vnir_fwhm, 2);
    CW_FWHM = cat(1, FWHM_VNIR, FWHM_SWIR); % CW:central wavelength 中心波长；FWHM:full width at half maxima 半宽
    %% 提取波段
    RGB_dst = [640, 550, 460]; % 红色640nm  绿色550nm   蓝色460nm
    % SMK_dst = [679, 412, 2130]; % 烟尘分布
    RGB_index = GetWavIndex(CW_FWHM, RGB_dst);
    WV_dst = [935, 870]; % 935nm / 870nm
    WV_index = GetWavIndex(CW_FWHM, WV_dst);

    series = xmlStruct.SatelliteID(1:2); % 卫星系列（GF/ZY）
    if(strcmp(series, 'GF'))
        Cloud_index = [11:20, 30:60, 192, 270:272];
    elseif(strcmp(series, 'ZY'))
        Cloud_index = [5:10, 15:30, 96, 135:136];
    else
        fprintf("请检查xmlStruct.SatelliteID是否准确\n");
    end

    Inv_index = GetBandIndex(Tabledata, SatelliteID, CW_FWHM);

    % 合并所需波段
    if(Process_Cloud)
        if(Process_WV)
            Index_all = [RGB_index, Inv_index, WV_index, Cloud_index];
        else
            Index_all = [RGB_index, Inv_index, Cloud_index];
        end
        Index_all = unique(Index_all); % 去除重复波段
        [~, rho_idx] = ismember(Cloud_index, Index_all); % Cloud索引在all索引中的索引
    else
        if(Process_WV)
            Index_all = [RGB_index, Inv_index, WV_index];
        else
            Index_all = [RGB_index, Inv_index];
        end
        Index_all = unique(Index_all);
    end
    [~, rgb_idx] = ismember(RGB_index, Index_all); % rgb索引在all索引中的索引
    [~, wv_idx] = ismember(WV_index, Index_all); % rgb索引在all索引中的索引
    [~, inv_idx] = ismember(Inv_index, Index_all); % inv索引在all索引中的索引
    %% 栅格数据
    mydisp('读取VNIR&SWIR数据...')
    tic;
    vnir = readgeoraster(flist.VNIR); % 1.3GB
    swir = readgeoraster(flist.SWIR);
    t = toc;
    mydisp('读取完成')
    fprintf('用时%.2fs\n', t);
    Radiance = MergeData(vnir, swir); % 2.6GB
    clear swir vnir % 清理不需要的变量
    Radiance = Radiance(:,:,Index_all); % 0.8GB
    vn_rad_cal = Read_Raw(flist.vnir_rad, 1); % 辐射定标系数
    sw_rad_cal = Read_Raw(flist.swir_rad, 1); % 辐射定标系数 radiance calibration
    rad_cal = [vn_rad_cal; sw_rad_cal];
    rad_cal = rad_cal(Index_all);
    Radiance = Rad_Calibration(Radiance, rad_cal); % 辐射定标  单位：W·m-2·sr-1·μm-1
    fprintf("内存占用：%.2fGB\n",sum([whos().bytes])/1e9);
    %% 表观反射率
    if(Process_Cloud)
        series = xmlStruct.SatelliteID(1:2); % 卫星系列（GF/ZY）
        if(strcmp(series, 'GF'))
            load('data/mat_data_GF5B/Esun.mat', 'Esun');
        elseif(strcmp(series, 'ZY'))
            load('data/mat_data_ZY1E/Esun.mat', 'Esun');
        else
            fprintf("请检查xmlStruct.SatelliteID是否准确\n");
        end
        rho = Radiance2Reflectance(Radiance(:, :, rho_idx), xmlStruct, Esun(Cloud_index));
    end
    %% 反演窗口
    swir_inv = Radiance(:, :, inv_idx);
    cw_inv = CW_FWHM(Inv_index, 1);
    %% RGB
    rgb = Radiance(:, :, rgb_idx);
    rgb_s = ImgStretch(rgb, 0.5, 99.5);
    rgb_t = color_AHSI(rgb_s);
    % rgb = ReadDataRGB(datapath, flist, FWHM_VNIR); % 读取RGB

    %% Water Vapor
    if(Process_WV)
        WV_data = Radiance(:, :, wv_idx);
        WV = WV_data(:,:,1) ./ WV_data(:,:,2);
        % WV = ReadDataWaterVapor(datapath, flist, FWHM_VNIR); % 计算水汽透过率
        WaterVapor = GetWaterVapor(WV, WaterVaporTable);
    end
    %% Cloud
    if(Process_Cloud)
        % 判断卫星是高分系列还是资源系列
        series = xmlStruct.SatelliteID(1:2); % 卫星系列（GF/ZY）
        if(strcmp(series, 'GF'))
            index = {11:20; 30:60; 192; 270:272};
            load('data/mat_data_GF5B/Esun.mat', 'Esun');
        elseif(strcmp(series, 'ZY'))
            index = {5:10; 15:30; 96; 135:136};
            load('data/mat_data_ZY1E/Esun.mat', 'Esun');
        else
            fprintf("请检查xmlStruct.SatelliteID是否准确\n");
        end
        % 合并为对应的一个宽波段数据
        % T1波段 433~471nm (第11~第20波段) % T2波段 510~640nm (第30~第60波段)
        % T3波段 1350~1360nm (第192波段)   卷云波段(1360~1390 nm) % T4波段 2000nm (第270~第272波段)

        T = cell(4, 1); % 初始化
        rho_CW_FWHM = CW_FWHM(Cloud_index, :);
        for idx_cloud = 1:4
            [~, T_idx] = ismember(index{idx_cloud}, Cloud_index); % Ti索引在Cloud索引中的索引
            T{idx_cloud} = mergebanddata(rho(:, :, T_idx), rho_CW_FWHM(T_idx, 2));

        end
        % 计算云掩膜并编码
        c1 = double((T{1} > 0.3) & (T{2} > 0.3)); % 厚云
        c2 = double((T{3} > 0.04) & (T{1} > 0.15)); % 薄云
        c3 = double((T{1} ./ T{4} > 7.5) & (T{4} ./ T{3} < 1)); % 云影
        Cloud = uint8(c1 * 4 + c2 * 2 + c3);
        Cloud = Cloud_post(Cloud, 15); % 后处理
        % Cloud = DetectCloud(datapath, flist, xmlStruct); % 计算云掩膜
    end

    %% 气体吸收系数
    % col = ChooseKtableCol(xmlStruct, WaterVapor); % 选择数据
    k_ = k_unit.(SatelliteID);
    k_ = k_(Inv_index - length(FWHM_VNIR), 3);
    
    %% 匹配滤波算法
    IsAC = false; % 不使用albedo factor矫正
    swir_inv = PreProcessSWIR(swir_inv);
    swir_unit = ConvertUnit(swir_inv, cw_inv); % 单位转换为 W/(cm2×Sr×cm-1)
    if(Process_Cloud)
        swir_final = RemoveCloudData(swir_unit, Cloud); % 去除云像元
    else
        swir_final = swir_unit;
    end
    %%
    if(strcmp(method, 'RMF'))
        result = RMF(swir_final, k_, IsAC);
    elseif(strcmp(method, 'ILMF'))
        result = ILMF(swir_final, k_, IsAC);
    end
    %%
    result_thr = Threshdata(result);
    %% TIFF 数据
    % [rgb, result] = CheckSize(rgb, result); % 检查尺寸是否一致
    if(~exist(tiffname_RGB, 'file'))
        ExportResultTIFF(rgb_t, flist.SWIR_rpb, tiffname_RGB); % 导出结果
    end
    ExportResultTIFF(result, flist.SWIR_rpb, tiffname_RMF);
    ExportResultTIFF(result_thr, flist.SWIR_rpb, tiffname_RMF_thr);
%%
end
