MERRA-2 ——M2T1NXAER——小时级气溶胶数据了解、下载、读取、重采样、碎碎念-优快云博客

本文链接：https://blog.youkuaiyun.com/qq_41359193/article/details/149078266

一、MERRA-2数据集信息查看——————————————————————————————————————

MERRA-2的主页网站
https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/

MERRA-2 包含的数据集介绍
https://gs6102dsc-goldsmr4.gesdisc.eosdis.nasa.gov/dods/
————————————————————————————————————————————————
本篇主要记录的数据集的信息查看，其他参数的信息也可以在此查看

M2T1NXAER——参数信息
https://gs6102dsc-goldsmr4.gesdisc.eosdis.nasa.gov/dods/M2T1NXAER.info

M2T1NXAER——详细的参数信息
https://gs6102dsc-goldsmr4.gesdisc.eosdis.nasa.gov/dods/M2T1NXAER.dds

https://gs6102dsc-goldsmr4.gesdisc.eosdis.nasa.gov/dods/M2T1NXAER.das

MERRA-2 M2T1NXAER下载地址：
https://disc.gsfc.nasa.gov/datasets?keywords=M2T1NXAER&page=1
————————————————————————————————————————————————

二、下载方式的选择——————————————————————————————————————-

主页网址进入https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/data_access/ 点击MDISC
在这里插入图片描述

下载选择方式一定要选择OpeNDAP，别问，一个字，快，两个字，超快，不到半小时下载一个参数的二十年，要是多台设备，可以一次性下载更多，不在赘述！！！记得参考使用IDM进行超高速下载！！！
在IDM中站点管理器设置网址为：https://opendap.earthdata.nasa.gov/ ///然后是nasa的账号和密码
在这里插入图片描述
————————————————————————————————————————————————

三、MATLAB查看数据集信息———————————————————————————————————————

其实就两行代码ncinfo() ncread() 我给复杂化了而已！！！

如下代码查看变量的信息，以及结构体的信息：

%% 读取MACA v2气候数据集中的单次散射反照率(SSA)及相关变量
clear;clc
%% 定义文件路径
nc_file = 'G:\MERRA-2\M2T1NXAER.5.12.4_MERRA2_200.tavg1_2d_aer_Nx.20000213.nc4.nc4';

%% 检查NC文件是否存在
if ~exist(nc_file, 'file')
    error('指定的NC文件不存在: %s', nc_file);
end

%% 获取并显示文件元数据
nc_info = ncinfo(nc_file);
fprintf('文件: %s\n', nc_file);
fprintf('格式: %s\n', nc_info.Format);
fprintf('变量数量: %d\n', length(nc_info.Variables));

%% 显示所有变量的原始信息
fprintf('\n==== 变量原始信息 ====\n');
for i = 1:length(nc_info.Variables)
    var = nc_info.Variables(i);
    fprintf('\n变量 #%d: %s\n', i, var.Name);
    disp(var);  % 直接显示变量结构体的完整内容
end

%% 显示所有变量的原始信息
fprintf('\n==== 变量原始变量结构体信息 ====\n');
for i = 1:length(nc_info.Variables)
    var = nc_info.Variables(i);
    fprintf('\n变量 #%d: %s\n', i, var.Name);
    disp(var);  % 直接显示变量结构体的完整内容
end

%% 显示全局属性
fprintf('\n==== 全局属性 ====\n');
for i = 1:length(nc_info.Attributes)
    att = nc_info.Attributes(i);
    fprintf('%s: %s\n', att.Name, mat2str(att.Value));
end

结果展示：可以不看，我为啥要贴上去，直接disp也可以哇，不更改了，直接跳过吧；；；

==== 变量原始信息 ====

变量 #1: lon
            Name: 'lon'
      Dimensions: [1×1 struct]
            Size: 576
        Datatype: 'double'
      Attributes: [1×5 struct]
       ChunkSize: 576
       FillValue: 'disable'
    DeflateLevel: 2
         Shuffle: 1


变量 #2: time
            Name: 'time'
      Dimensions: [1×1 struct]
            Size: 24
        Datatype: 'int32'
      Attributes: [1×8 struct]
       ChunkSize: 1
       FillValue: 'disable'
    DeflateLevel: 2
         Shuffle: 1


变量 #3: lat
            Name: 'lat'
      Dimensions: [1×1 struct]
            Size: 361
        Datatype: 'double'
      Attributes: [1×5 struct]
       ChunkSize: 361
       FillValue: 'disable'
    DeflateLevel: 2
         Shuffle: 1


变量 #4: TOTANGSTR
            Name: 'TOTANGSTR'
      Dimensions: [1×3 struct]
            Size: [576 361 24]
        Datatype: 'single'
      Attributes: [1×11 struct]
       ChunkSize: [144 91 1]
       FillValue: 'disable'
    DeflateLevel: 2
         Shuffle: 1

