%% 输入数据
load 'gp_chla_ISS' %%野外数据
% load 'yue_chlaISS' %%实验室混合数据
% load 'Totaluse' %%2010-2016年数据
%% 计算相关矩阵和协方差矩阵
cor=corrcoef(m_gp');
covv=cov(m_gp');
%% 主成分分析方法
for i=1:5
[corrr(i).loading,corrr(i).c]=c_pca(cor,3,i);
[covvv(i).loading,covvv(i).c]=c_pca(covv,3,i);
end
%% 余弦距离
for i=1:5
corrr(i).D=pdist2(corrr(i).loading,m_gp','cos');
covvv(i).D=pdist2(covvv(i).loading,m_gp','cos');
for j=1:3
corrr(i).minD(j)=find(corrr(i).D(j,:)==min(corrr(i).D(j,:)));
covvv(i).minD(j)=find(covvv(i).D(j,:)==min(covvv(i).D(j,:)));
end
end
%% 水体光谱曲线图
set(0,'defaultfigurecolor','w') ;
y1=plot(ubd,m_gp,'color','k','linewidth',1);
xlabel '波长（nm）';ylabel '遥感反射率（sr^-^1）';
set(gca,'FontSize',12);
set(gca,'LineWidth',2);
set(get(gca,'YLabel'),'FontSize',18);
set(get(gca,'XLabel'),'FontSize',18);
set(gca,'Xlim',[400 850]);
set(gca,'Ylim',[0 0.02]);
%% 协方差载荷系数曲线图
set(0,'defaultfigurecolor','w') ;
for i=1:3
figure;
[AX,h1,h2] = plotyy([ubd',ubd'],[covvv(1).loading(i,:)',covvv(3).loading(i,:)'],[ubd',ubd'],[covvv(2).loading(i,:)',covvv(4).loading(i,:)'],'plot');
xlabel '波长（nm）';ylabel '载荷系数';
set(get(AX(1),'YLabel'),'FontSize',18);
set(get(AX(2),'YLabel'),'FontSize',18);
set(AX(1),'FontSize',15);
set(AX(2),'FontSize',15);
set(AX(1),'Position',[.2 .2 .7 .65]);
set(h1(1),'LineStyle','--','Color','k','linewidth',2);
set(h1(2),'LineStyle','-','Color','k','linewidth',2);
set(h2(1),'LineStyle',':','Color','k','linewidth',2);
set(h2(2),'LineStyle','-.','Color','k','linewidth',2);
set(gca,'LineWidth',2);
set(gca,'Xlim',[400 850]);
legend([h1(1),h2(1),h1(2),h2(2)],'方法一','方法二','方法三','方法四')
end
%% 相关矩阵载荷系数曲线图
for i=1:3
figure;
[AX,h1,h2] = plotyy([ubd',ubd'],[corrr(1).loading(i,:)',corrr(3).loading(i,:)'],[ubd',ubd'],[corrr(2).loading(i,:)',corrr(4).loading(i,:)'],'plot');
xlabel '波长（nm）';ylabel '载荷系数';
set(get(AX(1),'YLabel'),'FontSize',18);
set(get(AX(2),'YLabel'),'FontSize',18);
set(AX(1),'FontSize',15);
set(AX(2),'FontSize',15);
set(AX(1),'Position',[.2 .2 .7 .65]);
set(h1(1),'LineStyle','--','Color','k','linewidth',2);
set(h1(2),'LineStyle','-','Color','k','linewidth',2);
set(h2(1),'LineStyle',':','Color','k','linewidth',2);
set(h2(2),'LineStyle','-.','Color','k','linewidth',2);
set(gca,'LineWidth',2);
set(gca,'Xlim',[400 850]);
legend([h1(1),h2(1),h1(2),h2(2)],'方法一','方法二','方法三','方法四')
end
%% NIPALS载荷系数曲线图
figure;
cycsm{1}='--';
cycsm{2}='-';
cycsm{3}=':';
for i=1:3
plot(ubd,covvv(5).loading(i,:),'color','k','linewidth',2,'linestyle',cycsm{i});hold on;
end
xlabel '波长（nm）';ylabel '载荷系数';
set(gca,'FontSize',12);
set(gca,'LineWidth',2);
set(get(gca,'YLabel'),'FontSize',18);
set(get(gca,'XLabel'),'FontSize',18);
set(gca,'Xlim',[400 850]);
legend('第一主成分','第二主成分','第三主成分')
figure;
for i=1:3
plot(ubd,corrr(5).loading(i,:),'color','k','linewidth',2,'linestyle',cycsm{i});hold on;
end
xlabel '波长（nm）';ylabel '载荷系数';
set(gca,'FontSize',12);
set(gca,'LineWidth',2);
set(get(gca,'YLabel'),'FontSize',18);
set(get(gca,'XLabel'),'FontSize',18);
set(gca,'Xlim',[400 850]);
legend('第一主成分','第二主成分','第三主成分')
