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