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	dimitri	1.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