One option would be to ‘decompose’ the angles into their ‘x’ and ‘y’ components, calculate the means of those, and then use atan2 to return the result —
D = [2020 4 1 0 77.0016556593828
2020 4 1 3 61.6846224380681
2020 4 1 6 31.4624383596130
2020 4 1 9 35.1994739111645
2020 4 1 12 56.7123189304487
2020 4 1 15 337.212660715222
2020 4 1 18 332.720394142200
2020 4 1 21 359.309990882483
2020 4 2 0 1.81332665386765
2020 4 2 3 316.084846935155
2020 4 2 6 297.362308097189
2020 4 2 9 286.447414863766
2020 4 2 12 253.377128455776
2020 4 2 15 272.755841963236
2020 4 2 18 283.300130880562
2020 4 2 21 305.805458428448
2020 4 3 0 303.435585960465
2020 4 3 3 305.408085546923
2020 4 3 6 314.945219141016
2020 4 3 9 301.363197206495
2020 4 3 12 297.344102342364
2020 4 3 15 294.845524112617
2020 4 3 18 301.544742721005
2020 4 3 21 320.320593041074];
dir = D(:,5);
meandir = @(a) atan2d(mean(sind(a)),mean(cosd(a)));
dir = accumarray(ceil((1:numel(dir))/8)',dir(:),[],meandir)
dir = wrapTo360(dir) % Mapping Toolbox Function
.
