基于matlab進(jìn)行NDVI對降水量和溫度的長時間柵格數(shù)據(jù)進(jìn)行偏相關(guān)分析
????????最近在做工作的時候，想到NDVI數(shù)據(jù)受到降水量、溫度的雙重影響，所以在進(jìn)行分析時，有些up主的工作只是對降水量和NDVI或者是溫度和NDVI的相關(guān)分析。但是其實在做的時候應(yīng)該考慮使用控制變量法來進(jìn)行相關(guān)的研究，一般的研究只是假設(shè)NDVI值受到三個因素的影響，即與溫度、降水量和人類活動因素有關(guān)，NDVI值、溫度、降水量都可以獲得相應(yīng)的數(shù)據(jù)，人類活動因素則是使用以上的三個數(shù)據(jù)做殘差分析后得到的。具體的過程大家可以在小破站上找到，此處不再贅述。
????????于是呢，我就在網(wǎng)上也學(xué)習(xí)了一下，加了一點自己的理解寫了一個程序。我自己的水平也不是很高，公開出來呢只是為了方便后面的人在做到和我類似的工作的時候能夠減少一點在程序上面花費的時間，更好的去分析解決自己所遇到的問題，如果幫到你了就給我一個點贊吧，有問題的話大家也可以在評論區(qū)里留言，大家一起交流進(jìn)步，我也會不定時的回復(fù)大家的。廢話不多說了，程序就放在下面了：
????????
clc
clear
%導(dǎo)入投影信息
[a,R]=geotiffread('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\年度數(shù)據(jù)\年NDVI值\2000.tif');
info=geotiffinfo('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\年度數(shù)據(jù)\年NDVI值\2000.tif');
[m,n]=size(a);
%導(dǎo)入NDVI數(shù)據(jù)
NDVI=zeros(m*n,21);
for year=1:21
? ?filename=strcat('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\2000-2020年歸一化數(shù)據(jù)按月份、季度排序\二月\','NDVI',int2str(year),'.tif');%此處要修改,將年份修改為總共的年份。
? ?data=importdata(filename);
? ? data=reshape(data,m*n,1);
? ?NDVI(:,year-0)=data;
end
%導(dǎo)入降水?dāng)?shù)據(jù)
PRE=zeros(m*n,21);
for year=1:21
? ?filename=strcat('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\2000-2020年降水量裁剪數(shù)據(jù)按月份、季節(jié)排序\2月\','CNPre',int2str(year),'.tif');%此處要修改,將年份修改為總共的年份。
? ?data=importdata(filename);
? ? data=reshape(data,m*n,1);
? ?PRE(:,year-0)=data;
end
%導(dǎo)入溫度數(shù)據(jù)
TEMP=zeros(m*n,21);
for year=2000:2020
? ?filename=strcat('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\2000-2020年溫度裁剪數(shù)據(jù)按月份、季節(jié)排序\2月\','CNtmp_',int2str(year),'_M2_1x1.tif');%此處要修改,將年份修改為總共的年份。
? ?data=importdata(filename);
? ? data=reshape(data,m*n,1);
? ?TEMP(:,year-1999)=data;
end
%溫度的相關(guān)性和顯著性
NDVI_TEMP_xgx=zeros(m*n,1);
NDVI_TEMP_p=zeros(m*n,1);
for i=1:length(NDVI)
? ?ndvi=NDVI(i,:);
? ?ndvi=ndvi';
? ?temp=TEMP(i,:)';
? ?pre=PRE(i,:)';
? ?[r,p]=partialcorr(ndvi,pre,temp);
? ?NDVI_TEMP_xgx(i)=r;
? ?NDVI_TEMP_p(i)=p;
end
NDVI_TEMP_xgx=reshape(NDVI_TEMP_xgx,m,n);
NDVI_TEMP_p=reshape(NDVI_TEMP_p,m,n);
filename1=('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\溫度降水量偏相關(guān)分析\降水量_NDVI\降水量_NDVI_2月數(shù)據(jù)偏相關(guān)分析.tif');
filename2=('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\溫度降水量偏相關(guān)分析\降水量_NDVI\降水量_NDVI_2月數(shù)據(jù)偏相關(guān)分析顯著性.tif');
geotiffwrite(filename1,NDVI_TEMP_xgx,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(filename2,NDVI_TEMP_p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
%顯著性檢驗
NDVI_TEMP_xgx=reshape(NDVI_TEMP_xgx,m*n,1);
NDVI_TEMP_p=reshape(NDVI_TEMP_p,m*n,1);
trend=zeros(m*n,1);
for i=1:length(NDVI)
? ?if NDVI_TEMP_xgx(i)>0
? ? ? ?if NDVI_TEMP_p(i)<0.05
? ? ? ? ? ?trend(i)=1;
? ? ? ?else NDVI_TEMP_p(i)>0.05;
? ? ? ? ? ?trend(i)=2;
? ? ? ?end
? ?elseif NDVI_TEMP_xgx(i)<0
? ? ? ? ? ?if NDVI_TEMP_p(i)<0.05
? ? ? ? ? ? ? ?trend(i)=4;
? ? ? ? ? ?else NDVI_TEMP_p(i)>0.05;
? ? ? ? ? ? ? ?trend(i)=5;
? ? ? ? ? ?end
? ?else NDVI_TEMP_xgx(i)=0;
? ? ? ?trend(i)=3;
? ?end
end
trend=reshape(trend,m,n);
filename3=('G:\畢業(yè)論文材料\裁剪數(shù)據(jù)\溫度降水量偏相關(guān)分析\降水量_NDVI\降水量_NDVI_2月數(shù)據(jù)通過顯著性檢驗的偏相關(guān)分析.tif');
geotiffwrite(filename3,trend,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
這個程序會輸出三個圖，第一個圖是相關(guān)系數(shù)的圖像，第二張圖是顯著性的圖像，第三張圖是通過顯著性檢查的出結(jié)論的圖像，這里重點說說第三張圖。
????????有些人在做完顯著性檢驗之后，會發(fā)現(xiàn)自己的圖一片黑或者是偶爾有幾個點，所以我就把顯著性檢驗的圖改成了5個數(shù)，分別為1、2、3、4、5，分別代表著顯著正向關(guān)、不顯著正相關(guān)、非線性相關(guān)、顯著負(fù)相關(guān)和不顯著負(fù)相關(guān)，這些需要在整個圖出來以后在ArcGIS里進(jìn)行重分類，重分類之后就是你想要的結(jié)果了。重分類的步驟大家可以在網(wǎng)上找相關(guān)的文章，這里也不進(jìn)行贅述。
????????最后說幾個比較重要的問題：
這個程序不能對1年的數(shù)據(jù)進(jìn)行相關(guān)分析，不然不會有結(jié)果的。
所有進(jìn)行分析的柵格數(shù)據(jù)的行列數(shù)必須一致，這很重要?。?！
如果有問題的話私聊我，我一般每天中午吃飯時會刷B站，我看到就會回復(fù)。
