/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_calc/gcmfaces_bindata.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/gcmfaces_bindata.m

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


Revision 1.7 - (hide annotations) (download)
Tue Jul 31 21:10:57 2012 UTC (12 years, 11 months ago) by gforget
Branch: MAIN
Changes since 1.6: +6 -14 lines
- fix neastest neighbor and linear interpolation
  in case of old griddata (DelaunayTri was fine).

1 gforget 1.1 function [varargout]=gcmfaces_bindata(varargin);
2 gforget 1.4 %object: compute delaunay triangulation, then bin averaging
3     %inputs: are all optional, triggering different sections of the code
4     % if none then generate the very triangulation (myTri)
5     % if lon,lat vectors, then compute closest neighbor index(ices)
6     % if lon,lat,obs vectors, then further do the bin average
7     %outputs: are all optional, triggering different sections of the code
8     % if OBS then return the bin aberage
9     % if OBS,NOBS then return the bin sum and count
10     %
11     %note: should mytri become part of mygrid ?
12 gforget 1.1
13 gforget 1.2 warning('off','MATLAB:dsearch:DeprecatedFunction');
14    
15 gforget 1.6 global mygrid mytri myenv;
16    
17     if ~isfield(myenv,'useDelaunayTri');
18     myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
19     end;
20 gforget 1.1
21     if nargin==0;
22 gforget 1.5 %generate delaunay triangulation:
23     XC=convert2array(mygrid.XC);
24     YC=convert2array(mygrid.YC);
25 gforget 1.7 mytri=[];
26     kk=find(~isnan(XC)); mytri.kk=kk;
27     [mytri.ii,mytri.jj]=find(~isnan(XC));
28 gforget 1.6 if myenv.useDelaunayTri;
29     %the new way, using DelaunayTri&nearestNeighbor
30     mytri.TRI=DelaunayTri(XC(kk),YC(kk));
31     else;
32     %the old way, using delaunay&dsearch
33 gforget 1.7 TRI=delaunay(XC(kk),YC(kk)); nxy = length(kk);
34 gforget 1.6 Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
35 gforget 1.7 mytri.XC=XC(mytri.kk); mytri.YC=YC(mytri.kk);
36     mytri.TRI=TRI; mytri.Stri=Stri;
37 gforget 1.6 %usage: kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
38     % where lon and lat are vector of position
39     end
40 gforget 1.1 elseif nargin==2;
41 gforget 1.5 %compute grid point vector associated with lon/lat vectors
42     lon=varargin{1}; lat=varargin{2};
43 gforget 1.6 if myenv.useDelaunayTri;
44     kk = mytri.TRI.nearestNeighbor(lon,lat);
45     else;
46     kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
47     end;
48 gforget 1.5 if nargout==1;
49     varargout={mytri.kk(kk)};
50     elseif nargout==2;
51     ii=mytri.ii(kk); jj=mytri.jj(kk);
52     varargout={ii}; varargout(2)={jj};
53     else;
54     error('wrong output choice');
55     end;
56 gforget 1.1 elseif nargin==3;
57 gforget 1.5 %do the bin average (if nargout==1) or the bin sum+count (if nargout==2)
58     lon=varargin{1}; lat=varargin{2}; obs=varargin{3};
59     ii=find(~isnan(obs)); lon=lon(ii); lat=lat(ii); obs=obs(ii);
60 gforget 1.6 if myenv.useDelaunayTri;
61     kk = mytri.TRI.nearestNeighbor(lon,lat);
62     else;
63     kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
64     end;
65 gforget 1.5 kk=mytri.kk(kk);
66    
67 gforget 1.6 OBS=convert2array(mygrid.XC);
68 gforget 1.5 OBS(:)=0; NOBS=OBS;
69     for k=1:length(kk)
70     NOBS(kk(k))=NOBS(kk(k))+1;
71     OBS(kk(k))=OBS(kk(k))+obs(k);
72     end % k=1:length(kk)
73    
74     if nargout==1;%output bin average
75     in=find(NOBS); OBS(in)=OBS(in)./NOBS(in);
76     in=find(~NOBS); OBS(in)=NaN; NOBS(in)=NaN;
77     varargout={convert2array(OBS)};
78     elseif nargout==2;%output bin sum+count
79     OBS=convert2array(OBS);
80     NOBS=convert2array(NOBS);
81     varargout={OBS}; varargout(2)={NOBS};
82     else;
83     error('wrong output choice');
84     end;
85    
86 gforget 1.1 else;
87 gforget 1.5 error('wrong input choice');
88 gforget 1.1 end;
89    
90 gforget 1.2 warning('on','MATLAB:dsearch:DeprecatedFunction');
91 gforget 1.1
92    
93    

  ViewVC Help
Powered by ViewVC 1.1.22