matlab_class/gcmfaces_calc/gcmfaces_bindata.m

function [varargout]=gcmfaces_bindata(varargin);
%object:    compute delaunay triangulation, then bin averaging
%inputs:    are all optional, triggering different sections of the code
%   if   none          then generate the very triangulation (myTri)
%   if   lon,lat       vectors, then compute closest neighbor index(ices)
%   if   lon,lat,obs   vectors, then further do the bin average
%outputs:   are all optional, triggering different sections of the code
%   if   OBS            then return the bin aberage
%   if   OBS,NOBS       then return the bin sum and count
%
%pre-requisite : the grid of interest must have been loaded with grid_load
%
%note to self: should mytri become part of mygrid ?

warning('off','MATLAB:dsearch:DeprecatedFunction');

gcmfaces_global; global mytri;

if ~isfield(myenv,'useDelaunayTri');
    myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
end;

%1) check that triangulation is up to date
if isempty(mytri);
  mytri.dirGrid=mygrid.dirGrid; mytri.nFaces=mygrid.nFaces;
  mytri.fileFormat=mygrid.fileFormat; mytri.ioSize=mygrid.ioSize;
end;

test1=strcmp(mytri.dirGrid,mygrid.dirGrid)&(mytri.nFaces==mygrid.nFaces)&...
      strcmp(mytri.fileFormat,mygrid.fileFormat)&(sum(mytri.ioSize~=mygrid.ioSize)==0);
test2=isfield(mytri,'TRI');

%2) update tile arrays if not up to date
if (~test1&~test2)|(nargin==0);
    %
    mytri=[];
    mytri.dirGrid=mygrid.dirGrid; mytri.nFaces=mygrid.nFaces;
    mytri.fileFormat=mygrid.fileFormat; mytri.ioSize=mygrid.ioSize;
    %generate delaunay triangulation:
    XC=convert2array(mygrid.XC);
    YC=convert2array(mygrid.YC);
    kk=find(~isnan(XC)); mytri.kk=kk;
    [mytri.ii,mytri.jj]=find(~isnan(XC));
    if myenv.useDelaunayTri;
        %the new way, using DelaunayTri&nearestNeighbor
        mytri.TRI=DelaunayTri(XC(kk),YC(kk));
    else;
        %the old way, using delaunay&dsearch
        TRI=delaunay(XC(kk),YC(kk)); nxy = length(kk);
        Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
        mytri.XC=XC(mytri.kk); mytri.YC=YC(mytri.kk); 
        mytri.TRI=TRI; mytri.Stri=Stri;
        %usage: kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
        %       where lon and lat are vector of position
    end
elseif nargin==2;
    %compute grid point vector associated with lon/lat vectors
    lon=varargin{1}; lat=varargin{2};
    if myenv.useDelaunayTri;
        kk = mytri.TRI.nearestNeighbor(lon,lat);
    else;
        kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
    end;
    if nargout==1;
        varargout={mytri.kk(kk)};
    elseif nargout==2;
        ii=mytri.ii(kk); jj=mytri.jj(kk);
        varargout={ii}; varargout(2)={jj};
    else;
        error('wrong output choice');
    end;
elseif nargin==3;
    %do the bin average (if nargout==1) or the bin sum+count (if nargout==2)
    lon=varargin{1}; lat=varargin{2}; obs=varargin{3};
    ii=find(~isnan(obs)); lon=lon(ii); lat=lat(ii); obs=obs(ii);
    if myenv.useDelaunayTri;
        kk = mytri.TRI.nearestNeighbor(lon,lat);
    else;
        kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
    end;
    kk=mytri.kk(kk);
    
    OBS=convert2array(mygrid.XC); 
    OBS(:)=0; NOBS=OBS;
    for k=1:length(kk)
        NOBS(kk(k))=NOBS(kk(k))+1;
        OBS(kk(k))=OBS(kk(k))+obs(k);
    end  % k=1:length(kk)
    
    if nargout==1;%output bin average
        in=find(NOBS); OBS(in)=OBS(in)./NOBS(in);
        in=find(~NOBS); OBS(in)=NaN; NOBS(in)=NaN;
        varargout={convert2array(OBS)};
    elseif nargout==2;%output bin sum+count
        OBS=convert2array(OBS);
        NOBS=convert2array(NOBS);
        varargout={OBS}; varargout(2)={NOBS};
    else;
        error('wrong output choice');
    end;
    
else;
    error('wrong input choice');
end;

