util1_jan/gridinterp.m

   1 function vi=gridinterp(x,y,v,xi,yi)
   2 % vi=gridinterp(x,y,v,xi,yi)
   3 %
   4 % Interpolate from 2D irregular rectangular grid
   5 % In:
   6 %   x   matrix of x coordinates of the grid nodes
   7 %   y   matrix of y coordinates of the grid nodes
   8 %   v   matrix of values at the grid nodes
   9 %   xi  vector of x coordinates of points to interpolate to
  10 %   yi  vector of y coordinates to interpolate to
  11 % Out:
  12 %   vi  values at points (xi, yi)
  13
  14 if any(size(xi)~=size(yi)),
  15     error('xi and yi must be same size')
  16 end
  17 [mi,ni]=size(xi);
  18
  19 [mm,nn]=size(x);
  20 [ii,jj]=ndgrid(1:mm,1:nn);
  21
  22 vi=zeros(mi,ni);
  23 disp('interpolating')
  24 for i=1:mi
  25     for j=1:ni
  26         xx=xi(i,j);
  27         yy=yi(i,j);
  28         if ~isnan(xx) & ~isnan(yy),   % skip any NaNs
  29             % find the nearest point
  30             d2=(x-xx).*(x-xx)+(y-yy).*(y-yy);
  31             k=find(d2(:)==min(d2(:)));k=k(1);
  32             ix=ii(k);
  33             iy=jj(k);
  34             if ix>1 & ix < mm & iy >1 & iy < nn, % if not on the boundary
  35                 x3=x(ix-1:ix+1,iy-1:iy+1);
  36                 y3=y(ix-1:ix+1,iy-1:iy+1);
  37                 v3=v(ix-1:ix+1,iy-1:iy+1);
  38                 v_interp = scatteredInterpolant(x3(:),y3(:),v3(:));
  39                 vv = v_interp(xx,yy);
  40                 % disp(x3);disp(xx);disp(y3),disp(yy),disp(v3),disp(vv)
  41                 % hold off; mesh(x3,y3,v3);hold on;plot3(xx,yy,vv,'*k');hold off;drawnow
  42             else
  43                 vv = NaN; % boundary points get NaN
  44             end
  45         else
  46             vv=NaN;
  47         end
  48         vi(i,j)=vv;
  49     end
  50 end
  51 end