==== 全局属性 ====
History: 'Original file generated: Wed Jan  7 20:58:03 2015 GMT'
Comment: 'GMAO filename: d5124_m2_jan91.tavg1_2d_aer_Nx.20000213.nc4'
Filename: 'MERRA2_200.tavg1_2d_aer_Nx.20000213.nc4'
Conventions: 'CF-1'
Institution: 'NASA Global Modeling and Assimilation Office'
References: 'http://gmao.gsfc.nasa.gov'
Format: 'NetCDF-4/HDF-5'
SpatialCoverage: 'global'
VersionID: '5.12.4'
TemporalRange: '1980-01-01 -> 2016-12-31'
identifier_product_doi_authority: 'http://dx.doi.org/'
ShortName: 'M2T1NXAER'
GranuleID: 'MERRA2_200.tavg1_2d_aer_Nx.20000213.nc4'
ProductionDateTime: 'Original file generated: Wed Jan  7 20:58:03 2015 GMT'
LongName: 'MERRA2 tavg1_2d_aer_Nx: 2d,1-Hourly,Time-averaged,Single-Level,Assimilation,Aerosol Diagnostics'
Title: 'MERRA2 tavg1_2d_aer_Nx: 2d,1-Hourly,Time-averaged,Single-Level,Assimilation,Aerosol Diagnostics'
SouthernmostLatitude: '-90.0'
NorthernmostLatitude: '90.0'
WesternmostLongitude: '-180.0'
EasternmostLongitude: '179.375'
LatitudeResolution: '0.5'
LongitudeResolution: '0.625'
DataResolution: '0.5 x 0.625'
Source: 'CVS tag: GEOSadas-5_12_4'
Contact: 'http://gmao.gsfc.nasa.gov'
identifier_product_doi: '10.5067/KLICLTZ8EM9D'
RangeBeginningDate: '2000-02-13'
RangeBeginningTime: '00:00:00.000000'
RangeEndingDate: '2000-02-13'
RangeEndingTime: '23:59:59.000000'
build_dmrpp_metadata.build_dmrpp: '3.21.0-272'
build_dmrpp_metadata.bes: '3.21.0-272'
build_dmrpp_metadata.libdap: 'libdap-3.21.0-70'
build_dmrpp_metadata.configuration: ['' char(10) '# TheBESKeys::get_as_config()' char(10) 'AllowedHosts=^https?:\/\/' char(10) 'BES.Catalog.catalog.FollowSymLinks=Yes' char(10) 'BES.Catalog.catalog.RootDirectory=/tmp/tmp0815rv8u/' char(10) 'BES.Catalog.catalog.TypeMatch=dmrpp:.*\.(dmrpp)$;' char(10) 'BES.Catalog.catalog.TypeMatch+=h5:.*(\.bz2|\.gz|\.Z)?$;' char(10) 'BES.Data.RootDirectory=/dev/null' char(10) 'BES.LogName=./bes.log' char(10) 'BES.UncompressCache.dir=/tmp/hyrax_ux' char(10) 'BES.UncompressCache.prefix=ux_' char(10) 'BES.UncompressCache.size=500' char(10) 'BES.module.cmd=/usr/lib64/bes/libdap_xml_module.so' char(10) 'BES.module.dap=/usr/lib64/bes/libdap_module.so' char(10) 'BES.module.dmrpp=/usr/lib64/bes/libdmrpp_module.so' char(10) 'BES.module.fonc=/usr/lib64/bes/libfonc_module.so' char(10) 'BES.module.h5=/usr/lib64/bes/libhdf5_module.so' char(10) 'BES.module.nc=/usr/lib64/bes/libnc_module.so' char(10) 'BES.modules=dap,cmd,h5,dmrpp,nc,fonc' char(10) 'FONc.ClassicModel=false' char(10) 'FONc.NoGlobalAttrs=true' char(10) 'H5.EnableCF=false' char(10) 'H5.EnableCheckNameClashing=true' char(10) '']
build_dmrpp_metadata.invocation: 'build_dmrpp -c /tmp/bes_conf_boLQ -f /tmp/tmp0815rv8u//MERRA2_200.tavg1_2d_aer_Nx.20000213.nc4 -r /tmp/dmr__YhGLGa -u OPeNDAP_DMRpp_DATA_ACCESS_URL -M'
history: ['2025-07-02 13:22:32 GMT hyrax-1.17.1-219 https://opendap.earthdata.nasa.gov/collections/C1276812830-GES_DISC/granules/M2T1NXAER.5.12.4%3AMERRA2_200.tavg1_2d_aer_Nx.20000213.nc4.dap.nc4?dap4.ce=/TOTANGSTR;/time;/lat;/lon' char(10) '']
history_json: '[{"$schema":"https:\/\/harmony.earthdata.nasa.gov\/schemas\/history\/0.1.0\/history-0.1.0.json","date_time":"2025-07-02T13:22:32.989+0000","program":"hyrax","version":"1.17.1-219","parameters":[{"request_url":"https:\/\/opendap.earthdata.nasa.gov\/collections\/C1276812830-GES_DISC\/granules\/M2T1NXAER.5.12.4%3AMERRA2_200.tavg1_2d_aer_Nx.20000213.nc4.dap.nc4?dap4.ce=\/TOTANGSTR;\/time;\/lat;\/lon"},{"decoded_constraint":"dap4.ce=\/TOTANGSTR;\/time;\/lat;\/lon"}]}]'
_NCProperties: 'version=2,netcdf=4.9.0,hdf5=1.10.10'

>> **以上的信息特别全面，用于了解该数据集的基本信息，然后对此进行重采样以及转换投影信息！！！**

—————————————————————————————————————————————————

四、读取以及重采样代码———————————————————————————————————————

先裁剪后进行的重采样——速度快

%% ==== 灵活自适应的重采样函数（已检查并修复纬度翻转问题）====
% 本函数用于对输入的二维栅格数据进行空间重采样，支持指定目标空间分辨率和经纬度范围。
% 输入参数：
%   data            - 待重采样的二维数据矩阵（行代表纬度，列代表经度）
%   R_source        - 输入数据的空间参考对象（georasterref）
%   target_resolution - 目标空间分辨率（单位：度）
%   method          - 插值方法，支持 'nearest','bilinear','bicubic'，默认线性插值
%   Latlim_target   - 目标纬度范围 [南 纬度, 北 纬度]，可选，默认使用原始数据范围
%   Lonlim_target   - 目标经度范围 [西 经度, 东 经度]，可选，默认使用原始数据范围
% 输出参数：
%   resampled_data  - 重采样后的数据矩阵
%   R_target        - 重采样数据对应的空间参考对象，起点为北边，列从北到南

function [resampled_data, R_target] = resample_geotiff_flexible(data, R_source, target_resolution, method, Latlim_target, Lonlim_target)
    % 如果未指定目标范围，则使用原始范围
    if nargin < 5 || isempty(Latlim_target)
        Latlim_target = R_source.Latlim;
    end
    if nargin < 6 || isempty(Lonlim_target)
        Lonlim_target = R_source.Lonlim;
    end

    % 确保经纬度范围升序，避免坐标范围颠倒
    if Latlim_target(1) > Latlim_target(2)
        Latlim_target = fliplr(Latlim_target);
    end
    if Lonlim_target(1) > Lonlim_target(2)
        Lonlim_target = fliplr(Lonlim_target);
    end

    % 计算目标网格中心点坐标，确保每个像元中心对齐目标分辨率网格
    lon_target = (Lonlim_target(1) + target_resolution/2) : target_resolution : (Lonlim_target(2) - target_resolution/2);
    lat_target = (Latlim_target(2) - target_resolution/2) : -target_resolution : (Latlim_target(1) + target_resolution/2);

    % 原始数据尺寸与经纬度坐标生成
    [nrows_src, ncols_src] = size(data);
    lon_orig = linspace(R_source.Lonlim(1), R_source.Lonlim(2), ncols_src);
    % 注意：原始纬度从北到南，需要反转顺序以匹配行号
    lat_orig = linspace(R_source.Latlim(2), R_source.Latlim(1), nrows_src);

    % 生成原始和目标网格坐标矩阵，供插值使用
    [LON_ORIG, LAT_ORIG] = meshgrid(lon_orig, lat_orig);
    [LON_TARGET, LAT_TARGET] = meshgrid(lon_target, lat_target);

    % 处理数据中的 NaN 值：插值前将 NaN 置零，插值后再还原 NaN 区域
    mask = ~isnan(data);
    data_filled = data;
    data_filled(~mask) = 0;

    % 映射用户指定插值方法到 interp2 支持的类型
    switch lower(method)
        case 'nearest'
            interp_method = 'nearest';
        case 'bilinear'
            interp_method = 'linear';
        case 'bicubic'
            interp_method = 'cubic';
        otherwise
            interp_method = 'linear';
    end

    % 对数据进行二维插值重采样
    resampled_data = interp2(LON_ORIG, LAT_ORIG, data_filled, LON_TARGET, LAT_TARGET, interp_method);

    % 对数据有效区域进行最近邻插值，用于还原原始数据的 NaN 区域
    resampled_mask = interp2(LON_ORIG, LAT_ORIG, double(mask), LON_TARGET, LAT_TARGET, 'nearest');
    resampled_data(resampled_mask < 0.5) = NaN;

    % 创建目标空间参考对象，定义栅格尺寸及经纬度范围
    % 注意列从北到南，即从上到下逐列排列
    R_target = georasterref('RasterSize', size(resampled_data), ...
                            'Latlim', Latlim_target, ...
                            'Lonlim', Lonlim_target, ...
                            'ColumnsStartFrom', 'north');
