function [hout,hcb,indices]=plot_raster_levels(xx,yy,map_data,binEdges,my_colormap,labels);

%my_colormap = [204  204  204
%               153  153  153
%               255  255  178
%               254  204   92
%               253  141   60
%               240   59   32
%               189    0   38
%                 0    0    0]/255 ;
%binEdges = [0 0.001 0.005 0.01 0.05 0.1 0.25 0.5 1] ;

proj = ismap(gca);    % true if the axes are mapped-- then will use pcolorm. 

if nargin <6; 
labels = textscan(num2str(binEdges),'%s') ;
labels = labels{1} ;
end

%labels{length(labels)} = [labels{length(labels)} '%'] ;

%[~,indices] = histc(map_data,binEdges);
[~,~,indices] = histcounts(map_data,binEdges);

indices(isnan(map_data)) = NaN ;
indices(isinf(map_data)) = NaN ;

%figure ;
if proj == 1 
hout=pcolorm(yy,xx,indices-1)   % Instead of image(), to display NaN values as white
else
hout=pcolor(xx,yy,indices-1)   % Instead of image(), to display NaN values as white
end
shading flat

colormap(gca,my_colormap);   % gca as first argument prevents 
                             % colormap from changing for all 
                             % subsequent plots
hcb = colorbar('location','southoutside');
caxis([0 length(binEdges)-1]);
hcb.Ticks=0:(length(binEdges)-1);
hcb.YTickLabel = labels ;
