jmquintana79
9/5/2017 - 6:37 AM

Wind Analysis Plot

Plot wind rose and wind distribution. Furthermore, return a Wind Speed Frequency Matrix.

"""
WIND ANALYSIS:

Plot the wind rose and wind distribution. It is required a DF with hwd and hws data.
"""
def analysis_wind(DATA,svar_hws,svar_hwd,stitle=''):
    
    ## 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 HWS FREQUENCY MATRIX """

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


    ## FINAL DF CREATION

    # create frequency matrix
    WSFREQ = pd.DataFrame.from_dict(dfreq, orient='columns')
    WSFREQ['index'] = [v[0] for k,v in dbin.items()]
    WSFREQ = WSFREQ.set_index('index')
    # add total per bins
    WSFREQ['ALL'] = WSFREQ.sum(axis=1)
    # create avg hws
    WSAVG = pd.DataFrame.from_dict(davg, orient='columns')
    WSAVG['ALL'] = WSAVG.mean(axis=1)
    WSAVG['index'] = ['avg_hws(m/s)']
    WSAVG = WSAVG.set_index('index')
    # create sum freq
    WSFREQSUM = pd.DataFrame(WSFREQ.sum(axis=0)).transpose()
    WSFREQSUM['index'] = ['sum_freq(%)']
    WSFREQSUM = WSFREQSUM.set_index('index')
    # append total
    WSFREQ = WSFREQ.append(WSAVG)
    WSFREQ = WSFREQ.append(WSFREQSUM)
    # clean
    del(WSAVG); del(WSFREQSUM)



    """ PLOT WIND ROSE """

    from math import pi, radians
    import matplotlib.pyplot as plt
    from matplotlib import rcParams
    from matplotlib.colors import ListedColormap
    import matplotlib as mpl

    # CREATE A TWIN POLAR AXIS
    def polar_twin(ax):
        ax2 = ax.figure.add_axes(ax.get_position(), projection='polar', 
                                 label='twin', frameon=False,
                                 theta_direction=ax.get_theta_direction(),
                                 theta_offset=ax.get_theta_offset())
        ax2.xaxis.set_visible(False)

        # There should be a method for this, but there isn't... Pull request?
        ax2._r_label_position._t = (22.5 + 180., 0.0)
        ax2._r_label_position.invalidate()

        # Bit of a hack to ensure that the original axes tick labels are on top of
        # whatever is plotted in the twinned axes. Tick labels will be drawn twice.
        for label in ax.get_yticklabels():
            ax.figure.texts.append(label)

        return ax2


    ## COLLECT DATA
    PLOT = WSFREQ.transpose()[['avg_hws(m/s)','sum_freq(%)']][:-1]
    PLOT.columns = ['avg','freq']
    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.avg.values)
    #freq = list(PLOT.freq.values)
    lsector = list(PLOT.index.values)
    # clean
    del(PLOT)

    # calculate total freq per sector for 4 defined bins
    # bin1: [0,4)  bin2: [4,8)  bin3: [8,12)  bin4: [12,inf)  
    ib1 = (0,3); ib2 = (4,7); ib3 = (8,11); ib4 = (12,15) 
    lfsector1 = list(); lfsector2 = list(); lfsector3 = list(); lfsector4 = list()
    for isector in lsector:
        lfsector1.append(np.nansum(WSFREQ[ib1[0]:ib1[1]+1][isector]))
        lfsector2.append(np.nansum(WSFREQ[ib2[0]:ib2[1]+1][isector]))
        lfsector3.append(np.nansum(WSFREQ[ib3[0]:ib3[1]+1][isector]))
        lfsector4.append(np.nansum(WSFREQ[ib4[0]:ib4[1]+1][isector]))

    # calculate ylim
    freq = np.array(lfsector1)+np.array(lfsector2)+np.array(lfsector3)+np.array(lfsector4)
    ylim = int(np.max(freq))+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,lfsector1,width, color='blue', alpha=0.8)
    ax.bar(thetaf,lfsector2,width,bottom=np.array(lfsector1), color='cyan', alpha=0.8)
    ax.bar(thetaf,lfsector3,width,bottom=np.array(lfsector1)+np.array(lfsector2), color='orange', alpha=0.8)
    ax.bar(thetaf,lfsector4,width,bottom=np.array(lfsector1)+np.array(lfsector2)+np.array(lfsector3), color='red', 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)))


    ## PLOT CHART OF Avg. WS

    # define second axis as twin
    ax2 = polar_twin(ax)
    # plot
    ax2.plot(theta, values, color='black', alpha=1.)
    # set limits
    ax2.set_ylim([0, ylim])
    # set axi tick labels
    plt.setp(ax2.get_yticklabels(), color='red')
    plt.yticks(range(0,ylim+4,4),fontsize=10,rotation='horizontal')
    ax2.set_yticklabels(list(range(0,ylim+4,4)))


    # set axis labels
    ax2.text(radians(22.5+7.5), ylim + 5, '(m/s)', size=10,color='blue', horizontalalignment="center", verticalalignment="center")
    ax.text(radians(202.5-7.5), ylim -0, '(%)', size=10,color='red', horizontalalignment="center", verticalalignment="center")

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



    ## LEGEND 

    # add new axis for legend object
    ax3 = fig.add_axes([0.5-0.25/2., 0.08, 0.25, 0.025])

    # value colors.
    cmap = mpl.colors.ListedColormap(['blue', 'cyan', 'orange','red'])
    cmap.set_over('0.25')
    # create boundary object
    bounds = [0, 4, 8, 12, 15]
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
    # create color bar
    cb = mpl.colorbar.ColorbarBase(ax3, cmap=cmap,
                                    norm=norm,
                                    boundaries=bounds + [13],
                                    extend='max',
                                    ticks=bounds,
                                    spacing='proportional',
                                    orientation='horizontal')
    # set label of color bar
    cb.set_label('(m/s)',fontsize=9,x=1.15,labelpad=-10)

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

    # display
    plt.show()



    """ PLOT BAR FREQ """

    import matplotlib.pyplot as plt

    ## COLLECT DATA
    BAR = WSFREQ[['ALL']].iloc[0:16]

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


    ## PLOT

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

    # 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])

    # 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("freq(%)",fontsize=10)

    # title
    if len(stitle)>0: title = ': %s'%stitle
    else: title = stitle
    plt.title('Wind Speed Distribution%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 WSFREQ