#!/usr/bin/python
# -*- coding: utf-8 -*-
#windrose - Computes wind frequencies and plots windroses
#Lionel Roubeyrie - 10/2006
#Licence : this code is released under the matplotlib license
from matplotlib.patches import Rectangle
from matplotlib.axes import PolarSubplot
import matplotlib.cm as cm
from scipy.stats import stats
from math import ceil
from pylab import *
from numpy import *
VERSION="0.3"
def __version__():
return VERSION
def windplot( speed, direction, ax=None, counts=True, speed_classes=[0,1,2,3,4,5], sectors = 18, style='bar', colormap=cm.jet ):
'''Plots a windrose. A count (or percent) table is computed, for each speed_classe
and wind sector.
Returns the count or percent table, and resulting axes.
speed: array of wind speeds
direction: array of wind directions, same order as speed array
ax: a valid polar axes instance. If None, one is created
counts: If True (default), the windrose (and the table) will be based on the
numbers of counts each sectors and speed_classes are reached, else all is in percent.
speed_classes: bins of wind speeds to compute
sectors: normally 18 or 36, the number of wind sectors (directions)
style: bar | fill
colormap: a matplotlib.cm colormap'''
if ax is None:
ax = axes( polar=True )
#Get the wind frequency table
freq = windfreq( speed, direction, counts, speed_classes, sectors )
nb_sc = len( speed_classes )
#Building the colormap according to the number of speed classes
cmapindex = linspace(0.0, 1.0, nb_sc )
cmap = {}
for i in range( nb_sc ):
cmap[i] = colormap( cmapindex[i] )
#The angle of each sector
dtheta = 2*pi/sectors
#Building the list of angles
angles = arange( 0, -2*pi, -dtheta )
angles = angles + pi/2
patches_list = []
if style == 'line':
#To have closed lines
angles = resize( angles, (angles.shape[0]+1,) )
angles[-1] = angles[0]
vals = freq.tolist()
for i in vals:
i.append(i[0])
vals = array( vals )
offset = 0
for row in range( nb_sc ):
val = vals[row,:] + offset
offset += vals[row,:]
patch = polar( angles, val, color=cmap[row], zorder=nb_sc-row )
patches_list.append( patch )
if style == 'fill':
offset = 0
for row in range( nb_sc ):
val = freq[row,:] + offset
offset += freq[row,:]
patch = fill( angles, val, facecolor=cmap[row], edgecolor=cmap[row], zorder=nb_sc-row )
patches_list.append( patch )
if style == 'bar':
for col in range( sectors ):
offset = 0
for row in range( nb_sc ):
if row > 0:
offset += freq[row-1, col]
patch = Rectangle( (angles[col]-dtheta/2, offset), dtheta, freq[row, col], facecolor=cmap[row] )
ax.add_patch( patch )
if col == 0:
patches_list.append( patch )
if style == 'bar2':
opening = linspace( 0.0, pi/16, nb_sc )
for col in range( sectors ):
offset = 0
for row in range( nb_sc ):
if row > 0:
offset += freq[row-1, col]
patch = Rectangle( (angles[col]-opening[row]/2, offset), opening[row], freq[row, col], facecolor=cmap[row] )
ax.add_patch( patch )
if col == 0:
patches_list.append( patch )
# construct the legend
speed_classes.append(inf)
speed_classes = map( str, speed_classes )
leg = [ "[%s : %s["%(speed_classes[i], speed_classes[i+1]) for i in range( len(speed_classes)-1 ) ]
figlegend(patches_list, leg, 'lower left' )
ax.autoscale_view()
# Set labels
radii = linspace( 0, ax.get_rmax(), 6 )
#radii = linspace( 0, ceil( freq.max() ), 6 )
if counts == True:
labels = [ "%.1f" %r for r in radii ]
else:
labels = [ "%.1f%%" %r for r in radii ]
ax.set_rgrids( radii, labels, angle=90 )
# Set directions
ax.set_thetagrids( angles=arange(0,360,45), labels=['E','N-E','N','N-W','W','S-W','S','S-E'], frac=1.2 )
return freq, ax
def windfreq( speed, direction, counts=True, speed_classes=[0,1,2,3,4,5], sectors = 18 ):
''' Returns a wind frequency table where for each sector of wind direction,
we have the count of time wind comes with a particular speed.
speed: array of wind speeds
direction: array of wind directions (same order as speed). Must be in degrees
from north
counts: if True, results are the numbers of events for each speed_classes and
each sectors. If False, results are in percentage
speed_classes: list of wind speed bins against we're going to calculate the
frequency table
sectors : number of sectors (18 or 36 generally). The result will be centred
on the north.
THe result is an array of shape=( sectors, len(speed_classes) ) which
contains the number of time each speed_classe is reached by direction bins'''
if len( speed ) != len( direction ):
raise ValueError, "speed array and direction array must have same length"
freq = zeros( (len(speed_classes), sectors+1) )
angle = 360./sectors
wind_classes = linspace( 0, 360, sectors+1 )
wind_classes[1:] = wind_classes[1:] - angle/2
#We start by selecting values by direction classes, i.e. all values where
#the direction is between wind_classes(n) and wind_classes(n+1)
last = []
for i in range( len( wind_classes ) ):
if i == len( wind_classes ) - 1:
values = select( [greater_equal( direction, wind_classes[i] )],[speed], default=-1.e20 )
else:
values = select( [greater_equal( direction, wind_classes[i] ) * less( direction, wind_classes[i+1] )],[speed], default=-1.e20 )
# Now, for each direction class, we compute the number of time each speed class is reached
res = stats.histogram2( values, speed_classes )
freq[:,i] = res
# add the last value to the first to have the freq of North winds
freq[:,0] = freq[:,0] + freq[:,-1]
# and remove the last col
freq = freq[:,:-1]
# To be directly in the format adapted to windplot
if counts == False:
freq = freq*100/freq.sum()
return freq