v4_3deg/matlab/compute_bin_average.m

function [lon lat bin_average] = compute_bin_average(XC,YC,RAC,LonIncr,LatIncr,saveDir,saveFilename)

XC = XC(:);
YC = YC(:);
RAC = RAC(:);

%% 
%Abhishek code for grid area

% specify latitude and longitude increment
%LonIncr = 4;
%LatIncr = 5;

% calculate lat-long boundaries
lon =-180+LonIncr/2:LonIncr:180;
lat =-90+LatIncr/2:LatIncr:90;

LatIncr = nanmean(diff(lat));
LonIncr = nanmean(diff(lon));

[xx yy] = meshgrid(lon,lat);

nlat = 180/LatIncr;
nlon = 360/LonIncr;

% Calculate gridwise areas
r=6375.*1000;        
fjep=0.5*(nlat+1);
dlat=pi/(nlat-1);
dd=2*r^2*2*pi/nlon*sin(0.5*dlat);

for j=2:nlat-1
    
    dxyp(j) = dd*cos(dlat*(j-fjep));

end

dxyp(1)=2*r^2*2*pi*sin(0.25*dlat)*cos(0.25*(2*pi-dlat))/nlon;
dxyp(nlat)=dxyp(1);

AREA = repmat(dxyp,nlon,1);  % final output

%% 

% define edges of output grid
lat1 = lat - LatIncr/2;
lat1(1) = -90;

lat2 = lat + LatIncr/2;
lat2(end) = 90;

lon1 = lon - LonIncr/2;
lon2 = lon + LonIncr/2;

[LAT1 LON1] = meshgrid(lat1,lon1);
[LAT2 LON2] = meshgrid(lat2,lon2);

%% 

% put XC in same range as LON1 and LON2
ix = find(XC < (min(LON1(:))));

if length(ix) > 0
    
    XC(ix) = XC(ix) + 360;
    
end

clear ix

ix = find(XC >= (max(LON2(:))));

if length(ix) > 0
    
    XC(ix) = XC(ix) - 360;
    
end

%% 

% Compute bin-averaging template
LON1v = LON1(:);
LAT1v = LAT1(:);

LON2v = LON2(:);
LAT2v = LAT2(:);

XCv = XC(:);
YCv = YC(:);

RACv = RAC(:);
AREAv = AREA(:);

bin_average = spalloc(length(LON1v),length(XCv),length(XCv));

for i=1:length(LON1v)
    
    ix = find(XCv >= LON1v(i) & XCv < LON2v(i) & YCv >= LAT1v(i) & YCv < LAT2v(i));
    
    if length(ix) > 0
        
        %bin_average(i,ix) = RACv(ix) / AREAv(i); %flux-conserving bin average
        bin_average(i,ix) = 1/length(ix); %normal bin average
        
    end
    
    disp(num2str(i));
    
end

save([saveDir saveFilename],'lon','lat','AREA','bin_average');

end
1	function [lon lat bin_average] = compute_bin_average(XC,YC,RAC,LonIncr,LatIncr,saveDir,saveFilename)
2
3	XC = XC(:);
4	YC = YC(:);
5	RAC = RAC(:);
6
7	%%
8	%Abhishek code for grid area
9
10	% specify latitude and longitude increment
11	%LonIncr = 4;
12	%LatIncr = 5;
13
14	% calculate lat-long boundaries
15	lon =-180+LonIncr/2:LonIncr:180;
16	lat =-90+LatIncr/2:LatIncr:90;
17
18	LatIncr = nanmean(diff(lat));
19	LonIncr = nanmean(diff(lon));
20
21	[xx yy] = meshgrid(lon,lat);
22
23	nlat = 180/LatIncr;
24	nlon = 360/LonIncr;
25
26	% Calculate gridwise areas
27	r=6375.*1000;
28	fjep=0.5*(nlat+1);
29	dlat=pi/(nlat-1);
30	dd=2r^22pi/nlonsin(0.5*dlat);
31
32	for j=2:nlat-1
33
34	dxyp(j) = ddcos(dlat(j-fjep));
35
36	end
37
38	dxyp(1)=2r^22pisin(0.25dlat)cos(0.25(2pi-dlat))/nlon;
39	dxyp(nlat)=dxyp(1);
40
41	AREA = repmat(dxyp,nlon,1); % final output
42
43	%%
44
45	% define edges of output grid
46	lat1 = lat - LatIncr/2;
47	lat1(1) = -90;
48
49	lat2 = lat + LatIncr/2;
50	lat2(end) = 90;
51
52	lon1 = lon - LonIncr/2;
53	lon2 = lon + LonIncr/2;
54
55	[LAT1 LON1] = meshgrid(lat1,lon1);
56	[LAT2 LON2] = meshgrid(lat2,lon2);
57
58	%%
59
60	% put XC in same range as LON1 and LON2
61	ix = find(XC < (min(LON1(:))));
62
63	if length(ix) > 0
64
65	XC(ix) = XC(ix) + 360;
66
67	end
68
69	clear ix
70
71	ix = find(XC >= (max(LON2(:))));
72
73	if length(ix) > 0
74
75	XC(ix) = XC(ix) - 360;
76
77	end
78
79	%%
80
81	% Compute bin-averaging template
82	LON1v = LON1(:);
83	LAT1v = LAT1(:);
84
85	LON2v = LON2(:);
86	LAT2v = LAT2(:);
87
88	XCv = XC(:);
89	YCv = YC(:);
90
91	RACv = RAC(:);
92	AREAv = AREA(:);
93
94	bin_average = spalloc(length(LON1v),length(XCv),length(XCv));
95
96	for i=1:length(LON1v)
97
98	ix = find(XCv >= LON1v(i) & XCv < LON2v(i) & YCv >= LAT1v(i) & YCv < LAT2v(i));
99
100	if length(ix) > 0
101
102	%bin_average(i,ix) = RACv(ix) / AREAv(i); %flux-conserving bin average
103	bin_average(i,ix) = 1/length(ix); %normal bin average
104
105	end
106
107	disp(num2str(i));
108
109	end
110
111	save([saveDir saveFilename],'lon','lat','AREA','bin_average');
112
113	end