function [A,X,Y] = binning_std(xin,yin,zin,xres,yres,westl,eastl,southl,northl)

% Bins data into a grid of equal resolution gres in x and y
% xin, yin and zin are column vector all with the same size
% global northl southl westl eastl xres yres
% Aug 2, 2011. Adapted for the entire Southern Ocean. Gets rid of profiles
% or data points outside the grid (such that we have nonzero values in
% xbin, ybin and hence accumarray doesn't give error)
% 8/28/2020: modified to do binning in depth coordinate (reversing axis to
% use negative values)

% xin=WINDEVE(:,1);
% yin=WINDEVE(:,2);
% zin=WINDEVE(:,3);

% xin=LON(ind);
% yin=PRES(ind);
% zin=SIGMA(ind);

% Define Grid
xedges=westl:xres:eastl;
yedges=southl:yres:northl;

if southl<northl
yedges=southl:yres:northl;
elseif southl>northl % reverse, e.g. if working with depth
southl=-1*southl;
northl=-1*northl; 
yedges=southl:yres:northl;
yin=-1*yin;
end

lnx=length(xedges); lny=length(yedges);
%xborder=westl*ones(lny,1);
%yborder=southl*ones(lnx,1);
gsX=abs(westl-eastl)/xres; % Grid Size will be 151 x 76
gsY=abs(northl-southl)/yres;


% Get rid of all profiles outside the grid boundaries
        % Define polygon
        RLat=[northl northl southl southl northl];
        RLon=[westl eastl eastl westl westl]; % W long

        % Extract region
        in=inpolygon(xin,yin,RLon,RLat);
        
         xin=xin(in);
         yin=yin(in);
         zin=zin(in);


% Create a border with nan such that in every bin there's a value (or nan),
% otherwise accumarray doesn't work (i.e. a dummy border that then we crop)
% maxgsX=360/0.25;
% maxgsY=180/0.25;
% fxedges=-180:xres:180;
% fyedges=-90:yres:90;

y = [yin; yedges']; % to add border with NaN values in z
x = [xin; xedges'];
    
if length(y)==length(x)
    z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif length(y)>length(x) % adjust x to the length of y
    nrep=ceil((length(y)-length(x))./length(xedges));
    xfaked=repmat(xedges,1,nrep);
    x = [x; xfaked']; 
    x = x(1:length(y));
    z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif length(x)>length(y) % adjust y to the length of x
    nrep=ceil((length(x)-length(y))./length(yedges));
    yfaked=repmat(yedges,1,nrep);
    y = [y; yfaked'];
    y = y(1:length(x));
    z = [zin; nan(lnx,1)]; % need to modify the faked brder accordingly if grid size changes!
end
% Note: 
% % x,y,z have to be all the same size
% % need to to add a border of nan values to the grid so that at list there's
% % 1 value (or nan) in each bin

% Use histc to get the index for each bin.
[nx, xbin]=histc(x,xedges);
[ny, ybin]=histc(y,yedges);

% accumarray computes nanmean on values corresponding to each bin.
A = accumarray([xbin ybin],z,[],@nanstd);
A=A';
A(A==0)=NaN;

% Remove faked edge
A=A(1:end-1,1:end-1);

% Could implement the use of hist3, instead of histc twice
% edges{1} = xedges;
% edges{2} = yedges;
% N = hist3([x y],'Edges',edges); %hist3 ignores NaN
% N=N';
% meanZ = A./N;


% Jan 27, 2016: Added grid as output:
% Grid to plot in a map
[X,Y]=meshgrid(westl+xres/2:xres:eastl-xres/2,southl+yres/2:yres:northl-yres/2);
X=X-xres/2;
Y=Y-yres/2;