% 清除工作区和命令窗口
clear; clc; close all;
%% 1. 读取NetCDF数据
ncFilePath = 'G:\DC\CERES_EBAF_Ed4.2_Subset_200003-202407 (1).nc'; % NC文件路径
varName = 'toa_net_all_mon'; % 待读取的变量名（TOA净辐射）
timeIndex = 1; % 时间索引（读取第1个时间点的数据）
% 获取NC文件信息
ncInfo = ncinfo(ncFilePath);
% 查找目标变量信息
varInfo = ncInfo.Variables(strcmp({ncInfo.Variables.Name}, varName));
if isempty(varInfo)
error('变量 %s 不存在于NC文件中', varName);
end
varDims = {varInfo.Dimensions.Name}; % 变量的维度名称
varSize = varInfo.Size; % 变量的尺寸
% 读取经纬度数据
lon = ncread(ncFilePath, 'lon'); % 经度（一维）
lat = ncread(ncFilePath, 'lat'); % 纬度（一维）
if ~isvector(lon) || ~isvector(lat)
error('经纬度数据必须为一维数组');
end
lonSize = length(lon); % 经度格点数
latSize = length(lat); % 纬度格点数
% 将经度范围转换为[-180, 180]
lon(lon > 180) = lon(lon > 180) - 360;
% 读取目标变量数据
data = ncread(ncFilePath, varName);
% 确定纬度、经度、时间在变量维度中的索引
latIdx = find(strcmp(varDims, 'lat'));
lonIdx = find(strcmp(varDims, 'lon'));
timeIdx = find(strcmp(varDims, 'time'));
% 构建索引（仅保留指定时间点，其他维度全取）
index = cell(1, length(varSize));
for i = 1:length(index)
if i == timeIdx
index{i} = timeIndex; % 时间维度取指定索引
else
index{i} = 1:varSize(i); % 其他维度全取
end
end
% 提取指定时间点的数据切片
dataSlice = data(index{:});
% 调整维度顺序为[纬度, 经度]（确保与经纬度网格匹配）
if latIdx ~= 1 || lonIdx ~= 2
dataSlice = permute(dataSlice, [latIdx, lonIdx]);
end
% 去除单维度（若有）
dataSlice = squeeze(dataSlice);
% 检查数据尺寸是否与经纬度匹配
if ~isequal(size(dataSlice), [latSize, lonSize])
error('数据切片尺寸异常: %s（预期[%d, %d]）', ...
num2str(size(dataSlice)), latSize, lonSize);
end
%% 2. 绘图准备：生成经纬度网格（二维）
% 将一维经纬度转换为二维网格，与dataSlice（[latSize, lonSize]）维度匹配
[LON_grid, LAT_grid] = meshgrid(lon, lat);
%% 3. 绘制原始数据（南半球立体投影）
% 设置立体投影，聚焦南半球（-90°至-60°纬度）
axesm('stereo', 'MapLatLimit', [-90 -60], 'MapLonLimit', [-180 180]); % 显示边框和网格
% 绘制数据（使用二维经纬度网格匹配数据维度）
geoshow(LAT_grid, LON_grid, dataSlice, 'DisplayType', 'surface');
title('原始数据（TOA净辐射）'); % 子图标题
colorbar; % 显示色标
% 调整色标范围（去除1%和99%极端值，增强可视化效果）
caxis([prctile(dataSlice(~isnan(dataSlice)),1), prctile(dataSlice(~isnan(dataSlice)),99)]);
