Cacher is the code snippet organizer for pro developers

"""
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')
    # add total per bins
    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 """

    from math import pi, radians
    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