clear
%ncdisp?HadISST_sst.nc
?%time?=?1730??（UNLIMITED）
?%latitude??=?180
?%longitude?=?360
?%nv????????=?2
%units?=?‘days?since?1870-1-1?0:0:0‘

lon=ncread（‘HadISST_sst.nc‘‘longitude‘）;%提取經度數據
lat=ncread（‘HadISST_sst.nc‘‘latitude‘）;?%提取緯度數據
sst=ncread（‘HadISST_sst.nc‘‘sst‘）;??????%提取海表面sst數據%?讀取nc格式數據

sst1=sst（1:9021:801201:1440）;??????????%選取所需要區域的數據
sst2=sst（311:36021:801201:1440）;???????%130E-90W20N-80N
sst3=zeros（14060240）;
sst3（90:-1:11:601:240）=sst1;
sst3（140:-1:911:601:240）=sst2;
sst=sst3;??????????????????????????????????????????

sst_area1=zeros（2408400）;????????????????%zeros全零數組
for?i=1:240;
?squ=squeeze（sst（::i））;?????????????????%執行該指令后sst數據轉換為二維數組?????????????????????????
?sst_area1（i:）=reshape（squ18400）;??????%將數據轉變為二維
end?

%?剔除與其它站點相關系數小的站點的數據
sst_area1（sst_area1==-1000）=NaN;%?簡單的認為剔除陸地和冬季結冰點的數據
sst_nan=isnan（sst_area1）;
i=0;
for?j=1:8400
????if?sum（sst_nan（:j））==0;
????????i=i+1;
????????sst_region（:i）=sst_area1（:j）;
????end
end
%**************************************************************************
%?求距平~注意季節的變換
X=zeros（size（sst_region））;??????????%?學者給的程序
for?i=1:12
X（i:12:228+i:）=sst_region（i:12:end:）?-?repmat（?mean（sst_region（i:12:end:）1）??size（sst_region（i:12:end:）1）?1）;
end

%?由于變量數（8400）>觀測樣本數~協方差矩陣的階數較大，因此可用時空轉換的方法提高計算速度
???????????????????????????????????????????????????????
R=X*X‘;?????????%?協方差矩陣R=X*X‘是8400*8400的方陣~現定義矩陣R=X‘*X是240*240的矩陣
[vd]=eig（R）;???%?進行EOF分解~因為X‘*X與X*X‘的秩相同所以特征值相同~d為x的特征值組成的對角陣~v為X*X‘的特征向量
v=fliplr（v）;????%?矩陣作左右翻轉
d=rot90（d2）;???%?矩陣上下翻轉后再左右翻轉（查看生成的對角陣是由小到大排列的，用命令doc?rot90可查詢
diagonal=diag（d）;
spacef=X‘*v;???????????????????????????????????????
for?i=1:240;
????spacef（:i）=spacef（:i）/sqrt（diagonal（i））;?????%?空間本征函數
end
timef=X*spacef;????????????????????????????????????%?時間本征函數
sum_d=sum（diagonal）;
count=0;
for?i=1:240;
????count=count+diagonal（i）;???????????
????G1（i）=count/sum_d;????????????????????????????%?G1（i）是累積方差貢獻率
end
for?i=1:240;
????G2（i）=diagonal（i）/sum_d;?????????????????????%?G2（i）是方差貢獻率
end
?????????????????????????????????????????????????%?將刪去的陸地與冰點的填充值補回
sst_area2=zeros（2408400）;
sst_area2（::）=NaN;
i=1;
spacef2=spacef‘;
for?j=1:8400
???if?sum（sst_nan（:j））==0;
??????sst_area2（:j）=spacef2（:i）;
??????i=i+1;
???end
end
sst_area3=sst_area2‘;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%????????????????????????????????????????????????
%畫圖

x=1:240;
plot（xtimef（:6）‘b‘）;
ylim（[?-80?80?]）;
%?xlabel（‘1971-1990年240個月‘‘fontsize‘15‘fontname‘‘隸書‘）?
ylabel（‘INDEX‘‘fontsize‘12‘fontname‘‘黑體‘）
set（gca‘xtick‘（1:6:162））
set（gca‘xticklabel‘{‘1971‘‘‘‘1972‘‘‘‘1973‘‘‘‘1974‘‘‘‘1975‘‘‘‘1976‘‘‘‘1977‘‘‘‘1978‘‘‘‘1979‘‘‘‘1980‘‘‘‘1981‘‘‘‘1982‘‘‘‘1983‘‘‘‘1984‘‘‘‘1985‘‘‘‘1986‘‘‘‘1987‘‘‘‘1988‘‘‘‘1989‘‘‘‘1990‘‘‘‘1991‘}）
title（‘北太平洋第6模態1971至1990年SST時間序列‘?‘color‘?‘k‘‘fontsize‘15‘fontna

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????337557??2014-04-24?23:32??010022011028藍方寶.docx
?????文件???????14363??2014-04-24?20:19??010022011028lfb\圖片\001.png
?????文件???????14334??2014-04-24?20:38??010022011028lfb\圖片\002.png
?????文件???????14316??2014-04-24?20:38??010022011028lfb\圖片\003.png
?????文件???????14305??2014-04-24?20:39??010022011028lfb\圖片\004.png
?????文件???????14305??2014-04-24?20:40??010022011028lfb\圖片\005.png
?????文件???????14310??2014-04-24?20:42??010022011028lfb\圖片\006.png
?????文件???????61709??2014-04-24?19:34??010022011028lfb\圖片\011.png
?????文件???????61474??2014-04-24?19:35??010022011028lfb\圖片\022.png
?????文件???????65711??2014-04-24?19:40??010022011028lfb\圖片\033.png
?????文件???????59179??2014-04-24?19:42??010022011028lfb\圖片\044.png
?????文件???????69114??2014-04-24?19:44??010022011028lfb\圖片\055.png
?????文件???????68440??2014-04-24?19:45??010022011028lfb\圖片\066.png
?????文件????????4372??2014-04-24?20:02??010022011028lfb\程序\newEOF.m
?????文件?????????918??2014-04-24?19:47??010022011028lfb\程序\new_draw.m
?????文件?????????445??2014-04-24?20:41??010022011028lfb\程序\timef_draw.m