Using my data and current code, how would I create polarhistograms for wave data?

I am struggling. All I want to do is create a polarhistograms for my wave data containing the information regarding wave energy and direction like this (got image from google images):

However for the life of me i cannot seem to get it to work. I am trying to do a FOR LOOP for five years of data. Here is the code i have written so far:

close all
 clear all
 clc
 
 %Load in the data

 fpath=uigetdir;
 files=dir([fpath '\*.csv']);
 files=char(files.name);

 %for loop to go through the variable files, creating variables for wave
 %height (Hmax [H])), wave power (P kW/m2), max period (Tmax [T]), and geographical direction in
 %degrees (Directiondeg [direction]). Using textscan function to read the
 %data files.

 for n=1:size(files,1)
     fid=fopen([fpath '\' files(n,:)]);
     data{n}=textscan(fid, ['%s' repmat('%f', 1,34)],'Delimiter',',','CollectOutput',1, 'Headerlines',1);
     fclose(fid);

     %readtable to set datetime array for first column in the data.
     T1=readtable([fpath '\' files(n,:)], 'VariableNamingRule','preserve');
     VN=T1.Properties.VariableNames;
     DT=T1.('Date/Time (GMT)');
     T1.('Date/Time (GMT)').TimeZone='UTC';
     T1.('Date/Time (GMT)').TimeZone;

     %Identify and generate wave parameter variables suitable to portray the wave climate of the
     %chosen sites. This identifies wave height (H), wave period (T), wave
     %direction (direction), wave power (power), and time (time).
     H=T1.("Hmax (m)");
     T=T1.("Tmax (s)");
     direction=T1.("Direction (deg)");
     time=T1.("Date/Time (GMT)");
     power=T1.("P (kW/m2)");

     %Identify and replace anomalous data with NaN.
     bad=9999;
     ind=find(H==bad);
     H(ind)=NaN;
     T(ind)=NaN;
     direction(ind)=NaN;
     power(ind)=NaN;
     ind2=find(H >10);
     H(ind2)=NaN;

    %calculate wave frequency and energy
    frequency=1./(T);
    energy=1/8.*1025.*9.81.*(H.^2);
    
    %Plot wave rose diagram with wave energy and wave direction
     %(direction).
     figure(n);
     pax1=polaraxes;
    polarhistogram(deg2rad(direction(H<10)),deg2rad(0:10:360),'DisplayStyle', 'bar')
    hold on
    polarhistogram(deg2rad(direction(H<10)),deg2rad(0:10:360),'Facecolor','#0000FF','FaceAlpha',.6,'DisplayStyle','bar')
    title(['n=' num2str(n) ' ' files(n,:)],'FontSize',12);

All that it gives me are more or less identical plots that has zero colour gradiant.

Please can someone help me, I am very new to MATLAB and struggle with the terminology and understanding the codes when using commands such as "doc polarhistogram".

I have attached two of the data years to this question as well.

Answer 1

[Mathieu NOE]

hello

well , I am not using polar histograms very often , but it seems to me the function polarhistogram does not allow any kind of stacked plot option (I amy be wrong but the doc does not show any example alike)

I spent some time to search what you are looking for , but for the time being, I haven't seen anything very close. need further ingress in that topic

otherwise I made a few corrections in the code

1 / This is not needed as we anyway load the data later on with readtable :

     fid=fopen([fpath '\' files(n,:)]);
     data{n}=textscan(fid, ['%s' repmat('%f', 1,34)],'Delimiter',',','CollectOutput',1, 'Headerlines',1);
     fclose(fid);

2 / The "9999" non valid data correction must be done differently as it appears not in the same rows for the 4 variables (H,T,direction, power).

3/ lastly I opted for this Fex submission,

histogram2Polar - File Exchange - MATLAB Central (mathworks.com)

which at least gives the possibility to create a bivariate histogram in polar coordinates ; you can opt either to plot the counts or a probability density function (pdf)

like that :

code :

     close all
 clear all
 clc
 
 %Load in the data
 fpath=uigetdir;
%   fpath=pwd;
files=dir([fpath '\*.csv']);
 files=char(files.name);
 %for loop to go through the variable files, creating variables for wave
 %height (Hmax [H])), wave power (P kW/m2), max period (Tmax [T]), and geographical direction in
 %degrees (Directiondeg [direction]). Using textscan function to read the
 %data files.
 for n=1:size(files,1)
%      fid=fopen([fpath '\' files(n,:)]);
%      data{n}=textscan(fid, ['%s' repmat('%f', 1,34)],'Delimiter',',','CollectOutput',1, 'Headerlines',1);
%      fclose(fid);
     %readtable to set datetime array for first column in the data.
     T1=readtable(fullfile(fpath,files(n,:)), 'VariableNamingRule','preserve');
     VN=T1.Properties.VariableNames;
     DT=T1.('Date/Time (GMT)');
     T1.('Date/Time (GMT)').TimeZone='UTC';
     %T1.('Date/Time (GMT)').TimeZone;
     %Identify and generate wave parameter variables suitable to portray the wave climate of the
     %chosen sites. This identifies wave height (H), wave period (T), wave
     %direction (direction), wave power (power), and time (time).
     H=T1.("Hmax (m)");
     T=T1.("Tmax (s)");
     direction=T1.("Direction (deg)");
     time=T1.("Date/Time (GMT)");
     power=T1.("P (kW/m2)");
     
     %Identify and replace anomalous data with NaN. Corrected
     bad=9000;
     H(H>bad)=NaN;
     T(T>bad)=NaN;
     direction(direction>bad)=NaN;
     power(power>bad)=NaN;
     
     % keep values only for H<=10  Corrected
      ind2=(H <=10);
      H=H(ind2);
      T=T(ind2);
      direction=direction(ind2);
      power=power(ind2);
     
      
    %calculate wave frequency and energy
    frequency=1./(T);
    energy=1/8.*1025.*9.81.*(H.^2);
    
    %Plot wave rose diagram with wave energy and wave direction
     %(direction).
     figure(n);
%      pax1=polaraxes;
%     polarhistogram(deg2rad(direction(H<10)),deg2rad(0:10:360),'DisplayStyle', 'bar')
%     hold on
%     polarhistogram(deg2rad(direction(H<10)),deg2rad(0:10:360),'Facecolor','#0000FF','FaceAlpha',.6,'DisplayStyle','bar')
%     title(['n=' num2str(n) ' ' files(n,:)],'FontSize',12);
    
    histogram2Polar(direction, H, 0.15, 'ThetaTicks', 0:45:360, 'thetaZeroLocation', 'right', 'edgeColor', 'k','RLim',[0 ceil(max(H))], 'Normalization', 'pdf');
%     histogram2Polar(direction, H, 0.15, 'ThetaTicks', 0:45:360, 'thetaZeroLocation', 'right', 'edgeColor', 'k','RLim',[0 ceil(max(H))]);
    title(['File : ' strrep(files(n,:),'_','-')]);
 end