/[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.4 - (show annotations) (download)
Thu Jul 11 12:41:06 2013 UTC (11 years ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64k, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.3: +28 -10 lines
if available use DelaunayTri instead of delaunay and tseach, code by
Yuan Lian

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.3 $ $Date: 2009/05/29 14:03:30 $
11
12 % $Header: /u/gcmpack/MITgcm/utils/matlab/griddata_preprocess.m,v 1.3 2009/05/29 14:03:30 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 % older version of matlab do not have DelaunayTri, later versions do not
58 % have tsearch
59 msg=which('DelaunayTri');
60 if length(msg) > 0 % version is 2012 or newer
61 % Triangularize the data
62 tri=DelaunayTri(x,y);
63 if isempty(tri),
64 warning('Data cannot be triangulated.');
65 return
66 end
67
68 % Find the nearest triangle (t)
69 [t, bcs] = pointLocation(tri, xi, yi);
70
71 else % use delaunay and tsearch (and hope it is available)
72
73 % Triangularize the data
74 tri = delaunayn([x y]);
75 if isempty(tri),
76 warning('Data cannot be triangulated.');
77 return
78 end
79
80 % Find the nearest triangle (t)
81 t = tsearch(x,y,tri,xi,yi);
82
83 end % end of selecting version
84
85 % Only keep the relevant triangles.
86 out = find(isnan(t));
87 if ~isempty(out), t(out) = ones(size(out)); end
88 tri = tri(t,:);
89
90 % Compute Barycentric coordinates (w). P. 78 in Watson.
91 del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
92 (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
93 w(:,3) = ((x(tri(:,1))-xi).*(y(tri(:,2))-yi) - ...
94 (x(tri(:,2))-xi).*(y(tri(:,1))-yi)) ./ del;
95 w(:,2) = ((x(tri(:,3))-xi).*(y(tri(:,1))-yi) - ...
96 (x(tri(:,1))-xi).*(y(tri(:,3))-yi)) ./ del;
97 w(:,1) = ((x(tri(:,2))-xi).*(y(tri(:,3))-yi) - ...
98 (x(tri(:,3))-xi).*(y(tri(:,2))-yi)) ./ del;
99 w(out,:) = zeros(length(out),3);
100
101 delau.tri=tri;
102 delau.w=w;
103 delau.siz=siz;
104 delau.out=out;
105
106 %------------------------------------------------------------

  ViewVC Help
Powered by ViewVC 1.1.22