/[MITgcm]/MITgcm/utils/matlab/griddata_preprocess.m
ViewVC logotype

Contents of /MITgcm/utils/matlab/griddata_preprocess.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (show annotations) (download)
Fri May 29 14:03:30 2009 UTC (15 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64i, checkpoint64h, checkpoint64j, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint62, checkpoint63, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.2: +8 -4 lines
fix interpolation to square grid (xi & yi of same size). Pb found by Martha.

1 function [del] = griddata_preprocess(x,y,xi,yi,method)
2 %GRIDDATA_PREPROCESS Pre-calculate Delaunay triangulation for use
3 % with GRIDDATA_FAST.
4 %
5 % DEL = GRIDDATA_PREPROCESS(X,Y,XI,YI)
6
7 % Based on
8 % Clay M. Thompson 8-21-95
9 % Copyright 1984-2001 The MathWorks, Inc.
10 % $Revision: 1.2 $ $Date: 2007/02/17 23:49:43 $
11
12 % $Header: /u/gcmpack/MITgcm/utils/matlab/griddata_preprocess.m,v 1.2 2007/02/17 23:49:43 jmc Exp $
13 % $Name: $
14
15 error(nargchk(4,5,nargin))
16
17 %- To interpolate to a section (with position vector xi,yi), skip the
18 % ndgrid call hereafter. However, since this is poorly documented,
19 % leave the unconditional call to ndgrid to avoid falling into this trap
20 % when doing 2-D interpolation on to a square grid (xi & yi of equal size)
21 %if prod(size(xi)) ~= prod(size(yi))
22 [yi,xi]=ndgrid(yi,xi);
23 %end
24
25 if nargin<6, method = 'linear'; end
26 if ~isstr(method),
27 error('METHOD must be one of ''linear'',''cubic'',''nearest'', or ''v4''.');
28 end
29
30
31 switch lower(method),
32 case 'linear'
33 del = linear(x,y,xi,yi);
34 % case 'cubic'
35 % zi = cubic(x,y,z,xi,yi);
36 % case 'nearest'
37 % zi = nearest(x,y,z,xi,yi);
38 % case {'invdist','v4'}
39 % zi = gdatav4(x,y,z,xi,yi);
40 otherwise
41 error('Unknown method.');
42 end
43
44
45
46 %------------------------------------------------------------
47 function delau = linear(x,y,xi,yi)
48 %LINEAR Triangle-based linear interpolation
49
50 % Reference: David F. Watson, "Contouring: A guide
51 % to the analysis and display of spacial data", Pergamon, 1994.
52
53 siz = size(xi);
54 xi = xi(:); yi = yi(:); % Treat these as columns
55 x = x(:); y = y(:); % Treat these as columns
56
57 % Triangularize the data
58 tri = delaunayn([x y]);
59 if isempty(tri),
60 warning('Data cannot be triangulated.');
61 return
62 end
63
64 % Find the nearest triangle (t)
65 t = tsearch(x,y,tri,xi,yi);
66
67 % Only keep the relevant triangles.
68 out = find(isnan(t));
69 if ~isempty(out), t(out) = ones(size(out)); end
70 tri = tri(t,:);
71
72 % Compute Barycentric coordinates (w). P. 78 in Watson.
73 del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
74 (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
75 w(:,3) = ((x(tri(:,1))-xi).*(y(tri(:,2))-yi) - ...
76 (x(tri(:,2))-xi).*(y(tri(:,1))-yi)) ./ del;
77 w(:,2) = ((x(tri(:,3))-xi).*(y(tri(:,1))-yi) - ...
78 (x(tri(:,1))-xi).*(y(tri(:,3))-yi)) ./ del;
79 w(:,1) = ((x(tri(:,2))-xi).*(y(tri(:,3))-yi) - ...
80 (x(tri(:,3))-xi).*(y(tri(:,2))-yi)) ./ del;
81 w(out,:) = zeros(length(out),3);
82
83 delau.tri=tri;
84 delau.w=w;
85 delau.siz=siz;
86 delau.out=out;
87
88 %------------------------------------------------------------

  ViewVC Help
Powered by ViewVC 1.1.22