Plotting Polar Plots using time series data

I have time series data for angle and speed (sample attached). I can plot the angel on a polarhistogram as follows.

load("sampleData.mat")
whos

  Name       Size                Bytes  Class     Attributes

  angle      1x15000            120000  double              
  ans        1x38                   76  char                
  speed      1x15000            120000  double              
  time       1x15000            120000  double              

figure
h=polarhistogram(angle);
set(gca,'ThetaTick',     [0 30 60 90 120 150 180 210 240 270 300 330], ...
        'ThetaTickLabel',{0 30 60 90 120 150 180 210 240 270 300 330})
set(gca,'ThetaZeroLocation','top','ThetaDir','clockwise');
h.DisplayStyle = 'stairs';

But I want to include the speed within the plot too, something similar to a countourf plot or Perfect Polar Plots as below. I tried using the custom function, but I couldn't make it work. How can I do this?

[Mathieu NOE]

hello

your angle range is very narrow (178 to 180 °) and I suspect that if you want to plot your density on a polar diagram this will be unreadable as the plot will appear as a very narrow area

instead you can still do a density plot in cartesian coordinates

I am using dscatteravailable here : Flow Cytometry Data Reader and Visualization - File Exchange - MATLAB Central

or here : Ratiometric FRET - File Exchange - MATLAB Central

load('sampleData.mat')

% whos
%   Name       Size                Bytes  Class     Attributes
% 
%   angle      1x15000            120000  double              
%   speed      1x15000            120000  double              
%   time       1x15000            120000  double 

figure(1)
[hAxes,col,ctrs1,ctrs2,F] = dscatter(speed(:),angle(:));
xlabel('speed');
ylabel('angle (°)');
grid on 

hold on
lb = min(F,[],'all');
ub = max(F,[],'all');
threshold = lb+0.1*(ub-lb); % threshold is 10% of F range above min value
[row,col]=find(F>threshold);

% k = boundary(___,s) specifies shrink factor s using any of the previous syntaxes.
% s is a scalar between 0 and 1. Setting s to 0 gives the convex hull,
% and setting s to 1 gives a compact boundary that envelops the points.
% The default shrink factor is 0.5.
s = 0.75;
k = boundary(row,col,s);
plot(ctrs1(col(k)),ctrs2(row(k)),'r')

hold off

[Jake]

Hi @Mathieu NOE, thank you for this suggestion. I agree that this is kind of a better option when the angle range is very small. I like this one and the answer you suggested :)

[Mathieu NOE]

hello Jake

glad I could be of some help !

Answer 1

[Mathieu NOE]

hello again

well, maybe I'm stupid but I coud not yet figure out how to use Perfect Polar Plots without having to spend some time to create a density map with the correct dimensions.

so the quick and dirty alternative I can propose is a modified version of dscatter (mentionned above in my comment ) , so let's call it dscatterpolar (see attached)

sorry to the community and to TMW to have frankensteined the original dscatter function !!!

result :

as i said before , the data does not cover a wide angular range so we're left with that narrow beam plot

the ticks in radial correspond to the wind speed

code : in one line

load('sampleData.mat')

% whos
%   Name       Size                Bytes  Class     Attributes
% 
%   angle      1x15000            120000  double              
%   speed      1x15000            120000  double              
%   time       1x15000            120000  double 

figure(1)
[hAxes,col,ctrs1,ctrs2,F] = dscatterpolar(speed(:),angle(:)*pi/180);  % make sure you pass the data speed first and angle in second (and in rad).

[Mathieu NOE]

EDIT

sligth modification : I finally prefer to use rlim and thetalim in the main code and not inside dscatterpolar function. So use the new (attached) version of it

now you can focus on your data if the angualr range is so small you would not see much in the full 0-360° polarplot

clc
clearvars

load('sampleData.mat')

% whos
%   Name       Size                Bytes  Class     Attributes
% 
%   angle      1x15000            120000  double              
%   speed      1x15000            120000  double              
%   time       1x15000            120000  double 

figure(1)
[hAxes,col,ctrs1,ctrs2,F] = dscatterpolar(speed(:),angle(:)*pi/180);  % make sure you pass the data speed first and angle in second (and in rad).

