Changeset 1438 in lmdz_wrf for trunk/tools/drawing_tools.py

Timestamp:

Feb 8, 2017, 7:54:43 PM (8 years ago)

Author:

lfita

Message:

adding: plot_WindRose: Function to plot a wind rose (from where the dinw blows)

File:

• trunk/tools/drawing_tools.py (modified) (2 diffs)

trunk/tools/drawing_tools.py

@@ -83 +83 @@
 # plot_2lines_time: Function to plot two time-lines in different axes (x/x2 or y/y2)
 # plot_lines: Function to plot a collection of lines
+# plot_WindRose: Function to plot a wind rose (from where the dinw blows)
 # plot_ZQradii: Function to plot following radial averages only at exact grid poins
 # transform: Function to transform the values and the axes
@@ -8042 +8043 @@
 
     return
+
+def plot_WindRose(ang, speed, dime, lpvals, thetalabs, windu, titimg, kfig, figname, \
+  close, origfilen, outputfile=False, ev=None, eunit=None):
+    """ Function to plot a wind rose (from where the dinw blows)
+      ang= angle of the winds
+      speed= speed of the winds
+      Nang= number of angles to split the total circumference
+      Nspeed= number of speed sections to split the wind distribution
+      maxspeed= maximum wind speed to be plotted
+      dime= number of winds to plot
+      lpvals= values to configure the kind of Wind Rose to plot
+        lpvals[0]: kind of Wind Rose
+        lpvals[1]: specie of kind of Wind Rose
+        lpvals[>=2]: specific values for the specie Wind Rose
+        'anglespeedfreq';[Nang];[Nspeed];[maxspeed];[cbar];[maxfreq] grid of frequencies of each angle and speed 
+          following a given discretization
+        'linepoint';[kind];[val1];[...;[valN]]: consecutive (time, height, level, ...) line-point angle and speed valules
+          'multicol';[extravarn];[line];[marker];[colbar]
+          'multicoltime';[extravarn];[line];[marker];[colbar];[timekind];[timefmt];[timelabel]
+          'singlecol';[line];[marker];[col]
+        'scatter';[kind];[val1];[...;[valN]]: a cross for each wind at different values (time, height, level, ...)
+          'multicol';[extravarn];[marker];[colbar]
+          'multicoltime';[extravarn];[line];[marker];[colbar];[timekind];[timefmt];[timelabel]
+          'singlecol';[marker];[col]
+      windu= units of the wind speed
+      thetalabs= kind of labels for the angles of the wind-rose
+        'cardianals': Following combinations of 'N', 'E', 'S', 'W' according to Nang
+      titimg= title of the image
+      kfig= kind of figure
+      figname= name of the figure
+      close= wether figure has to be closed or not
+      origfilen= Name of the file from which the data cames from
+      outputfile= wether file with wind and angle frequencies has to be created (only valid for 'anglespeedfreq')
+      ev= extra values (for certain Wind Roses: 'linepoint/scatter' 'multicol' & 'multicoltime')
+      eunits= units of the extravalues (for certain Wind Roses)
+    """
+    fname = 'plot_WindRose'
+
+    if lpvals[0] == 'anglespeedfreq':
+        # Computing statistics
+        Nang = gen.auto_val(lpvals[1], 8)
+        Nspeed = gen.auto_val(lpvals[2], 8)
+        maxspeed = gen.auto_val(lpvals[3], 40.)
+        freqcbar = gen.auto_val(lpvals[4], 'spectral_r')
+        maxfreq = lpvals[5]
+
+        dang = 2.*np.pi / (Nang)
+        dang2 = dang/2.
+
+        nspeed = 0.
+        xspeed = maxspeed
+        dspeed = maxspeed/(Nspeed)
+
+        # Adding one extra partition for the last range
+        #Nang = Nang + 1
+        Nspeed = Nspeed + 1
+
+        # Angle localized
+        angfound = np.zeros((dime), dtype=bool)
+        # Angles's distirbution
+        angdistrib = np.zeros((Nang), dtype=np.float)
+        # Wind speed distirbution following angle
+        ang_speeddistrib = {}
+        # Wind speed distirbution following angle & speed
+        speedang_speeddistrib = {}
+        for it in range(dime):
+            for iang in range(Nang):
+                # initial/ending of the angle section
+                iasec = iang*dang - dang/2.
+                easec = iasec + dang
+    
+                if ang[it] >= iasec and ang[it] < easec: 
+                    Siang = str(iang)
+                    angfound[it] = True
+                    angdistrib[iang] = angdistrib[iang] + 1
+                    if ang_speeddistrib.has_key(iang):
+                        speeds = ang_speeddistrib[iang]
+                        ang_speeddistrib[iang] = speeds + [speed[iang]]
+                    else:
+                        ang_speeddistrib[iang] = [speed[it]]
+                    # initial/ending of the speed section
+                    for ispd in range(Nspeed):
+                        issec = ispd*dspeed - dspeed/2.
+                        essec = issec + dspeed
+                        if speed[it] >= issec and speed[it] < essec:
+                            Siaspd = Siang + '_' + str(ispd)
+                            if speedang_speeddistrib.has_key(Siaspd):
+                                speedangs = speedang_speeddistrib[Siaspd]
+                                speedang_speeddistrib[Siaspd]= speedangs + [speed[it]]
+                            else:
+                                speedang_speeddistrib[Siaspd]= [speed[it]]
+                            break
+            if not angfound[it]:
+                print warnmsg
+                print 'angle:',ang[it],': not located within:', dang*np.arange(Nang)-\
+                  dang/2.
+    
+    lkind = lpvals[1]
+
+    # Getting time values and ticks
+    if lkind == 'multicoltime':
+        if lpvals[0] == 'linepoint':
+            timekind = lpvals[7]
+            timefmt = lpvals[8]
+        elif lpvals[0] == 'scatter':
+            timekind = lpvals[6]
+            timefmt = lpvals[7]
+        timepos, timelabels = CFtimes_plot(ev, eunit, timekind, timefmt)
+    
+    if lpvals[0] == 'anglespeedfreq':
+        # Distribution of Frequencies
+        TOTwinds = 0
+        freqs = np.zeros((Nspeed,Nang), dtype=np.float)
+        for ian in range(Nang):
+            for isp in range(Nspeed):
+                Siaspd = str(ian) + '_' + str(isp)
+                if speedang_speeddistrib.has_key(Siaspd):
+                    Nwind = len(speedang_speeddistrib[Siaspd])
+                    freqs[isp,ian] = Nwind
+                    TOTwinds = TOTwinds + Nwind
+    
+        if np.sum(freqs) != TOTwinds:
+            print errormsg
+            print 'number of located frequencies:', freqs, 'and number of values:',  \
+              TOTwinds, 'differ!!'
+            quit(-1)
+
+        # Normalizing
+        freqs = freqs*1./TOTwinds
+        freqs = ma.masked_equal(freqs, 0.)
+    
+        # Filling with zeros up to the highest speed
+        for iang in range(Nang):
+            xfreq = freqs[:,iang].max()
+            if not np.all(freqs.mask[:,iang]):
+                indicesc = gen.multi_index_mat(freqs.mask[:,iang], False)
+                xind = np.max(indicesc)
+                freqs[0:xind+1,iang]= np.where(freqs.mask[0:xind+1,iang], 0.,        \
+                  freqs[0:xind+1,iang])
+
+        if type(maxfreq) == type('s'):
+            if maxfreq == 'auto': maxfreq = freqs.max()
+
+    # Plot
+    plt.rc('text', usetex=True)
+    ax = plt.subplot(111, projection='polar')
+
+    # Rosekind
+    # Scatter
+    if lpvals[0] == 'fill':
+        plt.fill_between(ang, speed, 0.)
+    elif lpvals[0] == 'anglespeedfreq':
+        angs = np.arange(Nang,dtype=np.float)*dang - dang/2.
+        spds = np.arange(Nspeed,dtype=np.float)*dspeed - dspeed/2.
+        spds[0] = 0.
+        x, y = gen.lonlat2D(angs,spds)
+        plt.pcolormesh(x, y, freqs, cmap=freqcbar, vmin=0., vmax=maxfreq)
+        cbar = plt.colorbar()
+        cbar.set_label('Frequencies (1)')
+        plt.grid()
+        if type(thetalabs) != type('auto'):
+            ax.set_thetagrids((angs+dang/2.)*180./np.pi, thetalabs)
+
+    elif lpvals[0] == 'linepoint':
+        if lkind == 'multicol':
+            ltyp = gen.auto_val(lpvals[3],'-')
+            lmrk = gen.auto_val(lpvals[4],'x')
+            colbar = gen.auto_val(lpvals[5],'spectral_r')
+            Nang = gen.auto_val(lpvals[6],8)
+            # Setting up colors for each label
+            #   From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar
+            vcolmin = np.min(ev)
+            vcolmax = np.max(ev)
+            my_cmap = plt.cm.get_cmap(colbar)
+            norm = mpl.colors.Normalize(vcolmin, vcolmax)
+
+            for ie in range(dime-1):
+                labcol = my_cmap(norm(ev[ie]))
+                plt.plot([ang[ie], ang[ie+1]], [speed[ie], speed[ie+1]], ltyp, marker=lmrk, color=labcol)
+            # FROM: http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
+            # For matplotlib v1.3 or greater the code becomes:
+            sm = plt.cm.ScalarMappable(cmap=colbar, norm=norm)
+            # fake up the array of the scalar mappable. Urgh...
+            sm._A = []
+            cbar = plt.colorbar(sm)
+            cbar.set_label(eunit)
+        elif lkind == 'multicoltime':
+            ltyp = gen.auto_val(lpvals[3],'-')
+            lmrk = gen.auto_val(lpvals[4],'x')
+            colbar = gen.auto_val(lpvals[5],'spectral_r')
+            Nang = gen.auto_val(lpvals[6],8)
+            timekind = lpvals[7]
+            timefmt = lpvals[8]
+            timelabel = lpvals[9].replace('!',' ')
+            # Setting up colors for each label
+            #   From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar
+            vcolmin = np.min(ev)
+            vcolmax = np.max(ev)
+            my_cmap = plt.cm.get_cmap(colbar)
+            norm = mpl.colors.Normalize(vcolmin, vcolmax)
+
+            for ie in range(dime-1):
+                labcol = my_cmap(norm(ev[ie]))
+                lcolS = gen.coldec_hex(labcol[0:3])
+                plt.plot([ang[ie], ang[ie+1]], [speed[ie], speed[ie+1]], ltyp, marker=lmrk, color=lcolS)
+    
+            # FROM: http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
+            # For matplotlib v1.3 or greater the code becomes:
+            sm = plt.cm.ScalarMappable(cmap=colbar, norm=norm)
+            # fake up the array of the scalar mappable. Urgh...
+            sm._A = []
+            cbar = plt.colorbar(sm, ticks=timepos)
+            cbar.set_label(timelabel)
+            cbar.ax.set_yticklabels(timelabels)
+
+            # An attempt using text...
+            #maxt = np.max(timepos)
+            #for it in range(len(timelabels)):
+            #    daxistT = 0.6
+            #    iaxistT = (1. - daxistT)/2.
+            #    itt = (timepos[it]-timepos[0])*daxistT/(maxt-timepos[0])
+            #    labcol = my_cmap(norm(timepos[it]))
+            #    plt.annotate(timelabels[it], xy=(0.87,iaxistT+itt), color=labcol,        \
+            #      xycoords='figure fraction', multialignment='center')
+            #plt.annotate(timelabel, xy=(0.97,0.5), color='k', rotation=90,               \
+            #  rotation_mode="anchor", xycoords='figure fraction', horizontalalignment='center')
+        elif lkind == 'singlecol':
+            ltyp = gen.auto_val(lpvals[2],'-')
+            lmrk = gen.auto_val(lpvals[3],'x')
+            lcol = gen.auto_val(lpvals[4],'b')
+            Nang = gen.auto_val(lpvals[5],8)
+            plt.plot(ang, speed, ltyp, marker=lmrk, color=lcol)
+    
+    elif lpvals[0] == 'scatter':
+        if lkind == 'multicol':
+            lmrk = gen.auto_val(lpvals[3],'x')
+            colbar = gen.auto_val(lpvals[4],'spectral_r')
+            Nang = gen.auto_val(lpvals[5],8)
+            plt.scatter(ang, speed, c=ev, cmap=colbar, marker=lmrk)
+            cbar = plt.colorbar()
+            cbar.set_label(eunit)
+        elif lkind == 'multicoltime':
+            lmrk = gen.auto_val(lpvals[3],'x')
+            colbar = gen.auto_val(lpvals[4],'spectral_r')
+            Nang = gen.auto_val(lpvals[5],8)
+            timekind = lpvals[6]
+            timefmt = lpvals[7]
+            timelabel = lpvals[8].replace('!',' ')
+            plt.scatter(ang, speed, c=ev, cmap=colbar, marker=lmrk)
+            cbar = plt.colorbar(ticks=timepos)
+            cbar.set_label(timelabel)
+            cbar.ax.set_yticklabels(timelabels)
+        elif lkind == 'singlecol':
+            lmrk = gen.auto_val(lpvals[2],'x')
+            lcol = gen.auto_val(lpvals[3],'b')
+            Nang = gen.auto_val(lpvals[4],8)
+            plt.plot(ang, speed, ',', marker=lmrk, color=lcol)
+    
+    if thetalabs == 'cardinals':
+        if Nang == 4:
+            thetalabs = ['N', 'E', 'S', 'W']
+        elif Nang == 8:
+            thetalabs = ['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']
+        elif Nang == 16:
+            thetalabs= ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE', 'S', 'SSW',\
+              'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW']
+        else:
+            print errormsg
+            print '  ' + fname + ': number of angles:', Nang, "not ready for '" +    \
+              thetalabs + "' !!"
+            print '    available ones: 4, 8, 16'
+            quit(-1)
+
+        dang = 2.*np.pi / (Nang)        
+        angs = np.arange(Nang,dtype=np.float)*dang - dang/2.
+        ax.set_thetagrids((angs+dang/2.)*180./np.pi, thetalabs)
+    
+    ax.set_theta_zero_location('N')
+    ax.set_theta_direction(-1)
+    plt.annotate('wind (' + units_lunits(windu) + ')', xy=(0.75,0.06), color='k',    \
+      xycoords='figure fraction', horizontalalignment='center')
+
+    plt.title(gen.latex_text(titimg))
+
+    output_kind(kfig, figname, close)
+
+    # Creation of the output file
+    if lpvals[0] == 'anglespeedfreq' and outputfile:
+        print '  ' + fname + ': Writting netCDF of frequencies of wind angle and ' + \
+          'speeds...'
+        newnc = NetCDFFile(figname + '.nc', 'w')
+
+        # Dimension creation
+        newdim = newnc.createDimension('angle', Nang)
+        newdim = newnc.createDimension('speed', Nspeed)
+        newdim = newnc.createDimension('bounds', 2)
+
+        # Dimension variable
+        newvar = newnc.createVariable('angle', 'f8', ('angle'))
+        ncvar.basicvardef(newvar, 'angle', 'angle from which wind blows', 'degree')
+        newvar[:] = angs
+        newvar.setncattr('cell_method','angle_bounds')
+
+        newvar = newnc.createVariable('angle_bounds', 'f8', ('angle','bounds'))
+        ncvar.basicvardef(newvar,'angle_bounds', 'bounds of the angle sections',     \
+         'degree')
+        for iang in range(Nang):
+            # initial/ending of the angle section
+            iasec = iang*dang - dang/2.
+            if iasec < 0.: iasec = 2.*np.pi + iasec
+            easec = iasec + dang
+            newvar[iang,0] = iasec*180./np.pi
+            newvar[iang,1] = easec*180./np.pi
+        newnc.sync()
+
+        newvar = newnc.createVariable('speed', 'f8', ('speed'))
+        ncvar.basicvardef(newvar, 'speed', 'speed of wind', windu)
+        newvar[:] = spds
+        newvar.setncattr('cell_method','speed_bounds')
+
+        newvar = newnc.createVariable('speed_bounds', 'f8', ('speed', 'bounds'))
+        ncvar.basicvardef(newvar,'speed_bounds','bounds of the speed sections',windu)
+        for ispd in range(Nspeed):
+            # initial/ending of the speed section
+            issec = ispd*dspeed - dspeed/2.
+            if issec < 0.: issec = 0.
+            essec = issec + dspeed
+            newvar[ispd,0] = issec
+            newvar[ispd,1] = essec
+        newnc.sync()
+
+        # Variables
+        newvar = newnc.createVariable('frequency', 'f4', ('speed', 'angle'),         \
+          fill_value=gen.fillValueF)
+        ncvar.basicvardef(newvar, 'frequency', 'frequency of angle and speed', '1')
+        newvar[:] = freqs[:]
+        newnc.sync()
+
+        # Global values
+        newnc.setncattr('author', 'L. Fita')
+        newattr = ncvar.set_attributek(newnc, 'institution', unicode('Laboratoire ' +\
+          'de M' + unichr(233) + 't' + unichr(233) + 'orologie Dynamique'), 'U')
+        newnc.setncattr('university', 'Pierre Marie Curie - Jussieu')
+        newnc.setncattr('center', 'Centre National de Recherches Scientifiques')
+        newnc.setncattr('city', 'Paris')
+        newnc.setncattr('country', 'France')
+        newnc.setncattr('script', 'nc_var_tools.py')
+        newnc.setncattr('function', fname)
+        newnc.setncattr('version', '1.0')
+        newnc.setncattr('original_file', origfilen)
+
+        newnc.sync()
+        newnc.close()
+        print "Successful creation of file with frequencies '" + figname + ".nc' !!"    
+
+    return