ERA5数据下载，直接下载，笨办法；
数据集：ERA5 hourly data on single levels from 1940 to present
下载的变量名称为：Total column water vapour（tcwv）/ Total column ozone（tco3）
总柱水汽、臭氧——读取与重采样代码
代码以水汽柱为例，需要读取臭氧的话，直接替换变量名称就可以。
其中日均值的计算里面有详细的注释也有函数
1、小时直接读取
%% ERA5 TCWV数据处理：输出原始分辨率与重采样GeoTIFF（逐小时）
% 本脚本读取ERA5的总柱水汽（tcwv）数据，逐小时输出：
% - 原始分辨率GeoTIFF（0.25度）
% - 重采样至指定分辨率（如1度）并裁剪输出GeoTIFF
% - 自动识别时间戳，无需预设开始时间
% - 每小时生成两个文件：原始精度和重采样精度
% - 终端实时打印处理进度
clc; clear; close all;
%% =================== 输入文件配置 ===================
nc_file = 'G:/a0002exprocess/eERA5/tcwv/tcwv/data_stream-oper_stepType-instant.nc'; % NetCDF 文件路径
varname = 'tcwv'; % 变量名称（总柱水汽）
%% =================== 数据读取 ===================
lon = ncread(nc_file, 'longitude'); % 经度向量（递增）
lat = ncread(nc_file, 'latitude'); % 纬度向量（递减）
time_raw = ncread(nc_file, 'valid_time'); % 时间戳（单位：秒自1970-01-01）
base_time = datetime(1970,1,1,0,0,0); % 时间基准
time_all = base_time + seconds(time_raw); % 转换为datetime数组
tcwv = double(ncread(nc_file, varname)); % 读取变量数据，格式：(lon, lat, time)
tcwv = permute(tcwv, [2, 1, 3]); % 转换为 (lat, lon, time)
%% =================== 原始空间参考 ===================
R_source = georasterref('RasterSize', size(tcwv(:,:,1)), ...
'Latlim', [min(lat), max(lat)], ...
'Lonlim', [min(lon), max(lon)], ...
'ColumnsStartFrom', 'north');
%% =================== 重采样参数设置 ===================
target_res = 1; % 重采样目标分辨率（度）
method = 'bilinear'; % 插值方法：nearest | bilinear | bicubic
% Latlim_target = [20, 50]; % 目标裁剪纬度范围
% Lonlim_target = [70, 140]; % 目标裁剪经度范围
Latlim_target = [18, 45]; % 目标裁剪纬度范围
Lonlim_target = [100, 125]; % 目标裁剪经度范围
%% =================== 输出路径配置 ===================
out_path_raw = 'G:/a0002exprocess/eERA5/tcwv/output/hour/tiff_raw_hourly/';
out_path_resample = 'G:/a0002exprocess/eERA5/tcwv/output/hour/tiff_resampled_hourly/';
if ~exist(out_path_raw, 'dir'); mkdir(out_path_raw); end
if ~exist(out_path_resample, 'dir'); mkdir(out_path_resample); end
%% =================== 按小时输出GeoTIFF ===================
ntime = size(tcwv, 3);
for t = 1:ntime
time_str = datestr(time_all(t), 'yyyymmdd_HH');
tcwv_hour = tcwv(:, :, t);
% === 输出原始分辨率GeoTIFF ===
filename_raw = sprintf('%s%s_tcwv_raw.tif', out_path_raw, time_str);
geotiffwrite(filename_raw, single(tcwv_hour), R_source);
% === 重采样并裁剪 ===
[data_out, R_out] = resample_and_clip(tcwv_hour, R_source, target_res, method, Latlim_target, Lonlim_target);
filename_resample = sprintf('%s%s_tcwv.tif', out_path_resample, time_str);
geotiffwrite(filename_resample, single(data_out), R_out);
% === 打印进度 ===
fprintf('完成 %s：输出原始与重采样 TIF 文件（%d/%d）\n', time_str, t, ntime);
end
disp('全部逐小时GeoTIFF文件输出完毕。');
效果展示：图中重采样是1度
2、日均值计算
%% ERA5 TCWV数据处理：输出原始分辨率与重采样GeoTIFF
% 本脚本功能说明：
% 1. 读取ERA5卫星的总可降水量（TCWV）数据（NetCDF格式）
% 2. 计算每日平均值（将每24个时次数据取平均）
% 3. 输出两种GeoTIFF文件：
% - 原始分辨率（0.25度）数据
% - 重采样至指定分辨率并裁剪到目标区域的数据
% 4. 自动打印处理进度，并根据数据时间戳命名输出文件
% 清除工作区变量、命令行窗口内容，关闭所有图形窗口（确保环境干净）
clc; clear; close all;
%% 输入文件路径与变量名
% 定义输入NetCDF文件的路径（根据实际数据位置修改）
nc_file = 'G:/a0002exprocess/eERA5/tcwv/tcwv/data_stream-oper_stepType-instant.nc';
% 定义要读取的变量名（TCWV：Total Column Water Vapor，总可降水量）
varname = 'tcwv';
%% 读取变量与坐标
% 读取经度数据（单位：度，通常为0-360或-180-180）
lon = ncread(nc_file, 'longitude');
% 读取纬度数据（单位：度，注意ERA5数据通常为降序，即从90°N到90°S）
lat = ncread(nc_file, 'latitude');
% 读取时间数据（单位：秒，通常以1970-01-01 00:00:00为基准）
time_raw = ncread(nc_file, 'valid_time');
% 定义基准时间（Unix时间戳起点：1970年1月1日0时0分0秒）
base_time = datetime(1970,1,1,0,0,0);
% 将原始时间（秒）转换为datetime格式（便于日期处理）
time_all = base_time + seconds(time_raw);
% 读取TCWV数据：原始维度通常为（longitude, latitude, time）
tcwv = double(ncread(nc_file, varname));
% 调整数据维度顺序：从（lon, lat, time）转为（lat, lon, time）
% 原因：MATLAB中图像显示通常为（行, 列, 页），对应（纬度, 经度, 时间）
[Lon, Lat] = meshgrid(lon, lat); % 生成经纬度网格（验证维度用，后续未直接使用）
tcwv = permute(tcwv, [2, 1, 3]); % 维度置换：第2维（lat）移到第1位，第1维（lon）移到第2位
%% 构建原始空间参考对象（0.25°）
% georasterref对象用于存储栅格数据的地理信息，便于后续地理空间操作（如写入GeoTIFF）
R_source = georasterref('RasterSize', size(tcwv(:,:,1)), ... % 栅格尺寸：（纬度数, 经度数）
'Latlim', [min(lat), max(lat)], ... % 纬度范围：[最小纬度, 最大纬度]
'Lonlim', [min(lon), max(lon)], ... % 经度范围：[最小经度, 最大经度]
'ColumnsStartFrom', 'north'); % 列从北侧（高纬度）开始计数
%% 参数设置
target_res = 1; % 重采样目标分辨率（单位：度，此处设为1度）
method = 'nearest'; % 重采样插值方法：最近邻插值
Latlim_target = [20, 50]; % 目标区域纬度范围：20°N - 50°N
Lonlim_target = [70, 140]; % 目标区域经度范围：70°E - 140°E
%% 输出路径
% 原始分辨率TIFF文件输出路径
out_path_raw = 'G:/a0002exprocess/eERA5/tcwv/output/tiff_raw/';
% 重采样后TIFF文件输出路径
out_path_resample = 'G:/a0002exprocess/eERA5/tcwv/output/tiff_resampled/';
% 若输出目录不存在，则创建目录（避免写入文件时出错）
if ~exist(out_path_raw, 'dir'); mkdir(out_path_raw); end
if ~exist(out_path_resample, 'dir'); mkdir(out_path_resample); end
%% 每日平均并输出GeoTIFF（每24个小时平均）
% 获取总时次数（第三维长度）
ntime = size(tcwv, 3);
% 计算总天数：假设每24个时次为一天（ERA5通常为每小时一个时次）
daily_num = floor(ntime / 24);
% 循环处理每一天的数据
for day = 1:daily_num
% 计算当天数据在时次序列中的索引范围
idx1 = (day - 1) * 24 + 1; % 当天第一个时次索引
idx2 = day * 24; % 当天最后一个时次索引（第24个时次）
% 计算当天24个时次的平均值（沿第三维取平均）
tcwv_daily = mean(tcwv(:, :, idx1:idx2), 3);
% 获取当天日期字符串（以当天第一个时次的时间为准，格式：yyyymmdd）
date_str = datestr(time_all(idx1), 'yyyymmdd');
% === 输出原始分辨率GeoTIFF ===
% 构建原始分辨率文件名称（包含路径和日期）
filename_raw = sprintf('%s%s_tcwv_raw.tif', out_path_raw, date_str);
% 写入GeoTIFF文件：数据类型转为single（节省空间），附加原始地理参考
geotiffwrite(filename_raw, single(tcwv_daily), R_source);
% === 重采样并裁剪后输出GeoTIFF ===
% 调用重采样和裁剪函数，获取处理后的数据和新的地理参考
[data_out, R_out] = resample_and_clip(tcwv_daily, R_source, target_res, method, Latlim_target, Lonlim_target);
% 构建重采样后文件名称
filename_resample = sprintf('%s%s_tcwv.tif', out_path_resample, date_str);
% 写入重采样后的GeoTIFF文件
geotiffwrite(filename_resample, single(data_out), R_out);
% 打印处理进度：显示当前日期、完成的文件数/总文件数
fprintf('完成 %s：输出原始与重采样 TIF 文件（%d/%d）\n', date_str, day, daily_num);
end
% 所有文件处理完成后，在命令行显示提示信息
disp('全部GeoTIFF文件输出完毕。');
%% ==================================================================================
%% ==================================================================================
%% ========== 函数定义：重采样+裁剪 ==========
% 功能：对输入的地理栅格数据进行重采样（改变分辨率）和空间裁剪（提取指定区域）
% 输入参数：
% data - 输入的二维数据矩阵 (lat, lon)，包含待处理的栅格数据
% R_source - 原始数据的地理参考对象（georasterref类型），包含原始数据的空间信息
% target_res - 目标重采样分辨率（单位：度）
% method - 插值方法，可选：'nearest'（最近邻）、'bilinear'（双线性）、'bicubic'（双三次）
% Latlim_target - 目标区域的纬度范围 [最小纬度, 最大纬度]
% Lonlim_target - 目标区域的经度范围 [最小经度, 最大经度]
% 输出参数：
% data_out - 重采样并裁剪后的二维数据矩阵
% R_out - 处理后数据的地理参考对象（georasterref类型）
function [data_out, R_out] = resample_and_clip(data, R_source, target_res, method, Latlim_target, Lonlim_target)
% 获取输入数据的行列数（纬度方向行数，经度方向列数）
[nrow, ncol] = size(data);
% 生成原始数据的经度网格：从原始经度范围的起点到终点，等间隔生成ncol个点
lon_src = linspace(R_source.Lonlim(1), R_source.Lonlim(2), ncol);
% 生成原始数据的纬度网格：注意纬度是从高到低（降序），因为多数地理数据纬度从上到下递减
lat_src = linspace(R_source.Latlim(2), R_source.Latlim(1), nrow);
% 构建原始数据的经纬度网格矩阵（meshgrid将向量转换为网格矩阵）
[LON_SRC, LAT_SRC] = meshgrid(lon_src, lat_src);
% 生成重采样后的经度序列：
% 从原始经度范围起点+分辨率一半到终点-分辨率一半，间隔为目标分辨率
% 这样确保每个网格中心对准目标分辨率的网格点
lon_all = (R_source.Lonlim(1) + target_res/2) : target_res : (R_source.Lonlim(2) - target_res/2);
% 生成重采样后的纬度序列：
% 从原始纬度范围高值-分辨率一半到低值+分辨率一半，间隔为负的目标分辨率（确保纬度降序）
lat_all = (R_source.Latlim(2) - target_res/2) : -target_res : (R_source.Latlim(1) + target_res/2);
% 构建重采样后的经纬度网格矩阵
[LON_RES, LAT_RES] = meshgrid(lon_all, lat_all);
% 根据输入的插值方法字符串，映射到interp2函数支持的插值方法
switch lower(method) % lower()将输入转为小写，确保大小写不影响判断
case 'nearest'; interp_method = 'nearest'; % 最近邻插值（适用于分类数据，保持原值）
case 'bilinear'; interp_method = 'linear'; % 双线性插值（适用于连续数据，平滑过渡）
case 'bicubic'; interp_method = 'cubic'; % 双三次插值（更高平滑度，可能改变极值）
otherwise; interp_method = 'linear'; % 默认使用双线性插值
end
% 处理数据中的NaN值（缺失值）：
% 1. 创建掩码矩阵：非NaN值为1，NaN值为0
mask = ~isnan(data);
% 2. 将原始数据中的NaN替换为0（避免插值时受NaN影响）
data_filled = data;
data_filled(~mask) = 0;
% 3. 对填充后的数据进行插值，得到重采样结果
resampled_data = interp2(LON_SRC, LAT_SRC, data_filled, LON_RES, LAT_RES, interp_method);
% 4. 对掩码矩阵进行最近邻插值（确保掩码边界准确）
resampled_mask = interp2(LON_SRC, LAT_SRC, double(mask), LON_RES, LAT_RES, 'nearest');
% 5. 用重采样后的掩码恢复缺失值：掩码值<0.5的位置视为缺失，重新赋值为NaN
resampled_data(resampled_mask < 0.5) = NaN;
% 计算裁剪索引：
% 1. 筛选出在目标经度范围内的重采样经度索引
lon_idx = lon_all >= Lonlim_target(1) & lon_all <= Lonlim_target(2);
% 2. 筛选出在目标纬度范围内的重采样纬度索引
lat_idx = lat_all >= Latlim_target(1) & lat_all <= Latlim_target(2);
% 3. 根据索引裁剪重采样数据，得到目标区域数据
data_out = resampled_data(lat_idx, lon_idx);
% 提取裁剪后的经纬度序列
lon_clip = lon_all(lon_idx);
lat_clip = lat_all(lat_idx);
% 构建输出数据的地理参考对象：
% 包含裁剪后的数据尺寸、纬度范围、经度范围，以及列从北开始的标记
R_out = georasterref('RasterSize', size(data_out), ...
'Latlim', [min(lat_clip)-target_res/2, max(lat_clip)+target_res/2], ... % 纬度范围包含网格边缘
'Lonlim', [min(lon_clip)-target_res/2, max(lon_clip)+target_res/2], ... % 经度范围包含网格边缘
'ColumnsStartFrom', 'north'); % 列从北侧（高纬度）开始
end
图中重采样为0.1度
3、月均值计算
%% ERA5 TCWV数据处理：输出原始分辨率与重采样GeoTIFF（每月平均）
% 本脚本读取ERA5的总柱水汽（tcwv）数据，计算每月平均值后：
% - 输出原始分辨率GeoTIFF（0.25度）
% - 重采样至指定分辨率（如1度）并裁剪输出GeoTIFF
% - 自动识别时间戳，无需预设开始时间
% - 每月生成两个文件：原始精度和重采样精度
% - 终端实时打印处理进度
clc; clear; close all;
%% =================== 输入文件配置 ===================
nc_file = 'G:/a0002exprocess/eERA5/tcwv/tcwv/data_stream-oper_stepType-instant.nc'; % NetCDF 文件路径
varname = 'tcwv'; % 变量名称（总柱水汽）
%% =================== 数据读取 ===================
lon = ncread(nc_file, 'longitude'); % 经度向量（递增）
lat = ncread(nc_file, 'latitude'); % 纬度向量（递减）
time_raw = ncread(nc_file, 'valid_time'); % 时间戳（单位：秒自1970-01-01）
base_time = datetime(1970,1,1,0,0,0); % 时间基准
time_all = base_time + seconds(time_raw); % 转换为datetime数组
tcwv = double(ncread(nc_file, varname)); % 读取变量数据，格式：(lon, lat, time)
tcwv = permute(tcwv, [2, 1, 3]); % 转换为 (lat, lon, time)
%% =================== 原始空间参考 ===================
R_source = georasterref('RasterSize', size(tcwv(:,:,1)), ...
'Latlim', [min(lat), max(lat)], ...
'Lonlim', [min(lon), max(lon)], ...
'ColumnsStartFrom', 'north');
%% =================== 重采样参数设置 ===================
target_res = 1; % 重采样目标分辨率（度）
method = 'bilinear'; % 插值方法：nearest | bilinear | bicubic
% Latlim_target = [20, 50]; % 目标裁剪纬度范围
% Lonlim_target = [70, 140]; % 目标裁剪经度范围
Latlim_target = [18, 45]; % 目标裁剪纬度范围
Lonlim_target = [100, 125]; % 目标裁剪经度范围
%% =================== 输出路径配置 ===================
out_path_raw = 'G:/a0002exprocess/eERA5/tcwv/output/month/tiff_raw/';
out_path_resample = 'G:/a0002exprocess/eERA5/tcwv/output/month/tiff_resampled/';
if ~exist(out_path_raw, 'dir'); mkdir(out_path_raw); end
if ~exist(out_path_resample, 'dir'); mkdir(out_path_resample); end
%% =================== 月均值计算与输出 ===================
ntime = size(tcwv, 3);
[yy, mm, ~] = ymd(time_all); % 提取年月信息
months = unique(yy * 100 + mm); % 生成唯一年月索引（如200001）
for m = 1:length(months)
yyyymm = months(m);
y = floor(yyyymm / 100); % 年
mo = mod(yyyymm, 100); % 月
idx = find(yy == y & mm == mo); % 找出该月所有时间索引
if isempty(idx); continue; end
tcwv_month = mean(tcwv(:, :, idx), 3); % 计算该月平均
date_str = sprintf('%04d%02d', y, mo); % 构造文件名时间戳
% === 输出原始分辨率GeoTIFF ===
filename_raw = sprintf('%s%s_tcwv_raw.tif', out_path_raw, date_str);
geotiffwrite(filename_raw, single(tcwv_month), R_source);
% === 重采样并裁剪 ===
[data_out, R_out] = resample_and_clip(tcwv_month, R_source, target_res, method, Latlim_target, Lonlim_target);
filename_resample = sprintf('%s%s_tcwv.tif', out_path_resample, date_str);
geotiffwrite(filename_resample, single(data_out), R_out);
% === 打印进度 ===
fprintf('完成 %s：输出原始与重采样 TIF 文件（%d/%d）\n', date_str, m, length(months));
end
disp('全部每月平均GeoTIFF文件输出完毕。');
重采样函数代码：
%% ========== 函数定义：重采样+裁剪 ==========
% 功能：对输入的地理栅格数据进行重采样（改变分辨率）和空间裁剪（提取指定区域）
% 输入参数：
% data - 输入的二维数据矩阵 (lat, lon)，包含待处理的栅格数据
% R_source - 原始数据的地理参考对象（georasterref类型），包含原始数据的空间信息
% target_res - 目标重采样分辨率（单位：度）
% method - 插值方法，可选：'nearest'（最近邻）、'bilinear'（双线性）、'bicubic'（双三次）
% Latlim_target - 目标区域的纬度范围 [最小纬度, 最大纬度]
% Lonlim_target - 目标区域的经度范围 [最小经度, 最大经度]
% 输出参数：
% data_out - 重采样并裁剪后的二维数据矩阵
% R_out - 处理后数据的地理参考对象（georasterref类型）
function [data_out, R_out] = resample_and_clip(data, R_source, target_res, method, Latlim_target, Lonlim_target)
% 获取输入数据的行列数（纬度方向行数，经度方向列数）
[nrow, ncol] = size(data);
% 生成原始数据的经度网格：从原始经度范围的起点到终点，等间隔生成ncol个点
lon_src = linspace(R_source.Lonlim(1), R_source.Lonlim(2), ncol);
% 生成原始数据的纬度网格：注意纬度是从高到低（降序），因为多数地理数据纬度从上到下递减
lat_src = linspace(R_source.Latlim(2), R_source.Latlim(1), nrow);
% 构建原始数据的经纬度网格矩阵（meshgrid将向量转换为网格矩阵）
[LON_SRC, LAT_SRC] = meshgrid(lon_src, lat_src);
% 生成重采样后的经度序列：
% 从原始经度范围起点+分辨率一半到终点-分辨率一半，间隔为目标分辨率
% 这样确保每个网格中心对准目标分辨率的网格点
lon_all = (R_source.Lonlim(1) + target_res/2) : target_res : (R_source.Lonlim(2) - target_res/2);
% 生成重采样后的纬度序列：
% 从原始纬度范围高值-分辨率一半到低值+分辨率一半，间隔为负的目标分辨率（确保纬度降序）
lat_all = (R_source.Latlim(2) - target_res/2) : -target_res : (R_source.Latlim(1) + target_res/2);
% 构建重采样后的经纬度网格矩阵
[LON_RES, LAT_RES] = meshgrid(lon_all, lat_all);
% 根据输入的插值方法字符串，映射到interp2函数支持的插值方法
switch lower(method) % lower()将输入转为小写，确保大小写不影响判断
case 'nearest'; interp_method = 'nearest'; % 最近邻插值（适用于分类数据，保持原值）
case 'bilinear'; interp_method = 'linear'; % 双线性插值（适用于连续数据，平滑过渡）
case 'bicubic'; interp_method = 'cubic'; % 双三次插值（更高平滑度，可能改变极值）
otherwise; interp_method = 'linear'; % 默认使用双线性插值
end
% 处理数据中的NaN值（缺失值）：
% 1. 创建掩码矩阵：非NaN值为1，NaN值为0
mask = ~isnan(data);
% 2. 将原始数据中的NaN替换为0（避免插值时受NaN影响）
data_filled = data;
data_filled(~mask) = 0;
% 3. 对填充后的数据进行插值，得到重采样结果
resampled_data = interp2(LON_SRC, LAT_SRC, data_filled, LON_RES, LAT_RES, interp_method);
% 4. 对掩码矩阵进行最近邻插值（确保掩码边界准确）
resampled_mask = interp2(LON_SRC, LAT_SRC, double(mask), LON_RES, LAT_RES, 'nearest');
% 5. 用重采样后的掩码恢复缺失值：掩码值<0.5的位置视为缺失，重新赋值为NaN
resampled_data(resampled_mask < 0.5) = NaN;
% 计算裁剪索引：
% 1. 筛选出在目标经度范围内的重采样经度索引
lon_idx = lon_all >= Lonlim_target(1) & lon_all <= Lonlim_target(2);
% 2. 筛选出在目标纬度范围内的重采样纬度索引
lat_idx = lat_all >= Latlim_target(1) & lat_all <= Latlim_target(2);
% 3. 根据索引裁剪重采样数据，得到目标区域数据
data_out = resampled_data(lat_idx, lon_idx);
% 提取裁剪后的经纬度序列
lon_clip = lon_all(lon_idx);
lat_clip = lat_all(lat_idx);
% 构建输出数据的地理参考对象：
% 包含裁剪后的数据尺寸、纬度范围、经度范围，以及列从北开始的标记
R_out = georasterref('RasterSize', size(data_out), ...
'Latlim', [min(lat_clip)-target_res/2, max(lat_clip)+target_res/2], ... % 纬度范围包含网格边缘
'Lonlim', [min(lon_clip)-target_res/2, max(lon_clip)+target_res/2], ... % 经度范围包含网格边缘
'ColumnsStartFrom', 'north'); % 列从北侧（高纬度）开始
end