Changeset 652 in lmdz_wrf for trunk/tools

Timestamp:

Sep 23, 2015, 11:28:59 AM (10 years ago)

Author:

lfita

Message:

Adding `plot_ptZvals' to plot lon,lat stations at a given time
Slicing also dimensions on the '2D_shad' plot
Adding monotone removal on `slice_variable'

File:

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

trunk/tools/drawing_tools.py

@@ -203 +203 @@
           * -1: all along the dimension
           * [beg]:[end] slice from [beg] to [end]
+          * -9: last value of the dimension
+
     """
     fname = 'slice_variable'
@@ -223 +225 @@
     varvalsdim = []
     dimnslice = []
-
+    monodim = []
     for idd in range(Ndimvar):
         found = False
@@ -241 +243 @@
                         dimnslice.append(vardims[idd])
                     elif int(dimcutv) == -9:
-                        varvalsdim.append(int(varobj.shape[idd])-1)
+                        varvalsdim.append(varobj.shape[idd]-1)
+                        monodim.append(vardims[idd])
                     else:
                         varvalsdim.append(int(dimcutv))
+                        monodim.append(vardims[idd])
                 found = True
                 break
-        if not found and not searchInlist(dimnslice,vardims[idd]):
+        if not found and not searchInlist(dimnslice,vardims[idd]) and                \
+          not searchInlist(monodim,vardims[idd]):
             varvalsdim.append(slice(0,varobj.shape[idd]))
             dimnslice.append(vardims[idd])
     varvalues = varobj[tuple(varvalsdim)]
+
+    varvalues = np.squeeze(varobj[tuple(varvalsdim)])
 
     return varvalues, dimnslice
@@ -3507 +3514 @@
     if len(varsv.shape) != 2:
         print errormsg
-        print '  ' + fname + ': wrong variable shape:',varv.shape,'is has to be 2D!!'
+        print '  ' + fname + ': wrong variable shape:',varsv.shape,'is has to be 2D!!'
         quit(-1)
 
@@ -3547 +3554 @@
         ydims = '0'
 
-    lon0 = dimxv
-    lat0 = dimyv
-#    lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)
+#    lon0 = dimxv
+#    lat0 = dimyv
+    lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)
 
     if not mapv is None:
@@ -4711 +4718 @@
         [val]: value (-1 for all the range)
     """
-
     fname = 'dxdy_lonlatDIMS'
+
+    print 'Lluis dd:',dd
 
     slicex = []
@@ -4719 +4727 @@
         for idd in range(len(dd)):
             dname = dd[idd].split('|')[0]
-            dvalue = int(dd[idd].split('|')[1])
+            dvalue = dd[idd].split('|')[1]
             if dn == dname:
-                if dvalue == -1:
-                    slicex.append(slice(0,dxv.shape[ipos]))
+                if dvalue.find('@') != -1:
+                    slicex.append(slice(int(dvalue.split('@')[0]),                   \
+                      int(dvalue.split('@')[1])))
                 else:
-                    slicex.append(dvalue)
-                break
+                    if int(dvalue) == -1:
+                        slicex.append(slice(0,dxv.shape[ipos]))
+                    elif int(dvalue) == -9:
+                        slicex.append(dxv.shape[ipos]-1)
+                    else:
+                        slicex.append(int(dvalue))
+                    break
         ipos = ipos + 1
 
@@ -4733 +4747 @@
         for idd in range(len(dd)):
             dname = dd[idd].split('|')[0]
-            dvalue = int(dd[idd].split('|')[1])
+            dvalue = dd[idd].split('|')[1]
             if dn == dname:
-                if dvalue == -1:
-                    slicey.append(slice(0,dyv.shape[ipos]))
+                if dvalue.find('@') != -1:
+                    slicey.append(slice(int(dvalue.split('@')[0]),                   \
+                      int(dvalue.split('@')[1])))
                 else:
-                    slicey.append(dvalue)
-                break
+                    if int(dvalue) == -1:
+                        slicey.append(slice(0,dyv.shape[ipos]))
+                    elif int(dvalue) == -9:
+                        slicey.append(dyv.shape[ipos]-1)
+                    else:
+                        slicey.append(int(dvalue))
+                    break
         ipos = ipos + 1
-
 
     lonv = dxv[tuple(slicex)]
@@ -5493 +5512 @@
     return
 
