2

I am currently trying to plot in matlab a wind rose diagram with data wind velocities and directions for a given period. The main program is such that after plotting several plots on the Weibull distribution, it calls another matlab program to produce a wind rose.

The wind rose program is essentially here: https://www.mathworks.com/matlabcentral/fileexchange/47248-wind-rose

and the main program is mainly based on https://www.mathworks.com/matlabcentral/fileexchange/41996-computing-weibull-distribution-parameters-from-a-wind-speed-time-series?focused=3786165&tab=function

I've been working yesterday on a very old matlab edition and I had serious problems with making the code run. Today with Octave in a Ubuntu machine and after some efforts I managed to get a result with a minor problem: the wind rose did not had all the information it should have.

I run the program in a new version of matlab and then I got the following message:

Error using WindRose (line 244)
 is not a valid property for WindRose function.
Error in octavetestoforiginalprogram (line 184)
[figure_handle,count,speeds,directions,Table] = WindRose(dir,vel,Options); 

How can the program run in Octave and now produces such an error, I don't understand.

Does anyone have an idea of what this error means?

Note: I am posting the entire code below if anyone wants to read it:

%% EXTRACT AND PLOT RAW DATA

% Extract wind speed data from a file
v = xlsread('1981-1985_timeseries.xlsx');

% Plot the measured wind speed
plot(v)
title('Wind speed time series');
xlabel('Measurement #');
ylabel('Wind speed [m/s]');





%% PROCESS DATA

% Remove nil speed data (to avoid infeasible solutions in the following)
v(find(v==0)) = [];
% Extract the unique values occuring in the series
uniqueVals = unique(v);
uniqueVals(isnan(uniqueVals))=[];

% Get the number of unique values
nbUniqueVals = length(uniqueVals);

% Find the number of occurences of each unique wind speed value
for i=1:nbUniqueVals
    nbOcc = v(find(v==uniqueVals(i)));
    N(i) = length(nbOcc);
end

% Get the total number of measurements
nbMeas = sum(N);

% To take into account the measurement resolution
% (i.e., a measured wind speed of 2.0 m/s may actually correspond to a
% real wind speed of 2.05 or 1.98 m/s), compute the delta vector which
% contains the difference between two consecutive unique values
delta(1) = uniqueVals(1);
for i=2:nbUniqueVals
    delta(i) = uniqueVals(i) - uniqueVals(i-1);
end

% Get the frequency of occurence of each unique value
for i=1:nbUniqueVals
    prob(i) = N(i)/(nbMeas*delta(i));
end



% Get the cumulated frequency
freq = 0;
for i=1:nbUniqueVals
    freq = prob(i)*delta(i) + freq;
    cumFreq(i) = freq;
end


%% PLOT THE RESULTING DISTRIBUTION

% Plot the distribution
figure
subplot(2,1,1);
pp=plot(uniqueVals,prob)
title('Distribution extracted from the time series');
xlabel('Wind speed [m/s]');
ylabel('Probability');

% Plot the cumulative distribution
subplot(2,1,2);
plot(uniqueVals,cumFreq)
title('Cumulative distribution extracted from the time series');
xlabel('Wind speed [m/s]');
ylabel('Cumulative probability');


%% EXTRACT THE PARAMETERS USING A GRAPHICAL METHOD

% See the following references for more explanations:
% - Akdag, S.A. and Dinler, A., A new method to estimate Weibull parameters
% for wind energy applications, Energy Conversion and Management,
% 50 :7 1761�1766, 2009
% - Seguro, J.V. and Lambert, T.W., Modern estimation of the parameters of
% the Weibull wind speed distribution for wind energy analysis, Journal of
% Wind Engineering and Industrial Aerodynamics, 85 :1 75�84, 2000


% Linearize distributions (see papers)
ln = log(uniqueVals);
lnln = log(-log(1-cumFreq));

% Check wether the vectors contain inifinite values, if so, remove them
test = isinf(lnln);
for i=1:nbUniqueVals
    if (test(i)==1)
        ln(i)= [];
        lnln(i)= [];
    end
end

