Plotting Weibull distribution on Wind Speed

Question

I have a 720 hourly set of wind speed and wind direction data and I want to fit the Weibull Distribution on it. After searching for some time, I wrote the following code in Python to get my distribution, I will also share my image for clarification.

ws.hist(bins=np.arange(0, ws.max()), alpha=0.5, normed=True)
weibull_params = (sc.stats.exponweib.fit(ws,floc=0,f0=1))
x= ws
def weib(x,lamb,k):
   return (k / lamb) * (x / lamb)**(k-1) * np.exp(-(x/lamb)**k)
k_shape, lamb_scale = weibull_params[1], weibull_params[3]

plt.plot(x, weib(x, lamb_scale, k_shape), label='self-defined Weibull')

Does this seem alright? How do I make sense of this? I also found a package here but I don't understand how to implement it to my code. https://github.com/cqcn1991/Wind-Speed-Analysis

Thank you.

Answer 1

Tried a few things in the OP comments that didn't work. However, you may apply Wolfram Language to your project. There is a free Wolfram Engine for developers and with the Wolfram Client Library for Python you can use these functions in Python.

import datetime

from wolframclient.evaluation import WolframLanguageSession
from wolframclient.language import wl, wlexpr

Start the Wolfram session

wolfSession = WolframLanguageSession()

I will collect my own WindSpeedData from a WeatherData weather station closest to the "Country" Entity of Bermuda during October 2016.

weatherStation = wolfSession.evaluate(
    wl.WeatherData(wl.Entity('Country','Bermuda'))
);
print(weatherStation)

Entity['WeatherStation', 'TXKF']

windData = wolfSession.evaluate(
    wl.WindSpeedData(
        weatherStation,
        [
            datetime.datetime(2016,10,1),
            datetime.datetime(2016,10,31)
        ]
    )
)

windData is a TimeSeries object with observation count comparable to your problem.

print(wolfSession.evaluate(windData('PathLength')))

1058

The result can be visualised with DateListPlot by Exporting in one of the supported Raster Image Formats or Vector Graphics Formats.

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        wl.DateListPlot(windData, 
            PlotTheme='Detailed', 
            PlotRange=wl.All, 
            FrameLabel=wl.Automatic
        )
    )
)

A hurricane passed over in that month with is the cause of that spike in wind speed.

Now that there is data Histogram on the PDF scale gives.

hist=wolfSession.evaluate(
    wl.Histogram(windData, wl.Automatic, 'PDF',
        PlotTheme='Detailed', 
        ChartStyle=wl.ColorData('Crayola','Silver'),
        PlotRange=wl.All
    )
);

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        hist
    )
)

Use EstimatedDistribution for the WeibullDistribution.

windDistrbution=wolfSession.evaluate(
    wl.EstimatedDistribution(
        wl.QuantityMagnitude(windData('Values')),
        wl.WeibullDistribution(wl.Global.alpha, wl.Global.beta, wl.Global.mu)
    )
);
print(windDistrbution) 

WeibullDistribution[1.883495945177254, 28.34295076324276, -0.7675654467340361]

Or use FindDistribution restricted to the continuous distribution functions to auto find and fit a distribution.

autoDistribution=wolfSession.evaluate(
    wl.FindDistribution(
        wl.QuantityMagnitude(windData('Values')),
        TargetFunctions='Continuous'
    )
);

print(autoDistribution)

ExtremeValueDistribution[18.436141153779957, 10.204118744677338]

Check for goodness of with DistributionFitTest.

print(
    wolfSession.evaluate(
        wl.Map(
            wl.Function(wl.DistributionFitTest(windData,wl.Slot(1),'PValue')), 
            [windDistrbution,autoDistribution]
        )
    )
)

[0.0002862317707539308, 0.1802410209978217]

The Weibull fit is not very good and Extreme Value not that strong. In any case the Plot the PDF functions.

pdfs=wolfSession.evaluate(
    wl.Map(
        wl.Function(wl.PDF(wl.Slot(1), wl.Global.x)), 
        [windDistrbution,autoDistribution]
    )
);


pdfPlot=wolfSession.evaluate(
    wl.Plot(
        pdfs,
        [wl.Global.x, 0, wl.QuantityMagnitude(wl.Max(windData))],
        PlotTheme='Detailed',
        PlotLegends=wl.Placed(
            wl.LineLegend(wl.Automatic, [windDistrbution,autoDistribution]), 
            wl.Below)
    )
)

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        pdfPlot
    )
)

Combine the plots with Show.

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        wl.Show(hist, pdfPlot)
    )
)

Terminate the session

wolfSession.terminate()

Hope this helps.