warning('on','MATLAB:dsearch:DeprecatedFunction');


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	gforget	1.8	%pre-requisite : the grid of interest must have been loaded with grid_load
12			%
13			%note to self: should mytri become part of mygrid ?
14	gforget	1.1
15	gforget	1.2	warning('off','MATLAB:dsearch:DeprecatedFunction');
16
17	gforget	1.8	gcmfaces_global; global mytri;
18	gforget	1.6
19			if ~isfield(myenv,'useDelaunayTri');
20			myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
21			end;
22	gforget	1.1
23	gforget	1.8	%1) check that triangulation is up to date
24			if isempty(mytri);
25			mytri.dirGrid=mygrid.dirGrid; mytri.nFaces=mygrid.nFaces;
26			mytri.fileFormat=mygrid.fileFormat; mytri.ioSize=mygrid.ioSize;
27			end;
28
29			test1=strcmp(mytri.dirGrid,mygrid.dirGrid)&(mytri.nFaces==mygrid.nFaces)&...
30			strcmp(mytri.fileFormat,mygrid.fileFormat)&(sum(mytri.ioSize~=mygrid.ioSize)==0);
31			test2=isfield(mytri,'TRI');
32
33			%2) update tile arrays if not up to date
34			if (~test1&~test2)\|(nargin==0);
35			%
36			mytri=[];
37			mytri.dirGrid=mygrid.dirGrid; mytri.nFaces=mygrid.nFaces;
38			mytri.fileFormat=mygrid.fileFormat; mytri.ioSize=mygrid.ioSize;
39	gforget	1.5	%generate delaunay triangulation:
40			XC=convert2array(mygrid.XC);
41			YC=convert2array(mygrid.YC);
42	gforget	1.7	kk=find(~isnan(XC)); mytri.kk=kk;
43			[mytri.ii,mytri.jj]=find(~isnan(XC));
44	gforget	1.6	if myenv.useDelaunayTri;
45			%the new way, using DelaunayTri&nearestNeighbor
46			mytri.TRI=DelaunayTri(XC(kk),YC(kk));
47			else;
48			%the old way, using delaunay&dsearch
49	gforget	1.7	TRI=delaunay(XC(kk),YC(kk)); nxy = length(kk);
50	gforget	1.6	Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
51	gforget	1.7	mytri.XC=XC(mytri.kk); mytri.YC=YC(mytri.kk);
52			mytri.TRI=TRI; mytri.Stri=Stri;
53	gforget	1.6	%usage: kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
54			% where lon and lat are vector of position
55			end
56	gforget	1.1	elseif nargin==2;
57	gforget	1.5	%compute grid point vector associated with lon/lat vectors
58			lon=varargin{1}; lat=varargin{2};
59	gforget	1.6	if myenv.useDelaunayTri;
60			kk = mytri.TRI.nearestNeighbor(lon,lat);
61			else;
62			kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
63			end;
64	gforget	1.5	if nargout==1;
65			varargout={mytri.kk(kk)};
66			elseif nargout==2;
67			ii=mytri.ii(kk); jj=mytri.jj(kk);
68			varargout={ii}; varargout(2)={jj};
69			else;
70			error('wrong output choice');
71			end;
72	gforget	1.1	elseif nargin==3;
73	gforget	1.5	%do the bin average (if nargout==1) or the bin sum+count (if nargout==2)
74			lon=varargin{1}; lat=varargin{2}; obs=varargin{3};
75			ii=find(~isnan(obs)); lon=lon(ii); lat=lat(ii); obs=obs(ii);
76	gforget	1.6	if myenv.useDelaunayTri;
77			kk = mytri.TRI.nearestNeighbor(lon,lat);
78			else;
79			kk=dsearch(mytri.XC,mytri.YC,mytri.TRI,lon,lat,mytri.Stri);
80			end;
81	gforget	1.5	kk=mytri.kk(kk);
82
83	gforget	1.6	OBS=convert2array(mygrid.XC);
84	gforget	1.5	OBS(:)=0; NOBS=OBS;
85			for k=1:length(kk)
86			NOBS(kk(k))=NOBS(kk(k))+1;
87			OBS(kk(k))=OBS(kk(k))+obs(k);
88			end % k=1:length(kk)
89
90			if nargout==1;%output bin average
91			in=find(NOBS); OBS(in)=OBS(in)./NOBS(in);
92			in=find(~NOBS); OBS(in)=NaN; NOBS(in)=NaN;
93			varargout={convert2array(OBS)};
94			elseif nargout==2;%output bin sum+count
95			OBS=convert2array(OBS);
96			NOBS=convert2array(NOBS);
97			varargout={OBS}; varargout(2)={NOBS};
98			else;
99			error('wrong output choice');
100			end;
101
102	gforget	1.1	else;
103	gforget	1.5	error('wrong input choice');
104	gforget	1.1	end;
105
106	gforget	1.2	warning('on','MATLAB:dsearch:DeprecatedFunction');
107	gforget	1.1
108
109