Group wind data in clusters based on direction

Question

Ancalagon8 2023-1-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1890502-group-wind-data-in-clusters-based-on-direction

评论： Star Strider 2023-2-23

I have an annual table with wind direction values and I want to group them per month in 16 sectors based on the degrees.

Is it possible to do it automatically either than manually?

5 个评论
显示 3更早的评论隐藏 3更早的评论

Adam Danz 2023-1-9

编辑：Adam Danz 2023-1-9

在 MATLAB Online 中打开

Grouping angular values

This demo works with degrees but it's the same process for radians. It assumes angular data are between 0 and 360. If that's not the case in your data, wrap your data to 360 (or to 2pi).

Given the number of bins nBins this computes equispaced bin edges, with a bin centered at 0 deg, and then computes the bin number for each angular data point. The bin centered at 0 is delt with separately.

Grouping isn't too difficult. However, depending on what you plan to do with these group values, you could run into trouble that @Walter Roberson explained. Phase wrapping is a potential solution to this problem (I'm looking at you bin #1).

theta = [rand(1,30)*360]; % vector of angular data [0,360] deg

nBins = 12; % select number of bins

% Compute edges

w = 360/nBins;

edges = w/2:w:(360-w/2);

% Bin data

bin = discretize(theta, edges)+1;

bin(theta>=edges(end) | theta<edges(1)) = 1; % bin at 0

% Plot results, show theta values, bin edges, and bin ID

polaraxes('ThetaZeroLocation','top',...

'RTickLabel','', ...

'ThetaDir','ClockWise')

hold on

polarplot(theta*pi/180, 1, 'ro')

polarplot([1;1].*edges*pi/180, [0;1], 'k-')

text(theta*pi/180, ones(size(theta))-.1, string(bin), ...

'HorizontalAlignment','Center', 'VerticalAlignment', 'middle')

grid off

Walter Roberson 2023-1-9

See also the discussion at https://www.mathworks.com/matlabcentral/answers/128381-how-to-calculate-mean-wind-direction

Ancalagon8 2023-1-9

@Adam Danz, @Walter Roberson thank you for your detailed and well-explained answers.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-1-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1890502-group-wind-data-in-clusters-based-on-direction#answer_1144307

编辑：Star Strider 2023-1-8

在 MATLAB Online 中打开

One way to do this would be to use the histcounts function and get the first three outputs —

WindDir = rand(100,1)*360;
edgev = linspace(0, 360, 9);
[N,Edges,Octant] = histcounts(WindDir, edgev);
WindOctants = table(WindDir,Octant)
WindOctants = 100×2 table
WindDirOctant_____________

39       3   
97       7   
43       8   
679       1   
56       5   
31       5   
28       6   
39       8   
07       4   
63       4   
69       8   
55       8   
531       2   
01       8   
321       3   
95       6   

The first output are the number of counts in each octant (bonus information), and the third is the respective octant.

EDIT — (8 Jan 2023 at 16:30)

Added table.

.

20 个评论
显示 18更早的评论隐藏 18更早的评论

Star Strider 2023-1-9

在 MATLAB Online 中打开

My pleasure!

Perhaps something like this —

WindDir = rand(100,1)*360;

edgev = linspace(0, 360, 17);

edgem = [edgev(1:end-1); edgev(2:end)].';

ctrs = mean(diff(edgev))/2 + edgev(1:end-1);

[N,Edges,Decihexctant] = histcounts(WindDir, edgev);

WindDecihexctants = table(WindDir,Decihexctant)

WindDecihexctants = 100×2 table

WindDir Decihexctant _______ ____________ 252.73 12 103.54 5 131.08 6 82.584 4 313.13 14 138.65 7 213.2 10 271.18 13 193.12 9 350.1 16 277.9 13 14.329 1 275.39 13 237.54 11 81.801 4 324.91 15

WindDirCounts = table(edgem(:,1),edgem(:,2),N(:), 'VariableNames',{'Start','End','Count'})

WindDirCounts = 16×3 table

Start End Count _____ _____ _____ 0 22.5 9 22.5 45 7 45 67.5 5 67.5 90 10 90 112.5 4 112.5 135 6 135 157.5 9 157.5 180 1 180 202.5 8 202.5 225 5 225 247.5 5 247.5 270 6 270 292.5 8 292.5 315 4 315 337.5 6 337.5 360 7

figure

bar(ctrs, N)

xlabel('Direction (°)')

ylabel('Counts')

xticks(ctrs)

axis('padded')

That only required a couple tweaks to my previous code.

.

