阅读:2311回复:1
直接对SRTM90m进行栅格内插的程序代码
直接对SRTM90m进行栅格内插的程序代码
|
|
1楼#
发布于:2008-05-30 12:44
<P>没学过MATLAB谁能帮我改成vb</P>
<P 0cm 0cm 0pt"><STRONG>function srtmfil(file_name1,file_name2,num,method)<BR>% SRTM data processing..............................................<BR>% Function: srtmfil V1.0 <BR>% Fills the blank hole of SRTM3 data set<BR>% Use Matlab 'griddata' function<BR>% Usage:<BR>% srtmfil(input_file_name,output_file_name,size_of_srtm_matrix,interpolation_method)<BR>% size_of_srtm_matrix: 1201 for SRTM3, 3601 for SRTM1<BR>% interpolation_method: Switch from 1-4<BR>% 1= 'linear' - Triangle-based linear interpolation (default).<BR>% 2= 'cubic' - Triangle-based cubic interpolation.<BR>% 3= 'nearest' - Nearest neighbor interpolation.<BR>% <BR>% Programed by: <BR>% Frank Chen @ </STRONG><a href="http://www.lvye.org/" target="_blank" ><STRONG>www.LVYE.ORG</STRONG></A><BR><STRONG>% History:<BR>% V1.0 : 11-Jan-2004 on Matlab <st1:chsdate w:st="on" IsROCDate="False" IsLunarDate="False" Day="30" Month="12" Year="1899">6.1.0</st1:chsdate>, by Frank chen<BR>%<BR><BR>% Read the SRTM data into memory...<BR>display('Read SRTM data ....')<BR>fid=fopen(file_name1,'r','s');<BR>[Z0,n]=fread(fid,[num,num],'int16');<BR>fclose(fid);<BR><BR>% Process data set band by band to solve memory limitation problem...<BR>% band width=100 ; step = 50 for SRTM3 data<BR><BR>b_width=100;b_step=50;<BR><BR>for band=1:b_step:num-b_width<BR><BR>display(band)<BR><BR>% Cut a band for processing... <BR>Zi=Z0(band:band+b_width, 1:num);<BR>size(Zi);<BR><BR>%% meshgrid<BR>[Xi,Yi]=meshgrid(1:num,band:band+b_width);<BR>size(Xi);<BR><BR>%display('srtm data sorting....')<BR>a=[Xi(:) Yi(:) Zi(:)];<BR>b=sortrows(a,3);<BR>[m,n]=size(b);<BR><BR>%display('srtm data sorting finished....')<BR>for i=1:m<BR>if(b(i,3) > -2000) <BR>j=i;<BR>break;<BR>end<BR>end<BR><BR>% b1=b(1:j-1,:);<BR>% b1=sortrows(b1,[1 2]);<BR>% size(b1)<BR><BR>b2=b(j:m,:);<BR>b2=sortrows(b2,[1 2]);<BR>size(b2);<BR><BR>% SRTM data griding ....<BR>x=b2(:,1); y=b2(:,2); z=b2(:,3);<BR><BR>switch method<BR>case 1,<BR>Zi=griddata(x,y,z,Xi,Yi, 'linear');<BR>case 2,<BR>Zi=griddata(x,y,z,Xi,Yi, 'cubic'); <BR>case 3,<BR>Zi=griddata(x,y,z,Xi,Yi, 'nearest');<BR>otherwise,<BR>Zi=griddata(x,y,z,Xi,Yi, 'linear'); <BR>end<BR><BR>size(Zi);<BR>% merge the bands .....<BR>if band ==1 Z = Zi(1:b_width-(b_step/2)-1,: );<BR>elseif band==num-b_width Z=[Z ; Zi((b_step/2):b_width+1, :)];<BR>else Z=[Z ; Zi((b_step/2):b_width-(b_step/2)-1,:)];<BR>end<BR>size(Z)<BR><BR>% loop control of bands...<BR>end <BR><BR>% SRTM data saving ....<BR>fid=fopen(file_name2,'w','s')<BR>fwrite(fid,Z,'int16');<BR>fclose(fid);</STRONG></P> |
|