matlab SPEI干旱指数计算 nc tif各种 数据，多个时间尺度 2000到2023年 1/3/6/12 尺度

最近在折腾气象干旱指数，发现SPEI（标准化降水蒸散指数）这玩意儿挺有意思。它比SPI多了蒸散发因子，更适合分析气候变化下的干旱特征。今天咱们就用Matlab实操一把多时间尺度的SPEI计算，处理nc和tif格式的栅格数据。

先准备好2000-2023年的月降水数据和潜在蒸散数据（PET）。建议用ncread处理nc文件：

prec = ncread('input.nc','precipitation'); % 维度通常是[lon,lat,time]
pet = ncread('pet.nc','pet');
time = ncread('input.nc','time'); % 注意时间单位转换

时间维度处理有个坑——很多nc文件的时间轴是"days since 1900-1-1"这种格式，得转成实际日期：

date_vec = datetime(double(time),'ConvertFrom','excel','Format','yyyy-MM');

不同时间尺度（1/3/6/12月）的SPEI需要先做累积计算。举个3个月尺度的滑动平均例子：

window_size = 3;
prec_acc = movmean(prec, window_size, 3, 'omitnan', 'Endpoints','discard');
pet_acc = movmean(pet, window_size, 3, 'omitnan', 'Endpoints','discard');

重点来了——SPEI计算核心。推荐使用现成的MATLAB SPEI函数包：

spei_3 = SPEI_12(prec_acc - pet_acc, 3); % 第三个参数是时间尺度
% 注意这里输入的应该是水平衡数据（降水-PET）

处理多尺度时，可以用循环结构批量生成：

scales = [1,3,6,12];
for s = scales
    window = ones(1,1,s)/s;
    balance = convn(prec - pet, window, 'valid');
    spei = SPEI_12(balance, s);
    % 输出为GeoTIFF时记得保留地理信息
    geotiffwrite(sprintf('SPEI%d_%s.tif',s,datestr(time)), spei, R);
end

当遇到tif数据时，地理参考信息要用geotiffread获取：

[prec, R] = geotiffread('prec.tif');
xll = R.LongitudeLimits(1);
cellsize = R.CellExtentInLatitude;

最后提醒几个避坑点： 1. 时间轴对齐要检查闰年数据 2. 空间数据注意Nodata值处理 3. 大区域计算时优先用单精度存储 4. 推荐用parfor并行加速格点计算

计算结果用ArcGIS或QGIS可视化时，突然发现某区域SPEI-12出现持续负值——恭喜你，这可能就是一次重大干旱事件的量化证据。这种从代码到现实灾害的映射感，或许就是气象编程的魅力所在吧。

最近在搞气象数据的朋友肯定绕不开干旱指数，SPEI这玩意儿比SPI多了温度因子更贴近实际。今天咱们就手把手用MATLAB从nc文件里抠数据，批量处理多时间尺度的SPEI计算，最后生成带地理坐标的tif文件。

先看数据预处理部分。假设咱们手头有2000-2023年的月降水（precip）和气温（temp）nc文件，先得把它们读进MATLAB：

% 读取nc文件
precip = ncread('input_precip.nc', 'precipitation');
temp = ncread('input_temp.nc', 'temperature');
lat = ncread('input_precip.nc', 'latitude');
lon = ncread('input_precip.nc', 'longitude');

% 处理缺失值
precip(precip < 0) = NaN;
temp(temp < -50) = NaN;

这里要注意很多公开数据集会用-9999这类值表示缺失，得先替换成NaN。接着计算潜在蒸散发(PET)，用Thornthwaite方法就够用：

% 计算PET函数
function pet = thornthwaite(T, lat, month)
    % 日照小时数估算（简化版）
    doy = mod(month-1,12)+1;
    phi = deg2rad(lat);
    delta = 0.409*sin(2*pi*doy/365 - 1.39);
    ws = acos(-tan(phi).*tan(delta));
    N = 24*ws/pi;

% 热量指数计算 I = sum(T/5.0.^1.514, 'omitnan'); a = (6.75e-7)*I.^3 - (7.71e-5)*I.^2 + (1.79e-2)*I + 0.492;

pet = 16*(N/12).*(T>0).*(10*T./I).^a;
end

温度数据需要按月循环处理，注意这里传入的month参数需要对应数据的时间维度。计算完PET后，水分盈亏量D=降水-PET，这是SPEI计算的基础：

D = precip - PET; % 逐月水分盈亏

重点来了——多时间尺度累积。咱们需要实现类似滑动平均的效果，但要注意NaN值的处理：

scales = [1 3 6 12]; % 四种时间尺度
for s = scales
    % 累积计算
    D_cum = movmean(D, [s-1 0], 1, 'omitnan');

% 伽马分布拟合
    parfor i = 1:size(D,2) % 并行加速纬度循环
        for j = 1:size(D,3)
            [gamma_params, p] = gamfit(squeeze(D_cum(:,i,j)));
            % ...后续标准化处理
        end
    end
end

这里用movmean做滑动累积，参数[s-1 0]表示只看过去s个月的数据。注意第三维是时间维度，第四维是空间维度。伽马分布拟合部分用parfor并行处理可以节省大量时间，特别是处理高分辨率全球数据时。

生成结果文件时要注意保持地理坐标信息。推荐用geotiffwrite：

% 输出tif示例
R = georasterref('LatitudeLimits',[min(lat) max(lat)],...
                 'LongitudeLimits',[min(lon) max(lon)]);
geotiffwrite('SPEI_12_mon.tif', SPEI_12, R);

最后可视化可以用自带的地理坐标系作图：

figure('Position', [100 100 1200 600])
ax = axesm('MapProjection','eqdcylin');
geoimg = geoshow(SPEI_12, R, 'DisplayType','texture');
geolimits([min(lat) max(lat)], [min(lon) max(lon)])
colorbar('southoutside')
caxis([-2.5 2.5]) % 干旱等级标准范围
colormap(jet(10))

几个避坑提示：1）输入数据的时间连续性要检查，避免存在月份缺失；2）伽马分布拟合可能遇到全零数据，记得加try-catch；3）输出tif时注意矩阵方向，部分数据集是纬度倒序存储的。搞过类似项目的朋友应该都懂这种坐标系对齐的痛，经常出现数据漂移或者镜像翻转的灵异事件。

完整流程跑下来，单时间尺度在16核服务器上大概需要半小时（1度分辨率全球数据）。建议把中间结果存成.mat文件，方便后续分析调用。下次可以试试把结果接进机器学习模型做干旱预测，这才是真正有意思的部分。