How to properly make a circular histogram?

I have data points which are the angle and distance around a central point. I would like to make a histogram of these points so I can visualise their distribution and also perform some analyses on their density etc.

For now, I make a 2D histogram, treating the data points as if they are cartesian X,Y coordinates. I then transform the bin centers into polar coordinates and replot the data, like this:

% create data points
theta = normrnd(0,0.5,1000,1);
rad = 100.*rand(1000,1);

figure
subplot(1,3,1)
polarscatter(theta,rad,10,'k','filled','o')
title('data points')

% prepare bins
theta_bins = linspace(-pi,pi,120);    
rad_bins = 0:5:100;

% 'normal' histogram
f = histcounts2(rad,theta,rad_bins,theta_bins);

subplot(1,3,2)
imagesc(f)
daspect([1 1 1])
title('histogram in cartesian coordinates')

% transform bin coordinates
[theta,rad] = meshgrid(linspace(min(theta_bins),max(theta_bins),length(theta_bins)-1),linspace(min(rad_bins),max(rad_bins),length(rad_bins)-1));
[X,Y] = pol2cart(theta,rad); 

subplot(1,3,3)
s = surf(X,Y,f);
s.EdgeColor = 'none';
view(0,90);
daspect([1 1 1])
axis xy
title('histogram in polar coordinates')

However, the resulting histogram is not really ideal - the bins are not homogeneous because they get larger as the radius increases. This is due to the fact that I made the histogram in a cartesian reference frame.

I was wondering if there is a better way to go about this? Some sort of circular histogram with circular bands of triangular bins? Or hexagonal bins? Can anyone suggest anything?

Thanks for any help.

Answer 1

[Paul]

Hi Neuroesc,

"However, the resulting histogram is not really ideal - the bins are not homogeneous because they get larger as the radius increases."

Assuming that you want the anglular bins to be of the same angular width, one approach would be to define the radial bins nonuniformly such that each polar rectangle has the same area.

For example, assume we want 10 angle bins and 10 radial bins

nthetabin = 10;
nrhobin = 10;
nbin = nthetabin*nrhobin;

Assume the maximum radius is 10 for this example

rmax = 10;

The area of each polar rectangle should then be

binarea = pi*rmax^2/nbin;

Define the edges of the bins in angle

theta = (0:nthetabin)/nthetabin*360;

Define the radial edges such that each polar rectangle has area equal to binarea

rho = zeros(1,nrhobin+1);
for ii = 2:numel(rho)
    % binarea = pi*(rho(ii)^2 - rho(ii-1)^2)/nthetabin;
    rho(ii) = sqrt(binarea*nthetabin/pi + rho(ii-1)^2);
end

The resulting bins would look like this

figure
pax  = polaraxes;
pax.RTick = rho;
pax.ThetaTick = theta;
pax.RLim = [0 rho(end)];

If I did the math correctly (you should verify that), each polar rectangle has equal area.

Not sure if that's what you want. If it isn't, I think we'll need to know more about what the "ideal" binning in polar coordinates should be.

Answer 2

[William Rose]

@Neuroesc,

Check out these items on the file exchange:

https://www.mathworks.com/matlabcentral/fileexchange/65174-wind_rose-wind_direction-wind_speed

https://www.mathworks.com/matlabcentral/fileexchange/167981-wind-rose-with-matlab

The above packages generate hitogram like plots for a 2D distribution, where th e2 dimensions are angle and intensity (i.e. radial distance). If your data is only 1 dimensional (angle only), then the above packages will not be ideal.

[William Rose]

@Neuroesc,

Also see https://www.mathworks.com/matlabcentral/fileexchange/10676-circular-statistics-toolbox-directional-statistics

which includes a function circ_plot() which makes a circular histogram plot.

[William Rose]

@Neuroesc,

[Edit: Correct the formula for area of each bin, which was incorrect in my original version of this comment. The area or each bin is .]

You can normalize the histogram data (f) by the area of each bin:

theta = normrnd(0,0.5,1000,1);
rad = 100.*rand(1000,1);
% prepare bins
theta_bins = linspace(-pi,pi,120);    
rad_bins = 0:5:100;
% 'normal' histogram
f = histcounts2(rad,theta,rad_bins,theta_bins);
rmid=(rad_bins(1:end-1)+rad_bins(2:end))/2; % midpoint radii of bins
binArea=(rmid.*diff(rad_bins))'*diff(theta_bins);
disp([size(f),size(binArea)])  % show that arrays have same dimensions

    20   119    20   119

fNorm=f./binArea;  % hstogram density, normalized by bin area

Then you plot fNorm instead of f.