%by yinlichang3064@163.com
在进行长时间序列的栅格数据分析时,如NDVI,fvc,LAI,NPP,GPP,需要知道每个格点的长期趋势。
如果再arcgis中进行一元回归计算,需要将整个公式展开,并且容易出错,而matlab提供了强大的回归
计算功能,通过逐像元的迭代,可以得到每个格点的趋势及显著性栅格图像,话不多说,见以下代码
[a,R]=geotiffread('D:\日降水插值\年pre2000.tif');%先导入某个图像的投影信息,为后续图像输出做准确
info=geotiffinfo('D:\日降水插值\年pre2000.tif');
[m,n]=size(a);
years=16; %表示有多少年份需要做回归
data=zeros(m*n,years);
k=1;
for year=2000:2015 %起始年份
file=['D:\日降水插值\年prec',int2str(year),'.tif'];%注意自己的名字形式,这里使用的名字是年prec2000.tif,根据这个可修改
bz=importdata(file);
bz=reshape(bz,m*n,1);
data(:,k)=bz;
k=k+1;
end
xielv=zeros(m,n);p=zeros(m,n);
for i=1:length(data)
bz=data(i,:);
if max(bz)>0 %注意这是进行判断有效值范围,如果有效范围是-1到1,则改成max(bz)>-1即可
bz=bz';
X=[ones(size(bz)) bz];
X(:,2)=[1:years]';
[b,bint,r,rint,stats] = regress(bz,X);
pz=stats(3);
p(i)=pz;
xielv(i)=b(2);
end
end
name1='D:\simulation\黄河流域降水0.1度\year\趋势\年prec一元线性回归00-15趋势值.tif';
name2='D:\simulation\黄河流域降水0.1度\year\趋势\年prec一元线性回归00-15_P值.tif';
geotiffwrite(name1,xielv,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(name2,p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
%一般来说,只有通过显著性检验的趋势值才是可靠的
xielv(p>0.05)=NaN;
name1='D:\simulation\黄河流域降水0.1度\year\趋势\通过显著性0.05检验的年prec一元线性回归00-15趋势值.tif';
geotiffwrite(name1,xielv,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
在用上述代码的时候请引用以下文献:
Yin et al. 2020. Irrigation water consumption of irrigated cropland and its dominant factor in China from 1982 to 2015. Advance in water resources. 143. doi: 10.1016/j.advwatres.2020.103661
Yin, L.; Wang, X.; Feng, X.; Fu, B.; Chen, Y. A Comparison of SSEBop-Model-Based Evapotranspiration with Eight Evapotranspiration Products in the Yellow River Basin, China. Remote Sens. 2020, 12, 2528.
更多需求,请查看个人介绍