Star Strider 2023-1-9

在 MATLAB Online 中打开

Thank you!

Possibly something like this —

DT = datetime('01-Jan-2022') + caldays(0:365).';
WindDir = rand(size(DT))*360;
TT = timetable(DT,WindDir)
TT = 366×1 timetable
DTWindDir__________________01-Jan-2022    209.94 
    02-Jan-2022    347.86 
    03-Jan-2022    61.247 
    04-Jan-2022    101.24 
    05-Jan-2022    116.63 
    06-Jan-2022     100.6 
    07-Jan-2022    210.46 
    08-Jan-2022    134.27 
    09-Jan-2022    238.92 
    10-Jan-2022    315.72 
    11-Jan-2022    27.034 
    12-Jan-2022    306.89 
    13-Jan-2022    354.17 
    14-Jan-2022     291.4 
    15-Jan-2022    2.3855 
    16-Jan-2022    233.24 
MM = month(TT.DT);
T = timetable2table(TT);
edgev = linspace(0, 360, 17);
edgem = [edgev(1:end-1); edgev(2:end)].';
ctrs = mean(diff(edgev))/2  + edgev(1:end-1);
[N,Edges,Decihexctant] = histcounts(WindDir, edgev);
WindDirCounts = table(edgem(:,1),edgem(:,2),'VariableNames',{'Start','End'});
for k = 1:12
    Month = T(MM==k,:);
    [N,Edges,Decihexctant] = histcounts(T.WindDir(MM==k,:), edgev);
    Nc{k} = N;
    Dch{k+2} = Decihexctant;
    WindDirCounts{:,k+2} = N(:);
    MMM{k} = month(Month{1,1},'shortname');
end
WindDirCounts.Properties.VariableNames = {'Start','End','Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'}
WindDirCounts = 16×14 table
    Start     End     Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct    Nov    Dec
    _____    _____    ___    ___    ___    ___    ___    ___    ___    ___    ___    ___    ___    ___

        0     22.5     1      3      1      2      0      1      3      3      1      2      4      3 
     22.5       45     1      1      1      2      2      0      3      3      6      3      3      1 
       45     67.5     3      0      3      1      2      1      2      1      2      0      4      2 
     67.5       90     0      4      1      1      3      4      3      2      2      0      0      2 
       90    112.5     3      1      1      1      1      1      2      2      2      6      4      3 
    112.5      135     4      2      2      2      2      2      1      2      1      0      0      0 
      135    157.5     0      2      1      1      5      0      4      2      2      5      2      3 
    157.5      180     1      2      4      4      1      3      0      2      2      1      1      0 
      180    202.5     0      2      2      4      3      1      0      2      0      2      1      1 
    202.5      225     4      1      3      2      2      1      1      3      2      2      0      1 
      225    247.5     3      2      2      1      2      4      2      1      3      1      1      2 
    247.5      270     1      1      1      1      2      2      2      1      3      1      3      5 
      270    292.5     4      3      3      3      2      4      3      2      1      2      2      1 
    292.5      315     2      1      3      2      1      2      3      2      1      1      2      5 
      315    337.5     2      2      1      3      1      3      1      2      2      1      0      1 
    337.5      360     3      1      2      0      2      1      1      1      0      4      3      1 

I thought about using splitapply or groupsummary for this, however it appears to be too complicated for those approaches.

My initial effort was to use the ‘MMM’ cell array to define the months for ‘WindDirCounts.Properties.VariableNames’ however that wouldn’t work and I finally gave up on that (accounting for the time delay) and just typed everything out. (That should be easier!)

.

Star Strider 2023-2-7

在 MATLAB Online 中打开

Possibly something like this —

DT = datetime('01-Jan-2022') + caldays(0:365).';

WindVel = lognrnd(log(5),log(3),size(DT));

TT = timetable(DT,WindVel)

TT = 366×1 timetable

DT WindVel ___________ _______ 01-Jan-2022 16.075 02-Jan-2022 2.1799 03-Jan-2022 45.175 04-Jan-2022 21.834 05-Jan-2022 1.3866 06-Jan-2022 55.71 07-Jan-2022 9.0376 08-Jan-2022 8.2512 09-Jan-2022 20.02 10-Jan-2022 19.658 11-Jan-2022 3.333 12-Jan-2022 34.172 13-Jan-2022 15.938 14-Jan-2022 2.7594 15-Jan-2022 2.9275 16-Jan-2022 5.6391

MM = month(TT.DT);

T = timetable2table(TT);

edgev = linspace(0, 100, 101);

