matlab 距平,[转载]基于Matlab软件进行EOF分解、回归趋势分析,并
clear
close
clc
ncload('E:/data/mon/ndvi79-06.nc','time','lat','lon','ndvi');
%第一部分:数据处理,剔除缺失值,求距平,并修正为等权重
%下面均针对温娜的sst进行标注
%选取区域:纬度:-9.5-59.5 ;经度:60.5-149.5
ilat=find(lat>=-10 &
lat<60);
nlat=length(ilat);
slat=lat(ilat);%选取区域纬度
ilon=find(lon>=60 &
lon<=150);
nlon=length(ilon);
slon=lon(ilon);%选取区域经度
itim=1:264; %选取时间为:1982.01-2003.12
ntim0=length(itim);
%表示月份数264 nyr=ntim0/12; %表示年数22
nph=nlat*nlon; %表示总共格点数
nph
st1=ndvi(itim,ilat,ilon);&4*70*90,将变量ndvi赋值给st1,以便后面统一使用公共变量st
st2=st1;&4*11*48,这里本该进行重新插值为低网格,即两两合并
%对每个格点进行计算月度年序列均值,然后计算月度年序列距平值
ndstp2=length(st2(1,:,1)); p
ndstp3=length(st2(1,1,:)); �
nlatlon=ndstp2*ndstp3;%表示低网格化后格点数目
% to get rid of the defaut value
%剔除缺失值,重新生成相应序列
st3=reshape(st2,ntim0,nlatlon);&4*6300(=70*90)
n1=find(abs(st3(1,:))<999); %寻找非缺失值位置 npstp=length(n1); %剔除缺失值后网个数:3136个
st4=st3(:,n1); %将剔除缺失值后的值重新赋值给264*3136
% to get rid of the seasonal cycle
st5=reshape(st4,12,nyr,npstp); *22*3136
st6_1=st5(3:5,:,:); %选取每个格点每年3到5月份值,3*22*3136 %从下面开始由于数据汇总,序列长度缩短为22个
st6_2=sum(st6_1,1); %对数据计算3-5月份ndvi和,得1*22*3136
st6=squeeze(st6_2); "*3136
% st7=st6-mean(st6,1); clear
ntim0;
ntim=nyr; %由于数据汇总 ,序列长度缩短,这里重新定义ntim;
for i=1:ntim
st7(i,:)=st6(i,:)-mean(st6,1); "*3136
end
st8=reshape(st7,ntim,npstp); "*3136
st9=st8; "*3136,这里本该进行剔除趋势
% to get equal weight in every point
%首先将缺失值部分补全,然后赋予等权重
st10(1:ntim,1:nlatlon)=-999; "*6300
st10(1:ntim,n1)=st8; "*6300
st11=reshape(st10,ntim,ndstp2,ndstp3);"*70*90
for i=1:ndstp2
st12(:,i,:)=st11(:,i,:)*sqrt(cos(slat(i)*3.14/180));"*70*90
end
%下面再次剔除缺失值
st13=reshape(st12,ntim,ndstp2*ndstp3); "*6300
st14(1:ntim,1:npstp)=st13(:,n1); "*3136
stp=st14
; "*3136
%第二部分 : 下面进行 EOF
分析,数据为全年观测值距平,但由于格点数大于时间维,故采用时空转换求解
cst=stp*stp'./ntim;"*22
[eof0,lamda0,explain0]=pcacov(cst);%eof0: 22*22
eof1=stp'*eof0; %
时空转换后空间模态:3136*22
for i=1:ntim
sd=norm(eof1(:,i));
eof2(:,i)=eof1(:,i)/sd;%将空间模态归一化
end
pca0=stp*eof1; %主成分得分22*22
for i=1:ntim
sd=norm(pca0(:,i));
pca(:,i)=pca0(:,i)/sd;%将时间序列归一化
end
lamda=(npstp/ntim)*lamda0;
explain=lamda/sum(lamda)*100;
%下面将EOF空间模态结果coef补全,缺失位置赋值-999,
eof3=eof2'; % 为了后面存储方便,进行转置:22*3136
eof4(1:ntim,1:nlatlon)=-999; "*6300
eof4(1:ntim,n1)=eof3; %将上述转置后空间模态,利用有值位置你n1进行填充真值,得到:22*6300空间模态,一行表示一个模态
eof=reshape(eof4,ntim,ndstp2,ndstp3);R8*11*48 一行一个转换,这利于后面存储
%第三部分: 求解空间场的回归趋势
for i=1:npstp
% x=[ones(ntim,1),1:ntim]';
% y=stp(:,i);
% [b,bint,r,rint,stats] = regress(y,x);
% stpcoef=b(2);
% stpfstat=stats(3);
x=[1:ntim]';
y=stp(:,i);
nstats =
regstats(y,x);%给出回归常用统计量
stpcoef0(i)=nstats.beta(2);
stpstat0(i)=nstats.fstat.pval;
end
%首先我们不考虑显著性水平时回归系数
stpcoef1=stpcoef0;
stpcoef2(1:nlatlon)=-999; stpcoef2(n1)=stpcoef0; stpcoef=reshape(stpcoef2,ndstp2,ndstp3);
%其次我们考虑显著性水平,
stpstat1=stpcoef0;
n2=find((stpstat1>0.10));
stpstat1(n2)=-999;
stpstat2(1:nlatlon)=-999; stpstat2(n1)=stpstat1; stpstat=reshape(stpstat2,ndstp2,ndstp3);
%第四部分:将结果输出到grads使用的格式资料
%将上述回归分析结果输出为grads使用的二进制格式:reg1022.dat,并给出描述ctl文件:reg1022.ctl;
fid=fopen(['D:/lee0921/Leedata/reg1022.dat'],'w')
B1=stpcoef;
B1=B1';
fwrite(fid,B1,'float32');
B2=stpstat;
B2=B2';
fwrite(fid,B2,'float32'); fclose(fid);
%=====================================================
a1=['dset
D:/lee0921/Leedata/reg1022.dat'];
a2=['undef -999'];
a3=['title Regression trend of MAM NDVI
(NDVI/decade:1982-2003'];
a4=['xdef 90
linear 60.5 1'];
a5=['ydef 70
linear -9.5 1'];
a6=['zdef 1
linear 0 0'];
a7=['tdef 1 linear 00Z01JAN1982 12mo'];
a8=['vars 2']; a91=['stpcoef 0 -999 the regression trend of MAM
NDVI']; a92=['stpstat 0 -999 the regression trend of MAM
NDVI']; a10=['endvars'];
fid=fopen('D:/lee0921/Leepro/r-ndvi-pre/reg1022.ctl','w'); fprintf (fid,
'%sn',a1);
fprintf (fid,
'%sn',a2);
fprintf (fid,
'%sn',a3);
fprintf (fid,
'%sn',a4);
fprintf (fid,
'%sn',a5);
fprintf (fid,
'%sn',a6);
fprintf (fid,
'%sn',a7);
fprintf (fid,
'%sn',a8);
fprintf (fid,
'%sn',a91);
fprintf (fid,
'%sn',a92);
fprintf (fid,
'%sn',a10);
fclose(fid);
clear a1 a2 a3 a4 a5 a6 a7 a8 a91 a92 a10 B1
B2
%将上述EOF分解的时间系数项输出到grads使用的二进制格式:pca1022.dat,并给出描述ctl文件:pca1022.ctl;
fid=fopen(['D:/lee0921/Leedata/pca1022.dat'],'w');
for i=1:ntim
%这里ntim代表变量个数22个,即以变量作为外循环,以时间作为内循环
spca=pca(i,:);
B=spca;
fwrite(fid,B,'float32'); end
fclose(fid);
%=====================================================
a1=['dset
D:/lee0921/Leedata/pca1022.dat'];
a2=['undef -999'];
a3=['title PCA from 198201-200312'];
a4=['xdef 1
linear 1 1'];
a5=['ydef 1
linear 1 1'];
a6=['zdef 1
linear 0 0'];
a7=['tdef 22 linear 00Z01JAN1982 12mo'];
a8=['vars 22']; % for i=1:22
% a9=['pca' num2str(i+9-1), ' 0
-999 mon mean eof of ndvi'];
% fprintf (fid,
'%sn',a9);
% end
a10=['endvars'];
fid=fopen('D:/lee0921/Leepro/r-ndvi-pre/pca1022.ctl','w'); fprintf (fid,
'%sn',a1);
fprintf (fid,
'%sn',a2);
fprintf (fid,
'%sn',a3);
fprintf (fid,
'%sn',a4);
fprintf (fid,
'%sn',a5);
fprintf (fid,
'%sn',a6);
fprintf (fid,
'%sn',a7);
fprintf (fid,
'%sn',a8);
for i=1:22
a9=['pca' num2str(i), ' 0 -999 mon mean eof of
ndvi'];
fprintf (fid, '%sn',a9);
end
fprintf (fid, '%sn',a10);
fclose(fid);
clear a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 B
%将上述EOF分解的结果输出为grads使用的二进制格式:eof1022.dat,并给出描述ctl文件:eof1022.ctl;
fid=fopen(['D:/lee0921/Leedata/eof1022.dat'],'w');
for i=1:22
seof=eof(i,:,:);
B=squeeze(seof);
B=B';
fwrite(fid,B,'float32'); end
fclose(fid);
%=====================================================
a1=['dset
D:/lee0921/Leedata/eof1022.dat'];
a2=['undef -999'];
a3=['title EOF from
198201-200312 ','(',num2str(explain(1)),'%)'];
a4=['xdef 90
linear 60.5 1'];
a5=['ydef 70
linear -9.5 1'];
a6=['zdef 1
linear 0 0'];
a7=['tdef 22 linear 00Z01JAN1982 12mo'];
a8=['vars 1']; a9=['eof 0
-999 mon mean eof of ndvi']; a10=['endvars'];
fid=fopen('D:/lee0921/Leepro/r-ndvi-pre/eof1022.ctl','w'); fprintf (fid,
'%sn',a1);
fprintf (fid,
'%sn',a2);
fprintf (fid,
'%sn',a3);
fprintf (fid,
'%sn',a4);
fprintf (fid,
'%sn',a5);
fprintf (fid,
'%sn',a6);
fprintf (fid,
'%sn',a7);
fprintf (fid,
'%sn',a8);
fprintf (fid,
'%sn',a9);
fprintf (fid,
'%sn',a10);
fclose(fid);
clear a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 B
%
将原始ndvi资料转换为grads使用的二进制格式:ndvi79_06_mon_mean.dat,并给出描述ctl文件:ndvi1019.ctl;
fid=fopen(['D:/lee0921/Leedata/ndvi79_06_mon_mean.dat'],'w');
for t=1:300
sndvi=ndvi(t,:,:);
B=squeeze(sndvi);
B=B';
fwrite(fid,B,'float32'); end
fclose(fid);
%=====================================================
a1=['dset
D:/lee0921/Leedata/ndvi79_06_mon_mean.dat'];
a2=['undef -999'];
a3=['title NDVI from 197901-200612'];
a4=['xdef
360 linear -179.5 1'];
a5=['ydef
180 linear -89.5 1'];
a6=['zdef 1
linear 0 0'];
a7=['tdef 300 linear 00Z01JAN1979 1mo'];
a8=['vars 1']; a9=['ndvi 0
-999 mon mean of ndvi']; a10=['endvars'];
fid=fopen('D:/lee0921/Leepro/r-ndvi-pre/ndvi1019.ctl','w'); fprintf (fid,
'%sn',a1);
fprintf (fid,
'%sn',a2);
fprintf (fid,
'%sn',a3);
fprintf (fid,
'%sn',a4);
fprintf (fid,
'%sn',a5);
fprintf (fid,
'%sn',a6);
fprintf (fid,
'%sn',a7);
fprintf (fid,
'%sn',a8);
fprintf (fid,
'%sn',a9);
fprintf (fid,
'%sn',a10);
fclose(fid);
clear a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 B
% % clear eof0 eof1 eof2 eof3 eof4 st1 st2 st3 st4
st5 st6 st6_1 st6_2 st7 st8 st9 st10 st11 st12 st13 st14;
% clear stpcoef0 stpcoef1 stpstat0
stpstat1;
% clear i ans explain0 lamda0 fid nstats pca0 sd
seof spca x y;
matlab 距平,[转载]基于Matlab软件进行EOF分解、回归趋势分析,并相关推荐
- matlab设计译码器,基于MATLAB的循环码编译码器设计与仿真.doc
扳昂旨螺冈唉陨裤外狸尿恨铸伸隧刽搅必勒诚天腑皖漂豌鲁靳碑缆键兽峙棘陶宽槐撒层僧袁廖颤渐魄货鼎躬薛扬衍逮西兰迫依煤鲁虐渠惫平合啥昭并屿己笆坍痞庐披吏去凄嘛兄察突徊溅今箩直藩潦咙锨谓崇若制匹扮复淌颐糖嗅你 ...
- 「电子万年历matlab仿真」——基于Matlab的电子万年历仿真实现
「电子万年历matlab仿真」--基于Matlab的电子万年历仿真实现 作为一种具有时间显示.日期查询.闹钟提醒等功能的电子产品,电子万年历已经成为了人们日常生活中不可或缺的一部分.而在现代科技的发展 ...
- 怎么用matlab分析孔隙度,基于MATLAB软件的声波测井孔隙度求取
·99·2016年 第 24 期 基于MATLAB软件的声波测井孔隙度求取 向旻 (新疆工程学院,新疆 乌鲁木齐 830091) 摘要 :目前,声波测井是一种重要的孔隙度测井方法,被各大油田广泛的使用 ...
- matlab论文 关于高数,高数和matlab论文,关于基于MATLAB软件的轨道交通高职院校高等数学课程教学相关参考文献资料-免费论文范文...
导读:本文是关于高数和matlab类函授毕业论文范文跟轨道交通有关函授毕业论文范文. 吴玲玲 [摘 要] 在轨道交通高职院校高等数学课程教学中嵌入 MATLAB 软件,通过 MATLAB 的可视化 ...
- matlab仿真直流电机,[转载]基于Matlab/Simulink的无刷直流电机控制仿真研究
摘要: 基于Matlab/Simulink,本文设计了一个无刷直流电机的控制方案,详细阐述了无刷直流电机的运行原理,并用Matlab/Simulink对其进行了仿真.实验证明,用Matlab/Simu ...
- FLAC3D可视化后处理matlab,一种基于Matlab的由Midas导入Flac3D的模型识别方法与流程...
本发明涉及岩土工程的仿真模拟研究领域,具体涉及一种基于Matlab的由Midas导入Flac3D的模型识别方法. 背景技术: 随着我国近几年经济的快速发展和基础设施等的大力投资,涉及复杂地质环境下的岩 ...
- matlab演示系统,基于Matlab的通信原理演示系统的设计与应用
基于 Matlab的通信原理演示系统的设计与应用 李 强 , 明 艳 , 吴坤君 (重庆邮电大学 通信学院 , 重庆 400065) 摘 要 : 利用 Matlab图形用户界面的开发环境和强大的通信仿 ...
- 基于matlab的霍夫变换,基于matlab的霍夫变换
MATLAB 三维绘图功能 Plot3函数(三维曲线图) Mesh函数(网格图) Surf函数(曲面...步骤: 1.利用hough()函数执行霍夫变换,得到霍夫矩阵; 2.利用houghpeaks( ...
- 基于matlab的回波,基于MATLAB回波信号产生与消除.doc
基于MATLAB回波信号产生与消除 摘 要 MATLAB可以进行矩阵运算.绘制函数和数据.实现算法.创建用户界面.连接其他编程语言的程序等,主要应用于工程计算.控制设计.信号处理与通讯.图像处理.信号 ...
- matlab距离保护程序,基于MATLAB的距离保护仿真.doc
基于MATLAB的距离保护仿真 摘要:本文阐述了如何利用Matlab中的Simulink及SPS工具箱建立线路的距离保护仿真模型,并用S函数编制相间距离保护和接地距离保护算法程序,构建相应的保护模块, ...
最新文章
- 如何获取主机的IP址址
- 用python画猫咪怎么画-Turtle库画小猫咪
- 这八大互联网金融商业模式,你都知道吗?
- DCL 管理用户 mysql
- 常用 AT 命令手册
- Android实现登录
- Mysql 复制原理以及配置 简要分析
- linux g++开启C++11/14支持
- keil p0 0c语言不了,Keil C51对C语言的关键词扩展之十三: sfr
- 2012禁用ip隧道 win_IMCP协议的魅力——IMCP隧道
- 3种方法实现Android按钮的点击事件,建议收藏!
- 强烈推荐几个BAT大佬技术公众号~值得学习!
- Atitit 2019技术趋势与没落技术 目录 1.1. abcdAtitit 技术领域趋势 abcd研究总结AI(人工智能)BlockChain(区块链)、Cloud(云)、和Data(大数据)
- Vista暴力破解器只是一个玩笑 谁说破解谁撒谎
- 金蝶K3对接数据库相应语句大全
- tableau 日周月筛选器_【数据可视化】Tableau教程(六)日历热力图
- 更实用 批量解析 Sanger 测序.ab1文件 出图出文本
- tableau度量值计算_度量值与度量名称
- 12 个动画设计方法,帮助你快速实现炫酷的网页动画效果
- 把时间当作朋友 读书笔记