GOCI二级叶绿素浓度数据简单处理

前置信息

按照老师布置的作业要求，选择某一天的GOCI数据作日变化图。我选择的是叶绿素数据，选择一块高值区一块低值区进行绘图。把过程简单记录一下。

首先有一些基础的东西需要知道：1、这颗卫星是对地静止卫星，覆盖范围只有不大的一块区域

2、每天从9：00-16：00，每小时会有一张影像，每日八张

3、二级数据分辨率300m

数据下载

普通下载（主要是官网）

首先，下载数据可以从ocean color网站下载，或者GOCI官网下载。

NASA Ocean Colorhttps://oceancolor.gsfc.nasa.gov/

해양위성센터https://kosc.kiost.ac.kr/index.nm

前者下载到的文件为nc格式，直接使用matlab读取失败，可能是我读取的方式不对。当然，这个nc文件是可以通过seadas打开并导出的。后者下载到的数据是he5格式，除了要下载对应数据外还要另外下载经纬度文件（最下面两个就是）。

不管先下载什么文件，初次下载都会自动下载一个插件，后续会使用这个插件进行文件下载，会很快，可以直接设置下载到哪。

通过点击图中左下角的GOCI图标右上角的加号进入下载流程，一步一步勾选点击就好。

如果只是想要单日影像的话，选择第一个方法即可。

下载的时段较长，可以选择第二个方法。

注意事项

下载时需要看清所下文件到底是什么内容，文件在网页上的排列方法是9：00的CODM、CHL等文件，然后是10：00的CODM、CHL等文件，依次到16：00.

爬虫批量下载

这个有同学已经实现了，我自己没这么做，可以参考这个博客GOCI L2数据批量下载_goci影像的批量下载-CSDN博客

简单处理

读取

每个下载到的文件解压缩之后均含有两个文件，一个是预览图，一个是包含数据的he5文件。以下提供一段代码来查看该文件的数据结构。

close all;clear;clc;

% 指定 HDF5 文件的名称
filename = "D:\GOCI\CHL\COMS_GOCI_L2A_GA_20180823021643.CHL.he5";

% 获取文件的信息
info = h5info(filename);

% 检查每个组
for i = 1:length(info.Groups)
    groupInfo = info.Groups(i);
    groupName = groupInfo.Name;

    % 打印组名
    fprintf('组名: %s\n', groupName);

    % 检查组内的数据集
    for j = 1:length(groupInfo.Datasets)
        datasetName = groupInfo.Datasets(j).Name;

        % 构建完整的数据集路径
        fullPath = strcat(groupName, '/', datasetName);

        % 读取数据集
        data = h5read(filename, fullPath);

        % 显示数据集的名称、数据和维度
        fprintf('  数据集名称: %s\n', datasetName);
        disp('  数据内容:');
        disp(data);
        fprintf('  维度: %s\n', mat2str(size(data)));
    end
end

数据结构乍一看很庞大，其实配合接下来的数据读取代码中的读取路径来看也不太复杂。

clc; clear;

% 指定文件夹路径
folder_path = 'D:\GOCI\CHL\';
lon_file = 'COMS_GOCI_L2P_GA_20110524031644.LON.he5';
lat_file = 'COMS_GOCI_L2P_GA_20110524031644.LAT.he5';

% 读取经度和纬度数据
LON = h5read(fullfile(folder_path, lon_file), '/HDFEOS/GRIDS/Image Data/Data Fields/Longitude Image Pixel Values');
LAT = h5read(fullfile(folder_path, lat_file), '/HDFEOS/GRIDS/Image Data/Data Fields/Latitude Image Pixel Values');

% 保存经纬度数据为MAT格式
save(fullfile(folder_path, 'LON.mat'), 'LON');
save(fullfile(folder_path, 'LAT.mat'), 'LAT');

% 循环处理每个CHL文件
for i = 1:8
    % 构建文件名
    chl_file = sprintf('COMS_GOCI_L2A_GA_201808230%d1643.CHL.he5', i-1);

    % 读取CHL数据
    CHL_data = h5read(fullfile(folder_path, chl_file), '/HDFEOS/GRIDS/Image Data/Data Fields/CHL Image Pixel Values');
    CHL_data(CHL_data == -999) = NaN; % 将-999替换为NaN
    CHL_data(CHL_data > 10) = NaN;    % 将大于10的值替换为NaN
    
    % 转换为double类型
    CHL_data = double(CHL_data);
    
    % 保存CHL数据为MAT格式
    save(fullfile(folder_path, sprintf('CHL_data%d.mat', i)), 'CHL_data');
end

disp('所有CHL和经纬度数据已保存为MAT格式');

上述代码将从官网下载的he5数据文件和lat、lon文件都读取成mat文件并保存到指定文件夹内。

总览图

为了展示方便，我将八景图像均进行了绘制。

clc; clear;

% 设定文件夹路径
folder_path = 'E:\新建文件夹\';

% 确保文件夹路径存在
if ~exist(folder_path, 'dir')
    mkdir(folder_path);
end

% 读取经纬度数据
load(fullfile(folder_path, 'LON.mat'), 'LON');
load(fullfile(folder_path, 'LAT.mat'), 'LAT');

% 定义时间数组
observation_times = {'09:00', '10:00', '11:00', '12:00', '13:00', '14:00', '15:00', '16:00'};

% 定义边框的经纬度范围
lonlat_boxes = {
    [117, 122, 122, 117, 117], [35, 35, 40, 40, 35];  % 第一个框的经度和纬度
    [125, 130, 130, 125, 125], [25, 25, 30, 30, 25]   % 第二个框的经度和纬度
};

