clc;clear;
%% 加载路径
fun_path = genpath('functions'); addpath(fun_path)
%%
foldername = 'E:/2024.10 EMIT电厂数据';
%% 设置目标
facility = 'Compare';
files = Dirfiles(foldername, facility);
fpath = "E:/2024.10 EMIT电厂数据/Compare";

% filename_RAD = fullfile(fpath, 'EMIT_L1B_RAD_001_20240225T161125_2405611_012.nc'); % 无排放
% filename_RAD = fullfile(fpath, 'EMIT_L1B_RAD_001_20240326T201210_2408613_012.nc'); % 有一点排放
% filename_RAD = fullfile(fpath, 'EMIT_L1B_RAD_001_20240801T173015_2421412_018.nc'); % 有排放 range 30
% filename_RAD = fullfile(fpath, 'EMIT_L1B_RAD_001_20241018T184038_2429212_008.nc'); % 有排放 水面也有一点影响
%% 指定文件路径
len = length(files);
for j = 1:len
    close all;
    %% 读取数据
    filename = files(j).name;
    filename_RAD = fullfile(foldername, facility, filename);
    finfo_RAD = ncinfo(filename_RAD);
    Radiance = ncread(filename_RAD, 'radiance'); % 'uW/cm^2/SR/nm'
    Radiance = permute(Radiance, [3 2 1]); % 维度调换
    Radiance = Radiance * 10; % 'W/m^2/SR/um'

    % bandmean = mean(reshape(Radiance, [], 285));
    % plot(bandmean) % 绘制波段平均值
    %% 读取中心波长和半峰全宽
    ncid = netcdf.open(filename_RAD, 'NC_NOWRITE'); % 使用 netcdf.open 函数打开文件
    groupName = 'sensor_band_parameters'; % 查看文件中的组
    groupId = netcdf.inqNcid(ncid, groupName); % 获取组的 ID
    varName = 'wavelengths';
    varId = netcdf.inqVarID(groupId, varName); % 获取变量的 ID
    wavelength = netcdf.getVar(groupId, varId); % 读取变量的数据
    varName = 'fwhm';
    varId = netcdf.inqVarID(groupId, varName); % 获取变量的 ID
    fwhm = netcdf.getVar(groupId, varId); % 读取变量的数据
    FWHM_VNIR = [wavelength(1:84), fwhm(1:84)];
    FWHM_SWIR = [wavelength(85:end), fwhm(85:end)];

    %% 读取经纬度
    groupName = 'location'; % 查看文件中的组
    groupId = netcdf.inqNcid(ncid, groupName); % 获取组的 ID
    varName = 'lon';
    varId = netcdf.inqVarID(groupId, varName); % 获取变量的 ID
    lon = netcdf.getVar(groupId, varId); % 读取变量的数据
    lon = lon';
    varName = 'lat';
    varId = netcdf.inqVarID(groupId, varName); % 获取变量的 ID
    lat = netcdf.getVar(groupId, varId); % 读取变量的数据
    lat = lat';
    %%
    load("data/FWHM_SWIR/EMIT.mat");
    load("data/CO2_k/CO2_k_v6.mat", "k_EMIT");
    %% 读取反演波段
    datafile = 'data/Band_CO2_v1.xlsx';
    Tabledata = readtable(datafile, 'VariableNamingRule', 'preserve');
    bandindex = GetBandIndex(Tabledata, 'EMIT', FWHM_SWIR);
    %% 获取RGB影像
    RGB_dst = [640, 550, 460]; % 红色640nm  绿色550nm   蓝色460nm
    RGB_index = GetWavIndex(FWHM_VNIR, RGB_dst);
    rgb = Radiance(:, :, RGB_index);
    rgb_s = ImgStretch(rgb, 0.5, 99.5);

    %% 匹配滤波算法
    bandlength = FWHM_SWIR(bandindex, 1);
    IsAC = false; % 不使用albedo factor矫正
    swir = Radiance(:, :, 85:end);
    swir_ = swir(:, :, bandindex);
    swir_ = PreProcessSWIR(swir_);
    swir_2 = ConvertUnit(swir_, bandlength);
    %%
    k_ = k_EMIT(bandindex, 3);
    result = RMF(swir_2, k_, IsAC);
    %%
    % facility_lat = 32.2597;
    % facility_lon = -94.5703;
    % % 计算每个点与给定经纬度的距离
    % distances = sqrt((lat - facility_lat).^2 + (lon - facility_lon).^2);
    % % 找到最小距离的索引
    % [min_dist, linear_idx] = min(distances(:));
    % % 转换为行列坐标
    % [rowIndex, colIndex] = ind2sub(size(distances), linear_idx);
    % % 输出结果
    % fprintf('最近的点坐标在 (%d, %d)，经纬度为 (%.4f, %.4f)\n', rowIndex, colIndex, lat(rowIndex, colIndex), lon(rowIndex, colIndex));
    %%
    [ROW, COL, ~] = size(Radiance);
    range = 50;
    % rindexs = max(1, rowIndex - range + 1): min(rowIndex + range, ROW);
    % cindexs = max(1, colIndex - range + 1): min(colIndex + range, COL);
    rindexs = 1:ROW;
    cindexs = 1:COL;
    cutdata = result(rindexs, cindexs);
    cutbkg = rgb_s(rindexs, cindexs, :);
    figure
    ShowImage(cutbkg, cutdata, '')
    colormap(gca, mycolorbar(4));  % 使用颜色映射
    clim([0 max(cutdata(:)) * 0.8]);
    % clim([0 18])
    axis square;
    colorbar;
    toc;
end
%%
function files = Dirfiles(foldername, facility)
% 定义目录路径
directory = fullfile(foldername, facility);

% 获取该目录下的所有文件信息
allFiles = dir(fullfile(directory, '*RAD*'));

% 创建逻辑索引以排除 "." 和 ".."
validFilesIdx = ~ismember({allFiles.name}, {'.', '..'});

% 应用逻辑索引过滤文件列表
files = allFiles(validFilesIdx);
end