end

%% ---------------------------------------------------------------------------------------

%% ==== 处理指定的单个 MERRA-2 文件，输出为 .tif 和 .mat ====
clear; clc;

% === 指定文件路径 ===
nc_file = 'G:\MERRA-2\data\ae\M2T1NXAER.5.12.4-MERRA2_200.tavg1_2d_aer_Nx.20000101.nc4.nc4';  % 你的单个 nc 文件
base_output_dir = 'G:\MERRA-2\output\test\world\resampled_single_0.5';   % 输出目录

target_resolution = 0.5;
% interp_method = 'nearest';   % 插值方法，支持 'nearest','bilinear','bicubic'
% interp_method = 'bilinear';   % 插值方法，支持 'nearest','bilinear','bicubic'    %%  效果较差
interp_method = 'bicubic';   % 插值方法，支持 'nearest','bilinear','bicubic'
%% 目标经纬度区域裁剪示例，注释掉的可替换，留空默认全范围
% Latlim_target = [18 45];        % 目标纬度范围（南北）
% Lonlim_target = [100 125];      % 目标经度范围（东西）

%% 覆盖中国大部分区域示例
% Latlim_target = [10 55];
% Lonlim_target = [70 140];

% 长江流域区域裁剪示例（默认启用）
% Latlim_target = [23 37];          % 裁剪区域纬度范围
% Lonlim_target = [89 122];         % 裁剪区域经度范围

% 长江流域区域裁剪示例（默认启用）--更改格网大小
% Latlim_target = [22.9 37];          % 裁剪区域纬度范围
% Lonlim_target = [89 122.1];         % 裁剪区域经度范围
% 
%% 任意裁剪示例
% Latlim_target = [-60 60];
% Lonlim_target = [-140 140];
%% 全球显示示例
 % 原始格网给出的纬度、经度与这个有细微差异， 读出的格网范围是[-90 90]  [-180, 179.375] 
 % 这里没有搞懂他是按照中心点赋值经纬网数据还是严格按照一个经纬度赋值一个数值，从上述读出的经纬度数据来看，他就不是按照中心点读取的，是赋值来的，因为假设是中心点赋值，那么最后给出的纬度max应该是89.25度纬度min应该是-89.25，经度max是179.375度，经度min是-179.375才是哦！

Latlim_target = [-90 90];    
Lonlim_target = [-180 180];   % 这里裁剪给和重采样强制经纬度按照这个进行限制了，下面的计算空间参考也是类似

% 从文件名或默认设定年份
[~, fname, ~] = fileparts(nc_file);
match = regexp(fname, '(\d{4})', 'match');
if ~isempty(match)
    year = match{1};
else
    year = 'unknown';
end
ncinfo = ncinfo(nc_file);
% === 创建输出目录 ===
resample_dir = fullfile(base_output_dir, sprintf('%.2fdeg_%s', target_resolution, interp_method));
out_day_dir = fullfile(resample_dir, 'day');
mat_output_dir = fullfile(resample_dir, 'matfile');

mkdir(out_day_dir);
mkdir(mat_output_dir);

% === 读取数据 ===%% 读取经纬度和时间数据
nc_info = ncinfo(nc_file);
lon = double(ncread(nc_file, 'lon'));
lat = double(ncread(nc_file, 'lat'));
time = ncread(nc_file, 'time');
%% 读取气溶胶数据 (变量名可能为 TOTEXTTAU)  TOTEXTTAU
tot_data = ncread(nc_file, 'TOTEXTTAU');

%% 替换无效值为 NaN
fill_value = -9.9999999e+14;     % 无效值填充
tot_data(tot_data == fill_value) = NaN;

%% 沿时间维度求均值（第三维）
tot_mean = nanmean(tot_data, 3)';  % 转置
tot_mean = flipud(tot_mean);   % 翻转数据

% === 计算空间参考 ===
dlon = mode(diff(lon));
dlat = mode(diff(lat));    %% 这里应该有一个是要绝对值的，忘记了；
% lonlim = [lon(1)-dlon/2, lon(end)+dlon/2];
% latlim = [min(lat)-dlat/2, max(lat)+dlat/2];
lonlim = [-180,180];   %% 这里没有参考读取的经度；
latlim  = [-90,90];    %% 这是读取的纬度；

R_source = georasterref('RasterSize', size(tot_mean),'Latlim', latlim, 'Lonlim', lonlim, 'ColumnsStartFrom', 'north');