% 循环绘制每个图像
for i = 1:8
    % 读取CHL数据
    chl_data_file = sprintf('CHL_data%d.mat', i);
    load(fullfile(folder_path, chl_data_file), 'CHL_data');

    % 绘图
    figure;
    pcolor(LON, LAT, CHL_data); % 使用pcolor函数来绘制数据
    shading interp; % 插值显示，去除网格线并平滑
    colormap(jet); % 使用jet颜色映射
    caxis([0, 2]); % 设置颜色范围为0到2

    % 绘制边框
%     hold on;
%     for j = 1:size(lonlat_boxes, 1)
%         plot(lonlat_boxes{j, 1}, lonlat_boxes{j, 2}, '-k', 'LineWidth', 1);
%     end
%     hold off;

    % 添加颜色条
    hcb = colorbar;
    title(hcb, 'mg/m^3'); % 为颜色条设置单位
    
    % 设置图像标题和坐标轴标签
    title(sprintf('CHL Concentration - %s', observation_times{i}));
    xlabel('Longitude');
    ylabel('Latitude');

    % 调整图像大小
    axis tight;

    % 保存图像
    img_file_name = fullfile(folder_path, sprintf('CHL_Image_%s.png', observation_times{i}));
    print(gcf, img_file_name, '-dpng', '-r300');  % '-dpng'指定保存为PNG格式，'-r300'指定分辨率为300 dpi

    % 关闭当前图窗
    close(gcf);
end

disp('绘图完成并保存于指定文件夹。');

其中可以选择将感兴趣的区域框起来。效果如下

平均及折线图

我将两个区域内的有效值点数目进行了统计，并将有效值进行了平均处理，从而得到所需日变化数据，存储于execl文件中。以下是计算一个区域数值的代码，如果要计算其他的需要手动更改经纬度范围。

% 清空工作区
clear;
close all;
clc;

% 创建“picture”子文件夹（如果不存在）
folder_path = 'D:\GOCI\CHL\';  % 修改为您的文件夹路径
if ~exist(folder_path, 'dir')
    mkdir(folder_path);
end

% 初始化存储结果的数组
CHL_averages = zeros(8, 1);  % 存储平均值
CHL_valid_points = zeros(8, 1);  % 存储有效值点数

% 指定经纬度范围
lat_range = [35, 40];
lon_range = [117, 122];

% 循环处理每个文件
for i = 1:8
    % 读取数据
    CHL_data_filename = sprintf('CHL_data%d.mat', i);
    LAT_filename = 'LAT.mat';  % 经纬度文件可能对于所有数据文件都是相同的
    LON_filename = 'LON.mat';
    load(fullfile(folder_path, CHL_data_filename), 'CHL_data');
    load(fullfile(folder_path, LAT_filename), 'LAT');
    load(fullfile(folder_path, LON_filename), 'LON');

    % 找出给定经纬度范围内的索引
    valid_indices = LAT >= lat_range(1) & LAT <= lat_range(2) & ...
                    LON >= lon_range(1) & LON <= lon_range(2);

    % 提取对应的CHL数据
    CHL_data_region = CHL_data(valid_indices);

    % 计算平均值和有效值点数
    CHL_averages(i) = mean(CHL_data_region(~isnan(CHL_data_region)));
    CHL_valid_points(i) = sum(~isnan(CHL_data_region));
end

% 创建表格并写入数据
T = table((1:8)', CHL_averages, CHL_valid_points, ...
          'VariableNames', {'FileNumber', 'AverageCHL', 'ValidDataPoints'});

% 写入Excel文件
excel_filename = fullfile(folder_path, 'CHL_Data_Summary.xlsx');
writetable(T, excel_filename);

disp('数据已经写入文件：');
disp(excel_filename);

绘制折线图：

% 清空工作区
clear;
close all;
clc;

% 读取Excel文件
folder_path = 'D:\GOCI\CHL\';  % 指定文件夹路径
excel_filename = fullfile(folder_path, 'high.xlsx');
T = readtable(excel_filename);

% 定义时间数组
times = {'09:00', '10:00', '11:00', '12:00', '13:00', '14:00', '15:00', '16:00'};

% 绘制图像
figure;
[ax, h1, h2] = plotyy(1:8, T.AverageCHL, 1:8, T.ValidDataPoints);
hold(ax(1), 'on');
hold(ax(2), 'on');
plot(ax(1), 1:8, T.AverageCHL, 'o-b', 'LineWidth', 2, 'MarkerSize', 8, 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot(ax(2), 1:8, T.ValidDataPoints, '*--r', 'LineWidth', 2, 'MarkerSize', 8, 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'r');

% 设置横轴标签为时间
set(ax(1), 'XTick', 1:8, 'XTickLabel', times);
set(ax(2), 'XTick', 1:8, 'XTickLabel', times);

% 设置纵轴标签
ylabel(ax(1), 'Average CHL (mg/m^3)');
ylabel(ax(2), 'Valid Data Points');

% 设置图表标题和轴标签
title('CHL Concentration and Valid Data Points Over Time');
xlabel('Time');

% 设置图例
legend([h1, h2], 'Average CHL', 'Valid Data Points');

% 调整图像大小
axis(ax(1), 'tight');
axis(ax(2), 'tight');

% 关闭图形保持状态
hold(ax(1), 'off');
hold(ax(2), 'off');

效果如下