jmquintana79
9/5/2017 - 8:13 AM

## Turbulence Analysis Plot

Turbulence Intensity plot where furthermore it is returned the turbulence matrix (10*-3).

``````"""
TURBULENCE INTENSITY ANALYSIS:

Plot the wind rose and turbulence per wind speed bins. It is required a DF with hwd, hws and stdhws data.

Furthermore return the turbulence intensity matrix (*10***-3)
"""
def analysis_turbulence(DATA,svar_hws,svar_stdhws,svar_hwd,stitle=''):

## CALCULATE TURBULENCE
DATA['turbulence'] = DATA[[svar_hws,svar_stdhws]].apply(lambda x: x[1]//x[0], axis=1)

## CALCULATE SECTORS
def setsector(vdir):
for k, v in dsector.items():
if k==0:
if (vdir>=v[1]) or (vdir<v[2]): return v[0]
else:
if (vdir>=v[1]) and (vdir<v[2]): return v[0]

## SECTOR DEFINITION (16 sectors)

# initialize
dsector = dict(); disector = dict()
lsectors = ['N','NNE','NE','ENE','E','ESE','SE','SSE','S','SSW','SW','WSW','W','WNW','NW','NNW']
llimits =[348.75,11.25],[11.25,33.75],[33.75,56.25],[56.25,78.75],[78.75,101.25],[101.25,123.75],[123.75,146.25],[146.25,168.75],[168.75,191.25],[191.25,213.75],[213.75,236.25],[236.25,258.75],[258.75,281.25],[281.25,303.75],[303.75,326.25],[326.25,348.75]
# loop of angles
for i,isector in enumerate(lsectors[:]):
# build dictionary
dsector[i] = [isector,llimits[i][0],llimits[i][1]]
disector[isector] = i

# collect sector
lsector = [v[0] for k,v in dsector.items()]
# calculate sector
DATA['sector'] = DATA[svar_hwd].apply(lambda col: setsector(col))

## BIN DEFINITION

# set bin
def setbin(vspeed):
for k, v in dbin.items():
if (vspeed>=v[1]) and (vspeed<v[2]): return int(k)

ibin = 0
dbin = dict()
for il in range(16):
# set limits
il1 = il
if il==15: il2 = 9999
else: il2 = il+1
# define bin and store
dbin[ibin]=['%s<=hws<%s'%(il1,il2),il1,il2]
# define number of bin
ibin += 1

# calculate bin
DATA['bin'] = DATA[svar_hws].apply(lambda col: setbin(col))

""" CREATE TURBULENTECE INTENSITY MATRIX """

# initialize
lbin = list(range(16))
ntotal = len(DATA)
davg = dict()
# loop of sectors
for isector in lsector:
# initialize
lavg = list()
# loop of bins
for ibin in lbin:
# collect data
AUX = DATA[(DATA['sector']==isector) & (DATA['bin']==ibin)]
# calculate avg of hws
iavg = np.nanmean(AUX['turbulence'].values)*1000.
# store
lavg.append(iavg)
# clean
del(AUX)
# store by sector
davg[isector] = lavg

## FINAL DF CREATION

# create frequency matrix
WSTURB = pd.DataFrame.from_dict(davg, orient='columns')
WSTURB['index'] = [v[0] for k,v in dbin.items()]
WSTURB = WSTURB.set_index('index')
WSTURB['TOTAL'] = WSTURB.sum(axis=1)

# create sum freq
WSTURBSUM = pd.DataFrame(WSTURB.sum(axis=0)).transpose()
WSTURBSUM['index'] = ['total']
WSTURBSUM = WSTURBSUM.set_index('index')
# append total
WSTURB = WSTURB.append(WSTURBSUM)
# clean
del(WSTURBSUM)

""" PLOT TURBULENCE ROSE """

import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib.colors import ListedColormap
import matplotlib as mpl

## COLLECT DATA
PLOT = WSTURB.transpose()[['total']][:-1]
PLOT['isector'] = [disector[i] for i in PLOT.index.values]
PLOT.sort_values(['isector'], ascending=[1], inplace=True)

## DATA PREPARATION

# collect total data
values = list(PLOT.total.values); values = values[:-1]
#freq = list(PLOT.freq.values)
lsector = list(PLOT.index.values)
# clean
del(PLOT)
# calculate y lim
ylim = int(np.max(values))+2

# define angles of sectors for labels of rose chart
langle = list(np.arange(0,360,22.5))

# define x
numpoints = len(lsector)
theta = [n / float(numpoints) * 2 * pi for n in range(numpoints)]

# Because our chart will be circular we need to append a copy of the first
# value of each list at the end of each list with data
values += values[:1]
theta = list(theta)
thetaf = copy.deepcopy(theta)
theta += theta[:1]

## create figure object
params = dict(projection='polar', theta_direction=-1, theta_offset=np.pi/2)
fig, ax = plt.subplots(subplot_kw=params,figsize=(6,7))

## PLOT CHART OF FREQUENCIES

# plot
width = 0.25
ax.bar(thetaf,values,width, color='blue', alpha=0.8)
# define limits
ax.set_ylim([0, ylim])
# set axi tick labels
plt.setp(ax.get_yticklabels(), color='blue')
plt.yticks(range(0,ylim+4,4),fontsize=10,rotation='horizontal')
ax.set_yticklabels(list(range(0,ylim+4,4)))

# set axis labels
ax.text(radians(22.5+10.), ylim + 5, '(*10**-3)', size=10,color='blue', horizontalalignment="center", verticalalignment="center")

# set ytick labels of rose chart
ax.set_thetagrids(langle,labels=lsector,fontsize=16, weight="normal", color="black")

# title
if len(stitle)>0: title = ': %s'%stitle
else: title = stitle
ax.set_title('Turbulence Intensity Rose%s'%title,fontsize=15,y=1.1)

# display
plt.show()

""" PLOT BAR FREQ """

import matplotlib.pyplot as plt

## COLLECT DATA
BAR = WSTURB[['TOTAL']].iloc[0:16]

## PREPARE DATA
nrow = len(BAR)
x = np.array(list(range(nrow)))+0.5
y = BAR.TOTAL.values
ylim = int(np.max(y))+2
# clean
del(BAR)

## PLOT

# create objects
fig, ax = plt.subplots(figsize=(6,5))

# plot bars
width = 1
ax.bar(x,y,width, color='r', edgecolor='black')

# limits of axes
ax.set_xlim([0,16])
ax.set_ylim([0,ylim+5])

# resolution of axis ticks
plt.xticks(list(range(nrow)),fontsize='small',rotation='horizontal')
ax.set_xticklabels(list(range(nrow)))
plt.yticks(range(0,ylim+4,4),fontsize='small',rotation='horizontal')
ax.set_yticklabels(list(range(0,ylim+4,4)))

# set axis labels
plt.xlabel("WS(m/s)",fontsize=10)
ax.set_ylabel("Turbulence (*10**-3)",fontsize=10)

# title
if len(stitle)>0: title = ': %s'%stitle
else: title = stitle
plt.title('Turbulence per WS%s'%title,fontsize=16)

# set labels on points for scatter chart
lvalue = ['%.1f'%iy for iy in y]
for i, svalue in enumerate(lvalue):
if float(svalue)>=10.: ax.annotate(svalue, (x[i]-0.48,y[i]+0.25),fontsize=9,color='red',rotation='horizontal')
else: ax.annotate(svalue, (x[i]-0.3,y[i]+0.25),fontsize=9,color='red',rotation='horizontal')

# display
plt.show()

# close plot
plt.cla()   # Clear axis
plt.clf()   # Clear figure
plt.close() # Close a figure window

# return data matrix
return WSTURB``````