% set R and angle limits to the plot
rlim(gca,[5*floor(min(speed/5)) 5*ceil(max(speed/5))]); % rounded speed values with 5 increments
thetalim(gca,[165 195])

% add R axis label
rruler = hAxes.RAxis;
rruler.Label.String = 'Speed (m/s)';

% add angle axis label
truler = hAxes.ThetaAxis;
truler.Label.String = 'Direction (°)';

[Mathieu NOE]

hello again

I f you like , this can be a demo for you , but again, as the angular range of your data is very small , you see only a very narrow beam and it's hard to see the density variations

full circle rendering

code :

load('sampleData.mat')

% whos
%   Name       Size                Bytes  Class     Attributes
% 
%   angle      1x15000            120000  double              
%   speed      1x15000            120000  double              
%   time       1x15000            120000  double 

% create 2D histogram data
[n,xedges,yedges] = histcounts2(angle, speed, 'BinMethod', 'auto');
xcenters = xedges(1:end-1) + diff(xedges) / 2;
ycenters = yedges(1:end-1) + diff(yedges) / 2;

% Points between the outmost bin are tricky, using the
% outermost bin (similar to imfilter's replicate option) is
% a decent approximation
xi = [xedges(1) xcenters xedges(end)];
yi = [yedges(1) ycenters yedges(end)];
n = [n(:,1) n n(:,end)];
n = [n(1,:); n; n(end,:)];
if numel(ycenters) == 1
    d = interp1(xi,n',angle);
elseif numel(xcenters) == 1
    d = interp1(yi,n',speed);
else
    [xi,yi] = meshgrid(xi, yi);
    d = interp2(xi, yi, n', angle, speed);
end
   
% Then set up the grid, as with the previous examples:
figure(5); clf
ph = makePolarGrid(...
  'ALabelScheme', 'normal',...
  'RTicks',        15:1:25,...  % Radial ticks (inner circles)
  'RUnits',        '[m/s]',...   % Add units to outer-most radial labels
  'RLabelAngle',  170,...        % Angle where the radial tick labels will go
  'RLabelFormat', '%.0f');      % Format for radial labels     


[md,mf] = meshgrid(xcenters, ycenters);
I = scatteredInterpolant(angle(:), speed(:), d(:),'linear','none');
Z = I(md,mf);
% Now convert to cartesian coords
[px,py] = polgrid2cart(md, mf, ph);

% And then use contourf to plot. Optional formatting adjustments follow.
[~,hc]  = contourf(px, py, Z);
uistack(hc, 'bottom');  % this will put the grid on top of the contours
% hc.LevelList(1:5) = [];    % this removes the countours on lowest level(s) 
hc.LineStyle = 'none';  % enable/disble contour lines
colormap jet
title('Wave Directional Spectrum')

[Mathieu NOE]

tried to narrow down the circle grid with the hope that I would zoom into the area of interest , but I just get a small portion of the full circle rendering, so not very useful here

code :

the only difference in the code is this added parameters :

'ALim',       [173,187],...     % span of circle grid

load('sampleData.mat')

% whos
%   Name       Size                Bytes  Class     Attributes
% 
%   angle      1x15000            120000  double              
%   speed      1x15000            120000  double              
%   time       1x15000            120000  double 

% create 2D histogram data
[n,xedges,yedges] = histcounts2(angle, speed, 'BinMethod', 'auto');
xcenters = xedges(1:end-1) + diff(xedges) / 2;
ycenters = yedges(1:end-1) + diff(yedges) / 2;

% Points between the outmost bin are tricky, using the
% outermost bin (similar to imfilter's replicate option) is
% a decent approximation
xi = [xedges(1) xcenters xedges(end)];
yi = [yedges(1) ycenters yedges(end)];
n = [n(:,1) n n(:,end)];
n = [n(1,:); n; n(end,:)];
if numel(ycenters) == 1
    d = interp1(xi,n',angle);
elseif numel(xcenters) == 1
    d = interp1(yi,n',speed);
else
    [xi,yi] = meshgrid(xi, yi);
    d = interp2(xi, yi, n', angle, speed);
end
   
% Then set up the grid, as with the previous examples:
figure(5); clf
ph = makePolarGrid(...
  'ALabelScheme', 'normal',...
  'ALim',       [173,187],...     % span of circle grid
  'RTicks',        15:1:25,...  % Radial ticks (inner circles)
  'RUnits',        '[m/s]',...   % Add units to outer-most radial labels
  'RLabelAngle',  170,...        % Angle where the radial tick labels will go
  'RLabelFormat', '%.0f');      % Format for radial labels     


[md,mf] = meshgrid(xcenters, ycenters);
I = scatteredInterpolant(angle(:), speed(:), d(:),'linear','none');
Z = I(md,mf);
% Now convert to cartesian coords
[px,py] = polgrid2cart(md, mf, ph);

% And then use contourf to plot. Optional formatting adjustments follow.
[~,hc]  = contourf(px, py, Z);
uistack(hc, 'bottom');  % this will put the grid on top of the contours
% hc.LevelList(1:5) = [];    % this removes the countours on lowest level(s) 
hc.LineStyle = 'none';  % enable/disble contour lines
colormap jet
title('Wave Directional Spectrum')

[Jake]

This is a very nice approach! I have two questions though - one related to this specifically, and the other one slightly different.

• In this generated plot, can we fill the rest of the space too, something similar to the one below (taken from here) with a colorbar?

• I tried plotting two polarhistogram plots with the following code:

clear; close all; clc
load("sampleData_02.mat");
whos

  Name          Size               Bytes  Class     Attributes

  angle_01      1x10000            80000  double              
  angle_02      1x10000            80000  double              
  speed_01      1x10000            80000  double              
  speed_02      1x10000            80000  double              
  time          1x10000            80000  double              

figure
h1=polarhistogram(angle_01); hold on
h2=polarhistogram(angle_02);
set(gca,'ThetaTick',     [0 30 60 90 120 150 180 210 240 270 300 330], ...
        'ThetaTickLabel',{0 30 60 90 120 150 180 210 240 270 300 330})
set(gca,'ThetaZeroLocation','top','ThetaDir','clockwise');

[min(angle_01), max(angle_01); min(angle_02), max(angle_02)]

ans = 2×2

179.3354 180.0000 177.7802 180.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

As can be seen, the angle range is small, but the plot shows the angles very differently. What would be the reason for this, and can it be fixed?

[Mathieu NOE]

be carefull that polarhistogram awaits angle in radians and not in degrees

h1=polarhistogram(angle_01*pi/180); hold on
h2=polarhistogram(angle_02*pi/180);

[Mathieu NOE]

fyi

you can of course do the same dscatter plot with the two data sets

load('sampleData_02.mat')

% whos
%   angle_01      1x10000            80000  double              
%   angle_02      1x10000            80000  double              
%   speed_01      1x10000            80000  double              
%   speed_02      1x10000            80000  double 

figure(1)
[hAxes,col,ctrs1,ctrs2,F] = dscatter(speed_01(:),angle_01(:));
hold on
[hAxes,col,ctrs1,ctrs2,F] = dscatter(speed_02(:),angle_02(:));
xlabel('speed');
ylabel('angle (°)');
grid on 

[Mathieu NOE]

you can do it also in polar plot with

dscatterpolar

as shown already above :

load('sampleData_02.mat')

% whos
%   angle_01      1x10000            80000  double              
%   angle_02      1x10000            80000  double              
%   speed_01      1x10000            80000  double              
%   speed_02      1x10000            80000  double 

figure(1)
[hAxes,col,ctrs1,ctrs2,F] = dscatterpolar(speed_01(:),angle_01(:)*pi/180);  % make sure you pass the data speed first and angle in second (and in rad).
hold on
[hAxes,col,ctrs1,ctrs2,F] = dscatterpolar(speed_02(:),angle_02(:)*pi/180);
grid on 

% set R and angle limits to the plot
rlim(gca,[0 25]); % rounded speed values with N increment
thetalim(gca,[180-10 180+10])

% add R axis label
rruler = hAxes.RAxis;
rruler.Label.String = 'Speed (m/s)';

% add angle axis label
truler = hAxes.ThetaAxis;
truler.Label.String = 'Direction (°)';

[Jake]

Hi @Mathieu NOE, Thank you for all the help - this is good :)

You're right - I made the silly mistake of passing degrees instead of radians to polarhistogram!

[Mathieu NOE]

as always, my pleasure !

tx for accepting my suggestions !