ERA5---MATLAB处理水汽数据与臭氧数据的读取与重采样-重复性工作

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:ntimetime_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);
enddisp('全部逐小时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值为0mask = ~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的位置视为缺失，重新赋值为NaNresampled_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; endtcwv_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));
enddisp('全部每月平均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值为0mask = ~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的位置视为缺失，重新赋值为NaNresampled_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