edgem = [edgev(1:end-1); edgev(2:end)].';

ctrs = mean(diff(edgev))/2 + edgev(1:end-1);

[N,Edges,Decihexctant] = histcounts(WindVel, edgev);

figure

bar(ctrs, N)

grid

xlabel('Velocity')

ylabel('Counts')

WindVelCounts = table(edgem(:,1),edgem(:,2));

for k = 1:12

Month = T(MM==k,:);

[N,Edges,Decihexctant] = histcounts(T.WindVel(MM==k,:), edgev);

Nc{k} = N;

Dch{k+2} = Decihexctant;

WindVelCounts{:,k+2} = N(:);

MMM{k} = month(Month{1,1},'shortname');

end

MMMv = [MMM{:}];

WindVelCounts.Properties.VariableNames = {'Start','End', MMMv{:}}

WindVelCounts = 100×14 table

Start End Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec _____ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ 0 1 0 0 2 5 4 0 1 2 4 1 0 2 1 2 4 4 2 3 6 5 5 4 1 4 8 4 2 3 5 3 4 0 1 4 3 1 4 3 6 3 3 4 4 6 4 3 2 7 5 1 3 4 1 2 4 5 4 2 5 2 1 1 3 5 4 3 4 2 5 6 1 2 0 3 2 2 3 3 2 4 1 1 6 7 1 1 2 3 2 1 1 1 1 3 0 3 7 8 0 1 2 0 2 3 3 2 1 1 1 1 8 9 1 0 1 2 1 3 1 0 3 1 0 1 9 10 1 2 3 1 1 0 0 2 0 0 2 1 10 11 0 1 1 0 1 0 1 3 0 1 1 1 11 12 0 0 3 0 2 0 0 0 0 1 2 2 12 13 0 1 0 1 0 0 0 0 1 1 0 1 13 14 0 1 0 4 1 0 0 2 1 0 0 0 14 15 0 0 0 0 0 1 0 2 1 0 1 0 15 16 2 1 0 1 0 1 1 0 0 1 1 1

Change the second and third parameters in ‘edgeV’ to reflect the actual data, if necessary. For 1 m/s increments, the third parameter must be 1 greater than the second parameter, and both would be integers (the third must be an integer), for best results. (This uses ‘MMM’ to define the variable names in the table.)

.

Star Strider 2023-2-8

编辑：Star Strider 2023-2-8

在 MATLAB Online 中打开

As always, my pleasure!

I create the data with those lines. I do not remember any data files being posted with wind velocity (or direction), and I need data to test my code. Use the actual data instead of my synthetic data. (The only actual data I have seen thus far as been what I assume to be air temperature.)

If you post the data (with wind velocity, wind direction, temperature, and anything else necessary), I will update my code to import the files and process those data.

EDIT — (8 Feb 2023 at )

LD = load(websave('T','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1289030/T.mat'))

LD = struct with fields:

T: [513077×1 timetable]

% DT = datetime('01-Jan-2022') + caldays(0:365).';

% WindVel = lognrnd(log(5),log(3),size(DT));

TT = LD.T

TT = 513077×1 timetable

date_time Vel_avg __________________ _______ 01-Jan-19 00:00:00 1.337 01-Jan-19 00:03:00 1.609 01-Jan-19 00:04:00 1.389 01-Jan-19 00:05:00 1.466 01-Jan-19 00:06:00 2.345 01-Jan-19 00:07:00 2.074 01-Jan-19 00:08:00 1.505 01-Jan-19 00:09:00 1.253 01-Jan-19 00:10:00 1.15 01-Jan-19 00:11:00 1.958 01-Jan-19 00:12:00 2.339 01-Jan-19 00:13:00 1.641 01-Jan-19 00:14:00 1.557 01-Jan-19 00:15:00 1.053 01-Jan-19 00:16:00 1.486 01-Jan-19 00:17:00 2.436

% return

MM = month(TT.date_time);

T = timetable2table(TT);

Vmax = ceil(max(TT.Vel_avg)) % Maximum Velocity In File

Vmax = 21

edgev = linspace(0, Vmax, Vmax+1); % Edge Vector

edgem = [edgev(1:end-1); edgev(2:end)].'; % Bin Limits (For Output Table)

ctrs = mean(diff(edgev))/2 + edgev(1:end-1); % Bin Centers

[N,Edges,mps] = histcounts(TT.Vel_avg, edgev); % Call 'histcounts'

figure

bar(ctrs, N)

grid

xlabel('Velocity')

ylabel('Counts')