% === 年度重采样 ===
fprintf('重采样年度数据...\n');
[aod_res, R_res] = resample_geotiff_flexible(tot_mean, R_source, target_resolution, interp_method, Latlim_target, Lonlim_target);

% === 保存年度数据 ===
geotiffwrite(fullfile(out_day_dir, sprintf('%s_aod_%.3fdeg_%s.tif', year, target_resolution, interp_method)), aod_res, R_res, 'CoordRefSysCode', 'EPSG:4326');
save(fullfile(mat_output_dir, sprintf('%s_ann_resampled.mat', year)), 'aod_res', 'R_res');

fprintf('完成处理文件: %s\n', nc_file);

关于MERRA-2经纬度格网问题

存疑：
不确定MERRA-2是按照中心点进行读取的经纬度格网

验证：
以MOD08_D3 1度 X 1度 全球数据集作为参照，MOD08_D3是读取的经纬度的中心点，读取的经纬度范围分别是[-89.5, 89.5] [-179.5, 179.5] 行和列分别为180和360；
MERRA-2 0.5度 X 0.625度 读取的经纬度范围分别是[-90 90] [-180, 179.375] ，行和列分别为361和576；

从上述读出的经纬度数据来看，看起来像是赋值来的，因为假设是中心点赋值，那么最后给出的纬度应该是[-89.75, 89.75]，经度应该是[-179.6875, 179.6875] ，但是呢，MERRA2的赋值好像还不对，经度少了一列，难绷；就相当于MERRA-2是按照经纬线来构成网格，不是中心点构成网格；

算了，再次假设，他就是按照中心点进行计算的，读出的就是中心点，那么，原始图层经纬度应该是[-90.25, 90.25] [-180.3125, 179.6875]，按照这个经纬网，导出的图是比[-90,90] [-180,180] 大一点，除了东经180，少了半个像元的格网，而且像元大小也没有压缩，是0.5 X 0.625 ；；；

再次显示验证：

为了再次验证：使用ArcMap 10.2直接将原始的nc通过多维工具——创建NetCDF栅格图层；也将MAC v3 1度的全球网格也创建出来，同时为了避免数据集相同，出现同样的问题，又添加了MOD08_D3数据集，直接拖进ArcMap中；全图展示如下：
MOD08和MACv3是完全重叠的，这是两者的经纬度范围：
在这里插入图片描述
这是图像展示：确保了MOD08和MAC v3是完全重叠的

下面以MERRA-2和MOD08进行对比：经纬度对比

图像显示：

结果：ArcMap按照经纬度中心点读取MERRA-2的数据，所以才会导致上述的结果；只有这一个结论，再见吧，毁灭吧！

尝试代码直接读取原始数据：得到的结果是读取的经纬度范围，就这样吧，代码是参考博主：
https://blog.youkuaiyun.com/weixin_45626690/article/details/122373638

%% MERRA-2 AOD 数据处理与 GeoTIFF 输出（按日均值和逐小时导出）
%% merra2_aod_to_geotiff.m   
%% 注意——输出的像元是按照原始像元计算的，与0.5 和 0.625 稍微有点差异，因为读取经纬度所致；；；
%% 要想严格控制格网精度，就重采样吧，避免同一个像元下，数据集总是差一列或者多一列的；；；如此烦人
% -------------------------------------------------------------------------
% 功能：读取 MERRA-2 AOD（气溶胶光学厚度）NetCDF 数据，提取每天的平均值
%       以及逐小时数据，并将其分别导出为 GeoTIFF 格式，用于可视化与后续处理。
% 输入：
%   - 输入路径 InPath 中包含的 .nc4 文件（每个文件包含 24 小时的 TOTEXTTAU 数据）
% 输出：
%   - 每个文件生成一个日均值 GeoTIFF，保存在 OutPath 路径下；
%   - 每个小时生成一个 GeoTIFF，共 24 个小时，保存在 hour_output_path 路径下；
% 注意：
%   - 代码会自动旋转矩阵（rot90），确保行对应纬度；
%   - 填充值（-9.9999999e+14）将被替换为 NaN；
%   - 如果输出目录不存在，将自动创建；
%   - 支持批量处理多个 .nc4 文件，并显示进度信息。
% 作者：小小程序牛马猿
% 日期：2025年7月
% -------------------------------------------------------------------------

clear; clc;

%% ======================= 路径设置 ======================= %%
InPath = 'G:\MERRA-2\data\ae\';  % NetCDF 输入路径
OutPath = 'G:\MERRA-2\output\ae_tiff2\mean\';  % 输出路径（日均值）
hour_output_path = 'G:\MERRA-2\output\ae_tiff2\hour\';  % 输出路径（逐小时）

