interp_dem.m

来自「StaMps最新测试版」· M 代码 · 共 58 行

function []=interp_dem(width,infile,outfile)
%INTERP_DEM Interpolate DEM 5 times in azimuth
%   INTERP_DEM(WIDTH,INFILE,OUTFILE) WIDTH is width of dem in pixels
%
%   Andy Hooper, Mar 2005
%
% ======================================================
% 08/06 AH: generalized to allow variable filenames
% 11/06 AH: interp used instead of interp2 - less memory
% 05/07 AH: more changes to use less memory
% ======================================================

if nargin < 2
   infile='refdem_1l.raw';
end

if nargin < 3
   outfile='refdem_1li.raw';
end

fid=fopen(infile);
dem=fread(fid,[width,inf],'float');
fclose(fid);
dem=dem';
[ny,nx]=size(dem);
ix=find(isnan(dem));
dem_ix=logical(ones(size(dem)));
dem_ix(ix)=0;
demi1=dem;
[X_patch,Y_patch]=meshgrid([1:11],[1:11]);

for i=ix'
    i2=ceil(i/ny);
    i1=i-(i2-1)*ny;
    ix1=max(1,i1-5):min(i1+5,ny);
    ix2=max(1,i2-5):min(i2+5,nx);
    dem_patch=dem(ix1,ix2);
    ix_patch=dem_ix(ix1,ix2);
    if sum(ix_patch(:))<40
        ix_patch(:)=1;
    end
    [X_patch,Y_patch]=meshgrid([1:length(ix2)],[1:length(ix1)]);
    demi1(i)=griddata(X_patch(ix_patch),Y_patch(ix_patch),dem_patch(ix_patch),i2-ix2(1)+1,i1-ix1(1)+1,'linear',{'QJ'});
end
demi1(isnan(demi1))=0;
demi2=zeros(size(demi1,1)*5,size(demi1,2));

for i=1:size(demi1,2)
    demi2(:,i)=interp(demi1(:,i),5); % oversample
end

demi2=single(demi2);
demi2=[demi2(1,:);demi2(1,:);demi2(1:end-2,:)];

fid=fopen(outfile,'w');
fwrite(fid,demi2','float');
fclose(fid);