Changeset 545 in lmdz_wrf for trunk

Timestamp:

Jun 30, 2015, 6:24:47 PM (9 years ago)

Author:

lfita

Message:

Adding `draw_points' to draw points on a shadded map

Location:

trunk/tools

Files:

• drawing.py (modified) (4 diffs)
• drawing_tools.py (modified) (6 diffs)

trunk/tools/drawing.py

@@ -31 +31 @@
 namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time',                 \
   'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line',                \
-  'draw_2D_shad_line_time', 'draw_barbs', 'draw_timeSeries', 'draw_topo_geogrid',    \
+  'draw_2D_shad_line_time', 'draw_barbs', 'draw_points', 'draw_timeSeries',          \
+   'draw_topo_geogrid',                                                              \
   'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
   'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', 'list_graphics',      \
@@ -2182 +2183 @@
       timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
 
+def draw_points(filen, values):
+    """ Function to plot a series of points
+      [values]= [ptasciifile]:[gtit]:[mapvalues]:[pointcolor]:[pointlabels]:[locleg]:
+        [figurek]:[figuren]
+        [ptasciifile]:[file],[comchar],[collon],[collat],[lab]
+          [file]: column ASCII file with the location of the points
+          [comchar]: '|' list of characters for commentaries
+          [collon]: number of column with the longitude of the points
+          [collat]: number of column with the latitude of the points
+          [collab]: number of column with the labels of the points ('None', and will get
+            the values from the [pointlabels] variable
+        [gtit]: title of the figure ('|' for spaces)
+        [mapvalues]: drawing coastaline ([proj],[res]) or None
+          [proj]: projection
+             * 'cyl', cilindric
+             * 'lcc', lambert conformal
+          [res]: resolution:
+             * 'c', crude
+             * 'l', low
+             * 'i', intermediate
+             * 'h', high
+             * 'f', full
+        [pointcolor]: color for the points ('auto' for "red")
+        [pointlabels]: ',' of labels [only used if [collab]='None'] ('None' for no labels)
+        [locleg]: location of the legend (0, autmoatic)
+          1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
+          5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
+          9: 'upper center', 10: 'center'
+        [figurek]: kind of figure
+        [figuren]: name of the figure
+      [filen]= [ncfile],[lonvarn],[latvarn][,[varn],[dimvals],[vargn],[min],[max],[cbar],[varu]]
+        [ncfile]: netCDF to use to geolocalize the points
+        [lonvarn]: name of the variable with the longitudes
+        [latvarn]: name of the variable with the latitudes
+        [varn]: optional variable to add staff into the graph
+        [dimval]: '@' list of [dimn]|[dimval] to get the values for [varn]
+          [dimn]: name of the dimension
+          [dimval]: value of the dimension (no value all range)
+        [vargn]: name of the variable in the graph
+        [min]: minimum value for the extra variable
+        [max]: maximum value for the extra variable
+        [cbar]: color bar
+        [varu]: units of the variable
+    """
+    fname = 'draw_points'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print draw_points.__doc__
+        quit()
+
+    expectargs = '[ptasciifile]:[gtit]:[mapvalues]:[pointcolor]:[pointlabels]:' +    \
+      '[locleg]:[figurek]:[figuren]'
+ 
+    drw.check_arguments(fname,len(expectargs.split(':')),values,':',                 \
+      expectargs.split(':'))
+
+    ptasciifile = values.split(':')[0]
+    gtit = values.split(':')[1]
+    mapvalues = values.split(':')[2]
+    pointcolor = values.split(':')[3]
+    pointlabels = values.split(':')[4]
+    locleg = int(values.split(':')[5])
+    figurek = values.split(':')[6]
+    figuren = values.split(':')[7]
+
+# Getting station locations
+## 
+    filev = ptasciifile.split(',')[0]
+    comchar = ptasciifile.split(',')[1].split('|')
+    collon = int(ptasciifile.split(',')[2])
+    collat = int(ptasciifile.split(',')[3])
+    collab = ptasciifile.split(',')[4]
+
+    if not os.path.isfile(filev):
+        print errormsg
+        print '  ' + fname + ": file '" + filev + "' does not exist!!"
+        quit(-1)
+
+# Getting points position and labels
+    oascii = open(filev, 'r')
+    xptval = []
+    yptval = []
+    if collab != 'None':
+        ptlabels = []
+        for line in oascii:
+            if not drw.searchInlist(comchar, line[0:1]):
+                linevals = drw.reduce_spaces(line)
+                xptval.append(np.float(linevals[collon].replace('\n','')))
+                yptval.append(np.float(linevals[collat].replace('\n','')))
+                ptlabels.append(linevals[int(collab)].replace('\n',''))
+    else:
+        ptlabels = None
+        for line in oascii:
+            if  not drw.searchInlist(comchar, line[0:1]):
+                linevals = drw.reduce_spaces(line)
+                xptval.append(np.float(linevals[collon].replace('\n','')))
+                yptval.append(np.float(linevals[collat].replace('\n','')))
+
+    oascii.close()
+
+    if pointlabels != 'None' and collab == 'None':
+        ptlabels = pointlabels.split(',')
+
+# Getting localization of the points
+##
+    filev = filen.split(',')
+    Nvals = len(filev)
+
+    ncfile = filev[0]
+    lonvarn = filev[1]
+    latvarn = filev[2]
+    varn = None
+    varextrav = None
+    if Nvals == 10:
+        varn = filev[3]
+        dimvals = filev[4]
+        varextrav = [filev[5], np.float(filev[6]), np.float(filev[7]), filev[8],     \
+          filev[9]]
+    
+    if not os.path.isfile(ncfile):
+        print errormsg
+        print '  ' + fname + ": file '" + ncfile + "' does not exist!!"
+        quit(-1)
+
+    onc = NetCDFFile(ncfile, 'r')
+    
+    lonv, latv = drw.lonlat2D(onc.variables[lonvarn], onc.variables[latvarn])
+
+    print 'Lluis varn',varn
+    if varn is not None:
+        objextra = onc.variables[varn]
+        vard = objextra.dimensions
+        dd = {}
+        for dn in dimvals.split('@'):
+            ddn = dn.split('|')[0]
+            ddv = dn.split('|')[1]
+            dd[ddn] = ddv
+        slicevar = []
+        for dv in vard:
+            found= False
+            for dn in dd.keys():
+                if dn == dv:
+                    slicevar.append(int(dd[dn]))
+                    found = True
+                    break
+            if not found:
+                slicevar.append(slice(0,len(onc.dimensions[dv])))
+
+        varextra = np.squeeze(objextra[tuple(slicevar)])
+
+    drw.plot_points(xptval, yptval, lonv, latv, varextra, varextrav, gtit, mapvalues,\
+      pointcolor, ptlabels, locleg, figurek, figuren)
+
+    onc.close()
+
+    return
+
 def draw_timeSeries(filen, values, variables):
     """ Function to draw a time-series
@@ -2641 +2800 @@
 Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time',                    \
   'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time',    \
-  'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
-  'variable_values']
+  'draw_points', 'draw_topo_geogrid_boxes', 'draw_trajectories',                     \
+  'draw_vals_trajectories', 'variable_values']
 
 ####### ###### ##### #### ### ## #
@@ -2678 +2837 @@
 elif oper == 'draw_lines_time':
     draw_lines_time(opts.ncfile, opts.values, opts.varname)
+elif oper == 'draw_points':
+    draw_points(opts.ncfile, opts.values)
 elif oper == 'draw_timeSeries':
     draw_timeSeries(opts.ncfile, opts.values, opts.varname)

trunk/tools/drawing_tools.py

@@ -58 +58 @@
 
 # From nc_var_tools.py
+def reduce_spaces(string):
+    """ Function to give words of a line of text removing any extra space
+    """
+    values = string.replace('\n','').split(' ')
+    vals = []
+    for val in values:
+         if len(val) > 0:
+             vals.append(val)
+
+    return vals
+
 def searchInlist(listname, nameFind):
     """ Function to search a value within a list
@@ -1368 +1379 @@
     return varvals
 
+def lonlat2D(lon,lat):
+    """ Function to return lon, lat 2D matrices from any lon,lat matrix
+      lon= matrix with longitude values
+      lat= matrix with latitude values
+    """
+    fname = 'lonlat2D'
+
+    if len(lon.shape) != len(lat.shape):
+        print errormsg
+        print '  ' + fname + ': longitude values with shape:', lon.shape,            \
+          'is different that latitude values with shape:', lat.shape, '(dif. size) !!'
+        quit(-1)
+
+    if len(lon.shape) == 3:
+        lonvv = lon[0,:,:]
+        latvv = lat[0,:,:]
+    elif len(lon.shape) == 2:
+        lonvv = lon[:]
+        latvv = lat[:]
+    elif len(lon.shape) == 1:
+        lonlatv = np.meshgrid(lon[:],lat[:])
+        lonvv = lonlatv[0]
+        latvv = lonlatv[1]
+
+    return lonvv, latvv
+
 #######    #######    #######    #######    #######    #######    #######    #######    #######    #######
 
@@ -2500 +2537 @@
 #quit()
 
-def plot_points(xval, yval, ifile, vtit, kfig, mapv):
+def plot_points(xval, yval, vlon, vlat, extravals, extrapar, vtit, mapv, color,      \
+  labels, lloc, kfig, figname):
     """ plotting points
       [x/yval]: x,y values to plot
-      vn,vm= minmum and maximum values to plot
-      unit= units of the variable
-      ifile= name of the input file
-      vtit= title of the variable
-      kfig= kind of figure (jpg, pdf, png)
+      vlon= 2D-matrix with the longitudes
+      vlat= 2D-matrix with the latitudes
+      extravals= extra values to be added into the plot (None for nothing)
+      extrapar= [varname, min, max, cbar, varunits] of the extra variable
+      vtit= title of the graph ('|' for spaces)
       mapv= map characteristics: [proj],[res]
         see full documentation: http://matplotlib.org/basemap/
         [proj]: projection
           * 'cyl', cilindric
+          * 'lcc', lambert-conformal
         [res]: resolution:
           * 'c', crude
@@ -2518 +2557 @@
           * 'h', high
           * 'f', full
+      color= color for the points/labels ('auto', for "red")
+      labels= list of labels for the points (None, no labels)
+      lloc = localisation of the legend
+      kfig= kind of figure (jpg, pdf, png)
+      figname= name of the figure
+
     """
     fname = 'plot_points'
-
-    dx=xval.shape[0]
-    dy=yval.shape[0]
+# Canging line kinds every 7 pts (end of standard colors)
+    ptkinds=['.','x','o','*','+','8','>','D','h','p','s']
+
+    Npts = len(xval)
+    if Npts > len(ptkinds)*7:
+        print errormsg
+        print '  ' + fname + ': too many',Npts,'points!!'
+        print "    enlarge 'ptkinds' list"
+        quit(-1)
+
+    N7pts = 0
+
+    if color == 'auto':
+        ptcol = "red"
+    else:
+        ptcol = color
+
+    dx=vlon.shape[1]
+    dy=vlat.shape[0]
 
     plt.rc('text', usetex=True)
 
     if not mapv is None:
-        lon00 = np.where(xval[:] < 0., 360. + olon[:], olon[:])
-        lat00 = yval[:]
-        lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
-        lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
-
-        for iy in range(len(lat00)):
-            lon0[iy,:] = lon00
-        for ix in range(len(lon00)):
-            lat0[:,ix] = lat00
+        vlon = np.where(vlon[:] < 0., 360. + vlon[:], vlon[:])
 
         map_proj=mapv.split(',')[0]
         map_res=mapv.split(',')[1]
 
-        nlon = lon0[0,0]
-        xlon = lon0[dy-1,dx-1]
-        nlat = lat0[0,0]
-        xlat = lat0[dy-1,dx-1]
-
-        lon2 = lon0[dy/2,dx/2]
-        lat2 = lat0[dy/2,dx/2]
+        nlon = np.min(vlon)
+        xlon = np.max(vlon)
+        nlat = np.min(vlat)
+        xlat = np.max(vlat)
+
+        lon2 = vlon[dy/2,dx/2]
+        lat2 = vlat[dy/2,dx/2]
 
         print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
@@ -2558 +2611 @@
               llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
 
-        lons, lats = np.meshgrid(olon[:], olat[:])
-        lons = np.where(lons < 0., lons + 360., lons)
-
-        x,y = m(lons,lats)
-        plt.plot(x,y)
+#        lons, lats = np.meshgrid(vlon, vlat)
+#        lons = np.where(lons < 0., lons + 360., lons)
+
+        x,y = m(vlon,vlat)
 
         m.drawcoastlines()
@@ -2571 +2623 @@
         parallels = pretty_int(nlat,xlat,5)
         m.drawparallels(parallels,labels=[False,True,True,False])
-    else:
+#    else:
 #        plt.xlim(0,dx-1)
 #        plt.ylim(0,dy-1)
 
-        plt.plot(xval, yval, '.')
-
-    figname = ifile.replace('.','_') + '_' + vtit
+# Extra values
+    if extravals is not None:
+        plt.pcolormesh(x, y, extravals, cmap=plt.get_cmap(extrapar[3]),              \
+          vmin=extrapar[1], vmax=extrapar[2])
+        cbar = plt.colorbar()
+        cbar.set_label(extrapar[0].replace('_','\_') +'('+ units_lunits(extrapar[4])+\
+          ')')
+
+    if labels is not None:
+        for iv in range(len(xval)):
+#            plt.annotate(labels[iv], xy=(xval[iv],yval[iv]), xycoords='data')
+
+            if np.mod(iv,7) == 0: N7pts = N7pts + 1
+#            print iv,xval[iv],yval[iv],labels[iv],ptkinds[N7pts]
+
+            plt.plot(xval[iv],yval[iv], ptkinds[N7pts],label=labels[iv])
+        plt.legend(loc=lloc)
+
+    else: 
+        plt.plot(xval, yval, '.', color=ptcol)
+
     graphtit = vtit.replace('_','\_').replace('&','\&')
 
-    plt.title(graphtit)
+    plt.title(graphtit.replace('|', ' '))
     
     output_kind(kfig, figname, True)