+def plot_ptZvals(vname,vunits,points,ptype,ptsize,graphlims,minmax,figtitle,cbar,    \
+  mapv,kfig):
+    """ Function to plot a given list of points and values 
+      vname= name of the variable in the graph
+      vunits= units of the variable
+      points= [lon,lat,val] matrix of values
+      ptype= type of the point
+      ptsize= size of the point
+      graphlims= minLON,minLAT,maxLON,maxLAT limits of the graph, None for the full size
+      minmax= minimum and maximum type
+        'auto': values taken from the extrems of the data
+        [min],[max]: given minimum and maximum values
+      figtitle= title of the figure
+      cbar= color bar
+      mapv= map characteristics: [proj],[res]
+        see full documentation: http://matplotlib.org/basemap/
+        [proj]: projection
+          * 'cyl', cilindric
+          * 'lcc', lambert-conformal
+        [res]: resolution:
+          * 'c', crude
+          * 'l', low
+          * 'i', intermediate
+          * 'h', high
+          * 'f', full
+      kfig= kind of figure
+    """
+    fname = 'plot_ptZvals'
+
+    figname = 'pointsZval'
+
+    minlon = points[:,0].min()
+    maxlon = points[:,0].max() 
+
+    minlat = points[:,1].min()
+    maxlat = points[:,1].max()
+
+    minval = points[:,2].min()
+    maxval = points[:,2].max()
+
+#    print 'min/max val;',minval,maxval
+
+    lonrange = (points[:,0] - minlon)/(maxlon - minlon)
+    latrange = (points[:,1] - minlat)/(maxlat - minlat)
+    colorrange = (points[:,2] - minval)/(maxval - minval)
+
+    plt.rc('text', usetex=True)
+
+    if mapv is not None:
+        vlon = points[:,0]
+        vlat = points[:,1]
+        dx = len(vlon)
+        dy = len(vlat)
+
+#        vlon = np.where(vlon[:] < 0., 360. + vlon[:], vlon[:])
+#        xvala = np.array(xval)
+#        xvala = np.where(xvala < 0., 360. + xvala, xvala)
+#        xval = list(xvala)
+
+        map_proj=mapv.split(',')[0]
+        map_res=mapv.split(',')[1]
+
+        if graphlims is not None:
+            nlon = graphlims[0]
+            xlon = graphlims[2]
+            nlat = graphlims[1]
+            xlat = graphlims[3]
+        else:
+            nlon = np.min(vlon)
+            xlon = np.max(vlon)
+            nlat = np.min(vlat)
+            xlat = np.max(vlat)
+
+        lon2 = vlon[dy/2]
+        lat2 = vlat[dy/2]
+
+        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
+          xlon, ',', xlat
+
+        if map_proj == 'cyl':
+            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
+              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        elif map_proj == 'lcc':
+            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
+              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
+        else:
+            print errormsg
+            print '  ' + fname + ": map projecion '" + map_proj + "' not ready!!"
+            print '    available: cyl, lcc'
+            quit(-1)
+
+#        lons, lats = np.meshgrid(vlon, vlat)
+#        lons = np.where(lons < 0., lons + 360., lons)
+
+        x,y = m(vlon,vlat)
+
+        m.drawcoastlines()
+
+        meridians = pretty_int(nlon,xlon,5)
+        m.drawmeridians(meridians,labels=[True,False,False,True])
+
+        parallels = pretty_int(nlat,xlat,5)
+        m.drawparallels(parallels,labels=[False,True,True,False])
+#    else:
+#        x = vlon
+#        y = vlat
+#        plt.xlim(0,dx-1)
+#        plt.ylim(0,dy-1)
+ 
+    if minmax == 'auto':
+        plt.scatter(points[:,0], points[:,1], c=points[:,2], s=ptsize, cmap=cbar,    \
+          marker=ptype)
+    else:
+        minv = np.float(minmax.split(',')[0])
+        maxv = np.float(minmax.split(',')[1])
+
+        plt.scatter(points[:,0], points[:,1], c=points[:,2], s=ptsize, cmap=cbar,    \
+          marker=ptype, vmin=minv, vmax=maxv)
+
+    cbar = plt.colorbar()
+    cbar.set_label(vname.replace('_','\_') +' ('+ units_lunits(vunits) + ')')
+
+    plt.title(figtitle)
+    if graphlims is not None:
+        plt.xlim(graphlims[0], graphlims[2])
+        plt.ylim(graphlims[1], graphlims[3])
+
+    output_kind(kfig, figname, True)
+
+    return
+
+#pts = np.zeros((10,3), dtype=np.float)
+#pts[:,0] = np.arange(10,20)*1.
+#pts[:,1] = np.arange(30,40)*1.
+#pts[:,2] = np.arange(-5,5)*1.
+
+#plot_ptZvals('vals','kgm-2',pts,'.',300, 'values of values', 'seismic', 'cyl,l', 'pdf')
+
 def plot_ZQradii(Zmeans, graphtit, kfig, figname):
     """ Function to plot following radial averages only at exact grid poins