|
|
【答案】应助回帖
★ ★ ★ ★ ★
pingxie: 金币+5, 博学EPI+1, ★★★很有帮助 2015-07-21 10:06:31
不知道你要什么分析，这就是个简单的处理程序，读取数据，排序，计算之类的程序，跟excel表的功能差不多，特点就是批量处理一堆txt文件里的数据。具体专业的计算不懂，只能给你简单注释下通用文件操作。
function n=spi(R) %函数头
P=zeros(84,1);
file_data=zeros(84,1);
filedata=zeros(84,1);
g=zeros(84,1);
A=zeros(84,1);
C=zeros(30,84); %以上初始化
c0=2.515517;
c1=0.802853;
c2=0.010328;
d1=1.432788;
d2=0.189269;
d3=0.001308; % 以上专业参数
D=importdata('D:\Program Files\气象要素实时资料处理系统\text\ZX旬降水.txt');
data=D.data;
text=D.textdata; % 数据
[c,r]=sort(text(3:86,1)); % 排序
for m=1:84
x(m,=data(r(m),; % 排序后的数据写入一个矩阵
end
pathname=['D:\SPI\1971－2000逐站各旬降水量\'];
files=dir([pathname '*.txt' ]); % 获取文件夹中每个文件名
[file_num,s]=size(files); % 共有几个文件
for m=1:file_num
B=load([pathname files(m).name]);
n=R;
C(:,m)=B(:,n);
end % 以上依次读取每个文件，存入数组C
C=C.*0.1;
for m=1:84
c=C(:,m);
index=find(c~=32766);
index0=find(c==0);
[num,p]=size(index);
[num0,q]=size(index0);
if num0~=0
c(index0)=0.001;
end
P(m,1)=sum(c(index))/num;
file_data(m,1)=log(P(m,1));
d=prod(c(index));
filedata(m,1)=log(d)/num;
end
A=file_data-filedata;
X=x(:,2);
for m=1:84
if X(m,1)==0
X(m,1)=0.001;
end
a(m,1)=(1+sqrt(1+4*A(m,1)/3))/(4*A(m,1));
b(m,1)=P(m,1)/a(m,1);
G(m,1)=gammainc(X(m,1)/b(m,1),a(m,1));
end
for M=1:84
if G(M,1)>0.5
t=sqrt(log(1/((1-G(M,1))^2)));
SPI(M,1)=(t-(c1*t+c2*t^2-c0))/(1+d1*t+d2*t^2+d3*t^3);
else t=sqrt(log(1/(G(M,1)^2)));
SPI(M,1)=((c1*t+c2*t^2-c0)-t)/(1+d1*t+d2*t^2+d3*t^3);
end
end %以上就不懂了，应该是专业的计算
fid=fopen('D:\SPI\data.txt','w'); % 存入文件
fprintf(fid,'%.1f\r\n',SPI); % 写入数据
fclose(fid); |
|
回复此楼