WindVelCounts = table(edgem(:,1),edgem(:,2)); % Initialise Table

for k = 1:12

Month = T(MM==k,:);

[N,Edges,mps] = histcounts(T.Vel_avg(MM==k,:), edgev); % Call 'histcounts¹

Nc{k} = N;

Dch{k+2} = mps;

WindVelCounts{:,k+2} = N(:);

MMM{k} = month(Month{1,1},'shortname');

end

MMMv = [MMM{:}];

WindVelCounts.Properties.VariableNames = {'Start','End', MMMv{:}}

WindVelCounts = 21×14 table

Start End Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec _____ ___ _____ ____ _____ _____ _____ ____ ____ _____ _____ _____ _____ _____ 0 1 7795 5600 7057 6124 6305 7404 7028 7613 7985 11406 6710 8092 1 2 10956 8951 14743 11734 11255 9244 8714 10383 12751 15673 13858 11413 2 3 6025 4791 5595 6726 6466 4887 6072 4637 5098 5946 7968 5106 3 4 4810 3813 4157 5829 5610 4540 5179 3997 3517 3420 5417 3886 4 5 3538 3597 3318 4266 5170 4518 4086 3874 3756 3175 3212 2486 5 6 2797 2717 2628 2841 3850 3537 3915 3844 3360 2199 1838 1675 6 7 2193 2305 2013 1998 2550 2837 3858 3938 2932 1038 1002 1527 7 8 1645 1942 1538 1558 1641 2534 3239 3269 1953 611 372 1657 8 9 1171 1355 1066 1016 666 1538 1474 2109 1027 334 169 1549 9 10 787 1040 779 617 494 652 342 529 269 205 79 1248 10 11 609 1066 664 278 350 152 65 39 33 92 27 762 11 12 359 957 506 113 145 27 11 1 7 58 7 427 12 13 194 780 249 14 38 3 5 0 0 15 7 287 13 14 42 461 98 0 2 0 0 0 0 1 2 154 14 15 0 253 31 0 0 0 0 0 0 0 0 77 15 16 0 177 7 0 0 0 0 0 0 0 0 31

figure

TL = tiledlayout(4,3);

for k = 1:12

nexttile

bar(ctrs,WindVelCounts{:,k+2})

grid

ylim([0 1.8E+4])

title(MMM{k})

end

title(TL,year(TT.date_time(1,1)))

I added the tiledlayout plots out of my own curiosity.

.

Star Strider 2023-2-23

在 MATLAB Online 中打开

I thought I could do this with one of the table functions, however extracting and transposing proved to be easiest —

LD = load(websave('T','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1289030/T.mat'))

LD = struct with fields:

T: [513077×1 timetable]

% DT = datetime('01-Jan-2022') + caldays(0:365).';

% WindVel = lognrnd(log(5),log(3),size(DT));

TT = LD.T

TT = 513077×1 timetable

date_time Vel_avg __________________ _______ 01-Jan-19 00:00:00 1.337 01-Jan-19 00:03:00 1.609 01-Jan-19 00:04:00 1.389 01-Jan-19 00:05:00 1.466 01-Jan-19 00:06:00 2.345 01-Jan-19 00:07:00 2.074 01-Jan-19 00:08:00 1.505 01-Jan-19 00:09:00 1.253 01-Jan-19 00:10:00 1.15 01-Jan-19 00:11:00 1.958 01-Jan-19 00:12:00 2.339 01-Jan-19 00:13:00 1.641 01-Jan-19 00:14:00 1.557 01-Jan-19 00:15:00 1.053 01-Jan-19 00:16:00 1.486 01-Jan-19 00:17:00 2.436

% return

MM = month(TT.date_time);

T = timetable2table(TT);

Vmax = ceil(max(TT.Vel_avg)) % Maximum Velocity In File

Vmax = 21

edgev = linspace(0, Vmax, Vmax+1); % Edge Vector

edgem = [edgev(1:end-1); edgev(2:end)].'; % Bin Limits (For Output Table)

ctrs = mean(diff(edgev))/2 + edgev(1:end-1); % Bin Centers

[N,Edges,mps] = histcounts(TT.Vel_avg, edgev); % Call 'histcounts'

figure

bar(ctrs, N)

grid

xlabel('Velocity')

ylabel('Counts')

WindVelCounts = table(edgem(:,1),edgem(:,2)); % Initialise Table

for k = 1:12

Month = T(MM==k,:);

[N,Edges,mps] = histcounts(T.Vel_avg(MM==k,:), edgev); % Call 'histcounts¹

Nc{k} = N;