% 创建输出目录（如不存在）
if ~exist(OutPath, 'dir'), mkdir(OutPath); end
if ~exist(hour_output_path, 'dir'), mkdir(hour_output_path); end

% 获取所有 .nc4 文件列表
filelist = dir(fullfile(InPath, '*.nc4'));
k = length(filelist);  % 文件数量

fprintf('共检测到 %d 个 MERRA-2 文件，开始批量处理...\n\n', k);

%% ======================= 主循环处理 ======================= %%
for i = 1:k
    fprintf('>>> 正在处理第 %d/%d 个文件：%s\n', i, k, filelist(i).name);

    % 构建完整路径
    InFile = fullfile(InPath, filelist(i).name);

    % 提取日期（例如：20200101），若失败则用索引编号
    match = regexp(InFile, '(\d{8})', 'match');
    if ~isempty(match)
        ymd = match{1};
    else
        ymd = sprintf('F%03d', i);
    end

    %% ========== 数据读取与处理 ==========
    % 读取 AE 数据（纬度 × 经度 × 时间）
    ae_raw = ncread(InFile, 'TOTANGSTR');
    
    % 旋转数据（MERRA-2 格式通常需旋转90度）
    ae = rot90(ae_raw, 1);  % 顺时针旋转90度

    % 替换填充值为 NaN
    fill_value = -9.9999999e+14;
    ae(ae == fill_value) = NaN;

    % 计算每日平均值（忽略 NaN）
    ae_mean = nanmean(ae, 3);

    % 读取经纬度
    lon = double(ncread(InFile, 'lon'));
    lat = double(ncread(InFile, 'lat'));

    % 构建空间参考对象（仅构建一次）
    GeoRef = georasterref( ...
        'RasterSize', size(ae_mean), ...
        'Lonlim', [min(lon), max(lon)], ...
        'Latlim', [min(lat), max(lat)] ...
    );

    %% ========== 保存日均值 GeoTIFF ==========
    mean_tif_path = fullfile(OutPath, [ymd, '_mean.tif']);
    geotiffwrite(mean_tif_path, flip(ae_mean), GeoRef);
    fprintf('    √ 已保存日均值 GeoTIFF：%s\n', mean_tif_path);

    %% ========== 保存逐小时 GeoTIFF ==========
    for j = 1:24
        ae_hour = ae(:, :, j);  % 提取第 j 小时数据
        hour_tif_path = fullfile(hour_output_path, [ymd, '_hour_', sprintf('%02d', j), '.tif']);
        geotiffwrite(hour_tif_path, flip(ae_hour), GeoRef);
        fprintf('        --- 已完成并保存第 %02d 小时\n', j);
    end

    fprintf('>>> 完成 %d/%d 个文件处理。\n\n', i, k);
end

fprintf('✅ 全部文件处理完成！\n');

先记录到这里了，去学习一下师妹的Python代码重采样，学习探讨一波。
下次记录再会！！！
适当的休息一下，不然脾气会暴躁，没有耐心，会一直碎碎念，烦人，的！！！

这里是先重采样后裁剪的代码，与上述代码相比，这个代码得到的数据精度会更高一点；
再次记录防止忘记：
开整：

%% 特别说明
%% 该代码读取是先重采样然后进行裁剪的，得到的结果没有严格限制经纬度范围，
%% 是按照中心点进行裁剪的，所以得到的数据是包含裁剪范围的，但是会多出半个像元；
%% 这是因为没有在函数中进行/2pixel的设置，重新设置一下就可以；
%% 这样裁剪也不能说没有问题哈，行列是一样的，就是对不齐而已，（哈哈哈，小问题改大问题）；；；留作备份

%% ==== 修正重采样裁剪函数，确保纬度升序 ====
%% ==== 修正重采样裁剪函数，严格按经纬度裁剪 ====
%% ========================================================================
% 函数名称：resample_then_crop
% -------------------------------------------------------------------------
% 功能：
%   本函数用于对输入的二维地理栅格数据进行：
%     1. 空间重采样：将原始数据插值到指定的空间分辨率；
%     2. 严格裁剪：根据用户指定的经纬度范围裁剪输出结果，确保输出数据严格位于目标区域之内。
%
% 输入参数：
%   - data (二维矩阵)：原始地理数据矩阵，维度为 [纬度 × 经度]；
%   - R_source (georasterref)：输入数据的空间参考对象，描述原始数据的地理边界和分布；
%   - target_resolution (标量)：目标空间分辨率，单位为度（°）；
%   - method (字符串)：插值方法，支持如下几种：
%       * 'nearest'  - 最近邻插值（适用于分类数据）；
%       * 'bilinear' - 双线性插值（默认值，适用于连续变量）；
%       * 'bicubic'  - 三次卷积插值（精度更高但计算更慢）；
%   - Latlim_target (1×2 向量, 可选)：输出数据的纬度范围 [南界, 北界]；
%   - Lonlim_target (1×2 向量, 可选)：输出数据的经度范围 [西界, 东界]；
%
% 输出参数：
%   - resampled_data_crop (二维矩阵)：重采样并裁剪后的输出数据矩阵；
%   - R_target_crop (georasterref)：裁剪后数据的空间参考对象。
%
% 主要特性：
%   - 支持缺失值处理，插值前将 NaN 设为 0，插值后再掩膜恢复为 NaN；
%   - 裁剪逻辑基于插值后数据网格中心点是否完全落入用户定义的范围；
%   - 裁剪结果空间参考对象与目标区域严格匹配。
%
% 使用示例：
%   [outData, outRef] = resample_then_crop(data, R, 0.1, 'bicubic', [25 35], [105 120]);
%
% 更新记录：
%   - 2025-07-05：初版，支持严格裁剪与多种插值方式。
% -------------------------------------------------------------------------
%% ========================================================================

