/[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.1 - (show annotations) (download)
Fri Jun 4 15:50:52 2004 UTC (20 years, 1 month ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint58e_post, checkpoint57v_post, checkpoint53d_post, checkpoint58u_post, checkpoint54a_pre, checkpoint57m_post, checkpoint55c_post, checkpoint54e_post, checkpoint57s_post, checkpoint54a_post, checkpoint57k_post, checkpoint55d_pre, checkpoint57d_post, checkpoint57g_post, checkpoint57b_post, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, checkpoint56b_post, checkpoint57i_post, checkpoint57y_post, checkpoint57e_post, checkpoint55h_post, checkpoint58n_post, checkpoint57g_pre, checkpoint54b_post, checkpoint55b_post, checkpoint58t_post, checkpoint58h_post, checkpoint54d_post, checkpoint56c_post, checkpoint57y_pre, checkpoint55, checkpoint57f_pre, checkpoint57a_post, checkpoint54, checkpoint58q_post, checkpoint54f_post, checkpoint55g_post, checkpoint58j_post, checkpoint55f_post, checkpoint57r_post, checkpoint58, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, eckpoint57e_pre, checkpoint57h_done, checkpoint58f_post, checkpoint53g_post, checkpoint57x_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint57f_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint57c_post, checkpoint55e_post, checkpoint58k_post, checkpoint58v_post, checkpoint53f_post, checkpoint55a_post, checkpoint54c_post, checkpoint58s_post, checkpoint58p_post, checkpoint57j_post, checkpoint58b_post, checkpoint57h_pre, checkpoint58m_post, checkpoint57l_post, checkpoint57h_post, checkpoint56a_post, checkpoint55d_post
New scripts to help use cube and netcdf

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: 5.28 $ $Date: 2001/04/15 11:59:14 $
11
12 error(nargchk(4,5,nargin))
13
14 if prod(size(xi)) ~= prod(size(yi))
15 [yi,xi]=ndgrid(yi,xi);
16 end
17
18 if nargin<6, method = 'linear'; end
19 if ~isstr(method),
20 error('METHOD must be one of ''linear'',''cubic'',''nearest'', or ''v4''.');
21 end
22
23
24 switch lower(method),
25 case 'linear'
26 del = linear(x,y,xi,yi);
27 % case 'cubic'
28 % zi = cubic(x,y,z,xi,yi);
29 % case 'nearest'
30 % zi = nearest(x,y,z,xi,yi);
31 % case {'invdist','v4'}
32 % zi = gdatav4(x,y,z,xi,yi);
33 otherwise
34 error('Unknown method.');
35 end
36
37
38
39 %------------------------------------------------------------
40 function delau = linear(x,y,xi,yi)
41 %LINEAR Triangle-based linear interpolation
42
43 % Reference: David F. Watson, "Contouring: A guide
44 % to the analysis and display of spacial data", Pergamon, 1994.
45
46 siz = size(xi);
47 xi = xi(:); yi = yi(:); % Treat these as columns
48 x = x(:); y = y(:); % Treat these as columns
49
50 % Triangularize the data
51 tri = delaunayn([x y]);
52 if isempty(tri),
53 warning('Data cannot be triangulated.');
54 return
55 end
56
57 % Find the nearest triangle (t)
58 t = tsearch(x,y,tri,xi,yi);
59
60 % Only keep the relevant triangles.
61 out = find(isnan(t));
62 if ~isempty(out), t(out) = ones(size(out)); end
63 tri = tri(t,:);
64
65 % Compute Barycentric coordinates (w). P. 78 in Watson.
66 del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
67 (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
68 w(:,3) = ((x(tri(:,1))-xi).*(y(tri(:,2))-yi) - ...
69 (x(tri(:,2))-xi).*(y(tri(:,1))-yi)) ./ del;
70 w(:,2) = ((x(tri(:,3))-xi).*(y(tri(:,1))-yi) - ...
71 (x(tri(:,1))-xi).*(y(tri(:,3))-yi)) ./ del;
72 w(:,1) = ((x(tri(:,2))-xi).*(y(tri(:,3))-yi) - ...
73 (x(tri(:,3))-xi).*(y(tri(:,2))-yi)) ./ del;
74 w(out,:) = zeros(length(out),3);
75
76 delau.tri=tri;
77 delau.w=w;
78 delau.siz=siz;
79 delau.out=out;
80
81 %------------------------------------------------------------

  ViewVC Help
Powered by ViewVC 1.1.22