Dch{k+2} = mps;

WindVelCounts{:,k+2} = N(:);

MMM{k} = month(Month{1,1},'shortname');

end

MMMv = [MMM{:}];

WindVelCounts.Properties.VariableNames = {'Start','End', MMMv{:}}

WindVelCounts = 21×14 table

Start End Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec _____ ___ _____ ____ _____ _____ _____ ____ ____ _____ _____ _____ _____ _____ 0 1 7795 5600 7057 6124 6305 7404 7028 7613 7985 11406 6710 8092 1 2 10956 8951 14743 11734 11255 9244 8714 10383 12751 15673 13858 11413 2 3 6025 4791 5595 6726 6466 4887 6072 4637 5098 5946 7968 5106 3 4 4810 3813 4157 5829 5610 4540 5179 3997 3517 3420 5417 3886 4 5 3538 3597 3318 4266 5170 4518 4086 3874 3756 3175 3212 2486 5 6 2797 2717 2628 2841 3850 3537 3915 3844 3360 2199 1838 1675 6 7 2193 2305 2013 1998 2550 2837 3858 3938 2932 1038 1002 1527 7 8 1645 1942 1538 1558 1641 2534 3239 3269 1953 611 372 1657 8 9 1171 1355 1066 1016 666 1538 1474 2109 1027 334 169 1549 9 10 787 1040 779 617 494 652 342 529 269 205 79 1248 10 11 609 1066 664 278 350 152 65 39 33 92 27 762 11 12 359 957 506 113 145 27 11 1 7 58 7 427 12 13 194 780 249 14 38 3 5 0 0 15 7 287 13 14 42 461 98 0 2 0 0 0 0 1 2 154 14 15 0 253 31 0 0 0 0 0 0 0 0 77 15 16 0 177 7 0 0 0 0 0 0 0 0 31

VN = WindVelCounts.Properties.VariableNames;

WindVelCountsTransposed = array2table(table2array(WindVelCounts(:,3:end)).', 'RowNames',VN(3:end));

WindVelCountsTransposed.Properties.VariableNames = compose('[%d-%d]',edgem)

WindVelCountsTransposed = 12×21 table

[0-1] [1-2] [2-3] [3-4] [4-5] [5-6] [6-7] [7-8] [8-9] [9-10] [10-11] [11-12] [12-13] [13-14] [14-15] [15-16] [16-17] [17-18] [18-19] [19-20] [20-21] _____ _____ _____ _____ _____ _____ _____ _____ _____ ______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ _______ Jan 7795 10956 6025 4810 3538 2797 2193 1645 1171 787 609 359 194 42 0 0 0 0 0 0 0 Feb 5600 8951 4791 3813 3597 2717 2305 1942 1355 1040 1066 957 780 461 253 177 108 68 33 12 1 Mar 7057 14743 5595 4157 3318 2628 2013 1538 1066 779 664 506 249 98 31 7 0 0 0 0 0 Apr 6124 11734 6726 5829 4266 2841 1998 1558 1016 617 278 113 14 0 0 0 0 0 0 0 0 May 6305 11255 6466 5610 5170 3850 2550 1641 666 494 350 145 38 2 0 0 0 0 0 0 0 Jun 7404 9244 4887 4540 4518 3537 2837 2534 1538 652 152 27 3 0 0 0 0 0 0 0 0 Jul 7028 8714 6072 5179 4086 3915 3858 3239 1474 342 65 11 5 0 0 0 0 0 0 0 0 Aug 7613 10383 4637 3997 3874 3844 3938 3269 2109 529 39 1 0 0 0 0 0 0 0 0 0 Sep 7985 12751 5098 3517 3756 3360 2932 1953 1027 269 33 7 0 0 0 0 0 0 0 0 0 Oct 11406 15673 5946 3420 3175 2199 1038 611 334 205 92 58 15 1 0 0 0 0 0 0 0 Nov 6710 13858 7968 5417 3212 1838 1002 372 169 79 27 7 7 2 0 0 1 1 1 0 0 Dec 8092 11413 5106 3886 2486 1675 1527 1657 1549 1248 762 427 287 154 77 31 6 3 3 1 0

% figure

% TL = tiledlayout(4,3);

% for k = 1:12

% nexttile

% bar(ctrs,WindVelCounts{:,k+2})

% grid

% ylim([0 1.8E+4])

% title(MMM{k})

% end

% title(TL,year(TT.date_time(1,1)))

All the original information is still there in case you need it. Plotting it using bar3 might be another option, in addition to tiledlayout.

.