function [resampled_data_crop, R_target_crop] = resample_then_crop(data, R_source, target_resolution, method, Latlim_target, Lonlim_target)
    % === 1. 原始经纬度网格生成 ===
    [nrows_src, ncols_src] = size(data);
    lon_orig = linspace(R_source.Lonlim(1), R_source.Lonlim(2), ncols_src);
    lat_orig = linspace(R_source.Latlim(2), R_source.Latlim(1), nrows_src);  % 注意纬度方向：北->南
    [LON_ORIG, LAT_ORIG] = meshgrid(lon_orig, lat_orig);

    % === 2. 构建目标重采样网格中心点 ===
    lon_full = (R_source.Lonlim(1)+target_resolution/2):target_resolution:(R_source.Lonlim(2)-target_resolution/2);
    lat_full = (R_source.Latlim(2)-target_resolution/2):-target_resolution:(R_source.Latlim(1)+target_resolution/2);
    [LON_TARGET, LAT_TARGET] = meshgrid(lon_full, lat_full);

    % === 3. 选择插值方法 ===
    switch lower(method)
        case 'nearest'
            interp_method = 'nearest';
        case 'bilinear'
            interp_method = 'linear';
        case 'bicubic'
            interp_method = 'cubic';
        otherwise
            interp_method = 'linear';
    end

    % === 4. 处理缺失值（NaN 设为 0） ===
    mask = ~isnan(data);
    data_filled = data;
    data_filled(~mask) = 0;

    % === 5. 重采样数据（插值） ===
    resampled_full = interp2(LON_ORIG, LAT_ORIG, data_filled, LON_TARGET, LAT_TARGET, interp_method);

    % === 6. 用掩膜恢复 NaN 区域 ===
    resampled_mask = interp2(LON_ORIG, LAT_ORIG, double(mask), LON_TARGET, LAT_TARGET, 'nearest');
    resampled_full(resampled_mask < 0.5) = NaN;

    % === 7. 创建重采样后的空间参考对象 ===
    lat_sorted = sort([min(lat_full), max(lat_full)]);  % 确保升序
    R_full = georasterref('RasterSize', size(resampled_full), ...
                          'Latlim', lat_sorted, ...
                          'Lonlim', [min(lon_full), max(lon_full)], ...
                          'ColumnsStartFrom', 'north');

    % === 8. 设置裁剪范围（如未提供则默认用重采样全范围） ===
    if nargin < 5 || isempty(Latlim_target)
        Latlim_target = R_full.Latlim;
    end
    if nargin < 6 || isempty(Lonlim_target)
        Lonlim_target = R_full.Lonlim;
    end

    % === 9. 计算裁剪区域在重采样网格中的行列索引 ===
    lat_vec = linspace(R_full.Latlim(2) - target_resolution/2, R_full.Latlim(1) + target_resolution/2, R_full.RasterSize(1));  % 北到南
    lon_vec = linspace(R_full.Lonlim(1) + target_resolution/2, R_full.Lonlim(2) - target_resolution/2, R_full.RasterSize(2));  % 西到东

    row_idx = find(lat_vec >= Latlim_target(1) & lat_vec <= Latlim_target(2));
    col_idx = find(lon_vec >= Lonlim_target(1) & lon_vec <= Lonlim_target(2));

    % === 10. 执行裁剪操作 ===
    resampled_data_crop = resampled_full(row_idx, col_idx);

    % === 11. 创建裁剪后的空间参考对象 ===
    lat_crop = [Latlim_target(1), Latlim_target(2)];
    lon_crop = [Lonlim_target(1), Lonlim_target(2)];
    R_target_crop = georasterref('RasterSize', size(resampled_data_crop), ...
                                 'Latlim', lat_crop, ...
                                 'Lonlim', lon_crop, ...
                                 'ColumnsStartFrom', 'north');
