- 积分
- 234
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-4-21
- 最后登录
- 1970-1-1
|
发表于 2017-12-31 09:25:28
|
显示全部楼层
楼主提出的通过插值的方法把地形文件插值到其他层次,我在这里做了一个较为简单的方法。
思路是假设地形随高度线性变化的。
MATLAB程序如下
clc
clear
in=fopen(~\orog.bin','r')
orog0=fread(in,'float32');
fclose(in);
orog=reshape(orog0,144,73,17);
lev1=[1000;925;850;700;600;500;400;300;250;200;150;100;70;50;30;20;10];
lev2=[1000;975;950;925;900;875;850;825;800;775;750;700;650;600;550;500;450;400;350;300;250;225;200;175;150;125;100];
for i=1:12
for j=1:27
if lev1(i)==lev2(j)
orog_27(:,:,j)=orog(:,:,i);
end
end
end
for i=1:27
if orog_27(1,1,i)==0
for n=1:10
if orog_27(1,1,i+n)~=0
break
end
end
for j=0,n
orog_27(:,:,i+j)=orog_27(:,:,i-1)+(orog_27(:,:,i+n)-orog_27(:,:,i-1))/(lev2(i+n)-lev2(i-1))*(lev2(i+j)-lev2(i-1));
end
end
end
out=fopen('~\orog_27.dat','w')
for l=1:27
for y=1:73
fwrite(out,orog_27(:,y,l),'float32');
end
end
结果用30°N,60°E-130°E的剖面测试了一下,放不了图就很尴尬,对于2.5°×2.5°的精度来说,误差在可接受范围内。
存在的问题:地形随高度不是线性变化的,有可能有较大误差。
另一种思路,用相邻格点的值进行插值,还没有尝试
才疏学浅,在此仅抛砖引玉,希望有大佬能给出更好的方案 |
|