水蚀_RUSLE_K._

clear all in='D:\2019\20190402\SHUISHI\K\'; T_SAND=dlmread([in,'sand.txt'],' ',6,0); T_SAND=T_SAND(:,1:2338); T_SILT=dlmread([in,'silt.txt'],' ',6,0); T_SILT=T_SILT(:,1:2338); T_CLAY=dlmread([in,'clay.txt'],' ',6,0); T_CLAY=T_CLAY(:,1:2338); T_OC=dlmread([in,'oc.txt'],' ',6,0); T_OC=T_OC(:,1:2338); mask=dlmread('D:\DHSVM\studyregion\ASCII\mask01.txt'); SILT_NEW=nan(1857,2338); SAND_NEW=nan(1857,2338); K=nan(1857,2338); %判断输入数据的大小，根据输入数据的大小设置函数运行中涉及的相关变量的空矩阵 OM=1.724.*T_OC; %计算土壤有机质这一因子 x=[0.002,0.05,2]; xx=[0.02,2]; %将数据中的美国制转换为国际制，x是美国制的分割点，xx是国际制中的分割点 for i=1:1857 for j=1:2338 if mask(i,j)==1 y=[T_CLAY(i,j),T_CLAY(i,j)+T_SILT(i,j),T_CLAY(i,j)+T_SILT(i,j)+T_SAND(i,j)]; %计算每一个分割点对应了累计量 yy=spline(x,y,xx); %按照国际制的分割点，进行三次样条插值，得到国际制分割点对应的值 SILT_NEW(i,j)=yy(1)-T_CLAY(i,j); SAND_NEW(i,j)=yy(2)-SILT_NEW(i,j)-T_CLAY(i,j); %插值后分别计算分割点对应的数值 K(i,j)=(0.2+0.3*exp(-0.0256*SAND_NEW(i,j)*(1-SILT_NEW(i,j)/100)))*... (SILT_NEW(i,j)/(SILT_NEW(i,j)+T_CLAY(i,j)))^0.3*... (1-0.25*OM(i,j)/(OM(i,j)+exp(3.72-2.95*OM(i,j))))*... (1-(0.7*(1-SAND_NEW(i,j)/100))/((1-SAND_NEW(i,j)/100)+exp(-5.51+22.9*(1-SAND_NEW(i,j)/100)))); %将每一组计算好的数据带入公式中计算土壤可蚀性因子 else K(i,j)=-9999; end end end dlmwrite('D:\2019\20190402\SHUISHI\K\Kvalue.txt',K,'delimiter',' ','precision','%.3f')