end

%% ========================================================================
%% ========================================================================
%% 按照先重采样然后裁剪的顺序——缺点就是运行的时间慢，优点是低像元采样到高像元会比先裁剪后重采样效果好些；；；
% 
%  MERRA-2 气溶胶变量 NetCDF 批处理与严格裁剪输出脚本
%  功能：将指定目录下的 MERRA-2 NetCDF 文件中的气溶胶变量（如 TOTANGSTR）
%       读取后进行：
%       1. 沿时间维度平均；
%       2. 空间重采样（支持最近邻、线性、三次插值）；
%       3. 严格按照用户定义的经纬度范围裁剪（不扩展或偏移）；
%       4. 输出为 GeoTIFF（带空间参考）和 .mat 文件（含数据与空间参考）。
%
%  输入：路径中所有符合条件的 NetCDF 文件
%  输出：按变量名与年份分类的 .tif 与 .mat 文件
%  说明：默认裁剪区域为长江流域，可自行修改 Latlim_target 与 Lonlim_target
%
%  更新日志：
%  ------------------------------------------------------------------------
%  2025-07-05 初版（支持重采样 + 严格裁剪，经纬度范围严格控制）
%             修复纬度方向与空间参考定义一致性问题
%  ------------------------------------------------------------------------
%% ========================================================================
%% ==== 主程序部分：使用严格裁剪函数 ====
clear; clc;
input_pattern = 'G:\MERRA-2\data\ae\*.nc4';
base_output_dir = 'G:\MERRA-2\output\ae_test_crop_lat2';
varname = 'TOTANGSTR';
%% 重采样的像元
target_resolution = 0.1;
%% 插值方式
interp_method = 'bilinear';   %% bilinear   bicubic
%% 裁剪方式
Latlim_target = [22.9 37];
Lonlim_target = [89 122.1];
%% nc文件所在位置导向，以及总数据量列表
nc_files = dir(input_pattern);
total_files = length(nc_files);
fprintf('共检测到 %d 个 NetCDF 文件，准备处理变量 <%s> ...\n\n', total_files, varname);
%% 开始循环读取
for idx = 1:total_files
    try
        nc_file = fullfile(nc_files(idx).folder, nc_files(idx).name);
        [~, fname, ~] = fileparts(nc_file);
        %% 匹配时间 年月日
        match = regexp(fname, '(\d{8})', 'match');
        if ~isempty(match)
            ymd = match{1};
            year_str = ymd(1:4);
        else
            ymd = sprintf('F%03d', idx);
            year_str = 'unknown';
        end

        tif_year_dir = fullfile(base_output_dir, 'tif', varname, year_str);
        mat_year_dir = fullfile(base_output_dir, 'mat', varname, year_str);
        if ~exist(tif_year_dir, 'dir'), mkdir(tif_year_dir); end
        if ~exist(mat_year_dir, 'dir'), mkdir(mat_year_dir); end

        lon = double(ncread(nc_file, 'lon'));
        lat = double(ncread(nc_file, 'lat'));
        data_raw = ncread(nc_file, varname);
        fill_value = -9.9999999e+14;
        data_raw(data_raw == fill_value) = NaN;

        data_mean = nanmean(data_raw, 3)';
        data_mean = flipud(data_mean);

        lat_min = min(lat);
        lat_max = max(lat);
        lon_min = min(lon);
        lon_max = max(lon);
        %% 原始数据空间参考
        R_source = georasterref( ...
            'RasterSize', size(data_mean), ...
            'Latlim', [lat_min, lat_max], ...
            'Lonlim', [lon_min, lon_max], ...
            'ColumnsStartFrom', 'north');

        [data_res, R_res] = resample_then_crop( ...
            data_mean, R_source, ...
            target_resolution, interp_method, ...
            Latlim_target, Lonlim_target);

        tif_name = sprintf('%s_%s_%.1fdeg_%s.tif', varname, ymd, target_resolution, interp_method);
        mat_name = sprintf('%s_%s_%.1fdeg_%s.mat', varname, ymd, target_resolution, interp_method);

        geotiffwrite(fullfile(tif_year_dir, tif_name), data_res, R_res, ...
            'CoordRefSysCode', 'EPSG:4326');
        save(fullfile(mat_year_dir, mat_name), 'data_res', 'R_res');

        fprintf('[%3d/%3d] %s\n', idx, total_files, tif_name);

    catch ME
        fprintf('[%3d/%3d] 错误处理文件: %s\n  错误信息: %s\n', ...
            idx, total_files, nc_file, ME.message);
    end
end

fprintf('\n 全部完成，共成功处理 %d 个文件。\n', total_files);