| 1 |
gforget |
1.1 |
function [val]=gcmfaces_interp_2d(fld,lon,lat,varargin); |
| 2 |
|
|
%object: linearly interpolate field to given positions |
| 3 |
|
|
%inputs: fld is a 2D gcmfaces field |
| 4 |
|
|
% lon,lat are position vectors |
| 5 |
|
|
%optional: doNearFill (1 by default) to use the nearest |
| 6 |
|
|
% neighbor to complement linear interp |
| 7 |
|
|
%outputs: val is the vector of interpolated values |
| 8 |
|
|
% |
| 9 |
|
|
%pre-requisite: generate the delaunay triangulation using gcmfaces_bindata |
| 10 |
|
|
%assumption: fld should show NaN for missing values |
| 11 |
|
|
|
| 12 |
|
|
warning('off','MATLAB:dsearch:DeprecatedFunction'); |
| 13 |
|
|
|
| 14 |
|
|
global mygrid mytri myenv; |
| 15 |
|
|
|
| 16 |
|
|
%inputs and pre-requisites |
| 17 |
|
|
if ~isfield(myenv,'useDelaunayTri'); |
| 18 |
|
|
myenv.useDelaunayTri=~isempty(which('DelaunayTri')); |
| 19 |
|
|
end; |
| 20 |
|
|
|
| 21 |
|
|
if isempty(mytri); |
| 22 |
|
|
error('missing triangulation (mytri; from gcmfaces_bindata)'); |
| 23 |
|
|
end; |
| 24 |
|
|
|
| 25 |
|
|
if nargin>3; |
| 26 |
|
|
doNearFill=varargin{1}; |
| 27 |
|
|
else; |
| 28 |
|
|
doNearFill=1; |
| 29 |
|
|
end; |
| 30 |
|
|
|
| 31 |
|
|
%switch longitude range to -180+180 or 0-360 according to grid |
| 32 |
|
|
if max(mygrid.XC)<0; |
| 33 |
|
|
lon(find(lon>180))=lon(find(lon>180))-360; |
| 34 |
|
|
end; |
| 35 |
|
|
|
| 36 |
|
|
if max(mygrid.XC)>180; |
| 37 |
|
|
lon(find(lon<0))=lon(find(lon<0))+360; |
| 38 |
|
|
end; |
| 39 |
|
|
|
| 40 |
|
|
%do the actual interpolation |
| 41 |
|
|
if ~myenv.useDelaunayTri;%(code from old griddata.m) |
| 42 |
|
|
|
| 43 |
|
|
% Find the nearest triangle (t) |
| 44 |
|
|
x=convert2array(mygrid.XC); x=x(mytri.kk); |
| 45 |
|
|
y=convert2array(mygrid.YC); y=y(mytri.kk); |
| 46 |
|
|
VEC=convert2array(fld); VEC=VEC(mytri.kk); |
| 47 |
|
|
t = tsearch(x,y,mytri.TRI,lon',lat')';%the order of dims matters!! |
| 48 |
|
|
|
| 49 |
|
|
% Only keep the relevant triangles. |
| 50 |
|
|
out = find(isnan(t)); |
| 51 |
|
|
if ~isempty(out), t(out) = ones(size(out)); end |
| 52 |
|
|
tri = mytri.TRI(t(:),:); |
| 53 |
|
|
|
| 54 |
|
|
% Compute Barycentric coordinates (w). P. 78 in Watson. |
| 55 |
|
|
del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ... |
| 56 |
|
|
(x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1))); |
| 57 |
|
|
w(:,3) = ((x(tri(:,1))-lon(:)).*(y(tri(:,2))-lat(:)) - ... |
| 58 |
|
|
(x(tri(:,2))-lon(:)).*(y(tri(:,1))-lat(:))) ./ del; |
| 59 |
|
|
w(:,2) = ((x(tri(:,3))-lon(:)).*(y(tri(:,1))-lat(:)) - ... |
| 60 |
|
|
(x(tri(:,1))-lon(:)).*(y(tri(:,3))-lat(:))) ./ del; |
| 61 |
|
|
w(:,1) = ((x(tri(:,2))-lon(:)).*(y(tri(:,3))-lat(:)) - ... |
| 62 |
|
|
(x(tri(:,3))-lon(:)).*(y(tri(:,2))-lat(:))) ./ del; |
| 63 |
|
|
|
| 64 |
|
|
val = sum(VEC(tri) .* w,2); |
| 65 |
gforget |
1.2 |
val(out)=NaN; |
| 66 |
gforget |
1.1 |
val=reshape(val,size(lon)); |
| 67 |
|
|
|
| 68 |
|
|
else;%(use TriScatteredInterp) |
| 69 |
|
|
|
| 70 |
|
|
%get interpolant |
| 71 |
|
|
VEC=convert2array(fld); VEC=VEC(mytri.kk); |
| 72 |
|
|
F = TriScatteredInterp(mytri.TRI, VEC); |
| 73 |
|
|
%do interpolation |
| 74 |
|
|
val=F(lon,lat); |
| 75 |
|
|
|
| 76 |
|
|
end; |
| 77 |
|
|
|
| 78 |
|
|
if doNearFill;%use the nearest neighbor to complement linear interp |
| 79 |
|
|
kk=gcmfaces_bindata(lon(:),lat(:)); |
| 80 |
|
|
ARR=convert2array(fld); |
| 81 |
|
|
val2=ARR(kk); |
| 82 |
|
|
val(isnan(val))=val2(isnan(val)); |
| 83 |
|
|
end; |
| 84 |
|
|
|
| 85 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph