I needed to generate a wind speed time series for a simulation, and couldn’t find any pre-existing code in Python to achieve this.

However, I did find the following post on using a Markov autoregressive function to achieve the same in R:

http://procomun.wordpress.com/2010/11/02/how-to-simulate-wind-speed-time-series-with-r/

So I’ve copied the mathematics used there and ported it across to Python, but with a bit of OOP rearrangement to give an object which will compute values on demand. In retrospect, a full-on generator class might be more pythonic, but I’ll leave that for now.

First some imports and helper functions:

 1234567891011121314151617 #!/usr/bin/env python import datetime import math import numpy import random import operator from collections import Counter #homegrown weibull probability function #scipy.stats.exonweib.pdf gives different result from R dweibull function def weibullpdf(data,scale,shape):     return [(shape/scale)*((x/scale)**(shape-1))*math.exp(-1*(x/scale)**shape) for x in data]     #matrix-vector multiplier from http://code.activestate.com/recipes/121574-matrix-vector-multiplication/ #equivalent to %*% in R def matmult4(m, v):     return [reduce(operator.add,map(operator.mul,r,v)) for r in m]

And then a WindResource object which will return values at each step:

 18192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899 class WindResource(object):     def __init__(self,mean_speed=9.0,max_speed=30.0,n_states=30,start_time=datetime.datetime(2012,1,1),hist_length=24):         self.mean_speed=mean_speed         self.max_speed=max_speed         self.n_states=n_states         self.time=start_time         self.hist_length=hist_length                 #setup matrix         n_rows=n_states                                     n_columns=n_states                   self.l_categ=float(max_speed)/float(n_states)    #position of each state                 #weibull parameters         weib_shape=2.0         weib_scale=2.0*float(mean_speed)/math.sqrt(math.pi);                 #Vector of wind speeds         self.WindSpeed = numpy.arange(self.l_categ/2.0,float(max_speed)+self.l_categ/2.0,self.l_categ)                 #distribution of probabilities, normalised         fdpWind = weibullpdf(self.WindSpeed,weib_scale,weib_shape)         fdpWind = fdpWind/sum(fdpWind)                 #decreasing function         G = numpy.empty((n_rows,n_columns,))         for x in range(n_rows):             for y in range(n_columns):                 G[x][y]=2.0**float(-abs(x-y))                     #Initial value of the P matrix         P0 = numpy.diag(fdpWind)                 #Initital value of the p vector         p0 = fdpWind                 #below comment from R source code:         #"The iterative procedure should stop when reaching a predefined error.         #However, for simplicity I have only constructed a for loop. To be improved!"         P,p=P0,p0         for i in range(10):             r=matmult4(P,matmult4(G,p))             r=r/sum(r)             p=p+0.5*(p0-r)             P=numpy.diag(p)                     N=numpy.diag([1.0/i for i in matmult4(G,p)])         MTM=matmult4(N,matmult4(G,P))         self.MTMcum = numpy.cumsum(MTM,1)                 #initialise series         self.state=0         self.states_series=[]         self.speed_series=[]         self.power_series=[]         self.randoms1=[]         self.randoms2=[]                 #tick over to first value (decrement time accordingly)         self.time=self.time+datetime.timedelta(hours=-1)         self.getNext()             #show current value without incrementing         def getCurrent(self):         wind_counter = Counter([int(round(x,0)) for x in self.speed_series])         return {'time': self.time,                 'data':{                     'wind_speed': self.speed_series[-1],                     'wind_speed_av' : sum(self.speed_series)/float(len(self.speed_series)),                     'wind_hist': dict(wind_counter),                     }                 }             #increment time by one hour and return new value     def getNext(self):         self.randoms1.append(random.uniform(0,1))         self.randoms2.append(random.uniform(0,1))         self.state=next(j for j,v in enumerate(self.MTMcum[self.state]) if v > self.randoms1[-1])         self.states_series.append(self.state)         self.speed_series.append(self.WindSpeed[self.state]-0.5+(self.randoms2[-1]*int(self.l_categ)))         self.time=self.time+datetime.timedelta(hours=1)         return self.getCurrent()

So, to use, create the resource object, which will initialise and create the first values. getCurrent() gives the current time and speed, getNext() increments the time by an hour and returns the next set of values. Each returned dataset also includes binned data by rounded wind speed value for a histogram, as well as the average speed to date.

>>> myWind = WindResource()
>>> myWind.getCurrent()
{'data': {'wind_speed': 2.4038549908024818, 'wind_hist': {2: 1}, 'wind_speed_av': 2.4038549908024818}, 'time': datetime.datetime(2012, 1, 1, 0, 0)}
>>> myWind.getNext()
{'data': {'wind_speed': 2.1456530189886505, 'wind_hist': {2: 2}, 'wind_speed_av': 2.2747540048955663}, 'time': datetime.datetime(2012, 1, 1, 1, 0)}

The distribution is a Weibull distribution with the given average, and a shape parameter of 2 (i.e. a Rayleigh distribution).

A run of the default settings for a year gives the following: Issues/improvements:

• No seasonal component
• No diurnal component
• Values can be seen to occasionally return to zero at an unrealistic delta
• Shape parameter fixed at 2
• Randoms not seeded – unreproducable
• Hourly values only – use aggregation for longer periods, interpolation for shorter?

1. Hasan says:

Great Job. I was looking for something like this. Really nice

2. Excellent post. I definitely love this website. Keep it
up!

3. Hi there, of course this post is really good and I have learned lot of things from it concerning blogging.
thanks.

4. Good idea! :evil:

5. I am really searching for a create a breeze speed time arrangement for a reenactment, and couldn’t discover any prior code in Python to accomplish this.So thanks for sharing this programing code here.

6. virginia says:

good

7. This is very serious topic. It was nice to read the article.

8. Looking at the codes makes my head spin. I really admire people who can easily decipher the meaning of that. Actually I took IT as my major course before but I really couldn’t handle the things that they are doing, it led me to shifting to another course which is engineering. My brother had a really hard time getting his degree on IT and I didn’t want to share the same fate. This is a very wonderful post that you have shared to us.

9. Alice says:

Was searching for it earlier and couldn’t find it eather. Thanks for sharing! Cheers

10. interesting! where did you get this idea?