% Extract the line parameters (y=ax+b) using the polyfit function
params = polyfit(ln,lnln',1);
a = params(1);
b = params(2);
y=a*ln+b;

% Compare the linealized curve and its fitted line
figure
plot(ln,y,'b',ln,lnln,'r')
title('Linearized curve and fitted line comparison');
xlabel('x = ln(v)');
ylabel('y = ln(-ln(1-cumFreq(v)))');

% Extract the Weibull parameters c and k
k = a
c = exp(-b/a)


%% CHECK RESULTS

% Define the cumulative Weibull probability density function
% F(V) = 1-exp(-((v/c)^k)) = 1-exp(-a2), with a1 = v/c, a2 = (v/c)^k
a1 = uniqueVals/c;
a2 = a1.^k;
cumDensityFunc = 1-exp(-a2); 

% Define the Weibull probability density function
%f(v)=k/c*(v/c)^(k-1)*exp(-((v/c)^k))=k2*a3.*exp(-a2), 
% with  k2 = k/c, a3 = (v/c)^(k-1)
k1 = k-1;
a3 = a1.^k1;
k2 = k/c;
densityFunc = k2*a3.*exp(-a2);  

% Plot and compare the obtained Weibull distribution with the frequency plot
figure
subplot(2,2,1);
pp=plot(uniqueVals,prob,'.',uniqueVals,densityFunc, 'r')
title('Weibull probability density function');
xlabel('v');
ylabel('f(v)');

subplot(2,2,3)
h=hist(v);
title('Wind speed time series');
xlabel('Measurement #');
ylabel('Wind speed [m/s]');
h=h/(sum(h)*10);
bar(h)

% Same for the cumulative distribution
subplot(2,2,2);
plot(uniqueVals,cumFreq,'.',uniqueVals,cumDensityFunc, 'r')
title('Cumulative Weibull probability density function');
xlabel('v');
ylabel('F(V)');


%inner
figure
hold on
pp=plot(uniqueVals,prob,'.',uniqueVals,densityFunc, 'r')
title('Weibull probability density function');
xlabel('v');
ylabel('f(v)');

bar(h)

hold off

%inner



%rose

w=xlsread('rose.xlsx');
dir=w(:,2)*10;
vel=w(:,1);

Options = {'anglenorth','FreqLabelAngle',0,'angleeast','FreqLabelAngle',90,'labels',{'N (0)','S (180)','E (90)','W (270)'},'freqlabelangle',45,'nDirections',20,'nFreq',25,'LegendType',1};
[figure_handle,count,speeds,directions,Table] = WindRose(dir,vel,Options);

close all; clear Options;
B--rian
  • 5,578
  • 10
  • 38
  • 89

1 Answers1

2

After a quick read of the script documentation, here is what I found concerning the creation of the windrose plot:

% With options in a cell array:
Options = {'anglenorth',0,'angleeast',90,'labels',{'N (0°)','S (180°)','E (90°)','W (270°)'},'freqlabelangle',45};
[figure_handle,count,speeds,directions,Table] = WindRose(dir,spd,Options);

% With options in a structure:
Options.AngleNorth     = 0;
Options.AngleEast      = 90;
Options.Labels         = {'N (0°)','S (180°)','E (90°)','W (270°)'};
Options.FreqLabelAngle = 45;
[figure_handle,count,speeds,directions,Table] = WindRose(dir,spd,Options);
close all;

% Usual calling:
[figure_handle,count,speeds,directions,Table] = WindRose(dir,spd,'anglenorth',0,'angleeast',90,'labels',{'N (0°)','S (180°)','E (90°)','W (270°)'},'freqlabelangle',45);

Your error is:

Error using WindRose (line 244)
 is not a valid property for WindRose function.

Error in octavetestoforiginalprogram (line 184)
[figure_handle,count,speeds,directions,Table] = WindRose(dir,vel,Options);

And it is being produced within the routine that undertakes option arguments sanitization. Since options must be provided in the form of name-value pairs... it seems that the script is detecting a mismatching number of elements and one or more names with a missing value. And here they are:

Options = {'anglenorth','FreqLabelAngle',0,'angleeast','FreqLabelAngle',90,'labels',{'N (0)','S (180)','E (90)','W (270)'},'freqlabelangle',45,'nDirections',20,'nFreq',25,'LegendType',1};

The two properties marked with a bold font have no value associated to them (unlike the examples in the tutorial) and this probably messes up the whole parametrization process. Probably, the first option being extracted is anglenorth = FreqLabelAngle, which is not correct.

Tommaso Belluzzo
  • 23,232
  • 8
  • 74
  • 98
  • 2
    Well, if that's the problem( and somehow Octave got around it) I can't thank you enough. Unfortunately I don't have access to matlab right now, to verify your claim, but it seems reasonable to believe you are correct and that this is in fact the root of the error. I will check the code in about 9-10 hours, and I will get back at you. But again, thank you for all your effort here. – Constantine Black Feb 26 '18 at 22:35
  • 1
    It runs. Now I just have to figure how to produce the desired result. Again, thank you very much; you consideration has been great. – Constantine Black Feb 27 '18 at 08:34