matlab_class/sample_processing/example_bin_average.m

function []=example_bin_average(choiceV3orV4,doSmooth);
%object:    a bin average example within gcmfaces
%inputs:    choiceV3orV4 ('v3' or 'v4') selects the sample GRID
%           doSmooth; if set to 1, follow on with smoothing example

%%%%%%%%%%%%%%%%%
%load parameters:
%%%%%%%%%%%%%%%%%

gcmfaces_global;
dir0=[myenv.gcmfaces_dir '/sample_input/'];
dirGrid=[dir0 '/GRID' choiceV3orV4  '/'];
dirIn=[dir0 '/SAMPLE' choiceV3orV4  '/'];
if strcmp(choiceV3orV4,'v4'); 
    nF=5; fileFormat='compact'; 
else; 
    nF=1; fileFormat='straight'; 
end;
grid_load(dirGrid,nF,fileFormat);

%%%%%%%%%%%%%%%%%%%%%%%
%generate delaunay triangulation

XC=convert2array(mygrid.XC);
YC=convert2array(mygrid.YC);

%needed when 0-360 longitude convention
if mygrid.nFaces==1;
%   ii=max(find(XC(:,1)<=180)); 
%   XC=circshift(XC,[-ii 0]); 
  XC(XC>180)=XC(XC>180)-360; 
%   YC=circshift(YC,[-ii 0]);
end;

%needed so that we do not loose data
if mygrid.nFaces==5;
    XC(isnan(XC))=270;
    YC(isnan(YC))=180;
end;

TRI=delaunay(XC,YC);
nxy = prod(size(XC));
Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);

%%%%%%%%%%%%%%%%%%%%%%%
%generate random data

nobs=1e6;
lat=(rand(nobs,1)-0.5)*2*90;
lon=(rand(nobs,1)-0.5)*2*180;  xx=find(lon>180);lon(xx)=lon(xx)-360;
obs=(rand(nobs,1)-0.5)*2;

%%%%%%%%%%%%%%%%%%%%%%%%
%bin average random data

OBS=0*XC; NOBS=OBS; 
ik=dsearch(XC,YC,TRI,lon,lat,Stri);
for k=1:length(ik)
        NOBS(ik(k))=NOBS(ik(k))+1;
        OBS(ik(k))=OBS(ik(k))+obs(k);
end  % k=1:length(ik)

in=find(NOBS); OBS(in)=OBS(in)./NOBS(in);
in=find(~NOBS); OBS(in)=NaN; NOBS(in)=NaN;

obsmap=convert2array(OBS);%put in gcmfaces format
obsmapcompact=convert2gcmfaces(obsmap); %put in gcm input format

whos obs* OBS

figure; set(gcf,'Units','Normalized','Position',[0.1 0.3 0.4 0.6]);
imagescnan(OBS','nancolor',[1 1 1]*0.8); axis xy; caxis([-1 1]/2); colorbar;
figure; set(gcf,'Units','Normalized','Position',[0.5 0.3 0.4 0.6]);
imagescnan(NOBS','nancolor',[1 1 1]*0.8); axis xy; caxis([0 25]); colorbar;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%now illustrate the smoothing filter
if doSmooth; example_smooth; end;


1	gforget	1.8	function []=example_bin_average(choiceV3orV4,doSmooth);
2			%object: a bin average example within gcmfaces
3			%inputs: choiceV3orV4 ('v3' or 'v4') selects the sample GRID
4			% doSmooth; if set to 1, follow on with smoothing example
5	gforget	1.1
6			%%%%%%%%%%%%%%%%%
7			%load parameters:
8			%%%%%%%%%%%%%%%%%
9
10	gforget	1.10	gcmfaces_global;
11			dir0=[myenv.gcmfaces_dir '/sample_input/'];
12	gforget	1.8	dirGrid=[dir0 '/GRID' choiceV3orV4 '/'];
13			dirIn=[dir0 '/SAMPLE' choiceV3orV4 '/'];
14	gforget	1.10	if strcmp(choiceV3orV4,'v4');
15			nF=5; fileFormat='compact';
16			else;
17			nF=1; fileFormat='straight';
18			end;
19			grid_load(dirGrid,nF,fileFormat);
20	gforget	1.1
21			%%%%%%%%%%%%%%%%%%%%%%%
22			%generate delaunay triangulation
23
24			XC=convert2array(mygrid.XC);
25			YC=convert2array(mygrid.YC);
26
27	gforget	1.8	%needed when 0-360 longitude convention
28			if mygrid.nFaces==1;
29			% ii=max(find(XC(:,1)<=180));
30			% XC=circshift(XC,[-ii 0]);
31			XC(XC>180)=XC(XC>180)-360;
32			% YC=circshift(YC,[-ii 0]);
33			end;
34
35	gforget	1.1	%needed so that we do not loose data
36	gforget	1.8	if mygrid.nFaces==5;
37			XC(isnan(XC))=270;
38			YC(isnan(YC))=180;
39			end;
40	gforget	1.1
41			TRI=delaunay(XC,YC);
42			nxy = prod(size(XC));
43			Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
44
45			%%%%%%%%%%%%%%%%%%%%%%%
46			%generate random data
47
48			nobs=1e6;
49			lat=(rand(nobs,1)-0.5)290;
50			lon=(rand(nobs,1)-0.5)2180; xx=find(lon>180);lon(xx)=lon(xx)-360;
51			obs=(rand(nobs,1)-0.5)*2;
52
53			%%%%%%%%%%%%%%%%%%%%%%%%
54			%bin average random data
55
56			OBS=0*XC; NOBS=OBS;
57			ik=dsearch(XC,YC,TRI,lon,lat,Stri);
58			for k=1:length(ik)
59			NOBS(ik(k))=NOBS(ik(k))+1;
60			OBS(ik(k))=OBS(ik(k))+obs(k);
61			end % k=1:length(ik)
62
63			in=find(NOBS); OBS(in)=OBS(in)./NOBS(in);
64			in=find(~NOBS); OBS(in)=NaN; NOBS(in)=NaN;
65
66	gforget	1.9	obsmap=convert2array(OBS);%put in gcmfaces format
67	gforget	1.1	obsmapcompact=convert2gcmfaces(obsmap); %put in gcm input format
68
69			whos obs* OBS
70
71			figure; set(gcf,'Units','Normalized','Position',[0.1 0.3 0.4 0.6]);
72			imagescnan(OBS','nancolor',[1 1 1]*0.8); axis xy; caxis([-1 1]/2); colorbar;
73			figure; set(gcf,'Units','Normalized','Position',[0.5 0.3 0.4 0.6]);
74			imagescnan(NOBS','nancolor',[1 1 1]*0.8); axis xy; caxis([0 25]); colorbar;
75
76	gforget	1.8	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77			%now illustrate the smoothing filter
78			if doSmooth; example_smooth; end;
79	gforget	1.1
80